Micromechanical simulation of porous asphalt mixture compaction using discrete element method (DEM)

The paper aims to simulate the micromechanical behavior of asphalt mixtures during the compaction process using the Discrete Element Method (DEM). The interactions between the components of a Porous Asphalt (PA) mixture were represented using an Elastic Viscoelastic Contact Model (EVCM), which is a user-defined model implemented in EDEM software, developed based on linear elastic and Burger’s viscoelastic constitutive equations. The macroscale parameters of asphalt mortar were characterized using the nonlinear regression analysis of master curves obtained from Dynamic Shear Rheometer (DSR) tests. The verification process of EVCM successfully indicated that the computations trends fall within the range of expected values for the typical asphalt mixture material. Further, a Superpave Gyratory Compaction (SGC) test was carried out and the obtained sample was scanned using X-ray Computed Tomography (X-ray CT) to capture the air void distributions. The DEM was utilized where digital samples were established to simulate the overall process of laboratory and field compaction. The simulation results showed that the model provided a comparable prediction of responses and demonstrated the capability of SGC to fabricate a representative sample. The influence of temperature on the asphalt compaction process was explored and the results implied that temperature decreasing adversely affects the compactability and dramatically increases the demanded compaction efforts which are consistent with the law of viscoelasticity. On the contrary, when the temperature is high, the asphalt binder becomes too fluid and roller loads will simply displace, or “shove” the mat rather than compact it. Tracking the change in the air voids proportion indicates that the motion of aggregates is rather compound. The aggregates flowed vertically downwards in line with the compacting orientation while moved horizontally outwards away from the center. All in all, the findings confirm that the concept is technically practicable, affording the model great potential to help researchers understand the microstructural phases of asphalt mixture during the compaction.


Introduction
Over the past few years, clients in the pavement sector have become more aware of the importance of the accuracy and quality of the production and construction phases of projects. Production of the asphalt mix is controlled and can be checked on the plant, the temperature of the mix can be checked upon arrival. Both are important but can be done prior to placing the mix. Assuming that this is taken care of, even a good mixture can lead to the bad pavement if the compaction effort is not homogeneous or incorrect for the material or working conditions. In practical terms, compaction is an essential factor in determining the final quality of the pavement [1]. Insufficient compaction results in lower pavement strength and higher air voids into which moisture can penetrate [2,3]. In contrast, over-compaction does not attain adequate voids needed accommodate the asphalt binder that may bleed as the temperature rises.
Compaction quality control requires a thorough understanding of this mechanism that is not easy to gain visually. Several studies have been conducted using experiential methods to ensure compaction quality, yet many of the existing tools are still empirical [4]. Over the years, laboratory compaction methods have been developed to fabricate samples for performance evaluation. However, current laboratory tests for asphalt mixture compaction are time-consuming, labor-intensive and incapable to provide structural information about the mixture.
Moreover, the asphalt mixture is compacted in a close and harsh environment where the temperature is relatively high. Thanks to advances in computer technology, sophisticated analyses can be performed effectively enabling researchers to develop more realistic techniques for studying and evaluating the compaction process.
Many prior researchers tried to use Finite Element Methods (FEM) to simulate the compaction. For instance, Koneru et al. [5] developed a constitutive theory within a thermodynamic setting to study the compaction of asphalt mixtures where ABAQUS commercial FEM software was used to simulate the SGC. Zheng et al. [6] built a mechanical model to simulate the pavement vibratory compacting process utilizing the FEM ANSYS software. Huerne [1] used the FEM with code DiekA to simulate the compaction process of asphalt mixture under roller compaction conditions. Besides, researchers [7] have also employed computer-aided pavement analysis CAPA-3D finite element software [8] to explore factors and processes that influence field compaction.
Although modeling the compaction process based on continuum media FEM has been well researched, simulating the heterogeneity of the asphalt mixture remains problematic. Recently, the discrete element method (DEM) has been increasingly utilized for the microscopic understanding of macroscopic particulate material behavior [9]. With its inherent advantages, it can potentially overcome the shortcomings of FEM to simulate the structures consisting of particles. After DEM was first proposed by Cundall [10] for the study of rock mechanics and later expanded for modeling granular materials [11], DEM grew into the most popular numerical method to simulate the particulate system. These merits have allowed researchers to simulate some aspects of asphalt mixture behavior that cannot be seen or monitored in experiments. It has been reckoned that DEM could also be a promising and competent approach to modeling the micromechanical behavior of asphalt mixture [12]. DEM was successfully applied to predict creep strains of a laboratory sample of asphalt mixture subjected to diametric loads in the Superpave indirect tension test (IDT). The study indicated that micromechanical modeling has tremendous potential benefits in the field of asphalt technology for reducing or eliminating test costs [13]. Similarly, DEM was used to simulate a highly idealized asphalt mixture under uniaxial and triaxle compressive creep loading. The predicted results had a reasonable agreement with experimental data [14,15]. A clustered model was exploited to predict the complex modulus by using a 2D DEM simulation. This approach presented a low modulus compared to the experimental tests for some fine mixes due to insufficient aggregate--aggregate contact with the 2-D model [16,17].
Researchers [18] found that DEM can consistently describe the compaction phenomena with field observations and empirical experience. Past research studies [19,20] demonstrated that DEM models could simulate the compaction of asphalt mixture pretty well and concluded that simulation results match the experimental data and previous research results. Olsson [21] developed a new DEM modeling approach for studying the asphalt compaction process, incorporating contact and damage laws based on granular mechanics. The results showed that the proposed DEM is capable of capturing qualitatively and quantitatively responses in both cases, in addition to providing predictions of aggregate damage.
The literature review revealed that DEM offers a powerful tool for simulating the asphalt mixture; this is due to the discrete inhomogeneous nature of these substances. While a few models have been introduced, some limitations still persist. For instance, the influence of the ambient conditions on microstructural interaction, which in turn affects the macroscopic performance, has not been fully understood yet. On the other hand, determining the efficiency of a laboratory compaction device to accurately simulate the real field construction conditions still remains a challenge for pavement researchers and engineers.
This paper proposed a technique to develop a multi-physics model using DEM to simulate the compaction of a porous asphalt mixture in which the interaction between its components, the influence of temperature and the contact law were considered. To achieve this goal, a new bonded-contact model, hereafter referred to as an Elastic Viscoelastic Contact Model (EVCM), has been developed based on Burger's and elastic constitutive equations [20,[22][23][24]. EVCM is a materialdependent model as it can identify the components of an asphalt mixture and establish appropriate contact forces that suit their physical characteristics. The microscale parameters of asphalt mortar were characterized using the nonlinear regression analysis of master curves obtained from Dynamic Shear Rheometer (DSR) tests.
It is noteworthy that EVCM can be perceived as a straightforward tool for simulating the mechanical performance of asphalt mixtures, helping to understand their behavior and gain in-depth knowledge about their microstructure during the compaction process.

Research methodology
In order to compare laboratory and field compacted porous asphalt samples, SGC laboratory testing was carried out. The fabricated sample was then scanned using X-ray Computed Tomography (X-ray CT) to capture the air void distribution. DEM was then utilized to simulate the overall process of the SGC laboratory testing to assess the air voids distribution in the established digital sample. The developed model was further employed to simulate the behavior of asphalt mixture subjected to compaction in the field. The model's capability was checked by comparing the model's output against laboratory test data.

Contact model
In DEM, materials are modeled as particles and the interaction between these particles is governed by a contact law. Different materials can be distinguished by applying different force-displacement laws, namely, contact models. Contact models are physics laws that describe the behavior of the particles when they interact with each other and with equipment and effectively mimic the diverse behavior of particulate materials [25,26]. One major concern when using DEM is to suggest an appropriate constitutional contact law that includes microscopic parameters needed to represent the particulates on a macroscopic scale. The following sections discuss the development of an elastic and viscoelastic contact.

Theoretical principles
The ability of DEM to represent cementitious materials depends on the implementation of an inter-particle bonded-contact model [27]. One popular method for forming particle assemblies is by using a bonded contact. The bond has the ability to: (i) hold the adjacent particles together, (ii) prescribe their interactions, and (iii) transmit forces and moments.

The geometry of the bond
The force-displacement law for a bond between two particles requires the determination of the bond's stiffness. When particles are spherical, the bond is typically suggested to be a cylinder [28]. In this paper, the centers of the adjacent particles are connected by a cylindrical beam with a uniform cross-section. The beam length L b is the distance of the straight line between the centers; whereas, the beam radius R b is half the beam length. Fig. 1 shows the proposed bond.

The geometry of materials
The materials in DEM are typically represented by a set of particles that vary in shapes, sizes and properties. The maximum number of particles is limited by computational power. Due to this fact, it is common to simulate aggregates bigger than their actual sizes and to assume fine aggregates and asphalt as mortar which in turn is considered by the contact forces forces [19]. In this paper, all fine aggregates smaller than 2.0 mm, mineral fillers and asphalt binder were assumed as mortar. Consequently, the volume of the coarse aggregates is increased by multiplying each individual radius by increasing ratio λ to reflect the proportion of mortar volume. The aggregates in asphalt mixtures often have irregular shapes; this means that defining the shape of the DEM elements is not explicit and some approximations are required. In order to form more advanced shapes with spheres, particles could be formed from many spheres using a clustering technique [29]. Besides, EDEM allows introducing shape into a simulation by adding multiple surface spheres to a particle. For each surface, a user can change the radius of the sphere and also adjust its position so that shape of the particle (or cluster) will be changed [25]. In this paper, a cluster of particles has been shaped by joining composed of spherical surfaces and holding them together with relatively strong bonds.

EVCM process
The interaction between the constituents of an asphalt mixture could be seen through the role of each individual and how it contributes to determining the overall performance. As the aggregates and mortars have different physical and mechanical properties, the contact between them must have different tendencies that interpret those in the mixtures. Once the particles come into contact, EVCM is cable to identify each material and create corresponding elastic or viscoelastic contact forces antiparticles. Simple algorithms can be executed to detect the contact between spheres where the contact emerged only if the variable straightline distance between particle centers less than the sum of the radii of the two mortar-covered particles; i.e: Where: X A and X B are the positions of particles (A and B) in the global coordinate system, R [A] and R [B] are the physical radii of the two particles and λ is the ratio used to increase the actual aggregate sizes to reflect the volume of mortar. The difference between these two limitations is the overlap δ (See equation (2) and Fig. 2).
In this sense, the bonds act in a viscoelastic manner when the contact is between asphalt mortar, i.e.: A special case of the viscoelastic behavior occurs when the contact is between mortar and aggregates, that is, if one of the following two formulas is true: (4) EVLM will behave in an elastic manner when the contact is between aggregates, i.e.: EVCM's starting point is controlled by a nominated timestep T b0 , at which the bond is initiated. The flow chart in Fig. 3 summarizes the whole process of EVCM.

Burger's macroscale parameters
Asphalt mixtures are rheological and thermal materials that typically behave in a viscoelastic manner. It has been well documented that DEM combined with Burger's model could be utilized as a reliable tool to capture the rheological behavior of asphalt mixture [22]. Thereafter, Burger's model has been widely applied in DEM-based simulations to analyze the viscoelastic response of asphalt mixture, where the predictions compared favorably with the measurements [20,23,24]. In this paper, Burger's constitutive law was implemented to describe micromechanical behavior and to consider the time-temperature-dependent property of asphalt mixture. Dynamic Shear Rheometer (DSR) tests were carried out to determine the microscale parameters of asphalt mortar. The compositions of the mortars consisted of asphalt binder of 70/100 penetration grade, fine fillers (10-75 µm) type lime filler (Wigras 60) and sand (fractions less than 0.5 mm). The mortar composition is given in Table 1.
The tests were conducted under different conditions of five temperatures (− 10.0, 4.0, 21.0, 37.0 and 54.0) • C. For each one, different frequencies ranging from (0.001 to 10000.0) Hz were applied. Consequently, five different curves, represent the temperature-frequency relations, were produced. It has been demonstrated that viscoelastic data collected at one temperature can be superimposed upon data obtained at a different temperature simply by shifting one of the curves along the frequency (or time) axis [30]. Based on this principle, which is called Time-Temperature Superposition (TTS), the master curve was generated by setting a reference curve at 21.0 • C, then shifting the other data with respect to frequency using appropriate shift factors. The conventional method comprehended in the MEPDG [31] was applied in fitting the master curves of dynamic shear modulus |G*| and phase angle (φ). Fig. 4 shows the resulting curves.
For the purposes of inferring parameters of Burger's model at high temperatures, the curve at the reference temperature of 21.0 • C was shifted into temperatures (150.0, 125.0 and 100.0) • C. The resulting master curves for dynamic shear modulus and phase angle for the mortar  Typically, Burger's parameters can be obtained through nonlinear regression analysis of master curves [20,23,24,32]. In this paper, the four parameters of Maxwell and Kelvin models (E 1 , η 1 , E 2 and η 2 ), that form Burger's model, were determined by fitting the experimental test data and applying the complex compliance and phase angle equations [24]. Table 2 summarizes the results that correspond to temperatures (150.0, 125.0 and 100.0) • C.

The elastic contact
The bonds are subjected to internal forces caused by the relative displacements and rotations of the particles in contact. Mathematically, elastic contact forces are calculated based on Hooke's law, where the increment of the contact force linearly connects with the increment of the displacement. In each timestep, contact normal and shear forces of the elastic contact model are incrementally calculated using equations (7) and (8): [33] where: K n , K s are the normal and shear stiffness, δ U → n , δ U → s are normal and shear displacement. The stiffness of a cylindrical bond in the normal and shear directions are calculated by applying Equations (9) and (10): Viscoelastic Contact

Calculation of contact forces and moment
Summation of internal and external forces acting on each particle External forces Integration of equations of motion where the equivalent Young's modulus: E [eq] and equivalent shear modulus G [eq] are defined using Equations (11) and (12): [35,36] Table 4. Similarly, internal moments of the bond can be calculated incrementally using Equations (13) and (14): [33] δM In which, δM → n , δM → s are the incremental torsional moment, and flexural (bending) moment, respectively; δ θ → n , δ θ → s are the increment of the normal and shear rotational displacement, while K torq , K bend are torsional and bending stiffness of the bond respectively, and they are defined using Equations (15) and (16): Where: I, J: are area and polar moments of inertia for the cylindrical bond.

The viscoelastic contact
The same bond technique employed to simulate the elastic micromechanical interaction between the aggregate was applied to simulate the viscoelastic micromechanical interaction between mortars or between mortar and aggregates. Burger's model microscale parameters represent the stiffness contribution of the force-displacement law. Equations (17) to (24) express the correlations between micro-macro parameters of Burger's model when the contact is between the asphalt mortars: [24] Where: ν is Poisson's ratio for asphalt mortar, which in this paper was assumed to be: 0.35 (see Table 4). K mbend , K mtorq , K kbend , K ktorq , C mbend , C mtorq , C kbend and C ktorq are torsional and bending stiffness of Maxwell and Kelvin sections respectively.
When the contact, however, is between asphalt mortar and aggregates, i.e. one of the formulas (4) or (5) is true, K mn , K ms , K mbend and K mtorq are calculated using equations (25) and (26) instead of equations (17) and (21): [24] Where: ν is Poisson's ratio for asphalt mortar, E [eq] is the equivalent Young's modulus of the mortar and aggregate that are being in contact and can be calculated based on equation (11). The internal moments of the bond in the viscoelastic fashion can be calculated incrementally using Equations (27) and (28): [33] δM In which, δM → Bn , δM → Bs are the incremental torsional moment, and flexural (bending) moment, respectively; while K Btorq , K Bbend are torsional and bending stiffness of the viscoelastic bond respectively, and they are defined using the Equations (29) and (30):

The strength
The bond can resist the compressive, tensile and shear forces, up to the maximum strength, at which point the bond breaks. In this paper, the viscoelastic bond is broken when the normal (compressive, tensile) stress σ max or shear stress τ max exceed predefined values of compressive σ c , tensile σ t and shear strengths τ. The stresses value can be determined for each timestep using Equations (31) and (32): [33] The strength of asphalt mortar has been investigated and the obtained results demonstrate that the strength is exponentially related to the temperature [37]. The curve plotted in Fig. 7 shows the relationship between temperature and the logarithm of mortar strength:

Slip model
The slip model is necessary to resist the sliding of contacting particles. It serves to control the relationship between shear forces f t and normal forces f n . The model parameter for the slip model involves the coefficient of static friction μ. In this paper, the shear friction force follows Coulomb's law of friction model (f t ≤ μ s f n ). Besides, μ s equals 0.7, 0.3 for contacts between particles or particle and equipment respectively; and the coefficient of rolling friction is 0.01.

Damping forces
Realistically, both normal and shear forces contain damping components that (i) reflect the kinetic energy dissipated through microscopic actions such as the wave scattering and internal friction, and (ii) bring the simulation, if necessary, to a quasi-static state [33]. In this paper, in addition to viscous local damping acting at each contact, a bond damping has been proposed to consider the energy dissipation as where: β b is the bond damping coefficient that varies in the range (0-1) to determine the energy dissipation rate. (in this paper, β b = 0.1). m eq is an equivalent mass of particles A, B and can be calculated by applying equation (37) [39]:

Evaluating EVCM
Computationally, DEM simulation is a time-consuming process. For this reason, newly developed or modified DEM models are often validated using a simple and straightforward model with a small number of particles in a simple analytical problem. Practically, the simplified model can be relied upon to study the trends of the process and assure the reliability of the DEM model. In the simplified model, external forces were applied via two plates pushing horizontally against the two particles, both with a radius of 12.75 mm, passing through three stages: 1. The two particles were pressed together by plates that moved towards each other at a velocity of 1.0 mm/sec until the distance between them was reduced to 8.0 mm. During this stage, the plates push the particles together, causing them to overlap; 2. Then, they stopped moving for 1.0 sec (between 4.0 sec and 5.0 sec); 3. Finally, they returned in the opposite direction away of the particles at the same velocity 1.0 mm/sec, allowing the applied external forces to gradually decrease to the point the contacts between the particles and the plates disappeared.
To explore the behavior of EVCM with respect to materials in contact, three different scenarios have been simulated using the simplified model. Firstly, the entire material of the two particles was assumed mortar, meaning that the calculations will be carried out based on Burger's model. The resulting normal viscoelastic contact forces and the corresponding normal overlap recorded during the aforementioned three stages are plotted in Fig. 8.
Secondly, the entire material of the two particles was assumed aggregates (λ = 1.0), meaning that the calculations are carried out based on the elastic model. The resulting normal elastic contact forces and the corresponding normal overlap recorded during the aforementioned three stages are plotted in Fig. 9.
Finally, the materials of both particles were assumed aggregates covered with mortar including asphalt binder, filler and fine aggregates up to 2.0 mm. The volume is determined by λ, which is in this scenario corresponds to 1,085 (Fig. 2). In this case, bonding contact will initially perform in a viscoelastic manner until the point where the aggregates come into interaction and the contact starts acting with elastic behavior. Fig. 10 reflects the entire approach.
€The effects of the velocity of moving plates were also examined. The curves in Fig. 11 highlight the relationship between normal forces versus the normal overlap at different velocities (0.5, 1.0 and 2.0 mm/sec). The EVCM is obviously sensitive to the loading rate.
Implementing EVCM numerical analyses of the simplified model proved that the fundamental mathematical basis in the programmed code matches viscoelastic and elastic theories. The trends of DEM computations have been indeed found to fall within the range of values expected of typical asphalt mixture materials.

Simulating compaction process
Several factors can effect the overall compaction process, including temperature and thickness of the laid materials, asphalt binder content and type and grading of the aggregates in the mixture [40]. When attempting to simulate the compaction, particular attention should be paid to the roles of these factors. In the following texts, the simulation of the compaction process using DEM and analysis of the internal structure of the compacted mixture will be discussed.

Material specifications
Material specifications are usually determined through laboratory experiments. Among several laboratory compaction methods SGC, has popularity throughout asphalt production [41,42]. For instance, the ability of different compaction devices to simulate field compaction was evaluated and the results showed that SGC has an advantage in producing compacted samples with mechanical and engineering properties (resilient moduli, strengths, creep, stiffness, percent air voids, aggregate particle orientation, etc.) closer to those in the field [40,43].
To obtain the experimental data, an SGC test was conducted to  compact a PA mixture which was prepared according to the proportion mentioned in Table 3. The temperature during the mixing and compaction was 150 • C; and at the end of the compaction, the obtained sample was 55.0 cm in height and 150 cm in diameter. The vertical pressure was set at 600 kPa and the angle of gyration was set at 0.82 • ; while the gyration was applied at a rate of 30 rpm.

Formulating the digital samples
Based on the same distribution of coarse aggregate listed in Table 3, digital samples were generated to simulate the SGC test and field compaction. Particles were randomly generated and accommodated in a specific space within geometry elements, which are a cylindrical section represents the SGC mold and a cubic box represents the pavement domain. The geometry elements were used in the creation of boundary conditions for the particulate media. They are undeformable and the user controls their motion by setting prescribed velocity or acceleration. Parameters describing material characteristics were given based on the assumption that the geometry is made of steel. However, the bottom surface of the box was assumed to represent a base course of compacted crushed stones. The Hertz-Mindlin contact Model was employed to simulate the contact between the particles and geometry. Table 4 highlights the characterizations of the materials as well as the input parameters of the digital samples.
All fine aggregates (less than 2.0 mm), mineral fillers, and asphalt binder were assumed to be mortar and were represented by the viscoelastic contact forces. The volume fraction of the mortar phase has been counted by increasing the radius of each individual particle by 8.736% (λ = 1.08736) (see Fig. 2). Based on this assumption, Table 5 demonstrates the new particle volumes for both SGC and the field simulations. In this study, particle radii were randomly sized within a set size range that corresponds to the aggregate gradation. Particles were produced and introduced into the model by an efficient engine named "Particle Factory™" [25]. Particle factories are used to define where, when and how particles appear in a simulation. A static generation was set to produce the particles, which means the computational simulation is paused during the process. The particles were arranged in a random loosely packing where packing of spheres continued under gravity force until it is stable. Fig. 12 shows the particles arrangement for the digital samples.

SGC simulation
During the formulation of the digital sample, all laboratory conditions were asserted. More specifically, a loose asphalt mixture, with an initial height of 61.23 cm, was generated in a cylindrical compaction mold with a diameter of 150 cm. The load applied to the SGC digital sample consisted of a vertical pressure of 600 kPa, a gyration angle of 0.82 • , and a gyration rate of 30 rpm. Constant pressure was applied to particles by a moving steel plate using an EDEM coupling interface and stopped when the sample height was 55.0 cm.
The heights as a function of the number of gyrations for both the laboratory and DEM samples during compaction are shown in Fig. 13. In the simulation, this height was obtained by recording the centroid vertical position of the moving plate, in the test the set-up outputs the height at that position as well. Evidently, the EVCM model gives a reasonable agreement with the experimental data, showing the model's ability to predict the micromechanical behavior of an asphalt mixture during compaction.

Field compaction simulation
In the field, PA is typically compacted by a static-steel-wheel roller with several passes over the surface [37,44,45]. To assess the efficiency of the developed DEM model to simulate field compaction, a cylindrical slice of a steel drum with a width of 20 cm was set to simulate a 5-ton tandem roller with two drums each with a diameter of 122 cm and a width of 107 cm. During compaction, the roller exerts static forces on an asphalt mat laid inside a cubic box of dimensions 60.0 cm (length) × 40.0 cm (width) × 9.0 cm (thickness). The width and thickness of the sample were determined to minimize the influence of the boundary conditions. The center of the bottom face of the simulation box was taken as the center of the coordinate system (0,0,0). Every single pass, the roller traveled 60 cm distance at a constant speed of 1 m/sec (from  of the geometry. The center of rotation for the roller is equal to the center of mass and it moves as the position is updated by the coupling interface process. Compaction efficiency has been correlated with reducing air voids to quantify the compactability of asphalt mixtures [46]. In order to weigh the change in the air void content of the field DEM model during the compaction process, the model was divided into four sectors which are from left to right LL, CL, RL and RR. Each sector was split vertically into 9 slices. Only the central sectors (CL and CR) were compacted. The air void content of the compacted asphalt mixture in each slice was computed at the end of the pass. The chart in Fig. 15 shows the change in air void content across the nine slices of LL, CL, CR and RR sectors versus the number of passes. Based on these results, one can observe that the lowest slice at the bottom (slice 1) was not significantly affected by the compaction efforts whereas the compacted sectors (CL and CR) in the upper slice 9 disappeared because the compaction forces had pushed all aggregates towards the underlying slices. The decrease in air void content in CL and CR sectors of slice 8 during the first two passes attributes to the vertical movement of the aggregates from the upper slice 9. However, as the number of passes increases, a portion of the slice 8 aggregates continues to move towards the underlying slices causing an increase in the air void content.
The curves depicted in Fig. 16 record the change in the content of air voids over the height of the compacted sectors of the field model before starting and at the end of each pass.

Heterogeneity of compacted samples
Image processing techniques have been widely relied on to reconstruct the microstructure, identify the internal structure and capture the air void distribution of asphalt mixtures [47][48][49][50].
Moreover, the image processing results can be seen as physical benchmarks by which DEM results are interpreted and validated. To verify the reliability of DEM simulations, the SGC laboratory sample was scanned using CT X-ray. Image processing and segmentation techniques were implemented using Simpleware ScanIP version N-2018-3 commercial software to characterize the internal structure and quantify the air void distribution of the sample (Fig. 17). The process included the analyses of horizontal cascades of 106 X-ray images along the vertical axis with the aim of exploring the distribution of air void, which averaged 21.6%.
Similar to the procedure explained in the previous paragraph, the DEM sample was split into 30 slides perpendicular to the gyration axis and the air void distribution for each slide was computed and compared with those obtained from X-ray imaging processing. The average of the air void in the DEM-based SGC sample was about 24.4%, while it was 27.9% in the DEM-based field compacted asphalt mat. Fig. 18 illustrates a comparison between the air voids vertical distribution identified by DEM-based simulations of SGC and field compaction as well as by image processing of the SGC laboratory sample. The results show that the initial packing arrangement, initial void content and change in air voids distribution during the compaction are comparable.
Although the DEM simulations overestimated the air void content near the top, the results showed analogous distributions. This shows that both laboratory and field compactions apply comparable efforts to the asphalt mixtures where the prior is measured in terms of the number of gyrations and the latter in terms of the number of passes. The curves plotted in Fig. 19 show the relationships between decreasing air void content, in lab-based SGC, DEM-based SGC and DEM-based field samples, and the corresponding increase in compaction efforts.

Effect of temperature
As it is well known that the behavior of asphalt mixtures is temperature dependent [51,52]. Therefore under certain conditions changes in temperatures may have a greater influence on mechanical responses than changes in loading magnitude [53]. Thus, to better simulate compaction behavior, the temperature-dependency of the mixture should be reflected in the process. When the temperature of an asphalt mix drops during compaction, the ease of the process is adversely affected, and the compaction efforts lose their leverage. Typically, compactability is defined as the ease with which an asphalt mixture is compacted, i.e. reducing the air void content. There is a straightforward relationship between the mixture temperature and compactability. The compactability in this research has been represented by the compaction efforts exerted to decrease the proportion of the air void to a certain value. For this purpose, laboratory and field compaction modeling was conducted in order to explore the relationship between air void reduction and the corresponding compaction efforts. It has been argued that mechanical properties may decline by up to 30%, if an asphalt mixture is compacted outside the optimum temperature window [54]. Open-source software (PaveCool) has been developed to help users making quick decisions regarding coolweather paving [55]. The change in temperature of an asphalt mixture, as a function of time, can be predicted using this software. In this paper, PaveCool was implemented to determine the cooling rate of the asphalt mixture over time as well as the optimal compaction window in which the compaction process should be completed. (see Fig. 20).
To better understand the effect of temperature on compactability, three different temperatures (150, 125, and 100) • C were selected to   Table 2.

DEM-based SGC condition
As expected, simulation results imply that decreasing the temperature has an adverse consequence on the ease of compaction. Fig. 21 shows the relationships between air void content and the number of gyrations at different temperatures.

DEM-based field compaction condition
During the field compaction simulation, the materials were compacted at a constant temperature of 150 • C across all passes. It is obvious that it does not reflect the real field conditions, as the mixture temperature decreases over time. To gain a deeper insight into the effect of temperature on the field compaction, a different temperature was proposed during the second pass. This means two cases have been studied. In Case 1, the temperature of the mixture was at 150 • C during the first and second passes. While in Case 2, the temperatures during the first and second passes were set at 150 • C and 100 • C respectively. The 3D displacements of a particle, chosen at the top near the center of the asphalt mat, were tracked. Table 6 presents the coordinates of the particle center before starting compaction as well as at the end of each pass.
The curves in Fig. 22 show the vertical motion of the particle as the roller approaches the particle's position and moves away from it.
During the first pass, when the roller was close enough, the chosen particle moved up slightly before it turned sharply downwards. When the compaction roller was directly over the particle, the maximum depth was reached. Subsequently, the roller crossed the particle causing it to move upwards and settle at a greater depth than it was before the pass. In the second pass, the roller was moved in the opposite direction. As aforementioned, two different temperature values had been chosen to simulate the second pass. In case 1, the response showed approximately the same tendency as it did on the first pass. Initially, the particle raised slightly before dropping sharply to a greater depth than it did during the first pass. Likewise, while the roller was moving away, the particle was gradually pushed up. In case 2, maximum and final particle depths were less than they had been in case 1. However, the curve, which is relatively flatter (the blue in Fig. 22), reveals that the particle had more freedom to move, which is consistent with the fact that the adhesion between adjacent particles deteriorates as the temperature decreases. From this observation, it could be deduced, that the proposed model is successfully able to capture the influence of mix-temperature.

Conclusion
This paper presented micromechanical modeling to simulate the compaction of PA using DEM. The elastic and viscoelastic interaction between the components of an asphalt mixture during the compaction process was represented using EVCM. The macroscale inputs were determined using nonlinear regression analyses of master curves obtained from DSR tests. PA sample had been compacted and scanned using X-ray CT scans. Simulation of the entire process was carried out to generate a reliable DEM-based digital sample. The internal structure of the compacted mixtures was investigated. In order to study the influence of temperature on the asphalt compaction process, different compaction temperatures were considered. Based upon the results presented in this study, the following conclusions can be drawn: -The EVCM validation procedure demonstrated its capability to imitate a real viscoelastic behavior where the tendencies of DEM computations were consistent with the law of laboratory tests. -It was observed that the simulation outcomes were in good agreement with the experimental results in terms of reducing the sample height versus the gyrations number. Similarly, the X-ray CT scans profile showed a good match as obtained by the DEM models in terms of the air voids distribution. Overall, it could be concluded that the  proposed DEM model managed to simulate the laboratory compaction behavior reasonably well.
-Tracking the change in the air voids proportion of the field simulation indicates that the motion of aggregates is rather compound. In a way, the majority of aggregates moved in the vertical direction downwards in line with the compacting orientation. However, the slight decrease in the content of air voids across the non-compacted sectors (LL and RR) proves that the horizontal movement of aggregates from the adjacent compacted sectors (CL and CL).
-Although the DEM simulations overestimated the air void content near the top, in general, the air void distributions computed using both DEM analysis and image processing matched well. It is noted that the same set of inputs in SGC and field compactions did not produce the same level of volumetric reduction, which could be attributed to the fact that loading state and boundary conditions differ significantly between field and laboratory. -Simulation results confirmed that both laboratory and field compaction degree significantly depends on the compaction temperature. This could be associated with the fact that as the temperature decreases, the asphalt binder becomes more viscous and resistant to deformation, which results in a lesser reduction in air voids resulting from a given compaction effort. The number of gyrations required to achieve the same volume increased (5.3, 2.0 and 2.7) times with a decrease in temperature from (150 to100) • C, (150 to 120) • C and (125 to 100) • C respectively. Conversely, when the temperature is high, the asphalt binder becomes too fluid and roller loads will simply displace, or "shove" the mat rather than compact it, which is evident from the curve representing the vertical motion of the particle at a temperature of 150 • C (Fig. 22).

Table 6
The coordination of a particle.
Step  In summary, the developed model can be considered a useful tool that could help researchers gain insight into the compaction behavior of asphalt mixtures. In general, there is a good agreement between the experimental data and the analytical predictions. It is noted that further research studies could enhance the capabilities and accuracy of the proposed model. For example, simulating the aggregate is still challenging in DEM. The complex morphological characteristics, shapes and angularities of aggregates can be more realistically incorporated into the model and are recommended for future research studies.