Discrete element modelling of uniaxial constant strain rate tests on asphalt mixtures

Constant strain rate tests for a graded asphalt mixture under three constant strain rates have been undertaken in the laboratory. The Discrete Element Model has been used to simulate the laboratory tests with a numerical sample preparation procedure being developed to represent the physical specimen. The Burger’s model has been used to represent the time dependent behavior of the asphalt mixture. The Burger’s model was implemented to give bending and torsional resistance as well as in direct tension and compression. The stress-strain response for the laboratory tests and the simulations under three loading speeds were recorded. The results show reasonable agreement when the bond strengths in the model are made to be a function of strain rate. Both normal and Weibull distributions have been used for the bond strengths between the aggregate particles. The effects on the stress-strain response of bond strength variability and particle position are proved to be negligible. Bond breakage was recorded during the simulations to explain the internal damage within the sample. The modified Burger’s model has proved to be a useful tool in modeling the bending and torsional resistance at particle contacts in an asphalt mixture, in order to correctly predict observed behavior.


Introduction
The traditional approach to model asphaltic materials is to treat them at the macro-scale using continuum-based models which are implemented into a finite element program. In this approach, the micromechanical behavior of the mixture is not included. An alternative approach is to use the Discrete Element Method (DEM), which has been widely used to model the behavior of granular materials, in order to gain micro-mechanical insight.
Cundall [1] introduced a simple computer program to model the progressive failure of a system of discrete blocks. In that paper, the basic theories of DEM were described. Since then, DEM has been widely used to model the mechanical behavior of soil and rock. However, to the knowledge of the authors, it has not been applied to model the mechanical behavior of asphalt mixture until 1992. Rothenburg et al. [2] proved that it is possible to simulate the rutting problems in the pavement by using a discrete element model. In their model, the particles are treated as elastic elements and the binder as a linearly visco-elastic material.
Later, Chang and Meegoda [3] developed the ASBAL program by modifying the TRUBAL [4,5] DEM based program which is limited to modelling dry contacts of particles. They added combinations of springs and dashpots to simulate the behavior of asphalt concrete. The program has been used to simulate simple triaxial tests under simple loading [3] and rutting under repeated loading [6]. They have shown that the discrete element model is a useful tool to study the fundamental properties of asphalt concrete. They further used ASBAL to simulate triaxial tests on Hot Mix Asphalt under different loading combinations where the Burger's model was added in the normal and tangential directions at the contact between the particles [7]. Buttlar and You [8] introduced microfabric discrete element modelling (MDEM) by using the commercial DEM package PFC2D. In this model, various material phases were modelled as clusters of very small discrete elements and the microstructure was obtained by optically scanning smoothly sawn test specimens. A high resolution scanner was used to obtain gray-scale images of the sections. The MDEM has all the benefits of traditional DEM and can also model complex aggregate shapes and propagation of cracks around or through aggregates during a test. They used MDEM to predict the asphalt mixture complex modulus in extension/compression across a range of test temperatures and loading frequencies [9], the complex modulus of asphalt aggregate hollow cylinders subjected to internal pressure [10], the compressive dynamic moduli of asphalt mixtures [11] and the creep stiffness of asphalt mixtures [12]. The results proved to be reasonable and useful when compared with laboratory results. Except, for some fine mixes, the DEM approach provided a lower prediction of the mixture moduli compared with experimental tests. This is due to insufficient aggregate to aggregate contacts in the 2D model.
You et al. [13] developed a three-dimensional microstructure based discrete element model of asphalt to study the dynamic modulus from the stress-strain response under compressive load. The 3D microstructure of the asphalt mixture was obtained from a number of 2D images. They demonstrated that the 3D discrete element models were able to predict the mixture moduli across a range of temperatures and loading frequencies. Liu and You [14] introduced a randomly created polyhedron method where the aggregates are simulated with randomly created polyhedral specimens comprising a large number of spherical particles. The sample preparation method did not require 2D images from experiments; however, it proved to be time consuming due to the large number of particles involved.
Collop et al. [15] investigated the use of DEM to simulate the behavior of a highly idealized bituminous mixture comprising single-sized sand particles in a uniaxial compressive creep test. The effect of the bitumen was represented by time-dependent shear and normal contact stiffnesses. They used the model to predict the dilation in uniaxial compression tests for an axial stress of 400 kPa and triaxial tests for deviator stresses of 400 and 600 kPa at stress ratios of 0.6 and 0.8 [16]. The numerical results were validated with laboratory data. The comparison between the predicted and measured results has shown that the model is able to predict the effect of stress ratio on dilation. Later, they used the model to simulate the viscoelastic deformation behav-ior of an idealized asphalt mixture [17]. An elastic contact was used for the compressive normal contact stiffness and a viscoelastic contact was used for shear and tensile normal contact stiffness.
Carmona et al. [18] extended the two-dimensional discrete element model to capture microscopic failure mechanisms relevant for the process of fatigue of asphalt mixtures. They used polygons and beams to represent aggregates and binder respectively. The effect of healing has also been included. Later, Kun et al. [19] used DEM and a fiber bundle model to study the fatigue fracture of heterogeneous materials based on Basquin's law of fatigue.
Wu et al. [20] used DEM to simulate constant strain rate compressive tests for an idealized asphalt mixture. In their simulations, the bond strength was set as a power-law function of temperature compensated strain rate. It was proven that by using the strain rate dependent bond strength in the DEM simulations, a good agreement with the experimental data in both the pre-peak and post-peak regions of behavior was obtained.
Cai et al. [21] used PFC3D to model the asphalt mixtures with graded aggregates. A parallel bond comprising a set of elastic springs with constant normal and shear stiffnesses distributed over a circular disk with pre-defined radius was used to represent the elastic properties of asphalt mixtures. A minimum of 6,000 particles was needed to obtain a reliable Young's modulus and Poisson's ratio; a maximum of 0.02 m/s loading speed could be used to avoid dynamic stress wave propagation within the 6,000 graded particles sample. The Poisson's ratio was found to increase with the ratio of contact normal to shear stiffnesses. The Young's modulus was found to depend on both normal and shear contact stiffnesses. The Burger's model comprising a spring and dashpot in parallel, connected in series with a spring and a dashpot, were used to capture the time dependent properties of bitumen ( Fig. 1). In the viscoelastic simulation, the pre-peak slope of the stress-strain curve in the uniaxial compression tests at a constant strain rate of 0.005 s −1 was found to increase as the Maxwell stiffness increases and was independent of the Maxwell viscosity, Kelvin stiffness and Kelvin viscosity. A wider post peak softening was predicted as the ratio of normal to shear Burger's model parameters increases while peak stress increases as the inter-particle friction coefficient increases and a bond strength distribution (uniform, normal or Weibull) leads to a wider post peak softening behavior while reducing the peak stress compared with a single value of bond strength.
This paper provides a new contribution regarding the use of DEM to simulate the micromechanical behavior of a graded asphalt mixture based on the research mentioned above. Constant strain rate tests were performed in the laboratory and used to compare with the simulations obtained from DEM. The simulations were carried out by using Particle Flow Code in Three Dimensions (PFC3D [22]). A numerical sample preparation procedure has been developed to prepare cylindrical samples with graded particles. To prepare a sample which is initially isotropic with low residual stress and density packing, the following procedure was adopted.
At first, the cylindrical boundaries (walls in PFC3D) were generated according to the predefined volume. The wall in PFC3D has arbitrarily defined contact properties and the contact forces between the particle and wall are calculated by a force displacement law.
The number of particles for each diameter was calculated according to the grading curve and mixture design and given by: where range of particle radii is divided into n sub-ranges with lower and upper limits of r n and r n+1 respectively; N is the number of particles (with radius (r n+1 + r n ) /2) to be generated; P is the volume percentage of aggregate particles in the asphalt mixture; V is the total volume of the asphalt sample; p n+1 and p n are the percentage passing (by volume) of aggregate through sieve size 2r n+1 and 2r n . Particles are generated randomly to obtain an irregular packing. Particles are first generated with a small radius such that no particles overlap, then the radii are increased gradually to the target values. During the process, particles are allowed to re-arrange until the residual stress within the sample is approximately isotropic.
A sample prepared under the above conditions has a high level of isotropic stress due to the large amount of overlap at contact points between particles. ITASCA [22] suggested that the initial isotropic stress should be less than one percent of the uniaxial compressive strength (peak stress obtained from experimental data). This can be achieved by slightly reducing the radii for all the particles by 1-2 % of the radius value until the isotropic stress is less than the pre-set value which is determined from the experimental data (one percent of the peak stress). By doing this, the magnitude of overlap between particles decreases which leads to a reduction in contact forces in the sample. The particles are also allowed to re-arrange each time after decreasing the radii.
Typically, there are about 10-15 % particles with less than four contacts in the sample at this stage. Contact bonds are used to simulate the bonding effect of bitumen to the nearby aggregate. Rothenburg et al. [2] suggested that the particle assembly needs at least four contacts per particle on average to carry the load. Collop et al. [15] also showed that a minimum of four contacts per particle on average is needed to model the internal contact structure of a sand asphalt mixture. This can be achieved by the following procedure: (i) Fix the position of all the particles; (ii) Release the particles with less than four contacts (one by one); each released particle is taken in turn and expanded in increments of 0.1 % until four contacts are obtained and then the next released particle follows the same procedure. During this procedure, these particles are allowed to re-arrange to reach equilibrium.
The sample is then bonded by contact bonds. To simulate the uniaxial compression test, the lateral boundary is removed while the top and bottom boundaries act as loading platens. To reduce the unbalanced force due to particle radius enlargement as well as the removal of the lateral wall, the bonded sample is allowed to reach equilibrium. Figure 2a shows a typical sample. The sample containing 6,000 particles is 120.3 mm in height and 30.1 mm in radius; the numbers of particles for the different sizes are calculated by Eq. (1) and shown in Table 1. The average number of contacts per particle for the sample is 5.5 and particles occupy 70 % of the total volume. The bonded sample is shown in Fig. 2b, used to simulate the uniaxial compression tests.

Laboratory work
Stone mastic asphalt (SMA) is a dense gap-graded asphaltic mixture. It has been selected for this research because of its high coarse aggregate content and low fines content. The standard mean grading curve for 0-10 mm SMA [23] is show in Fig. 3 as the solid line. The aggregate sizes range from less than 0.063 to 14 mm. However, due to the limitation of computer processing power and simulation time, it is impossible to model all the fines in DEM. Therefore, a modified SMA has been produced by taking out particles which would pass a sieve size of 2 mm. The grading curve of the aggregate for the new asphalt mixture is shown by the dashed line in Fig. 3. Granite and limestone were chosen as aggregate and filler; a 40/60 penetration grade bitumen with a softening point of 51 • C was used. Several Schellenberg tests [24] were undertaken to obtain the required binder content for different combinations of graded aggregate, bitumen and filler. The Schellenberg test describes a method for determining binder drainage of bituminous mixtures. It is applicable to asphalt materials that are not porous asphalt or for those porous asphalts incorporating fibers. It can be used either for determining the binder drainage for different binder contents, or with a single binder content, eliminating successive repetitions. It also enables the effects of varying fine aggregate types or including any anti-draining additive to be quantified. A suitable design for the new asphalt mixture was determined ( Table 2).
The mixing of the aggregate, bitumen and filler was achieved by using the hot mix method [25]. The mixture was then compacted into 305 × 305 × 85 mm slabs using the Nottingham roller compactor [26]. The specimens were made by coring through the width of the slabs using a wet diamondtipped core. To avoid confinement effects in the specimen due to the friction between the specimen and the loading platens, the aspect ratio (ratio of height to diameter) of the specimen is usually chosen to be equal to 2 [27]. Cai et al. [21] found that the minimum specimen dimensions to give reasonably accurate estimates of bulk material properties in computer modeling are 60.16 mm in diameter and 120.32 mm in height (6,000 particles). Based on the above conditions and the diameter of available coring drills in the laboratory, the dimension of the test specimen in the laboratory was determined as 70 mm in diameter and 140 mm height. An INSTRON device with a temperature controllable cabinet was used for the uniaxial tests undertaken in this project.   The equipment consists of a 100 kN servo hydraulic actuator with ± 50 mm movement. The loading frame is operated by a 'Rubicon' digital servo control system. For a uniaxial constant strain rate test, a monotonic compression displacement is applied to the specimen through the hydraulic actuator. The axial and radial deformations were measured at pre-set intervals during testing by axial and radial Linear Variable Differential Transformers (LVDTs) mounted on the specimen vertically and horizontally (at mid height of the specimens). LVDTs have high resolution and can measure movements as small as a few millionths of a millimeter up to a few centimeters by suitable conditioning hardware and high resolution data acquisition. The constant strain rate tests were undertaken at three different strain rates of 0.005, 0.0005 and 0.00005 s −1 at 20 • C.
A typical result from the constant strain rate of 0.005 s −1 and 20 • C is shown in Fig. 4; the axial stress is plotted against the axial and radial strains. As can be seen for both curves, the axial stress versus axial and radial strain relationships are similar consisting of a hardening portion up to a peak stress followed by a softening portion. The onset of dilation occurs  Figure 5 show the axial stress against axial and radial strain for the three strain rates tests at 20 • C. The following general observation can be made: the peak stress increases when strain rate increases and the axial strain level corresponding to the peak stress decreases with increasing strain rate. Figure 6 summaries the peak stress as a function of strain rate on double logarithmic scales for these three tests. It can be seen that the peak stress is approximately proportional to strain rate when plotted on double logarithmic scales with a gradient of approximately 0.325 (n). This can be expressed by the power-law relation as where σ 0 is a reference stress (3.19 MPa at 20 • C) andε 0 is a reference strain rate (0.0005 s −1 at 20 • C). The gradient observed from Fig. 6 will be used in the further simulations to scale the contact normal and shear bond strengths.

Simulation
The ratio of sample height to diameter was set to be 2:1 (as for the laboratory specimen). The time dependent normal contact stiffness of the Burger's model (Fig. 1) is given by where k m is Maxwell normal stiffness; C m is Maxwell normal viscosity; k k is Kelvin normal stiffness; C k is the Kelvin normal viscosity; t is the loading time; τ = C k /K k is a relaxation time. The contact stiffness reduces as a function of loading time. The default Burger's model in PFC3D can transmit only force via a contact bond which allows the particles to roll over each other. However, in reality aggregate particles are not allowed to roll over each other in bonded aggregates. For this purpose, an alternative approach was developed.
The Burger's model was activated as the stiffness of the contact bond in tension and compression, and a "virtual" parallel bond was introduced to give a Burger's model in bending resistance and torsional resistance at a contact. The term "virtual" is used because the bond doesn't sustain direct tension or compression (this is given by the contact bond), nor can the virtual bond break; the contact bond governs breakage. If the current stiffness of the contact bond is k c , and the radius of the virtual parallel bond is R, and distributed stiffness k p per unit area, then the stiffness of the virtual parallel bond was taken as k p π R 2 = k c through all simulations.
For the moment resistance, each subsequent relative rotation increment at contacts will produce an increment of moment and adds to the current value. This can be imagined as a virtual circular disk (with no strength characteristics) lying on the contact plane and centered at the contact point which has bending and torsion Burger's stiffness with a chosen radius. The total moment at a contact can be resolved into respective normal (torsionM n ) and shear ("moment" M s ) components with respect to the contact plane as M =M n n i +M s t i (4) where n i is the unit normal to the contact plane and t i is a unit vector in the contact plane. The bending moment acts along the contact plane and the torsion moment acts normal to the contact plane. The increment of normal (torsion) and shear (bending) moments are given by where k n , k s are the Burger's stiffnesses for bending and torsion respectively; θ n , θ s are rotation increments due to the torque and moment, respectively; I and J , are the moment of inertia and polar moment of inertia of the image cross section, respectively, and are given by where R is the radius of the imagined circular disk; λ is known as the radius multiplier; R A and R B are the radii of particles in contact. It should be noted that the relative angular velocity which is used to calculate the rotation increment is stored as a vector in global coordinates in PFC3D and must be updated to account for the contact motion before the calculation. Finally the calculated moment needs to be updated in terms of global coordinates. A typical three-dimensional microstructural based discrete element viscoelastic modelling process is extremely Fig. 7 Bond strength normal distribution and Weibull distribution (modulus is the shape factor relating to the coefficient of variation and alpha is the scale factor (proportional to mean), see Eq. (10) for the probability density function) time-consuming. McDowell et al. [28] found a scaling method by using dimensional analysis for viscosity and velocity in discrete element modelling of strain-controlled tests on asphalt. They have proved that the effect of scaling applied velocity is the same as that of scaling both viscosities by the same factor. This scaling method, which offers great benefit in saving computer simulation time, has been used for the strain-controlled simulations in this paper.

Comparison with experimental results
In the constant strain rate simulations, normal and Weibull distributions have been used for the contact bond strengths. The mean and standard deviation for the normal distribution were set as 50 and 50 N, respectively. In PFC3D, when the bond strength value is set to be less than zero, no contact bond will be created; in this case the bond strength is not in a normal distribution any more. After consideration of all the bond strength values less than zero and replacing these strengths with zero values, the mean and standard deviation of this final bond strengths distribution are calculated as 54 and 44 N, respectively (Fig. 7). This distribution will still be referred to as the normal distribution even through the nonzero part has effectively been replaced by an impulse function at zero strength.
The modulus (shape parameter) and alpha (scale parameter) for the Weibull distribution were set as 1.1 and 60 N, respectively (Fig. 7). The probability density function of a Weibull distribution is given by [29] f (x; α, m) = m α where m > 0 is the shape factor and α > 0 is the scale factor.    Table 3. These values were chosen by trial and error to produce reasonable numerical results similar to previously presented experimental data at a constant strain rate of 0.005 s −1 (Fig. 5). The Burger's parameters could have been calibrated at a different strain rate, and then applied to the other two strain rates; however it is reassuring that in this case the chosen parameters have captured the essential features of the behavior across the three strain rates, even though the simulation do not exactly match the experimental data over a large range of strain.
Asphalt is a viscoelastic material for which the deformation depends on loading speed. Previous laboratory constant strain rate test results show that the peak stress is a function of strain rate and can be expressed by a power-law equation. In the constant strain rate test simulations, the previously justified power-law exponent (n) of 0.325 has been used in scaling both normal and shear bond strengths in the modelling for different strain rates. A numerical sample containing 6,000 particles was generated according to the sample preparation procedure mentioned above. The stress-strain response for   the simulations under three strain rates are shown in Fig. 8 with the laboratory results. In the simulations, for lower strain rates (0.0005 and 0.00005 s −1 ), both viscosities have been scaled instead of loading rate, as if both are scaled by the same factor as the strain rate, then this is equivalent as shown by dimensional analysis by McDowell et al. [28]. The results show reasonable agreement between laboratory and simulations with bond strengths having normal and Weibull distributions (Fig. 7). The quantitative differences between the numerical samples and the experimental result for the softening behavior at a strain rate 0.005 s −1 may result from the fact that in the numerical samples, when a bond breaks the viscoelastic behavior at the contact vanishes and visco-elastic behavior is only possible in compression and shear if the contact reforms; however, in the laboratory samples bitumen is present after bond breakage and thus visco-elastic behavior may still occur in compression, shear and torsion (not tension) if the adjacent particles should come into contact. The percentages of the normal, shear and total bond breakage can be recorded during the simulations. Figures 9 and  10 show the bond breakage for normal and Weibull bond strength distributions (Fig. 7) at three strain rates. As can be seen, for all the strain rates, the bond breakage increases as the axial strain increases and the maximum bond break-

Fig. 12
Axial stress versus a axial strain and b radial strain for three normal distributions of bond strengths and laboratory at 20 • C age rate corresponds to the peak stress. At lower strain rates (0.00005 s −1 ) fewer bond breakages are observed. Figure 11 shows the damage modes and locations of broken bonds at different axial strain levels for a constant strain rate of 0.005 s −1 with a normal distribution of bond strengths (Fig. 7). A grey circle represents a bond breaking in the normal direction while a black circle represents a bond breaking in the shear direction. The figure shows damage throughout the sample. Individual bond breakage can be tracked and the figure shows extensive internal damage.
As can be seen, the micro-cracks in both normal and shear damage modes are approximately uniformly distributed throughout the sample. An increase in the number of micro-cracks (individual bond breakages) can be observed as the axial strain increases. Few micro-cracks are observed at an axial strain of 0.5 % indicating that the damage to the sample at this level of strain is largely internal. However, at an axial strain of 1 % a dramatic increase in the number of micro-cracks is observed indicating that macroscopic cracks (a number of adjacent bond breakages) are forming in the sample. This is associated with the onset of the post peak softening region in Fig. 8.

Effect of random bond strength distribution
The bond strengths follow normal and Weibull distributions (Fig. 7). During the simulations, values of bond strengths are randomly chosen from the normal or Weibull distributions using a random number seed. Therefore, even for the same distribution, different bond strength values could be applied As can be seen, the axial stress is approximately independent of the variability of the bond strengths. The difference in behaviour appearing after the peak stress is due to contact bonds breaking at different times and positions due to the variability of the bond strengths.

Effect of particles position
The numerical sample preparation method randomly places the particles in a pre-defined area using a random number seed. Two separate constant strain rate simulations have been performed with 6,000 particles in different random positions at a strain rate of 0.005 s −1 for both the normal and Weibull distributions of bond strengths in Fig. 7. The sample Burger's model parameters were set to be the same as in the previous sections ( Table 3). The results of axial stress as a function of axial strain are shown in Fig. 14. It can be seen that the variability in the stress-strain behavior caused by the random positions of the particles in constant strain rate simulations is negligible.

Conclusions
This paper presents the use of a modified Burger's model with bending and torsional resistance to simulate constant strain rate tests on graded asphalt mixtures. The standard Burger's model in direct tension and compression has also been included at the contact bonds. Laboratory tests were performed on the same asphalt mixture under the same loading conditions as the simulations. Both normal and Weibull distributions have been used for the bond strengths. In the modelling, the bond strength has been made to be a powerlaw function of strain rate. Reasonable agreement is observed between simulation and laboratory results. The constant strain rate tests results are approximately independent of the variability of the bond strength and the positions of the particles for each of the constant strain rate tests. The simulations performed in this project show the potential of DEM as a useful tool to study the fundamental properties of asphalt mixtures.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.