Discrete element modelling of creep of asphalt mixtures

Creep tests on asphalt mixtures have been undertaken under four stress levels in the laboratory while the discrete element model (DEM) has been used to simulate the laboratory tests. A modified Burger’s model has been used to represent the time-dependent behaviour of an asphalt mixture by adding time-dependent moment and torsional resistance at contacts. Parameters were chosen to give the correct stress-strain response for constant strain rate tests in Cai et al. (2013). The stress-strain response for the laboratory creep tests and the simulations were recorded. The DEM results show reasonable agreement with the experiments. The creep simulation results proved to be dependent on both bond strength variability and positions of the particles. Bond breakage was recorded during the simulations and used to investigate the micro-mechanical deformation behaviour of the asphalt mixtures. An approach based on dimensional analysis is also presented in this paper to reduce the computational time during the creep simulation, and this analysis is also a new contribution.


Introduction
Asphaltic materials are traditionally modelled using continuumbased approaches and finite element programs, however this approach does not focus on the micromechanical behaviour of the mixture. An alternative approach is to use the discrete element method (DEM), which is commonly used to model the behaviour of granular materials in order to gain micromechanical insight. The first notable use of DEM to model asphaltic materials (to the authors' knowledge) was by Rothenburg et al. (1992), followed by work by authors such as Chang and Meegoda (1993), Meegoda and Chang (1994) and Buttlar (2004, 2006), all of which were 2D models. Various three-dimension models have since been developed, for example You et al. (2007), Collop et al. (2004Collop et al. ( , 2006Collop et al. ( , 2007, Carmona et al. (2007), and Liu and You (2009); these generally modelled idealised asphalt mixtures using various techniques to either represent or capture the behaviour of aggregates and binder.
Recently, the authors used a Burgers model to simulate three-dimensional constant strain rate tests using DEM (Cai et al., 2013). The samples were graded with a large number of particles, and good agreement was shown between numerical and laboratory results. This paper examines the creep behaviour of the same samples used by Cai et al. (2013); by using the same micro parameters and permitting bond breakage at the contacts, the simulated creep behaviour is compared with equivalent experimental data. As this paper is essentially a sequel to Cai et al. (2013), only key details will be repeated.

Numerical sample preparation
The material modelled in the current and previous work is Stone Mastic Asphalt (SMA), a dense gap-graded asphaltic mixture with a high content of course aggregate and binder rich mortar. The grading of the numerical material was obtained from the standard mean grading curve for a 0-10 mm SMA (BS EN 13108-5:2006), while neglecting particles with diameter less than 2 mm (Cai et al. 2014). It should be noted that this modified grading was also used in the laboratory experiments. The bitumen that is modelled has a penetration grade of 40/60 and a softening point of 51°; the modelled aggregate is granite and the filler limestone. A suitable mix design for this asphalt mixture was determined by several Schellenberg tests (BS EN 12697-18:2004) in the laboratory and is given in Table 1.
The DEM software used is PFC3D (Itasca 2008). The cylindrical samples are created within temporary rigid boundaries, using the radius expansion method (Itasca 2008, Cai et al. 2014, i.e. they are expanded from an initially small size to achieve a random packing and a specific porosity. The samples are created such that the residual stress within the sample is isotropic and less that 1% of the uniaxial compressive strength (obtained from experiments). Final adjustments are made to ensure that no particles have less than four contacts-both Rothenburg et al. (1992) and Collop et al. (2004) suggested that a minimum of four contacts per particle is required to accurately model the internal contact structure of a sand asphalt mixture. The final prepared numerical samples have a height of 120.3 mm, a diameter of 60.2 mm, and contain 6000 particles. The average number of contacts per particle is 5.5, and particles occupy 70% of the total volume. Full details of the preparation procedure are given in Cai et al. (2014), along with figures showing a typical sample.
The Burgers model ( Figure 1) can be used to capture the time-dependent properties of bitumen. The time-dependent normal contact stiffness of the Burgers model is given by: where K m is the Maxwell normal stiffness, C m is the Maxwell normal viscosity, K k is the Kelvin normal stiffness, C k is the Kelvin normal viscosity, t is the loading time, and τ is relaxation time. The contact stiffness reduces as a function of loading time. The default Burgers model in PFC3D can only transmit force via a contact bond, meaning the particles can roll relative to each other. In reality, aggregate particles cannot roll over each other when bonded, hence a modified Burgers model was developed including moment and torsional resistance in Cai et al. (2013). This Burgers model was activated as the stiffness of a contact bond in tension and compression, and a 'virtual' parallel bond was introduced to give the model moment and torsional resistance at the contact. This virtual bond does not sustain direct tension or compression (this is given by the contact bond), neither can the virtual bond breakthe contact bond governs breakage. The virtual parallel bond has a radius R = λ * min{R a , R b }, where λ is the radius multiplier, and R a , R b are the radii of the two particles, and a stiffness per unit area of k p . The stiffness of the parallel bond is obtained from the contact bond stiffness in all simulations using k p πR 2 = k c . For moment resistance, each subsequent relative rotation at a contact will result in an increment of moment which is added to the current value. This can be visualised as a virtual disk (with no strength characteristics) with a predefined radius lying on the contact plane, centred at the contact point which has bending and torsion Burgers stiffness. The total moment M is resolved into axial ( M n ) and shear ( M s ) directed moments: where n i and t i are unit vectors normal to the contact plane and lying in the contact plane respectively, and the increments of each component are given by: where k s and k n are the Burgers stiffnesses used in torsion and bending, Δθ n and Δθ s are the rotation increments due to torsion and bending, and J and I are the polar moment of inertia and moment of inertia of the disk cross section respectively (Cai et al. 2013). Cai et al. (2013) used this modified Burgers model to simulate uniaxial constant strain rate tests, and compared simulations with laboratory tests under three strain rates. They simulated both normal and Weibull distributions of bond strengths, where the bond strengths were assumed to be a function of the macroscopic strain rate (and were scaled accordingly): where σ 0 and _ ε 0 are reference values of stress and strain rate, and n = 0.325 (Cai et al. 2013). The effects of the bond strengths variability, as well as of different random samples were shown to be negligible for each strain rate investigated.
The sample parameters used in this work are set the same as in Cai et al. (2013), and are given in Table 2. Wei et al.
(accepted) presented a calibration of these parameters, including the ratio between normal and shear contact parameters (N/S = 11). The same bond strength distributions are also used, i.e. the mean and standard deviation of the normal distribution are 54 N and 44 N, and the modulus and characteristic scale factor for the Weibull distribution are 1.1 and 60 N respectively for a strain rate of _ ε = 0.005 s −1 . Hence, if the numerical behaviour demonstrates agreement with the  experimental data, it shows that it is possible to capture the essential features of asphalt behaviour using DEM and a single set of parameters.

Experimental materials
The same samples as used in laboratory tests by Cai et al. (2013) are used here. A modified grading curve to match that of the simulations was produced by removing particles that passed a sieve size of 2 mm. Granite and limestone were chosen as aggregate and filler, and a 40/60 penetration grade bitumen with a softening point of 51°was used; see Table 1 for the mix design. Based on both the numerical sample sizes and the available laboratory drills, samples 140 mm in height and 70 mm in diameter were produced for the laboratory tests. An aspect ratio of 2 was chosen to avoid confinement effects in the specimen due to friction of the loading platens. Full details of the specimen preparation are documented in Cai et al. (2013). An INSTRON device with a temperature-controlled cabinet was used for the creep tests presented in this paper. The axial deformations and applied axial force on the tested specimens were recorded using Rubicon software, and the axial deformations were measured by linear variable differential transformers.

Experimental results
The results of creep tests performed at 400, 500, 600, 1000 kPa at 20°C are shown in Figure 2. The axial strain is plotted against the time elapsed following the application of the load. As can be seen, the creep curve can be divided into three parts: a primary creep region where the strain rate decreases, a secondary creep region where the strain rate is almost constant, also known as the steady-state strain rate and a tertiary creep region where the strain rate increases rapidly as the specimen approaches failure. At the same loading time, the axial strain increases as the stress level increases. The sample tends to rupture earlier and more catastrophically as the stress increases.
Two more creep tests were undertaken under 400 and 500 kPa with a different type of granite (porphyritic andesite) at 20°C. The results are also shown in Figure 2. It can be seen that, even for the same type of material (granite), different sources lead to variability in the test results. In addition, the variability can also be caused by the variations in the air voids content in the specimens. Figure 3 shows four of the specimens used in this study. As can be seen, in some places, the sample is porous, while in others, it is quite dense. This variability can be reduced by using a gyratory compactor to make the specimens, but, in this case, the specimens would need to be made one by one which would mean that more material would be wasted and take longer to make the specimens. The new SMA only comprises aggregate larger than 2 mm in diameter, filler

Dimensional analysis
In DEM, it is highly desirable to develop new methods which are sufficient to describe the 3D mechanical behaviour but which keep computational time to a minimum. McDowell et al. (2009) derived a scaling method by using dimensional analysis for viscosities and velocity in DEM simulations of strain control tests on asphalt. They proved that the effect of scaling the applied velocity is the same as that of scaling both viscosities in the Burgers model by the same factor. This dimensional analysis is also used to accelerate the creep test simulations. The analytical expression for the time-deformation response captured by the Burger's model to a force F 0 applied at time t ¼ 0 is given by: More fundamentally, the deformation can be written as a function of the following six physical parameters: As indicated, the six parameters contain three physical dimensions: Length L, mass M, and time T. Based on the dimensional analysis, the physics of the problem can be described in terms of one dependent dimensionless group and three independent dimensionless groups. One of the possibilities is: It is now evident that if C k and C m are both increased/ decreased by the same factor x, then this is equivalent to scaling time by 1/x. The displacement calculated from Equation (6) is plotted against time in Figure 4a using the parameters given in Table 3. Figure 4b shows the analytical solution for the case which scales down both viscosities by a factor of 10. From  Figures 4a and 4b, it can be seen that, for the same displacement (0.3 m), the unscaled case takes 2000 s while the scaled case only needs 200 s. Figure 4c shows the analytical solution for unscaled parameters (Table 3) compared with increasing time by a factor of 10 and reducing both viscosities by a factor of 0.1. The two curves are identical, demonstrating that Equation (8) is correct. This scaling method is used in the creep simulations in this project to reduce the computational time. It is worth noting that You et al. (2011) also developed an approach to reduce the computational time with a time-temperature superposition principle. They also proved by a different analysis that the same creep compliance can be determined by scaling the Burger's viscosities and by scaling time. However, Equation (8) presents the argument in its fundamental and elegant form and is considered a new contribution within this paper.

Qualitative comparison with experimental results
For comparison with the experiments, creep simulations under 400, 500, 600, and 1000 kPa were performed, using the same parameters as the previous constant strain rate tests and given in Table 2. The bond strengths were set to be power-law function of strain rate (Cai et al. 2013). Six simulations were modelled for each loading condition with identical initial samples, i.e. for each load, two alternative bond strength distributions were used: Weibull and normal; and for each of the two distributions, three alternative random number seeds were considered for the same mean and standard deviation. The results are given in Figure 5, and as can be seen, at higher stresses (500, 600 and 1000 kPa) the bond strength variability can cause rupture (axial strain increases rapidly) at different times. Hence, the precise bond strength distribution within the sample plays an important role in the rupture of the sample. However at lower stress (400 kPa) the effect of bond strength variability is negligible. This is because at higher loadings, the    Figure 5 show similar patterns of behaviour compared with the laboratory results in Figure 2. The effect of the positions of the particles was also investigated. Four further simulations were performed, also with 6000 particles but each with different random positions (different random number seed for particle generation) for each loading condition (with same parameters given in Table 2). For each of these two additional configurations, the bond strengths followed alternative normal and Weibull distributions; for each normal distribution of bond strengths the seed was constant and for each Weibull distribution the seed was constant, so that only the particle configuration and type of distribution had any effect. The results are shown in Figure 6. As can be seen, the positions of the particles have an effect at higher stresses (500, 600 and 1000 kPa). Rupture can occur at different times and the axial strains can be different at the same loading time. This could also partly explain the variability of the creep results in the laboratory. The Burger's model is loading speed dependent.
In the creep simulations, bond strengths change as the strain rate changes (Equation (5)) which can lead to a different bond breakage pattern. At a lower stress (400 kPa) the effect of the positions of the particles is negligible. Bond breakages were also recorded during the creep simulations, shown in Figure 7. As can be seen, at a lower stress level (400 kPa), only a few bonds break after application of the load. This agrees well with laboratory results: no rupture was observed at lower stress levels. At higher stress levels, rupture corresponds to the highest bond breakage rate; see for example, Figure 7d. This behaviour is consistent with that observed in Cai et al. (2013) where a higher strain rate resulted in more brittle behaviour; in the creep simulations, a higher stress results in well-defined creep rupture when there is a rapid onset of bond breakage. Based on the bond breakage curves (Figure 7), the following general trend can be observed; bond breakage rate increase at the primary creep stage, and remains almost at constant at the secondary creep stage and increases again rapidly at the tertiary creep stage.

Conclusions
This paper presents the use of the modified Burger's model with moment and torsional resistance, for samples subjected to constant strain rate tests in the previously published paper (Cai et al. 2013). However, the purpose of this paper has been to ascertain whether those samples would also exhibit realistic behaviour under creep, using the same micro parameters. An approach was presented based on a new dimensional analysis to reduce the computation time during the creep simulation. The creep test simulations have shown reasonable agreement with laboratory results, using the previously determined parameters to model constant strain rate tests. In the modelling, the bond strength has been set to be a power-law function of strain rate as in the previous publication (Cai et al. 2013). The behaviour under creep is dependent on the variability of the bond strengths and the positions of the particles. This offers some explanation of the variability in experimental results, which follow a similar pattern. Bond breakage was recorded during the simulation and used to investigate the micro-mechanical deformation of the asphalt mixture. Creep rupture was observed to correspond to the maximum rate of bond breakage. Given that the micro parameters were chosen to model constant strain rate tests, it is reassuring that the same properties have resulted in creep behaviour for the samples which closely resembles that observed in the laboratory. In other words, the previous paper Cai et al. (2013) showed that constant strain rate tests could be modelled using DEM, but this paper has shown that, given that the correct creep behaviour is also exhibited by the same samples, that DEM is a very useful tool for capturing the time-dependent behaviour of asphalt mixtures under the most common loading conditions, namely constant strain rate and constant load. Future work can use the model described in this paper to examine the effect of temperature. It is well known that for asphalt there is a time-temperature equivalence. Here we have shown that for parameters calibrated to model constant strain rate tests, those parameters have also captured the creep behaviour of the same modelled material under different loads; therefore, given that the time-dependence of the material behaviour has been captured, it should be relatively straightforward to model temperature effects as a result. In addition the particle size distribution and particle shape will have an effect on performance; this can be modelled quite simply, and eventually, DEM is expected to help with the development of constitutive models and performance based models to solve highway industrial problems like rutting and fatigue.