Investigation of asphalt mixture internal structure consistency in accelerated discrete element models

Construction and Building Materials 244 (2020) 118272

Contents lists available at ScienceDirect

Construction and Building Materials journal homepage: www.elsevier.com/locate/conbuildmat

Investigation of asphalt mixture internal structure consistency in accelerated discrete element models Xiaodong Zhou, Siyu Chen, Dongdong Ge, Dongzhao Jin, Zhanping You ⇑ Department of Civil and Environmental Engineering, Michigan Technological University, Houghton, MI 49931, USA

h i g h l i g h t s
This study proposed a set of equations to calculate a material modulus reduction range, that is used to accelerate the discrete element method (DEM)

simulation of asphalt mixture compaction and meanwhile maintain simulation accuracy.
A three-step compaction process was proposed to compare the specimens. The results exhibited the DEM simulation agreed with the laboratory results

for specimens of three gradation designs under three compaction states.
Internal structure indexes in respect of specimen’s height, rotation angles of coarse aggregates, and spatial distribution of mastic particles were proposed

to evaluate the effects of the reduced material modulus.

a r t i c l e

i n f o

Article history: Received 14 October 2019 Received in revised form 21 January 2020 Accepted 22 January 2020

Keywords: Discrete element method Asphalt mixture compaction Speedup Material modulus reduction Internal-structure

a b s t r a c t The discrete element method (DEM) requires an enormous amount of computational resources when applied to asphalt mixture simulation. Reducing material modulus is recognized as an efficient method to cut down the computational cost. However, over-reduced material modulus would case unacceptable calculation errors for DEM models, and the allowable modulus reduction range is not clear. Based on the small displacement assumption in the DEM modeling, the overlap ratio between particles should be small enough to ensure efficient force transition and discrete system stability. Combining this assumption with time step determination and force-displacement law, a theoretical inference that specified a material modulus reduction range without causing an exceeded overlap ratio can be derived. To verify this theoretical inference, a three-step DEM compaction model, including gravity fall, static compaction, and gyratory compaction processes was established for asphalt mixtures. Three groups of asphalt mixtures with different mixture designs were tested both in the laboratory and in simulation. To evaluate the effects of the reduced material modulus on the internal-structure of asphalt mixture specimen under three compaction states, this study proposed 4 categories of internal-structure indexes in respect of specimen’s height, average coordination number, rotation angles of coarse aggregates, and spatial distribution of mastic particles. In the presented case, when the modulus reduction was 1/100 times, the maximum simulation error of internal-structure indexes (except the average coordinate number) was around 4% compared to the control group. Meanwhile, the model’s calculation efficiency was increased by 8–10 times depending on the mixture design. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction The image process is a major approach to evaluate the internal structure of the asphalt mixture in the past decades. Researchers utilized this approach to identify coarse aggregates, asphalt mastic, and voids by analyzing the images captured from sliced specimens and X-ray tomography. Internal-structure indexes were proposed ⇑ Corresponding author. E-mail address: [email protected] (Z. You). https://doi.org/10.1016/j.conbuildmat.2020.118272 0950-0618/Ó 2020 Elsevier Ltd. All rights reserved.

to evaluate and quantify the internal structure of the asphalt mixture [1]. Aggregate orientation, gradation change, and air voids were used to compare the specimens from gyratory compaction with field cores [2]. The analysis results of internal-structure indexes were used to improve the simulation of gyratory compaction to field compaction [3]. Three internal-structure indexes, including contact points, aggregate orientation, and spatial distribution, were proposed by Coenen, Kutay et al. to characterize the aggregate structure and volumetric properties of asphalt mixture [4]. The proposed indexes gave insight investigation to the aggre-

2

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

gate skeleton and mixture compaction. The aggregate gradation is proved to have direct impacts on the structure evolution of asphalt mixture. For further investigation, Sefidmazgi, Tashman et al. utilized the image process to characterize the rutting performance in the view of internal-structure indexes [5]. Xu, Guo et al. evaluated the structure evolution during the froze-tharw cycle with internal-structure indexes [6]. Numerical simulation was also used to characterize the internal structure of the asphalt mixture [7,8]. Compared with the image process, the numerical simulation has fewer needs on the costly test procedures, for example, the X-ray tomography. Meanwhile, the numerical simulation can provide more comprehensive data of the internal-structure. The discrete element method (DEM) was developed to solve granular matter problems [9] and has become one of the most effective approaches to simulate asphalt mixtures [10–12]. Mahmoud, Masad et al. employ the DEM to analyze the impacts of aggregate properties and inter-structure on the fracture behavior of asphalt mixture [13]. In this study, different gradation types and aggregates were used, and the results exhibited that the mixture fracture strength depends on aggregate strength and blending limit of the mixture design. The simulation of compacted asphalt mixture specimens received a huge amount of attention from researchers [14,15]. The mechanical behavior of compacted asphalt mixture specimens had been investigated by studies using various contact models [16–19]. Moreover, its modeling methods had been extended from 2-dimensional (2D) [20– 22] to 3-dimensional (3D) [23,24]. In addition, to represent the aggregate interlock effect, the user-defined aggregate shapes [25–27] and realistic coarse aggregate shapes modeling methods [28–30] were developed. However, the compaction of asphalt mixture, which is a critical procedure for the formation of asphalt pavement structures, received less attention due to its complex dynamic movement and a huge amount of computational needs. Some researchers established compaction models based on the methods mentioned above. For example, the movement of coarse aggregates was investigated through static compaction which applied a constant force on top of specimens [25]. These models did not change the basic model geometries and should be recognized as the extent applications of the compacted models. New generation methods and model geometries are needed to fully simulate the compaction process. A new method called flow test was proposed to describe the workability of the asphalt mixture, and the DEM was used to investigate the aggregate movement [31,32]. This test provides an approach to quantitatively measure the workability of asphalt mixture, however, the DEM simulation in these studies contained a small number of elements and employed ball as the basic elements. Another early study investigated the application of DEM to the SUPERPAVE Gyratory Compaction (SGC) test and proposed a parameter determination method for asphalt mixtures at high temperatures [33]. Following this study, more complex gyratory compaction models with realistic aggregate shapes were carried out to study the effects of aggregate shapes on the compactibility of the asphalt mixture [34,35]. One of the biggest barriers lingering in these new models from wide applications is their need for a large number of computational resources. To increase computational efficiency, researchers proposed to reduce the material modulus in DEM models to increase the time interval of calculation (timestep) [36]. The timestep in the DEM is the discrete time between two calculation cycles. It is set to limit the maxim relative movement of particles. Inappropriate large timestep leads to extreme overlap and incongruent movement between particles. The timestep is positive related to the contact stiffness and negative correlated with mass of particles. The larger timestep enables models to solve asphalt mixture flow with fewer calculation cycles, leading to less computational costs. Moreover,

when the materials are applied to large moduli, the discrete system could be disturbed by small stimulation [37,38]. Thus, the smaller material moduli also benefit to the stability of internal structure and mechanical response. Studies focused on the collisions between discrete particles in large space supported this speedup method [39,40]. However, the effects of reduced modulus on the compaction remain controversial. Some researchers reported that the stiffness of particles in both normal and shear directions would greatly influence the simulation results [41,42], while some others reported that the effects were not significant [36,43]. Three simulation cases on bulk properties of granular material were investigated by Lommen and Schott et al. [44]. In this study, two cases were reported as not being affected by the reduced modulus, while the other exhibited unaccepted error. The effects of this speedup method for the DEM simulation of asphalt mixture compaction are not fully revealed. It is not clear what is the exact modulus reduction range without causing unacceptable calculation error in respect of material structure, contact coordinate number, and aggregate arrangement. In the gyratory compaction, the specimens are exposed to gyratory rotation at a small angle all while being constrained by constant pressure on the top. This loading condition ensures the particles are well-contacted, meaning there is hardly any collisions between the particles. Therefore, the timestep of the simulated compaction system can be calculated by using contact stiffness, external pressure, and particle radius. On the other hand, to form a stable and well-connected discrete system, the DEM models need to ensure the particle overlap ratio below a maximum value [45,46]. Furthermore, the timestep determination and particle overlap ratio can be correlated by utilizing the forcedisplacement law. Then, a theoretical inference that represents the material modulus reduction range without causing instability of DEM models can be derived. 2. Objectives This study aims to theoretically derive a material modulus reduction range that serves to accelerate the DEM simulation of asphalt mixture compaction. Then, employing laboratory tests and simulations to verify its ability of computational speedup while maintaining simulation accuracy. Moreover, the contributions of this study are not only the verification of material modulus reduction method, but also the development of asphalt mixture internal-structures indexes in respect of specimen’s height, rotation angles of coarse aggregates, and spatial distribution of mastic particles. 3. Theories of contact models 3.1. Determination of timestep To restrict the DEM models within a stable state, the timestep needs to exceed a critical value that is related to the minimum intrinsic frequency of the whole system. According to Bathe and Wilson [47], the notion of a point mass m, and spring with stiffness k in a one-dimensional mass-spring system (Fig. 1(a)) is determined by Eq. (1). The critical timestep t crit then can be derived by a second-order finite-difference scheme as given in Eq. (2).

kx ¼ m€x t crit ¼ pT ;

ð1Þ pffiffiffi T ¼ 2p mk

ð2Þ

where T is the minimum of a period of the system; x is coordination in one-dimensional direction.

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

3

Fig. 1. Spring systems at contacts: (a) single mass-spring system; (b) multiple mass-spring systems.

The discrete particles are considered as a series of point mass and springs to simplify the determination procedures as shown in Fig. 1(b). Two equivalent systems of the single point mass system are shown in the same figure. Then, the critical timestep is determined by Eq. (3):

rffiffiffiffiffiffi rffiffiffiffiffi m m ¼ 4k k

ð3Þ

t crit ¼ 2

where k is the stiffness of the springs in Fig. 1(b). Considering the translational motion and rotational motion in this system, the critical timestep for the multiple mass-spring systems is given by the minimum value. The actual timestep (t step ) used in the calculation is taken as a fraction of this estimated critical value. In this case, the fraction equals 0.8. Therefore, the actual timestep is:

ffiffiffi 8 pm > k > > ffiffiffiffiffiffi ffi q < m ¼ 0:8  min ktran > ffi > > qffiffiffiffiffi : I

tstep

ð4Þ

rot

where k and k are the translational and rotational stiffness, respectively; I is the moment of inertia of the particle. Therefore, according to Eq. (4) the timestep is mainly determined by the mass of the particle and the contact stiffness. The relation between contact stiffness and material modulus can be expressed by:

kn ¼

AE ; L

ks ¼

kn  k

ð5Þ

A ¼ pr 2  r¼

R1 ; 

L¼

ð6Þ

min ðR1 ; R2 Þ;

ball  facet

R1 þ R2 ; R1 ;

ball  ball

ball  ball

ball  facet

d¼

ð7Þ

ð8Þ

where, kn and ks are the normal and shear stiffness of the contact; A is the area of the contact plane (shadow area in Fig. 2); E* is the effective modulus; L is the distance between the two contact balls; R1 and R2 are the radius of the contact balls (Fig. 2).

F ; kn

F ¼ PA

ð9Þ

where F is the contact force; P is the instant pressure. Then the displacement can be expressed by substituting Eqs. (5–8) into Eq. (9):

d¼

krot

tran

Fig. 3(a). The overlap divided by the particle diameter is then defined as the particle overlap ratio. Cleary recommended that the maximum particle overlaps ratio should be smaller than 1.0% [46]. Large overlaps between particles would significantly influence the model’s structure, force distribution, and void ratio, see Fig. 3(b, c). For a simple contact of two balls with different radii, the displacement at the contact under an instant pressure can be given by:

P ðR1 þ R2 Þ E

ð10Þ

Assuming R2 is the smaller radius, then the maximum overlap ratio (Rov erlap ) over the particle diameter at the contact is:

Rov erlap ¼

  d P R1 1 ¼  þ 2R2 E 2R2 2

ð11Þ

Then the maximum overlap ratio can be predicted by substituting the worst condition into Eq. (11). In this study, the P was set to 600 kPa during static compaction and gyratory compaction. The largest grain size (R1 ) was 19 mm while the smallest (R2 ) was 1.75 mm. The E of aggregate was set to 5.5  1010 Pa. Substituting the above values into Eq. (11), the Rov erlap at the contact equals 0.0065%, which is much lower than the required 1% overlap ratio for a stable system. Therefore, by decreasing the E to 5.5  108 Pa (1/100 times of its original value), the Rov erlap equals 0.65% which is still in the acceptable range. On the other hand, according to Eq. (4), the t crit is 10 times larger than the original value. Theoretically, the simulation error of asphalt mixture flow DEM models is in the acceptable range by reducing the material modulus to 1/100 of its actual value; this procedure shall increase the computational efficiency by 10 times. 4. Methodology for DEM simulations 4.1. Specimens preparation

3.2. Effects of the modulus values on particle contact Adequate contact modulus ensures the overlaps between particles are small enough to form a good contact network as shown in

Three mixture designs were adopted in the laboratory tests and DEM simulations. The mixes #1, #2, and #3 (see Table 1 for details) was designed with increasing nominal maximum aggregate sizes

4

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

Fig. 2. Contact schema: (a) ball-ball contact; (b) ball-wall contact.

Fig. 3. Overlaps between particles: (a) a good contact network with adequate contact modulus; (b) a squeezed contact network with low contact modulus; (c) large overlaps observed between particles.

Table 1 Mixture design. Sieve Size (mm)

Mix#1

Mix#2

Mix#3

100 100 97 75 54 36 25 15 7 4.8 5.8

100 94 86 71 54 38 26 16 8 4.4 5.4

Passing (%) 19 12.5 9.5 4.75 2.36 1.18 0.595 0.297 0.149 0.074 Asphalt Content (%)

100 100 100 94 69 46 32 20 13 8.5 7

(NMAS) to evaluate the effects of grain size on computational efficiency. In the DEM simulation, coarse aggregates were represented by clump elements with realistic shapes obtained through a 3D scanner. The bubble pack method was used to create clump elements. Non-overlap particles were generated in a cylinder container

according to the mixture designs, see Fig. 4 [28]. The grain sizes of realistic aggregate were determined through virtual sieving analysis [30]. Fine aggregates (grain size 0–2.36 mm) and asphalt binder were combined as asphalt mastic and represented by ball elements with diameters ranging from 1.75 mm to 2.36 mm. The numbers of elements are listed in Table 2. The Burger’s model is commonly used to describe the viscous behavior of the compacted asphalt mixtures. However, the particle movement in the compaction process of asphalt mixtures is far more dramatic than the compacted specimens. The coarse aggregates and asphalt mastic are also undergoing large displacement and rotation at high speed during the compaction. Thus, to simplify the parameter variables, the linear contact model was assigned to the DEM models. The asphalt mastic was treated as a noncompressible fluid with damping units. The parameters used in the control group are listed in Table 3. The parameters of coarse aggregate were selected based on previous studies [48,49]. The asphalt mastic at high temperature was treated as noncompressible fluid. To ensure the stability of the discrete system, the modulus of asphalt mastic was set to the same of coarse aggregate, and its normal to shear stiffness ratio is set to 1 due to its flowability.

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

5

Fig. 4. The DEM models of asphalt mixture with increased nominal maximum aggregate size (NMAS).

plate was recorded as the specimen’s height under loose compaction state (Fig. 6(b)).

Table 2 Element number and particle number in the DEM models.

Total Element Number Total Particle Number Asphalt Mastic Coarse Aggregate

Mix#1

Mix#2

Mix#3

154,115 94,168 91,050 3118

142,260 70,926 67,904 3022

135,801 78,573 76,174 2399

Table 3 Model parameters used in the DEM model control group. Parameters

Coarse Aggregate

Asphalt Mastic

Particle Diameter Range (mm) Normal Contact Modulus (Pa) Normal to Shear Stiffness Ratio Friction Coefficient Density (kg/m3) Damp Ratio

2.36–19 5.5  109 2.4 0.3 2650 0.7

1.75–2.36 5.5  109 1 0.5 2140 0.7

4.2. Gravity fall: loose stacking state In the laboratory, the aggregates and asphalt were mixed and poured into a steel mold used for gyratory compaction. The diameter of the mold was 100 mm. The upper surface of the asphalt mixture was carefully dealt with a scoop to form a flat plane. This state of the asphalt mixtures was called loose stacking of laboratory specimens. The height of the empty mold and mold with asphalt mixtures were recoded, separately (Fig. 5). Five groups of material modulus (Table 4) were assigned to specimens in the DEM models. Then, gravity was applied to the models. The particles fell to the bottom of the cylinder container and formed a loose stacking state (Fig. 6(a)).

4.3. Static compaction: loose compaction state The mold containing asphalt mixture was assembled to a SUPERPAVE gyratory compactor (SGC) test machine. A compression pressure of 600 kPa was applied to the top of the specimen. Once the applied pressure was stable, the specimen’s height was recorded as the height of the loose compaction state. In the DEM simulation, a horizontal plate was generated at the upper surface of the specimen. The force-velocity servo mechanism was employed to apply a 600 kPa pressure to the specimen via the top plate. As the model reached equilibrium, the position of the top

4.4. Gyratory compaction: dense compaction state Once the applied pressure was stable, the mold was tilted to 1.16°. Then, the specimen was compacted under a gyratory movement. 75 gyrations were carried out. The specimen’s height after each gyration was recorded. The final compacted state of the specimen after 75 gyrations was defined as the dense compaction state. In the DEM simulation, a complex algorithm that controls the mold dynamic movement was utilized to control the gyratory compaction. The mold was tilted to 1.16°, meanwhile the top wall kept at 600 kPa pressure (Fig. 6(c)). After the model was brought to equilibrium, the mold was assigned to a spin velocity of 30 rounds/minute around the middle axis. The position of the top wall was monitored during the gyratory compaction and recorded as the specimen’s height of dense compaction state. More details of the gyratory compaction model can be found in previous studies [34,35]. 5. Results 5.1. Specimen’s height in laboratory and in simulation For each mixture design, four replicate samples were tested in the laboratory. The Specimen’s height under loose stacking, loose compaction, and dense compaction states are listed in Fig. 7(a). In general, the specimen’s height showed a decreasing trend with a compaction degree and NMAS increased. Mix#1 achieved the highest specimen height in the loose stacking and the loose compaction, while Mix#3 was the lowest. After gyratory compaction, the specimen height of Mix#2 reached the highest under the dense compaction state. The results also indicated the standard deviation of the loose stacking height was greater than the other two states, which was caused by human error and the randomness in particle packing. The Fig. 7(b) compares the specimen’s height decrease due to the compaction processes. The first group of columns describe the height decrease from loose stacking state to static compaction; while the second group of columns are the height decrease from static compaction to gyratory compaction. The mixtures with smaller NMAS showed larger variation, especially under the static compaction states. Three replicated specimens of each mixture design were created by changing random seed values in DEM simulation. The spec-

6

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

Fig. 5. The laboratory tests measurement and the compacted asphalt mixture specimens: (a) measuring the height of an asphalt mixture specimen under the loose stacking state; (b) compacted asphalt mixture specimens of Mix#1, Mix#2, and Mix#3.

Table 4 Material modulus setting of the tested specimens.

Reduction Ratio Contact Modulus (Pa)

#Control

Group #2

Group #3

Group #4

Group #5

1 5.5  1010

1/10 5.5  109

1/100 5.5  108

1/100 5.5  107

1/10000 5.5  106

Fig. 6. The DEM simulation for asphalt mixture compaction: (a) loose stacking state; (b) loose compaction state; (c) dense compaction state.

imen’s height with three compaction states and their percentage error of simulation compared to laboratory tests are shown in Fig. 8. Similar to the laboratory results, the specimen’s height trended to decrease as the compaction degree and the NMAS increased (Fig. 8(a)). Mix#1 had the highest specimen’s height for all compaction states due to its NMAS being the lowest (4.75 mm) among the three mixture types. The specimen’s height maintained the decreasing trend of specimen’s height as the NMAS increased. This is different from the laboratory results in which Mix#2 had a reverse increase in specimen’s height under the dense compaction state. The standard deviation of the simulation results was smaller compared with laboratory results. Fig. 8(b) shows the percentage error of the DEM simulation in comparison to laboratory tests. Negative values represent smaller results in the simulation, while

positive values represent the opposite. The results indicated the simulated results were smaller than the lab results, except the samples of Mix#1 under dense compaction state. The percentage error at loose stacking state was larger than that of loose compaction and dense compaction state. In summary, the percentage error of the DEM simulations compared to the laboratory tests were within 5.0%. The control group of the asphalt DEM compaction model had reasonable reliability in simulating the compaction process of asphalt mixtures. In comparison with the laboratory tests, the simulation error under the loose compaction and dense compaction states were below 4%. The maximum error under the loose stacking state was 5.8%. Besides human factors and material randomness, the error was caused by the application of a linear contact model

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

7

Fig. 7. Laboratory volumetric properties of the three mixtures: (a) specimen height of loose stacking, loose compaction, and dense compaction; (b) specimen height variation after 600 kPa static compaction and 75 gyrations.

Fig. 8. Volumetric properties of the DEM simulation and percentage error against the laboratory tests: (a) specimen height under loose stacking, loose compaction, and dense compaction; (b) percentage error of the simulated results compared with laboratory results.

in the simulation. In the linear contact model, no adhesive forces are applied to the particles. However, the adhesive forces play an essential role in formatting the asphalt mixture structure under loose stacking state. The adhesive forces between connected particles help to maintain a higher stacking height during compaction. Moreover, the compaction forces of static compaction and gyratory compaction were significantly more significant than the adhesive forces. Particles were forced to move and rearrange by the compaction forces. Thus, the simulation error was much smaller for the loose compaction and dense compaction states.

5.2. Effect of reduced material modulus on calculation efficiency The average calculation time of the DEM models increased significantly as the material modulus increased (Fig. 9). The results agreed with the interpretation of Eq. (4). It took 340–450 h (on a workstation with 8  4.20 GHz cores) to complete the simulation of the compaction process when the material modulus was set to 5.5  1010 Pa. As the material modulus dropped to 5.5  106 Pa, the calculation only lasted 3–5 h. In respect of mixture design,

Fig. 9. Calculation time of the DEM models as a function of the material modulus.

8

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

Mix#3 had the largest NMAS and required the longest calculation time, while Mix#1 needed the shortest. The difference in calculation time increased along with the increase of material modulus for the three mixture designs. However, the effects of mixture design on calculation time were less evident compared with the effects of changing the material modulus. The reduced material modulus was observed profoundly increased the calculation efficiency of asphalt mixture DEM compaction models. The calculation time of the control group that used the normal material modulus required 300–450 h to solve the compaction process on a workstation with 8  4.20 GHz cores. By reducing the material modulus, the calculation time was cut down to one day or even several hours. The simulation results agreed with the interpretation of Eq. (4). The limited calculation efficiency is one of the most significant barriers that block developing the simulation of asphalt mixture compaction so that the higher calculation efficiency has pronounced importance. Meanwhile, the effects of reduced modulus on the structures of the simulated mixtures are the biggest concerns. Therefore, the key is to find a balance between efficiency and accuracy. 5.3. Effect of reduced material modulus on internal structure 5.3.1. Specimen’s height Fig. 10(a) shows the percentage error of the specimen’s height of Mix#3 as a function of the material modulus setting. For specimens under loose stacking state, the errors in height trended to be decreasing linearly as the modulus decreased exponentially. For the other two compaction states, the percentage error of the specimen’s height was smaller than 0.5% at 5.5  108 Pa and 5.5  109 Pa. However, the percentage error increased dramatically when the modulus dropped down to 5.5  107 Pa. Fig. 10(b) indicates the results from the three mixture designs have similar decreasing trends. Mix#3 had the largest error, while Mix#1 had the fewest error. This indicated the effects of reduced material modulus increased as the NMAS increased. In summary, when the material modulus reduction was 1/100 times, the percentage error for three groups of specimens under compacted states was within 2% compared with the control group. The specimen’s height is the most direct index to evaluate the internal structure of the specimens. The reduced material modulus caused a larger overlap ratio at contacts between the particles and then reduced the specimen’s height. Exceed overlap ratio affects the model’s internal-structure, force distribution, and void. It was found that the simulation error in the specimen’s height was less

than 0.5% compared to the control group for all three mixture types, except the maximum 2.07% error exhibited by the specimens under loose stacking states. 5.3.2. Average coordination number The average coordination number (Nco ) is defined as the total contact numbers on a particle. In the simulation, the total contact numbers were recorded. In the DEM, contact is the connection of two elements; thus, the average coordination number can be calculated by Eq. (12).

Nco ¼

2  Nc Np

ð12Þ



where, Nco is the average coordination number of DEM models; Nc is the total number of contacts; Np is the total number of particles. 

The result of N co as a function of material modulus in Mix#3 is 

shown in Fig. 11(a). The N co for the control group was 2.86, 3.87, and 4.58 under the dense compaction, loose compaction, and loose stacking states, respectively. The Nco increased as the material modulus decreased for all three compaction states. The dense com

paction state received the highest N co , while the loose stacking state received the least. Moreover, under the loose stacking state, 

the N co decreased slower than the other two states. Under the 

two compacted states, especially the dense compaction, the N co was more sensitive to the reduced material modulus. According 

to Fig. 11(b), the percentage error of N co compared with the control group increased significantly as the material modulus decreased. With 1/10 times material reduction, there were 5.55%, 8.86%, and 16.08% error under the dense compaction, loose compaction, and loose stacking state, respectively. The reduced material modulus had larger impacts on the dense compaction states. The dense compaction states exhibited the smallest error (5.55%) when the material modulus reduction was within 1/100 times. However, it exhibited the largest error (131.69%) when the reduction reached 1/10000 times. The average coordinate number reflects the generally contact the state of particles in a DEM model. A DEM model needs the 

N co to be large enough to form a stable and well-connected force 

transmit structure. In general, the N co needs to be larger than 4. 

In the presented case, the N co values of the control group were 2.86, 3.87, and 4.58 under the loose stacking, loose compaction,

Fig. 10. Percentage error of specimen’s height as a function of material modulus: (a) Mix#3 under three compaction states; (b) three mixtures under loose compaction state; (c) three mixtures under dense compaction state.

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

9

Fig. 11. Mean coordination number of elements in Mix#3 DEM models under different compaction states as a function of material modulus.

and dense compaction states, respectively. The results agreed with 

the general requirement of N co . On the other hand, the reduced material modulus causes a larger overlap ratio and create more 

contacts. According to the simulation results, the N co increased sig

nificantly as the material modulus decreased. The N co under dense compaction state of Mix#3 increased 5.55% and 23.93% as the material modulus reduction reached 1/10 and 1/100 times, respectively. These results indicated the reduced material modulus had 

relatively larger impacts on the N co values than the other internal-structure indexes. However, the material modulus is not 



the only control factor of N co in the DEM values. The N co can be adjusted to the general values by altering the contact gaps of the contact model. Thus, as long as the particles’ arrangement keep

Euler angle in the xy plane (a). Pebble A and Pebble B are the components of a clump that represents a coarse aggregate particle. The relative position of Pebble A and Pebble B are fixed. Thus, the vector (v) that connects the center of Pebble A and Pebble B is defined as the orientation of this clump. Based on the orientation in 2D, the 3D Euler angles (Fig. 13(b)) can be defined similarly using the vector that connects the central points of the two pebbles. The Euler angles in 3D (a, b, and c) of coarse aggregates with grain sizes larger than 4.75 mm were calculated for all five material modulus groups. The average absolute difference in the Euler angles of xy plane was compared with the control group (Da ). The calculation was defined by Eqs. (13) and (14).



a ¼ arccos



consistent, the differences in the N co caused by reduced material modulus would not affect the internal structure of a DEM model. In the DEM models, a contact is generated once their contact gaps value is smaller than a setting value. By extending or narrowing the contact gaps, one can decrease or increase the existing contact 

number. Therefore, the N co can be adjusted, without changing the model structures, by altering the contact gaps of at the contacts.

5.3.3. Arrangement and rotation angles of coarse aggregates The coarse aggregates are the main components that form the load-carrying structure within compacted asphalt mixtures. Thus, the arrangement and orientation of coarse aggregates are the key factors to evaluate the internal-structure of compacted asphalt mixtures. Five DEM models with increased material modulus under dense compaction state are shown in Fig. 12. In general, the coarse aggregates maintained their relative position to each other through the simulations. Coarse aggregates with small grain sizes received larger displacement and rotation. The relative position of the two particles (2.36–4.75 mm, in the white circles in Fig. 12) had significant differences; however, their position change in respect of the whole model was not evident. The effects of material modulus were less on the coarse aggregates with larger grain sizes. For example, the relative position of the two particles (4.75–9.5 mm) inside the white square had slight differences, especially when the material modulus was larger than 5.5  107 Pa. The Euler angles of clump orientation were used to evaluate the effects of the reduced material modulus on the spatial rotation angles of coarse aggregates. Fig. 13(a) shows the definition of the

Da ¼ N1p

vb n jv jb n



ð13Þ

 PNp  ai  a0  i¼1

i

ð14Þ

where a is the 3D Euler angle in the xy plane; v is the vector of b is the positive unit vector in the x-axis; Da is clump orientation; n the average absolute difference of a; N p is the total number of particles involved in the calculation; ai is the ith value of the experi0 mental group; ai is the ith value of the control group. The difference in the xz and yz plane (Db and Dc ) was calculated with a similar method. The average percentage error of Euler angles was divided by 360°. The calculated results are listed in Table 5. The differences of the Euler angles, which indicates the spatial rotation of coarse aggregates, increased as the material modulus decreased. The average difference in xy plane was smaller than the xz plane and yz plane for the three compaction states. Moreover, the reduced material modulus enlarged the differences caused by compaction, and the models under dense compaction state exhibited the most substantial differences. In summary, the percentage error in Euler angles were less than 3.16% when the modulus reduction was within 1/100 times. The arrangement of coarse aggregates directly impacts the force-carrying structure in the asphalt mixture specimens. It can be observed from the figures that most of the coarse aggregates maintained their relative positions through the material modulus reduced. To quantify the differences in the arrangement of coarse aggregates, the spatial rotation based on Euler angles were employed. The spatial rotation in the xz plane and yz plane were observed 20% larger than the zy plane. This was caused by the

10

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

Fig. 12. Arrangement of coarse aggregates under dense compaction state of Mix#3.

Fig. 13. The orientation and Euler angles of a clump: (a) Clump orientation and the Euler angle in 2D; (b) the Euler angles coordination in 3D.

Table 5 Average difference (D) and percentage error in rotation angels compared to the control group. Specimen States

Euler Angle

Material Modulus 5.50E+06

Loose Stacking

Loose Compaction

Dense Compaction

xy plane xz plane yz plane xy plane xz plane yz plane xy plane xz plane yz plane

5.50E+07

D(°)

Percentage Error (%)

6.65 9.14 7.70 8.47 11.66 10.34 13.79 14.62 14.65

1.85 2.54 2.14 2.35 3.24 2.87 3.83 4.06 4.07

5.50E+08

D(°)

Percentage Error (%)

5.63 8.01 6.86 7.27 9.11 7.90 12.11 12.69 11.93

1.57 2.22 1.90 2.02 2.53 2.19 3.36 3.53 3.31

5.50E+09

D(°)

Percentage Error (%)

D(°)

Percentage Error (%)

5.04 7.55 6.28 7.20 9.45 7.90 9.90 11.34 11.38

1.40 2.10 1.74 2.00 2.63 2.19 2.75 3.15 3.16

4.58 6.61 5.59 6.53 8.65 7.03 8.12 8.26 8.18

1.27 1.84 1.55 1.81 2.40 1.95 2.26 2.30 2.27

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

11

Fig. 14. Calculation of distribution density of asphalt mastic particles in the vertical direction and radical direction.

coarse aggregates received more displacement and rotation during gyratory compaction in the vertical direction compared to the horizontal direction [50]. In summary, the percentage error of spatial rotation of coarse aggregates were less than 3.16% when the modulus reduction was within 1/100 times. 5.3.4. Spatial distribution of asphalt mastic particles The distribution density of asphalt mastic particles in the radical direction and vertical direction was used to evaluate the influence of the reduced material modulus on the spatial distribution of asphalt mastic particles. The distribution density of asphalt mastic particles is defined as the area occupation ratio, which is calculated as the area of the cross-section of a cutting plane in radical or vertical direction divided the cutting surface area of mastic particles, see Fig. 14. In the vertical direction, the projections of the crosssections are cycles of the same radius R (R is the radius of the specimen), while in the radical direction the projections of the cross-sections are cylinders of radius from 0 to R. The density functions in the radical direction (dðrÞ) and vertical direction (dðhÞ) are denoted by Eqs. (15) and (16). The average error (ed ) compared to the control group is calculated by Eq. (17)

dðrÞ ¼ 2SpmrrH

ð15Þ

dðhÞ ¼ pSmh R2

ð16Þ

ed ¼ 1n

Pn d1 ðiÞd0 ðiÞ i¼1  d0 ðiÞ 

ð17Þ

where r is the distance variable to the specimen center; Smr is the cutting surface area of mastic particles cut by a cylinder of radius r; H is the total height of the specimen; h is the height variable; Smr is the cutting surface area of mastic particles cut by a cycle of radius R at height h; d1 ðiÞ is the density at i of the experimental group; d0 ðiÞ is the density at i of the control group. As shown in Fig. 15(a), the density of asphalt mastic for all the models with different material modulus had similar trends. The density curves of the largest four groups were close to the control group. 5.5  106 Pa material modulus group had a higher density for all the measurements; however, its density curve shape kept the same trend with the control group. The density curves can be divided into three sections by the distance to the specimen center: S1 (0–20 mm), S2 (20–40 mm), S3(40–50 mm). The density of S1 was lower than the other two sections and increased as the distance increased. The density of S2 was stable at around 0.35. The density of S3 had a higher increase in speed as the distance increased. Sharp peaks were observed close to the specimen boundary. Thus, the simulation results indicate the center of the specimens had a smaller density of asphalt mastic in the radical direction. Caused by the differences in specimen’s height under dense compaction state, the asphalt mastic spatial distribution in the vertical direction had significant differences (Fig. 15(b)). The

Fig. 15. Spatial distribution of asphalt mastic in Mix#3: (a) density curves in the radical direction as a function of radius r; (b) density curves in the vertical direction as a function of height h.

12

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272

5.5  108 Pa and 5.5  109 Pa groups were close to the control group. However, the models with lower material modulus had greater differences from the control group. All specimens showed a U-shaped distribution that had higher density near the bottom and top surface. The average error (ed ) in the radical direction and vertical direction dramatically increased as the material modulus decreased. When the material modulus reduction was within 1/100 times, the average error was smaller than 4.08%. A few paths were available for the fine particles to travel during the compaction process in a container with confined boundaries. Therefore, the distribution density in the radical and vertical directions of specimens using the reduced material modulus had similar curve shapes compared with the control group. However, the smaller material modulus did lead to larger overlaps and shorter specimens, this caused the simulation error in the vertical direction to be larger than the radical direction. In general, the maximum simulation error in respect of the spatial distribution of asphalt mastic was smaller than 4.08% when the material modulus reduction was within 1/100 times. Moreover, the segregation of fine particles was observed in both horizontal and vertical directions, which agreed with the observation of the laboratory specimens [51].

6. Conclusions The conclusions are summarized as follows: (1) An inference was derived by combining the timestep determination equations, the force-displacement law, and small displacement assumption in respect of overlap ratio: exists a material modulus reduction range that serves to accelerate the DEM simulation of asphalt mixture compaction while maintaining its simulation accuracy. In the case of asphalt mixture compaction, the modulus reduction range was determined as 1/100 times; the increase of calculation efficiency was determined as ten times. (2) A three-step compaction process was carried out in the laboratory and in simulation to verify this inference. The simulation results were observed to have less than 5% error with respect to the specimen’s height compared with laboratory test results. The models exhibited a significant increase of calculation efficiency and agreed with the prediction of the inference. Depending on the mixture designs, the compaction models now can be solved by 35–47 h instead of 340–450 h on a workstation with 8  4.20 GHz cores when the material modulus was reduced by 1/100 times. (3) Four categories of internal-structure indexes were proposed to evaluate the effects of the reduced material modulus, namely the specimen’s height, average coordinate number, rotation angles of coarse aggregates, and spatial distribution of asphalt mastic particles. Compared with the control group, the models exhibited errors less than 4.08% in respect of the internal-structure indexes. The reduced material modulus had relatively larger impacts on the Nco compared with the other internal-structures indexes. However, as the particles’ arrangement maintain consistently, the average coordinate number can be adjusted through altering the contact gaps and without affecting the internal structure. Based on the proposed conclusions, two aspects of the research can be further developed. Firstly, the speedup models can be applied to other compaction processes. DEM models for the field compaction with a large number of elements can be established since the model calculation time can be significantly reduced by using the speedup method. Those compaction DEM models can

advance the understanding of the asphalt mixture compaction mechanism. Secondly, the asphalt mixture samples that compacted in the accelerated DEM models can be used to develop a series of performance tests. For example, the authors currently are working on the indirect tensile (IDT) test, disk-shaped compact tension (DCT) test, and rutting test. These gyratory compacted samples have more realistic aggregate structures compared with those randomly generated models and cost significantly less amount of resources than the conventional X-ray CT image-based models. CRediT authorship contribution statement Xiaodong Zhou: Conceptualization, Methodology, Software, Writing - original draft. Siyu Chen: Methodology, Resources, Data curation. Dongdong Ge: Resources, Data curation. Dongzhao Jin: Resources, Data curation. Zhanping You: Conceptualization, Methodology, Writing - review & editing, Project administration, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment We thank Julia Barnes for assistance with grammar checking that significantly improved the manuscript. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.conbuildmat.2020.118272. References [1] P. Liu et al., Primary investigation on the relationship between microstructural characteristics and the mechanical performance of asphalt mixtures with different compaction degrees, Constr. Build. Mater. 223 (2019) 784–793. [2] E. Masad et al., Internal structure characterization of asphalt concrete using image analysis, J. Comput. Civil Eng. 13 (2) (1999) 88–95. [3] L.S. Tashman et al., Internal Structure Analysis of Asphalt Mixes to Improve the Simulation of Superpave Gyratory Compaction to Field Conditions, Washington State University, 2000. [4] A.R. Coenen et al., Aggregate structure characterisation of asphalt mixtures using two-dimensional image analysis, Road Mater. Pavement Des. 13 (3) (2012) 433–454. [5] N.R. Sefidmazgi, L. Tashman, H. Bahia, Internal structure characterization of asphalt mixtures for rutting performance using imaging analysis, Road Mater. Pavement Des. 13 (sup1) (2012) 21–37. [6] H. Xu, W. Guo, Y. Tan, Internal structure evolution of asphalt mixtures during freeze–thaw cycles, Mater. Des. 86 (2015) 436–446. [7] P. Liu et al., Modelling and evaluation of aggregate morphology on asphalt compression behavior, Constr. Build. Mater. 133 (2017) 196–208. [8] G. Lu et al., Numerical analysis for the influence of saturation on the base course of permeable pavement with a novel polyurethane binder, Constr. Build. Mater. 240 (2020.) 117930. [9] P.A. Cundall, O.D. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65. [10] K.-N.G. Chang, J.N. Meegoda, Micromechanical simulation of hot mix asphalt, J. Eng. Mech. 123 (5) (1997) 495–503. [11] W. Buttlar, Z. You, Discrete element modeling of asphalt concrete: microfabric approach, Trans. Res. Rec.: J. Trans. Res. Board 1757 (2001) 111–118. [12] A.C. Collop, G.R. McDowell, Y. Lee, Use of the distinct element method to model the deformation behavior of an idealized asphalt mixture, Int. J. Pavement Eng. 5 (1) (2004) 1–7. [13] E. Mahmoud, E. Masad, S. Nazarian, Discrete element analysis of the influences of aggregate properties and internal structure on fracture in asphalt mixtures, J. Mater. Civ. Eng. 22 (1) (2009) 10–20. [14] T. Ma et al., Heterogeneity effect of mechanical property on creep behavior of asphalt mixture based on micromechanical modeling and virtual creep test, Mech. Mater. 104 (2017) 49–59.

X. Zhou et al. / Construction and Building Materials 244 (2020) 118272 [15] T. Ma et al., Simulation of wheel tracking test for asphalt mixture using discrete element modelling, Road Mater. Pavement Des. 19 (2) (2018) 367–384. [16] Z. You, W.G. Buttlar, Application of discrete element modeling techniques to predict the complex modulus of asphalt–aggregate hollow cylinders subjected to internal pressure, Transp. Res. Rec. 1929 (1) (2005) 218–226. [17] A. Abbas et al., Micromechanical modeling of the viscoelastic behavior of asphalt mixtures using the discrete-element method, Int. J. Geomech. 7 (2) (2007) 131–139. [18] Y. Liu, Q. Dai, Z. You, Viscoelastic model for discrete element simulation of asphalt mixtures, J. Eng. Mech. 135 (4) (2009) 324–333. [19] X. Li et al., Discrete element analysis of indirect tensile fatigue test of asphalt mixture, Appl. Sci. 9 (2) (2019). [20] H. Kim, M.P. Wagoner, W.G. Buttlar, Micromechanical fracture modeling of asphalt concrete using a single-edge notched beam test, Mater. Struct. 42 (5) (2008) 677. [21] Y. Liu, Z. You, Simulation of cyclic loading tests for asphalt mixtures using user defined models within discrete element method, in GeoCongress, Characterization, Monit., Modeling GeoSyst. 2008 (2008) 742–749. [22] X. Ding, T. Ma, X. Huang, Discrete-element contour-filling modeling method for micromechanical and macromechanical analysis of aggregate skeleton of asphalt mixture, J. Trans. Eng., Part B: Pavements 145 (1) (2018) 04018056. [23] S. Adhikari, Z. You, 3D discrete element models of the hollow cylindrical asphalt concrete specimens subject to the internal pressure, Int. J. Pavement Eng. 11 (5) (2010) 429–439. [24] Y. Peng, J.-X. Bao, Z. Wang, A comparison of two-dimensional and threedimensional micromechanical discrete element modeling of the splitting tests for asphalt mixtures, DEStech Trans. Eng. Technol. Res. (2016). [25] W. Liu et al., Investigation of motion of coarse aggregates in asphalt mixture based on virtual simulation of compaction test, Int. J. Pavement Eng. (2018) 1– 13. [26] Y. Liu, Z. You, Visualization and simulation of asphalt concrete with randomly generated three-dimensional models, J. Comput. Civil Eng. 23 (6) (2009) 340– 347. [27] S. Chen et al., Innovation of aggregate angularity characterization using gradient approach based upon the traditional and modified Sobel operation, Constr. Build. Mater. 120 (2016) 442–449. [28] Y. Liu et al., Discrete element modeling of realistic particle shapes in stonebased mixtures through MATLAB-based imaging process, Constr. Build. Mater. 143 (Supplement C) (2017) 169–178. [29] H.-C. Dan et al., Numerical simulation of an indirect tensile test for asphalt mixtures using discrete element method software, J. Mater. Civ. Eng. 30 (5) (2018) 04018067. [30] Y. Liu et al. Determining Aggregate Grain Size Using Discrete-Element Models of Sieve Analysis. 2019. 19(4): p. 04019014. [31] E. Ghafoori Roozbahany, M.N. Partl, A new test to study the flow of mixtures at early stages of compaction, Mater. Struct. 49 (9) (2016) 3547–3558. [32] E. Ghafoori Roozbahany, Modelling the flow behavior of asphalt under simulated compaction using discrete element, Mater. Des. (2017).

13

[33] J. Chen et al., Application of discrete element method to Superpave gyratory compaction, Road Mater. Pavement Des. 13 (3) (2012) 480–500. [34] F. Gong et al., Lab assessment and discrete element modeling of asphalt mixture during compaction with elongated and flat coarse aggregates, Constr. Build. Mater. 182 (2018) 573–579. [35] F. Gong et al., Using discrete element models to track movement of coarse aggregates during compaction of asphalt mixture, Constr. Build. Mater. 189 (2018) 338–351. [36] K.F. Malone, B.H. Xu, Determination of contact parameters for discrete element method simulations of granular systems, Particuology 6 (6) (2008) 521–528. [37] Y. Zhang et al., Predicting dynamic shear modulus of asphalt mastics using discretized-element simulation and reinforcement mechanisms, J. Mater. Civ. Eng. 31 (8) (2019) 04019163. [38] Y. Zhang et al., Prediction of dynamic shear modulus of fine aggregate matrix using discrete element method and modified Hirsch model, Mech. Mater. 138 (2019) 103148. [39] Y. Tsuji, T. Tanaka, T. Ishida, Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol. 71 (3) (1992) 239– 250. [40] X. Zhang, L. Vu-Quoc, Simulation of chute flow of soybeans using an improved tangential force–displacement model, Mech. Mater. 32 (2) (2000) 115–129. [41] H. Huang et al. Contact stiffness affecting discrete element modelling of unbound aggregate granular assemblies, 2008. 167–172. [42] S. Masson, J. Martinez, Effect of particle mechanical properties on silo flow and stresses from distinct element simulations, Powder Technol. 109 (1) (2000) 164–178. [43] J. Härtl, J.Y. Ooi, Experiments and simulations of direct shear tests: porosity, contact friction and bulk friction, Granular Matter 10 (4) (2008) 263. [44] S. Lommen, D. Schott, G. Lodewijks, DEM speedup: stiffness effects on behavior of bulk material, Particuology 12 (2014) 107–112. [45] P.W. Cleary DEM, prediction of industrial and geophysical particle flows Particuology 8 2 2010 106 118 [46] P.W. Cleary, DEM simulation of industrial particle flows: case studies of dragline excavators, mixing in tumblers and centrifugal mills, Powder Technol. 109 (1-3) (2000) 83–104. [47] K.-J. Bathe, E.L. Wilson, Numerical methods in finite element analysis, 1976. [48] J. Chen et al. DEM Simulation of Laboratory Compaction of Asphalt Mixtures Using Open Source Code. 2015. [49] X. Yang et al., Integrated experimental-numerical approach for estimating asphalt mixture induction healing level through discrete element modeling of a single-edge notched beam test, J. Mater. Civ. Eng. 27 (9) (2014) 04014259. [50] X. Wang et al., Characterization of particle movement in Superpave gyratory compactor at meso-scale using SmartRock sensors, Constr. Build. Mater. 175 (2018) 206–214. [51] E. Masad et al., Quantifying laboratory compaction effects on the internal structure of asphalt concrete, Trans. Res. Record: J. Trans. Res. Board 1681 (1999) 179–185.