Generalised Kelvin contact models for DEM modelling of asphalt mixtures

ABSTRACT

Rigid particle models based on the discrete element method (DEM) have been adopted to model creep, fracture, and the viscoelastic behaviour of asphalt mixtures considering an irregular micro-structure and particle contacts. Within a DEM framework, the Burgers contact model, which is known to have a narrow frequency and temperature range, is usually adopted to model viscoelastic properties. In this study, a new explicit three-dimensional generalised Kelvin (GK) contact model formulation for the DEM model is proposed for asphalt materials. The model is implemented following two different methodologies (GK${}_1$ and GK${}_2$). The models are validated in uniaxial tension-compression sinusoidal tests for predicting the dynamic modulus ($\left|E^{\ast }\right|$) and phase angle (ϕ) of these composites at a frequency range of 1–10 Hz at $20{}^{\circ }$C. Four mixtures are investigated based on the modelling of their mastic. The influence of the GK contact parameters on the dynamic response of mastics is assessed. Results show that ${\kappa }_m$ has an important influence on both rheological properties and that ${\eta }_m$ can be used for small adjustments focussing on the predicted phase angle. Moreover, it is shown that a viscoelastic contact model should be adopted to simulate aggregate-to-mastic contacts in mixtures. As expected, the response obtained for both GK models is the same but the simulations with the GK${}_1$ are 14% faster. In addition, the response predicted with the proposed GK contact model is compared with the response predicted with a Burgers contact model. The DEM predictions obtained for the GK model are more accurate. For mastics, the average errors for $\left|E^{\ast }\right|$ and ϕ when adopting the GK model (Burgers) are 2.40% (3.46%) and 3.64% (4.17%), respectively. For mixtures, the average errors for $\left|E^{\ast }\right|$ for the GK model (Burgers) are 7.00% (7.92%). The proposed contact model greatly enhances the DEM ability to simulate the viscoelasticity of asphalt materials.

KEYWORDS:

• Asphalt mixture
• discrete element method
• three-dimensional model
• viscoelasticity
• dynamic behaviour

1. Introduction

Asphalt mixture is a heterogeneous material generally containing fine and coarse aggregates, air voids, and bitumen. During the manufacturing process, fine aggregates and bitumen mix together originating the asphalt mastic, which provides the time-temperature dependency to the asphalt matrix and influences important mechanical properties, such as the dynamic modulus, cohesion, permanent deformation resistance, and fatigue cracking (Zhang et al. 2019).

The finite element model (FEM) and the discrete element model (DEM) are two common numerical simulation techniques that are adopted by researchers to better understand, on a micro-level, the mechanical performance of asphalt materials. Compared to FEM, DEM has the following advantages: (1) large deformation and fracture are easily modelled (Pei et al. 2022), (2) the discrete and heterogeneous nature of asphalt materials is taken into account, and (3) adopts simple constitutive laws (Xie et al. 2022). Thus, the DEM is adopted in this study. The DEM method was first introduced in rock mechanics (Cundall, 1971). Regarding asphalt mixtures, one of the first applications was carried out by Rothenburg et al. (1992), who applied a two-dimensional (2D) DEM micro-mechanical model to analyse pavement rutting.

In initial studies to investigate bituminous materials with DEM models, elastic models have been adopted to simulate the interaction of the different constituents (Khattak and Roussel, 2009), not taking into account the time-temperature dependency nature. Nevertheless, when compared with more complex models, elastic models have lower computational costs. Viscoelastic models have later been developed and implemented in order to account for the time dependent behaviour of asphalt mixtures (Feng et al. 2014, 2015). With this purpose, the Burgers model has been widely adopted in the numerical modelling of asphalt materials. Recently, the Burgers contact model has been used to investigate the rheological behaviour (Zhang et al. 2019), creep stiffness (Wang et al. 2021), compaction (Gong et al. 2018), fracture (Wang and Buttlar, 2019), and indirect tensile (IDT) strength (Li et al. 2019). By investigating the influences of mechanical properties of aggregates and mastic on the creep behaviour of mixtures, Ma et al. (2017) showed that the heterogeneity in the coarse aggregate matrix had a great impact on the creep deformation. Cai et al. (2016) studied the creep behaviour in mixtures subjected to different stress levels, showing a reasonable agreement between the DEM predictions and the available experimental data.

Ma et al. (2018) reproduced the wheel tracking test and the high-temperature rutting deformation in asphalt mixtures in 2D, highlighting that the aggregate skeleton must have good stability and high stiffness to avoid large deformation. Following a similar investigation, Zhang et al. (2019) analysed the characteristics of aggregates on the permanent deformation of asphalt mixtures and mastics at high temperatures under uniaxial static creep test, where it is shown that a non-uniform distribution of aggregates has a negative impact on the rutting deformation. Peng et al. (2020) analysed the impact of properties, e.g. aggregate size, temperature, binder content, and air void, on the shear fatigue performance of mixtures using a three-dimensional (3D) DEM model, showing that all parameters exhibited an important influence on the shear fatigue life. Yi et al. (2021) proposed a method to obtain the parameters for a Burgers model based on a nano-indentation process. Peng et al. (2022) investigated the permanent deformation, shear stress, and transverse strain of crumb rubber-modified asphalt pavements under traffic loads, and the results presented showed that, at higher temperatures, larger shear stress is developed and that the mechanical response is also affected by the number of repeated loads and crumb rubber content.

In an effort to improve the DEM viscoelastic modelling simulations performance, Ren and Sun (2016) proposed a generalised Maxwell model in 2D to simulate the viscous response of mixtures at low temperatures under different testing conditions. The generalised Maxwell model has also been applied in the fracture analysis of asphalt mixtures (Ren and Sun, 2017, Sun, 2018) under different fracture loading modes in semi-circular bending tests at low temperatures. More recently, Ren et al. (2020) proposed a 2D generalised Kelvin contact model and assessed its performance in the dynamic behaviour of mixtures, comparing its numerical predictions with those obtained with a Burgers model and with a generalised Maxwell model. The results obtained show that the generalised Kelvin model predicts well the dynamic response at high temperatures, while the generalised Maxwell model is more suitable to predict the dynamic behaviour at low temperatures.

In summary, research previously carried out using the DEM provide meaningful insights into the behaviour of asphalt mixtures, that are difficult to be tracked or monitored experimentally. However, DEM research has been mainly focussed on the use of the Burgers model, which is known to have a narrow frequency and temperature range, and many have been developed in two dimensions. Given the above-mentioned limitations, the development of 3D DEM models, especially for larger structures, and the implementation of more complex viscoelastic models to better represent asphalt materials are still areas that require further research.

2. Objectives and scope

The objective of this study is to develop and calibrate a generalised Kelvin (GK) contact model with the use of a 3D DEM to investigate the viscoelastic behaviour of asphalt materials. Therefore, two different approaches are developed. In the first (GK${}_1$), a direct integration of the constitutive equations, using a time-centred difference scheme, similar to the adopted for the Burgers model in the study of margarines (Shimizu et al. 2013), is devised. The second method (GK${}_2$) is based on an incremental force-displacement formulation, used previously in the context of the FEM (Collop et al. 2003). The influence of the GK contact parameters on the dynamic response is also investigated. Moreover, the aggregate-to-mastic contacts are analysed under elastic and viscoelastic perspectives to assess which contact approach leads to a better agreement with known experimental data. Finally, the results predicted with the GK contact models are compared with those obtained using the traditionally adopted Burgers contact model.

3. Generalised Kelvin viscoelastic contact model

With the advance in computational resources over the years, modelling the behaviour of complex materials, such as asphalt materials, has become more accessible. Nevertheless, the precision of contact models describing these composites differs. Therefore, a generalised Kelvin contact model, following two different methodologies (GK${}_1$ and GK${}_2$), is developed to better represent the viscoelasticity in asphalt materials.

The GK contact model comprises a spring, a dashpot, and j-chains of Kelvin elements (springs and dashpots connected in parallel) placed in series, representing an elastic portion, a viscoelastic portion, and a viscoplastic portion, as illustrated in Figure 1. The spring and dashpot elements placed in series comprise a Maxwell unit.

Figure 1. (a) DEM ball-to-ball contact and (b) Generalised Kelvin contact model.

The governing equation of the contact model's displacement (u) results from the sum of the elastic displacement of the free spring ($u_{el}$), the irreversible displacement of the dashpot unit ($u_{vp}$), and the summation of the delayed elastic displacement from the Kelvin elements ($u_{ve}^i$). On the other hand, the contact model's force (f), the free spring element's force ($f_{el}$), the dashpot element's force ($f_{vp}$), and the Kelvin chain elements' force ($f_{ve}^i$) are equal to one another.

3.1. Model GK${}_1$

In this approach, the constitutive equations are integrated using a time-centred difference scheme, similar to the adopted for the Burgers model in the study of margarines (Shimizu et al. 2013). The force-displacement relationships of the composing elements of the contact model are expressed by: (1) $\{\begin{array}{l}{\displaystyle u_{el}=\frac{f}{{\kappa }_m}}\\ {\displaystyle {\dot{u}}_{vp}=\frac{f}{{\eta }_m}}\\ {\displaystyle {\dot{u}}_{ve}^i=\frac{f-{\kappa }_i{u}_{ve}^i}{{\eta }_i}}\end{array}$(1) Taking into consideration the force-displacement relationship of the viscoelastic portion in Equation (1(1) $\{\begin{array}{l}{\displaystyle u_{el}=\frac{f}{{\kappa }_m}}\\ {\displaystyle {\dot{u}}_{vp}=\frac{f}{{\eta }_m}}\\ {\displaystyle {\dot{u}}_{ve}^i=\frac{f-{\kappa }_i{u}_{ve}^i}{{\eta }_i}}\end{array}$(1) ), applying the time-centred difference approximation for the derivative element ${\dot{u}}_{ve}^i$ and average values for the force f and displacement $u_{ve}^i$ elements, the displacement corresponding to the ith-element (${}^i{u}_{ve}^{t+\mathrm{\Delta }t}$) can be calculated: (2) ${}^i{u}_{ve}^{t+\mathrm{\Delta }t}=\frac{1}{A_i}[{}^i{u}_{ve}^t{B}_i+\frac{\mathrm{\Delta }t(f^{t+\mathrm{\Delta }t}+f^t)}{2{\eta }_i}]$(2) where $\mathrm{\Delta }t$ is the time step, $f^{t+\mathrm{\Delta }t}$ is the contact force of the model, and $f^t$ and ${}^i{u}_{ve}^t$ are the contact force and viscoelastic displacement at the previous step. The constants $A_i$ and $B_i$ are given by: (3) $\{\begin{array}{l}{\displaystyle A_i=1+\frac{{\kappa }_i\mathrm{\Delta }t}{2{\eta }_i}}\\ {\displaystyle B_i=1-\frac{{\kappa }_i\mathrm{\Delta }t}{2{\eta }_i}}\end{array}$(3) The displacement of the Maxwell unit ($u_m$) corresponds to the sum of the elastic ($u_{el}$) and viscoplastic ($u_{vp}$) displacements of the generalised model. Accordingly, its first derivative versus time is defined in the following expression: (4) ${\dot{u}}_m={\dot{u}}_{el}+{\dot{u}}_{vp}$(4) Calculating the first derivative for the elastic component and taking the first derivative for the viscoplastic component, Equation (4(4) ${\dot{u}}_m={\dot{u}}_{el}+{\dot{u}}_{vp}$(4) ) can be simplified to: (5) ${\dot{u}}_m=\frac{\dot{f}}{{\kappa }_m}+\frac{f}{{\eta }_m}$(5) Applying the time-centred difference approximation for the time derivatives ${\dot{u}}_m$ and $\dot{f}$, calculating the average value for f and rearranging the terms, the displacement of the Maxwell component is defined as: (6) $u_m^{t+\mathrm{\Delta }t}=\frac{f^{t+\mathrm{\Delta }t}-f^t}{{\kappa }_m}+\frac{\mathrm{\Delta }t(f^{t+\mathrm{\Delta }t}+f^t)}{2{\eta }_m}+u_m^t$(6) The result of calculating the first derivative for the total displacement of the contact model and applying the time-centred difference approximation for the time derivatives is expressed as follows: (7) $u^{t+\mathrm{\Delta }t}-u^t=u_m^{t+\mathrm{\Delta }t}-u_m^t+\sum _{i=1}^j({}^i{u}_{ve}^{t+\mathrm{\Delta }t}-{}^i{u}_{ve}^t)$(7) Therefore, the contact force $f^{t+\mathrm{\Delta }t}$ can be determined after substituting Equations (2(2) ${}^i{u}_{ve}^{t+\mathrm{\Delta }t}=\frac{1}{A_i}[{}^i{u}_{ve}^t{B}_i+\frac{\mathrm{\Delta }t(f^{t+\mathrm{\Delta }t}+f^t)}{2{\eta }_i}]$(2) ) and (6(6) $u_m^{t+\mathrm{\Delta }t}=\frac{f^{t+\mathrm{\Delta }t}-f^t}{{\kappa }_m}+\frac{\mathrm{\Delta }t(f^{t+\mathrm{\Delta }t}+f^t)}{2{\eta }_m}+u_m^t$(6) ) into Equation (7(7) $u^{t+\mathrm{\Delta }t}-u^t=u_m^{t+\mathrm{\Delta }t}-u_m^t+\sum _{i=1}^j({}^i{u}_{ve}^{t+\mathrm{\Delta }t}-{}^i{u}_{ve}^t)$(7) ): (8) $f^{t+\mathrm{\Delta }t}=\frac{1}{C}[u^{t+\mathrm{\Delta }t}-u^t+\sum _{i=1}^j({}^i{u}_{ve}^t-\frac{B_i{}^i{u}_{ve}^t}{A_i})-D{f}^t]$(8) where constants C and D are subsequently expressed as: (9) $\{\begin{array}{cc}& {\displaystyle C=\sum _{i=1}^j\left(\frac{\mathrm{\Delta }t}{2A_i{\eta }_i}\right)+\frac{1}{{\kappa }_m}+\frac{\mathrm{\Delta }t}{2{\eta }_m}}\\ & {\displaystyle D=\sum _{i=1}^j\left(\frac{\mathrm{\Delta }t}{2A_i{\eta }_i}\right)-\frac{1}{{\kappa }_m}+\frac{\mathrm{\Delta }t}{2{\eta }_m}}\end{array}$(9)

3.2. Model GK${}_2$

In this approach, the generalised Kelvin contact model is based on an incremental formulation, in which the increment in displacement results in an increment in the contact force of the model, following closely the finite element formulation proposed by Collop et al. (2003). For this, the delayed elastic displacement and the irreversible displacement derive from the relationship between the force and the rate of creep compliance (J), which is the ratio between the respectively calculated displacement and the applied force in the numerical simulations. The equation governing the total displacement is defined with the use of an integral hereditary formulation that admits an initial value followed by an integration over time: (10) $u_{\xi }\left(t\right)=J_{\xi }\left(0\right)f\left(t\right)+{\int }_0^t\frac{\mathrm{d}{J}_{\xi }(t-t^{\prime })}{\mathrm{d}(t-t^{\prime })}f\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }$(10) where subscript ξ refers to ve and vp representing the viscoelastic and viscoplastic elements of the generalised model, respectively, and $t^{\prime }$ is an integration variable. The initial creep compliances for the viscoelastic and viscoplastic portions assume zero.

Therefore, it is possible to find the increment in the displacement of the viscoelastic and viscoplastic elements of the contact model by deriving Equation (10(10) $u_{\xi }\left(t\right)=J_{\xi }\left(0\right)f\left(t\right)+{\int }_0^t\frac{\mathrm{d}{J}_{\xi }(t-t^{\prime })}{\mathrm{d}(t-t^{\prime })}f\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }$(10) ), and consequently, the displacement governing the two components at time $t+\mathrm{\Delta }t$ is defined. On the other hand, the displacement corresponding to the elastic portion derives from the direct relationship between the contact force and the stiffness ${\kappa }_m$. The respective elastic, viscoelastic, and viscoplastic displacements can be expressed by: (11) $\{\begin{array}{cc}& u_{el}^{t+\mathrm{\Delta }t}=u_{el}^t+C_{el}{f}^t\\ & {\displaystyle {}^i{u}_{ve}^{t+\mathrm{\Delta }t}={}^i{u}_{ve}^t+{}^i{u}_{ve}^t({\mathrm{e}}^{(-\mathrm{\Delta }t/{\tau }_i)}-1)}\\ & +\mathrm{\Delta }t{\mathrm{e}}^{(-\mathrm{\Delta }t/2{\tau }_i)}{C}_{ve}^i(f^t+\frac{\mathrm{\Delta }f}{2})\\ & {\displaystyle u_{vp}^{t+\mathrm{\Delta }t}=u_{vp}^t+\mathrm{\Delta }t{C}_{vp}(f^t+\frac{\mathrm{\Delta }f}{2})}\end{array}$(11) where $C_{el}$, $C_{ve}^i$, and $C_{vp}$ correspond to the inverse of contact parameters ${\kappa }_m$, ${\eta }_i$, and ${\eta }_m$, respectively, and ${\tau }_i$ corresponds to the retardation time defined as the relationship between viscosity and stiffness for the ith-element of the Kelvin chain.

Accordingly, the total displacement increment for the GK${}_2$ contact model can be expressed by: (12) $\mathrm{\Delta }u=u_{el}^{t+\mathrm{\Delta }t}-u_{el}^t+\sum _{i=1}^j({}^i{u}_{ve}^{t+\mathrm{\Delta }t}-{}^i{u}_{ve}^t)+u_{vp}^{t+\mathrm{\Delta }t}-u_{vp}^t$(12) Substituting the elastic, delayed elastic, and viscoplastic displacement components from Equation (11(11) $\{\begin{array}{cc}& u_{el}^{t+\mathrm{\Delta }t}=u_{el}^t+C_{el}{f}^t\\ & {\displaystyle {}^i{u}_{ve}^{t+\mathrm{\Delta }t}={}^i{u}_{ve}^t+{}^i{u}_{ve}^t({\mathrm{e}}^{(-\mathrm{\Delta }t/{\tau }_i)}-1)}\\ & +\mathrm{\Delta }t{\mathrm{e}}^{(-\mathrm{\Delta }t/2{\tau }_i)}{C}_{ve}^i(f^t+\frac{\mathrm{\Delta }f}{2})\\ & {\displaystyle u_{vp}^{t+\mathrm{\Delta }t}=u_{vp}^t+\mathrm{\Delta }t{C}_{vp}(f^t+\frac{\mathrm{\Delta }f}{2})}\end{array}$(11) ) into Equation (12(12) $\mathrm{\Delta }u=u_{el}^{t+\mathrm{\Delta }t}-u_{el}^t+\sum _{i=1}^j({}^i{u}_{ve}^{t+\mathrm{\Delta }t}-{}^i{u}_{ve}^t)+u_{vp}^{t+\mathrm{\Delta }t}-u_{vp}^t$(12) ) and rearranging the terms, the displacement increment is determined. However, the resultant equation can be simplified into a pseudo-elastic expression relating the increment in the displacement to the increment in the contact force, as described in the following: (13) $\mathrm{\Delta }f=C_t^{-1}\mathrm{\Delta }\overline{u}$(13) where $C_t$, the total matrix of flexibility of the generalised model, is expressed by: (14) $C_t=C_{el}+\frac{\mathrm{\Delta }t}{2}{C}_1$(14) Furthermore, the matrix of flexibility of the time-dependent portion of the model, corresponding to the delayed elastic ($C_{ve}^i$) and viscoplastic ($C_{vp}$) elements, can be defined as: (15) $C_1=C_{vp}+\sum _{i=1}^j\left[C_{ve}^i{\mathrm{e}}^{-\left(\mathrm{\Delta }t/2{\tau }_i\right)}\right]$(15) Accordingly, the pseudo-elastic displacement increment $\left(\mathrm{\Delta }\overline{u}\right)$ for the normal and tangential directions is: (16) $\mathrm{\Delta }\overline{u}=\mathrm{\Delta }{u}_1-\mathrm{\Delta }t{C}_1{f}^t$(16) where the displacement increment $\mathrm{\Delta }{u}_1$ can be expressed by: (17) $\mathrm{\Delta }{u}_1=\mathrm{\Delta }u-\sum _{i=1}^j\left[{}^i{u}_{ve}^t({\mathrm{e}}^{(-\mathrm{\Delta }t/{\tau }_i)}-1)\right]$(17) Finally, the contact force $f^{t+\mathrm{\Delta }t}$ results from the sum of the contact force at time t and the increment computed in Equation (13(13) $\mathrm{\Delta }f=C_t^{-1}\mathrm{\Delta }\overline{u}$(13) ).

4. Experimental data

The experimental results reported by Silva (2005), who conducted uniaxial tension-compression cyclic tests on prismatic asphalt mastic and asphalt mixture specimens to investigate the effect of the mastic behaviour on the mechanical properties of asphalt mixtures, are used for the calibration and validation of the proposed DEM model. The experiments were selected because the adopted test (load configuration, amplitude, and frequencies) and specimen geometry for the asphalt mixtures and mastics are the same, and each mastic is representative of an asphalt mixture's type. In the study, four asphalt mixtures (AM${}_1$, AM${}_2$, AM${}_3$, and AM${}_4$) with different compositions were defined, and four asphalt mastics (M${}_1$, M${}_2$, M${}_3$, and M${}_4$) were designed in terms of aggregate gradation and bitumen content, aiming to represent the true mastic in these asphalt mixtures. Therefore, each mastic specimen corresponded to an asphalt mixture sample.

The reference mixture AM${}_1$ is a dense-graded (0/14 mm) mixture for surface courses, produced with paving grade bitumen 35/50, granite aggregates, and limestone filler. The optimum binder content was defined with the Marshall method. A paving grade bitumen 50/70 was used for the mixture AM${}_2$ to analyse the impact of a softer bitumen on the mechanical properties. The aggregate gradation passing the #10 sieve was made finer in mixture AM${}_3$ to evaluate the impact of a finer mastic gradation. Finally, the impact of the type of filler was studied in mixture AM${}_4$, where granite filler (baghouse dust collected in plant) was adopted in the composition. The four mastics were designed with a 0/2 mm aggregate gradation and a bitumen content adjusted to the specific surface of the aggregate following a specific method described in Silva (2005). Table 1 shows the material composition of asphalt mastics and asphalt mixtures.

Table 1. Aggregate gradation for asphalt mastics and asphalt mixtures.

	Aggregate size (percentage passing)
	Asphalt mastic			Asphalt mixture
Sieve size	M${}_1$/M${}_2$	M${}_3$	M${}_4$	AM${}_1$/AM${}_2$	AM${}_3$	AM${}_4$
#3/4”	–	–	–	100.0	100.0	100.0
#1/2”	–	–	–	83.7	83.7	83.7
#4	–	–	–	47.8	47.8	47.8
#10	100.0	100.0	100.0	32.5	32.5	32.5
#20	73.8	83.0	73.8	–	–	–
#40	52.1	68.8	52.1	16.9	22.4	16.9
#80	28.9	40.9	28.9	9.4	13.3	9.4
#200	16.2	23.3	16.2	5.3	7.6	5.3
Bitumen content (%)	15.9	19.9	15.9	5.2	6.5	5.2

The aggregates were preheated in an oven for 4 h at $170{}^{\circ }$C and the bitumen for 1 h at $155{}^{\circ }$C to fabricate the asphalt mixtures and mastics. The materials were mixed in a standard mixer and then transferred to a rectangular mould to be roller compacted. Finally, after cooling to room temperature and de-moulding, the slabs were sawn to obtain 50×50×80 mm${}^3$ prismatic specimens.

The specimens were subjected to uniaxial tension-compression sinusoidal dynamic loading with an imposed strain amplitude of $1\times {10}^{-4}$ m/m to characterise the stiffness properties of the materials. The dynamic stiffness modulus and phase angle were measured across 1, 2, 5, and 10 Hz loading frequencies (from the highest to the lowest frequency) at a testing temperature of 20 ${}^{\circ }$C. The prismatic specimen was positioned standing on the smaller side. The bottom and upper surfaces of the specimen were glued to the plates of the testing machine to allow the development of tension stresses. The same test conditions were applied for both asphalt materials. The test configuration is illustrated in Figure 2. The test was repeated with a minimum of four replicates for each material. Tables 2 and 3 show the dynamic response (average values) of the investigated composites. Further details on the sample preparation and experimental procedure techniques are available in Silva (2005).

Figure 2. Experimental analysis: (a) test configuration and (b) sinusoidal signal.

Table 2. Experimental results of asphalt mastics.

	M${}_1$		M${}_2$		M${}_3$		M${}_4$
Frequency (Hz)	$|$E*$|$ (MPa)	ϕ (degrees)	$|$E*$|$	ϕ	$|$E*$|$	ϕ	$|$E*$|$	ϕ
1	1825.4	35.7	1309.4	33.9	952.7	32.6	1503.3	35.9
2	2334.8	33.7	1613.3	34.6	1203.9	36.3	1914.3	32.3
5	3335.6	33.3	2308.6	35.0	1534.4	40.6	2776.7	31.8
10	4606.1	26.8	3187.5	30.6	2257.7	36.2	3820.4	26.5

Table 3. Experimental results of asphalt mixtures.

	AM${}_1$		AM${}_2$		AM${}_3$		AM${}_4$
Frequency (Hz)	$|$E*$|$ (MPa)	ϕ (degrees)	$|$E*$|$	ϕ	$|$E*$|$	ϕ	$|$E*$|$	ϕ
1	6870.8	25.5	5697.4	38.1	–	–	6062.6	30.7
2	8741.2	21.2	7376.9	28.5	7333.3	28.7	7112.6	25.1
5	11062.6	17.5	9905.2	19.4	10328.2	20.2	9413.0	15.5
10	12603.2	14.6	12231.7	17.0	12299.7	17.5	11373.9	13.3

5. Asphalt mastic modelling

5.1. Calibration of input parameters

5.1.1. Macroscopic numerical formulation

The DEM contact model parameters cannot be measured in laboratory tests. For this reason, a fitting procedure needs to be adopted to determine these contact parameters. This fitting procedure is more difficult for generalised contact models because more unknown variables exist. In order to ease the calibration process, an analytical solution is first established to fit experimental results of dynamic modulus and phase angle of mastics based on the macroscopic generalised Kelvin model (Figure 3). The set of parameters is calibrated for each mastic sample.

Figure 3. Macroscopic generalised Kelvin viscoelastic model representation.

The stress-strain relationships of the macroscopic model for the elastic, viscoplastic, and delayed elastic portions are expressed through the following: (18) $\{\begin{array}{l}{\sigma }_{el}=E_m{ϵ}_{el}\\ {\sigma }_{vp}=C_m{\dot{ϵ}}_{vp}\\ {\sigma }_{ve}^i=E_i{ϵ}_{ve}^i+C_i{\dot{ϵ}}_{ve}^i\end{array}$(18) where ${\sigma }_{el}$, ${\sigma }_{vp}$, and ${\sigma }_{ve}^i$ are the stress corresponding to the elastic, viscoplastic, and viscoelastic elements, respectively.

The applied sinusoidal loading and the corresponding dynamic strain are: (19) $\sigma \left(t\right)={\sigma }_0{\mathrm{e}}^{i\omega t}$(19) (20) $ϵ\left(t\right)={ϵ}_0{\mathrm{e}}^{i(\omega t-\varphi )}$(20) where ${\sigma }_0$ and ${ϵ}_0$ are the stress and strain amplitudes, ω is the angular frequency, and ϕ is the phase angle.

The result after substituting Equations (19(19) $\sigma \left(t\right)={\sigma }_0{\mathrm{e}}^{i\omega t}$(19) ) and (20(20) $ϵ\left(t\right)={ϵ}_0{\mathrm{e}}^{i(\omega t-\varphi )}$(20) ) into Equation (18(18) $\{\begin{array}{l}{\sigma }_{el}=E_m{ϵ}_{el}\\ {\sigma }_{vp}=C_m{\dot{ϵ}}_{vp}\\ {\sigma }_{ve}^i=E_i{ϵ}_{ve}^i+C_i{\dot{ϵ}}_{ve}^i\end{array}$(18) ) is: (21) $\{\begin{array}{l}{\sigma }_0{\mathrm{e}}^{i\omega t}=E_m{ϵ}_{el}{\mathrm{e}}^{i(\omega t-\varphi )}\\ {\sigma }_0{\mathrm{e}}^{i\omega t}=i\omega C_m{ϵ}_{vp}{\mathrm{e}}^{i(\omega t-\varphi )}\\ {\sigma }_0^i{\mathrm{e}}^{i\omega t}=E_i{ϵ}_{ve}^i{\mathrm{e}}^{i(\omega t-\varphi )}+i\omega C_i{ϵ}_{ve}^i{\mathrm{e}}^{i(\omega t-\varphi )}\end{array}$(21) The complex modulus ($E^{\ast }$) can be calculated by dividing the applied stress (Equation (19(19) $\sigma \left(t\right)={\sigma }_0{\mathrm{e}}^{i\omega t}$(19) )) by the resultant strain (Equation (20(20) $ϵ\left(t\right)={ϵ}_0{\mathrm{e}}^{i(\omega t-\varphi )}$(20) )). It accounts for a real component ($E^{\prime }$), corresponding to the storage modulus, and an imaginary component ($E^{″}$), corresponding to the loss modulus in a loading cycle. Therefore, it can be expressed as: (22) $E^{\ast }=E^{\prime }+i{E}^{″}=\frac{\sigma \left(t\right)}{ϵ\left(t\right)}=\frac{{\sigma }_0}{{ϵ}_0}{\mathrm{e}}^{i\varphi }=(\frac{{\sigma }_0}{{ϵ}_{el}}+\frac{{\sigma }_0}{{ϵ}_{vp}}+\sum _{i=1}^j\frac{{\sigma }_0^i}{{ϵ}_{ve}^i}){\mathrm{e}}^{i\varphi }$(22) Thus, the complex modulus is determined by substituting Equation (21(21) $\{\begin{array}{l}{\sigma }_0{\mathrm{e}}^{i\omega t}=E_m{ϵ}_{el}{\mathrm{e}}^{i(\omega t-\varphi )}\\ {\sigma }_0{\mathrm{e}}^{i\omega t}=i\omega C_m{ϵ}_{vp}{\mathrm{e}}^{i(\omega t-\varphi )}\\ {\sigma }_0^i{\mathrm{e}}^{i\omega t}=E_i{ϵ}_{ve}^i{\mathrm{e}}^{i(\omega t-\varphi )}+i\omega C_i{ϵ}_{ve}^i{\mathrm{e}}^{i(\omega t-\varphi )}\end{array}$(21) ) into Equation (22(22) $E^{\ast }=E^{\prime }+i{E}^{″}=\frac{\sigma \left(t\right)}{ϵ\left(t\right)}=\frac{{\sigma }_0}{{ϵ}_0}{\mathrm{e}}^{i\varphi }=(\frac{{\sigma }_0}{{ϵ}_{el}}+\frac{{\sigma }_0}{{ϵ}_{vp}}+\sum _{i=1}^j\frac{{\sigma }_0^i}{{ϵ}_{ve}^i}){\mathrm{e}}^{i\varphi }$(22) ). Considering its reciprocal value, the complex compliance ($D^{\ast }$) is given by: (23) $D^{\ast }=D^{\prime }+i{D}^{″}=\frac{1}{E_m}+\frac{1}{i\omega C_m}+\sum _{i=1}^j\left(\frac{1}{E_i+i\omega C_i}\right)$(23) Moreover, the real and imaginary parts of the complex compliance can be measured as: (24) $\{\begin{array}{l}{\displaystyle D^{\prime }=\frac{1}{E_m}+\sum _{i=1}^j\frac{E_i}{E_i^2+{\omega }^2{C_i}^2}}\\ {\displaystyle D^{″}=\frac{1}{\omega C_m}+\sum _{i=1}^j\frac{\omega C_i}{E_i^2+{\omega }^2{C_i}^2}}\end{array}$(24) Furthermore, the norms of the complex modulus and the phase angle are expressed as: (25) $\{\begin{array}{l}{\displaystyle \left|E^{\ast }\right|=\frac{1}{\sqrt{D^{\prime 2}+D^{″2}}}}\\ {\displaystyle \varphi ={\tan }^{-1}\left(\frac{D^{″}}{D^{\prime }}\right)}\end{array}$(25) The fitting process is based on the minimisation function (F) relating the predicted values of the real and imaginary components of the dynamic modulus and the results obtained in laboratory tests, defined as: (26) $F=\sum _{z=1}^n[{(\frac{E_z^{\prime }}{E_z^{{}^{\prime }0}}-1)}^2+{(\frac{E_z^{″}}{E_z^{{}^{″}0}}-1)}^2]$(26) where $E_z^{\prime }$ and $E_z^{″}$ are the real and imaginary values of dynamic modulus predicted at the zth frequency, $E_z^{{}^{\prime }0}$ and $E_z^{{}^{″}0}$ are the laboratory values of the real and imaginary components of the dynamic modulus at the zth frequency, and n is the data points.

5.1.2. Resultant calibration of input parameters

The macroscopic parameters for the GK model and the Burgers model are calibrated for loading frequencies of 1, 2, 5, and 10 Hz according to the equation group (24(24) $\{\begin{array}{l}{\displaystyle D^{\prime }=\frac{1}{E_m}+\sum _{i=1}^j\frac{E_i}{E_i^2+{\omega }^2{C_i}^2}}\\ {\displaystyle D^{″}=\frac{1}{\omega C_m}+\sum _{i=1}^j\frac{\omega C_i}{E_i^2+{\omega }^2{C_i}^2}}\end{array}$(24) )–(26(26) $F=\sum _{z=1}^n[{(\frac{E_z^{\prime }}{E_z^{{}^{\prime }0}}-1)}^2+{(\frac{E_z^{″}}{E_z^{{}^{″}0}}-1)}^2]$(26) ), taking as reference the rheological properties of mastics shown in Table 2. The number of Kelvin elements for the generalised model is selected according to the analysis of the error predictions in the previous equations. In this study, three Kelvin elements are sufficient to achieve the least error. In other asphalt studies that adopted generalised models, such as Dai (2010), Ren and Sun (2016), and Sun (2018), the number of elements varied from two to four. The set of macroscopic parameters for asphalt mastics M${}_1$, M${}_2$, M${}_3$, and M${}_4$ are shown in Tables 4 and 5. The influences of the different characteristics among the mastics are indirectly described in the different parameters.

Table 4. Fitted parameters for macroscopic generalised Kelvin model.

Mastic	M${}_1$	M${}_2$	M${}_3$	M${}_4$
E${}_1$ (kN/m)	$4.76\times {10}^6$	$5.60\times {10}^6$	$2.56\times {10}^6$	$3.76\times {10}^6$
E${}_2$	$3.03\times {10}^6$	$3.33\times {10}^6$	$2.67\times {10}^6$	$7.23\times {10}^7$
E${}_3$	$5.93\times {10}^5$	$3.04\times {10}^7$	$2.04\times {10}^7$	$8.50\times {10}^5$
C${}_1$ (kN·s/m)	$2.02\times {10}^5$	$1.50\times {10}^5$	$6.76\times {10}^4$	$1.60\times {10}^5$
C${}_2$	$8.51\times {10}^5$	$3.09\times {10}^5$	$3.08\times {10}^5$	$3.94\times {10}^7$
C${}_3$	$2.29\times {10}^7$	$2.82\times {10}^6$	$2.36\times {10}^6$	$4.75\times {10}^5$
E${}_{\mathrm{m}}$	$5.93\times {10}^6$	$4.69\times {10}^6$	$4.26\times {10}^6$	$4.88\times {10}^6$
C${}_{\mathrm{m}}$	$1.34\times {10}^6$	$6.13\times {10}^5$	$4.54\times {10}^5$	$2.58\times {10}^9$

Table 5. Fitted parameters for macroscopic Burgers model.

Mastic	M${}_1$	M${}_2$	M${}_3$	M${}_4$
E${}_1$ (kN/m)	$3.44\times {10}^6$	$2.31\times {10}^6$	$2.86\times {10}^6$	$1.56\times {10}^6$
C${}_1$ (kN·s/m)	$1.93\times {10}^5$	$1.20\times {10}^5$	$1.57\times {10}^5$	$6.85\times {10}^4$
E${}_{\mathrm{m}}$	$5.73\times {10}^6$	$4.34\times {10}^6$	$4.69\times {10}^6$	$3.70\times {10}^6$
C${}_{\mathrm{m}}$	$7.25\times {10}^5$	$5.25\times {10}^5$	$6.17\times {10}^5$	$3.61\times {10}^5$

Figure 4 illustrates the analytical result of the dynamic modulus and phase angle for the group of mastics. As shown, a good agreement between the fitting and laboratory data is achieved for both macroscopic models. Nevertheless, the generalised Kelvin model is able to generate a better calibration prediction in comparison with the response predicted with a Burgers model for both rheological properties. Therefore, similar correspondences are expected to occur in the numerical simulations.

Figure 4. Calibration of macroscopic parameters for the generalised Kelvin model and Burgers model for asphalt mastics: (a) M${}_1$, (b) M${}_2$, (c) M${}_3$, and (d) M${}_4$.

Subsequently, the macro-scale parameters listed in Tables 4 and 5 must be converted into contact parameters before their use in the DEM simulations. A viscoelastic beam composed of two particles in contact is used to relate the macroscopic response to the contact model parameters based on stress and strain relationships for the normal direction, as given by: (27) $\{\begin{array}{l}{\displaystyle {\kappa }_{\xi }^n=\frac{E_{\xi }A}{L}\delta }\\ {\displaystyle {\eta }_{\xi }^n=\frac{C_{\xi }A}{L}\delta }\end{array}$(27) where L is the sum of the two neighbouring particles' radii, A is the cross-sectional area of the spheres, δ is a coefficient of adjustment between the macroscopic and the DEM parameters related to the particle assembly and contact properties, and the subscript ξ assumes i for the ith viscoelastic element in the Kelvin chain and m for the elastic and viscoplastic elements of the model. The parameters for the tangent direction are obtained by the product between the contact stiffness ratio α (${\kappa }_s/{\kappa }_n$) and the values estimated for the normal direction.

5.2. Numerical specimen preparation and modelling results

The generation particle technique used in this study adopts a random distribution of particles in the domain and the Laguerre–Voronoi diagrams of the spherical particles to establish the particle contacts and the contact areas (Monteiro Azevedo et al. 2015).

After specifying the dimensions of the virtual specimen, the maximum and minimum particle diameters, the radius distribution, the porosity, and the density of particles, a uniform probability distribution is imposed to generate the centres of gravity of the particles. The particles are usually inserted with half of their final radius in order to speed up the generation process. Once the particles placement is complete, a radius expansion procedure is set to enlarge the particles' dimension to the real desired value. In the next step, a DEM cohesionless algorithm is carried out in order to reach a stable system, where the initial particle overlaps due to the radius enlargement procedure are greatly reduced and evenly distributed in the particle assembly.

Later, the centres of gravity of the created particles are triangulated using a 3D weighted Delaunay scheme, followed by the construction of the dual 3D Voronoi diagram (Monteiro Azevedo and Lemos 2013) and consequently, the polygonal-shaped particles are created. The contact area and the corresponding contact location are determined by the Voronoi tessellation. The contact area corresponds to the area of the associated Voronoi cell facet. The contact location is also defined by the Voronoi cell facet. In this particle generation scheme, particles are considered to interact with their neighbouring ones through a common polygonal facet and, when compared with traditional generation approaches, increases the contacts per particle. A similar procedure has been used for 2D and 3D simulations of rock fracture (Monteiro Azevedo et al. 2015, Monteiro Azevedo and Lemos, 2013), where it is shown that the adopted particle procedure enhances the particle model predictions in biaxial and triaxial tests when compared with traditional particle models.

Accordingly, a three-dimensional particle assembly is generated to represent the mastics M${}_1$, M${}_2$, M${}_3$, and M${}_4$ in the numerical simulations, as shown in Figure 5. The virtual sample, with dimensions of $4\times 2.5\times 2.5$ cm${}^3$, comprises 8873 randomly distributed particles with diameters ranging from 1.5 mm to 2.0 mm with a total of 53080 interactions. The mastic particles represent the fine aggregates ($<2.0$ mm) mixed with bitumen. The assembly is designed with half of the laboratory specimen's dimensions in order to reduce the computational cost and accelerate the calibration procedure of the 3D DEM mastic model.

Figure 5. Asphalt mastic sample: (a) experimental and (b) virtual.

In the numerical mastic sample, two contact types exist: contacts within mastic particles and mastic-to-wall contacts. For the purpose of describing the different contacts, a linear elastic model and the developed generalised Kelvin contact model are defined in the DEM model. The elastic model is adopted to describe the mastic-to-wall contacts, while the contacts within the mastic are represented by the viscoelastic model.

In order to determine the corresponding coefficient of adjustment δ for each mastic and calibrate the contact parameters for the input in the DEM model, laboratory dynamic loading tests are considered for numerical simulations. The same loading conditions adopted in the laboratory are applied. An imposed strain amplitude of $1\times {10}^{-4}$ m/m is adopted. The simulations are carried out for testing frequencies of 1, 2, 5, and 10 Hz. The load is applied to the upper wall in the vertical direction, while the bottom wall remains fixed in all directions. Based on the applied strain and stress response for each loading frequency, the dynamic modulus and phase angle can be obtained according to: (28) $\left|{\mathrm{E}}^{\ast }\right|=\frac{{\sigma }_{max}-{\sigma }_{min}}{{ϵ}_{max}-{ϵ}_{min}}$(28) (29) $\varphi =\frac{\mathrm{\Delta }t}{T}\times 360$(29) where ${\sigma }_{max}$ and ${\sigma }_{min}$, refer to the maximum and minimum stress in a loading cycle, ${ϵ}_{max}$ and ${ϵ}_{min}$ refer to the maximum and minimum strain in a loading cycle, $\mathrm{\Delta }t$ is the time difference between two adjacent peak stress and strain, and T is the loading period.

The predicted response of the dynamic modulus and phase angle for each asphalt mastic is compared to the laboratory testing results in Figure 6. A good correspondence to laboratory testing results is obtained with the coefficients of adjustment δ of 1.88, 1.90, 1.86, and 1.90 for specimens M${}_1$, M${}_2$, M${}_3$, and M${}_4$, respectively. The contact stiffness ratio α was set to 0.10. The response obtained with the two proposed generalised Kelvin contact models (GK${}_1$ and GK${}_2$) is as expected similar (difference in modulus less than ${10}^{-6}$ MPa). Additionally, the result achieved with the GK contact model is compared to the response obtained with a Burgers model and the relative errors between the predicted and experimental dynamic behaviours are measured.

Figure 6. Laboratory data and simulation results with the use of the generalised Kelvin contact model and Burgers contact model of dynamic modulus and phase angle for mastics: (a) M${}_1$, (b) M${}_2$, (c) M${}_3$, and (d) M${}_4$.

The GK contact model provided the best agreement with laboratory data when compared to the results obtained with a Burgers contact model. This is evidenced in Table 6, where the average predicted errors in the numerical simulations for both models are presented.

Table 6. Average errors between numerical simulations and laboratory data for dynamic modulus and phase angle.

	Dynamic modulus (%)		Phase angle (%)
Mastic	GK model	Burgers model	GK model	Burgers model
M${}_1$	1.27	3.65	2.59	2.54
M${}_2$	0.98	2.19	2.72	3.46
M${}_3$	3.51	4.02	5.67	6.64
M${}_4$	3.84	4.00	3.59	4.06

The average errors considering the group of mastics for the dynamic modulus in comparison with the experimental values are $2.40\text{\%}$ and $3.46\text{\%}$ for the GK contact model and Burgers contact model, respectively. In addition, the maximum errors in the simulations are in favour of the GK contact model, which for samples M${}_1$, M${}_2$, M${}_3$, and M${}_4$ are 3.83%, 1.92%, 5.69%, and 6.51%, respectively, while higher errors of 4.89%, 4.18%, 9.92%, and 6.82% are obtained with a Burgers contact model. Regarding the precision of the contact models on the testing frequencies, an improvement is noticed when using the GK contact model in most cases. The average errors are 2.79%, 2.78%, 1.86%, and 2.17% for the GK contact model, unlike 4.30%, 5.36%, 1.29%, and 2.91% obtained with a Burgers contact model at 1, 2, 5, and 10 Hz, respectively.

Similarly to the dynamic modulus, the phase angle is better predicted when adopting the GK contact model, which relative errors are on average $3.64\text{\%}$ (GK model) and $4.17\text{\%}$ (Burgers model). The maximum errors are in favour of the GK model for mastics M${}_1$, M${}_2$, and M${}_4$, with percentages of 4.77%, 4.97%, and 4.66%, while 7.24%, 9.37%, and 6.73% resulted from the use of a Burgers model. The only exception is verified for mastic M${}_3$, which in this case the maximum errors are on average 10.79% and 8.87% for the GK and Burgers contact models, respectively. Furthermore, the precision on most of the testing frequencies is improved when using the GK contact model, resulting in average errors of 2.87%, 2.91%, 5.00%, and 3.79%, in contrast with 3.03%, 3.72%, 2.61%, and 7.34% for the Burgers contact model at 1, 2, 5, and 10 Hz.

Three-dimensional numerical models are usually time-consuming and computationally demanding. Generally, the number of particles and contacts take an important role in the time of simulations. The numerical simulations computational runtimes varied from approximately 1 day to 8 days for loading frequencies of 10 and 1 Hz, respectively. The simulations were carried out with an Intel [email protected] GHz multi-core processor. When compared with the GK${}_1$ contact model, the GK${}_2$ contact model's computational runtimes were on average approximately 14% higher. The GK${}_1$ contact model computational runtimes were similar to the Burgers contact model computational runtimes.

5.3. Parametric study on the influence of contact parameters

In the calibration procedure previously presented, the same coefficient of adjustment δ is assigned to all the contact model parameters. In this section, a parametric study is carried out to evaluate the contribution of each of the elements composing the GK contact model on the main complex stiffness properties, which can significantly ease the DEM model calibration. Mastic sample M${}_1$ is chosen for the parametric analysis.

The following scenarios for the coefficient of adjustment ($\delta =1.88$) previously determined are studied: (1) δ only multiplying the free elastic spring, (2) δ only multiplying the set of all springs in the contact model, (3) δ only multiplying the free viscosity, and (4) δ only multiplying the set of all viscosities in the contact model. The results are compared with the results obtained with a non-calibrated system ($\delta =1.0$ for all properties) and with the results obtained with the calibrated response presented in Figure 6(a) for mastic M${}_1$ ($\delta =1.88$ multiplying all GK contact parameters). The loading conditions are equal to those adopted in Section 5.2. The predicted dynamic modulus and phase angle are shown in Figure 7 for each studied scenario.

Figure 7. Influence of the coefficient of adjustment for the contact parameters of asphalt mastic M${}_1$ for the analysis of the: (a) dynamic modulus and (b) phase angle.

Figure 7(a) shows that in all studied scenarios, with the exception of scenario (3) (δ multiplying viscosity ${\eta }_m$), there is an increase in the predicted dynamic modulus of the mastic when compared to the results obtained with the non-calibrated system. As expected, in scenario (2) (δ multiplying the set of all springs in the contact model) the variation in the predicted dynamic modulus is higher than in scenario (1) (δ only multiplying the free elastic spring). This is due to the fact that the predicted dynamic modulus is proportional to the global stiffness of the adopted GK contact model. Also, in scenarios (1) and (2) the increase in dynamic modulus is proportional to the frequency, which does not occur in scenario (4).

Regarding the predicted phase angle values, Figure 7(b) shows that the calibrated system predicts a phase angle similar to that predicted using a non-calibrated system. Figure 7(b) also shows that when compared with the non-calibrated response, scenarios (1) and (4) are the ones that lead to a higher variation in the predicted phase angle, more pronounced in scenario (4). For the latter, the effect of the frequency on the phase angle is notorious.

The results obtained indicate that a fine tune of the predicted phase angle numerical response can be obtained by adjusting the δ multiplier associated with ${\eta }_m$, as its influence on the dynamic modulus is negligible. The obtained numerical results also indicate that ${\kappa }_m$ is an influential contact parameter when assessing both the dynamic modulus and phase angle and that the set of all viscosities of the contact model has a significant influence on the predicted phase angle.

6. Asphalt mixture modelling

6.1. Numerical specimen preparation

Discrete element simulations of uniaxial dynamic tests following the same laboratory test conditions are undertaken for asphalt mixtures. A micro-structure asphalt mixture is generated using a random generation method. This methodology has the advantages of being laboratory non-dependent and easy to implement, different from the image-based method that requires laboratory analysis of every sample for representation, special equipment, and it is costly and time-consuming (Ding et al. 2019), particularly for 3D models with a considerable amount of particles. Besides, the random generation method can represent the dispersion and heterogeneity of mixtures (Xie et al. 2022). The adopted DEM model assumes aggregates as larger particles and mastics as smaller ones. For these reasons, the numerical modelling of large structures is possible.

The asphalt mixture particle generation, like the mastic generation (see Section 5.2), has two different stages: (1) the generation of the particles in the domain, and (2) building a 3D Voronoi polyhedral structure based on the Laguerre-Voronoi method. In the first stage of the process, the half-diameter generation is only adopted for the particles smaller than 2 mm, which in this study represent the asphalt mastic (mixture of bitumen and aggregate $<2$ mm). The coarse aggregate structure ($>2$ mm) is initially created following a sieve approach which considers the aggregate gradation (see Table 1). The coarse aggregate particles are inserted from the largest to the smallest diameter without any particle overlap. Additionally, the particles representing the aggregate structure are inserted with an enlarged radius (around 5%) to further guarantee that no overlap occurs at this stage. Later, the void space is filled with the particles representing the mastic ($<2$ mm), which are initially inserted with a half-diameter size in order to ease the process of particle insertion. When all the mastic particles are introduced, their radius is enlarged to its normal size. Following, a DEM modelling is performed (pure friction contacts only) in order to reach a stable system, where the initial particle overlaps are greatly reduced and evenly distributed in the particle assembly.

The second stage is similar to the one adopted for the asphalt mastic generation. A 3D Voronoi polyhedral structure is created based on the particle assembly that has been generated, taking into account the particles' centres of gravity and the particles' radius. The contacts of a given particle are set based on the adjacent 3D Voronoi Cells. This approach greatly increases the particle coordination number and directly sets the contact areas.

A prismatic virtual sample with dimensions of $8\times 5\times 5$ cm${}^3$ comprising 27646 particles and 160252 contacts is generated to represent the asphalt mixtures AM${}_1$, AM${}_2$, AM${}_3$, and AM${}_4$ in the DEM simulations. The numerical sample is created with the aggregate gradation based on the laboratory specimen characteristics shown in Table 1.

Table 7 shows the number of particles per aggregate size and volumetric characteristics of the numerical sample. The properties of the asphalt mastic portion derive from the previous sections and, correspondingly, particles are generated with diameters ranging from 1.5 mm to 2.0 mm. Moreover, the aggregate and asphalt mastic fractions, as well as the laboratory and numerical asphalt mixture samples, are shown in Figure 8. In addition, Figure 9 illustrates the internal structure of the virtual mixture.

Figure 8. Asphalt mixture sample: (a) coarse aggregates, (b) asphalt mastic, (c) experimental sample, and (d) virtual sample.

Figure 9. Virtual asphalt mixture specimen and its internal cross sections.

Table 7. Number and the corresponding volume of aggregates generated in the numerical sample.

	Sieve size (mm)
Aggregates	19	12.5	9.5	4.75	2
Number of particles retained on the sieve	–	13	29	309	2388
Volume (mm${}^3$)	–	27794.5	20957.6	39365.8	25788.0

In this study, aggregates are considered a pure elastic material. The aggregate modulus (E${}_{\mathrm{a}}$) is taken as 61 GPa, and the relationship between the elastic modulus and the contact stiffness (${\kappa }_a$) representing the aggregates assumes ${\kappa }_a=2E_a$.

In the asphalt mixture numerical structure, five different contact types exist: aggregate-to-aggregate, mastic-to-mastic, aggregate-to-mastic, aggregate-to-wall, and mastic-to-wall contacts. Two contact model types are chosen to represent the different interactions: an elastic contact model and a viscoelastic contact model (GK and Burgers). The elastic model is represented by a spring and defines the contact between two aggregate particles and the contact between any particle element to the wall. The viscoelastic contact model represents the behaviour within the mastic. In this analysis, only the GK${}_1$ contact model is adopted given that it has lower computational runtimes, but as previously shown, both GK contact models predict the same rheological response. The discussed contact models are illustrated in Figure 10.

Figure 10. Contact models in the asphalt mixture particle assembly: (a) aggregate-to-aggregate, (b) aggregate-to-wall, (c) mastic-to-wall, (d) viscoelastic approach to aggregate-to-mastic, and (e) elastic approach to aggregate-to-mastic.

For the purpose of analysing which model type better represents the contact between coarse aggregates and mastic particles, two possibilities are considered. In the first approach, the stiffness of the aggregate (${\kappa }_a$) is connected in series with the viscoelastic contact model representing the mastic, resulting in a viscoelastic model (Figure 10(d)). In the second approach, ${\kappa }_a$ is connected in series with the stiffness ${\kappa }_m$ of the viscoelastic contact model of the mastic, resulting in an elastic model (Figure 10(e)). In both cases, an equivalent stiffness component is set as the result of the combination in series between ${\kappa }_a$ and the stiffness ${\kappa }_m$ in the viscoelastic model.

6.2. Modelling results

The viscoelastic parameters of the mastics (Section 5.2) and the elastic properties of the aggregates are used as input parameters in the mixture modelling. Also, the same loading frequencies applied in the laboratory measurements are adopted. The stress and strain are monitored to determine the dynamic behaviour of the asphalt mixture numerical model. Similarly to the mastics, a coefficient of adjustment is required to properly calibrate the DEM results to the experimental values. When the aggregate-to-mastic contacts assume an elastic model, the coefficients of adjustment for AM${}_1$, AM${}_2$, AM${}_3$, and AM${}_4$ are 1.10, 1.75, 3.16, and 1.12, respectively. When adopting a viscoelastic approach for the aggregate-to-mastic contacts, coefficients of adjustment of 3.15, 6.75, 16.00, and 3.60 are required.

Experimental analysis showed that asphalt mixtures AM${}_2$ and AM${}_3$ presented lower stiffness due to a higher bitumen content and the use of a softer bitumen in comparison with the reference sample AM${}_1$. As a consequence, this behaviour is reflected in higher values of coefficients of adjustment for these specimens. Moreover, asphalt mixture AM${}_4$ produced similar rheological values to the reference mixture during experiments, and this fact contributed to an approximate value for δ with regard to AM${}_1$.

The predicted results for the DEM simulations are shown in Figure 11. The results predicted with the GK contact model are compared with the results predicted with a Burgers contact model and the experimental results. When comparing the response predicted with the elastic aggregate-to-mastic contact and with the response predicted with a viscoelastic aggregate-to-mastic contact, it is clear that a better agreement with the behaviour observed experimentally for the group of mixtures is obtained with a viscoelastic model. As expected, when compared with the response predicted with a Burgers contact model, a better agreement with the experimental results is obtained with the proposed GK contact model. Ren et al. (2020) reported a similar result when comparing the two viscoelastic contacts in a 2D DEM model of asphalt mixtures.

Figure 11. Numerical results for the rheological behaviour of asphalt mixtures: (a) and (b) AM${}_1$, (c) and (d) AM${}_2$, (e) and (f) AM${}_3$, and (g) and (h) AM${}_4$.

When adopting a viscoelastic aggregate-to-mastic GK contact, the average error regarding the dynamic modulus is reduced in comparison with the response predicted with a viscoelastic aggregate-to-mastic Burgers contact. Error values of 7.00% and 7.92% are predicted with each viscoelastic approach. The obtained average errors are in accordance with the study presented by Ren et al. (2020). Additionally, the maximum average error is registered for sample AM${}_1$ for both viscoelastic contact models, with values of $10.81\text{\%}$ (GK model) and 11.20% (Burgers model). Also, the precision of the GK contact model on the testing frequencies accounted an improvement in most cases. The average errors using the GK model (Burgers model) are 12.94% (13.04%), 7.07% (8.86%), 6.62% (8.45%), and 1.38% (1.31%) at 1, 2, 5, and 10 Hz. In general, the difference in the dynamic modulus decreased for higher loading frequencies. This can be attributed to the fact that the calibration procedure was initially performed for the highest loading frequency in order to achieve faster calibrations.

The numerical simulations with the elastic aggregate-to-mastic contact approach predicted higher average errors for the dynamic modulus. The average errors are $27.28\text{\%}$ and $23.66\text{\%}$ for ${\kappa }_m$ associated with the GK contact model and Burgers contact model, respectively. The maximum relative errors are verified for samples AM${}_3$ and AM${}_2$, with values reaching 33.80% (${\kappa }_m$ of the GK contact model) and 27.67% (${\kappa }_m$ of the Burgers contact model), respectively. Similarly to the observed when adopting a viscoelastic aggregate-to-mastic approach, the average errors increased with the decrease in the testing frequencies, which is also related to the adopted calibration procedure.

When adopting a viscoelastic aggregate-to-mastic contact model, the phase angle is over-predicted at higher frequencies, with the exception of asphalt mixture AM${}_3$, which generated a response very close to the lab-based data. Overall, the average error is reduced when adopting the GK contact model with respect to the use of a Burgers contact model, reaching values of $32.02\text{\%}$ and $33.29\text{\%}$, respectively. For an elastic aggregate-to-mastic contact model approach, the phase angle is considerably under-predicted for all mixtures, resulting in average errors for the ${\kappa }_m$ from the GK contact model and for the ${\kappa }_m$ from the Burgers contact model of $69.59\text{\%}$ and $69.22\text{\%}$, respectively. These results indicate the importance of taking into account the viscoelasticity at the aggregate-to-mastic contact when modelling asphalt materials.

The asphalt mixtures numerical simulations predicted a phase angle almost constant within the adopted frequency range. This may be due to the fact that the mastic-to-mastic contacts were calibrated based on the mastic experimental results, that as shown in Figure 4, have an almost constant phase angle. Note that, the aggregate-to-mastic viscoelastic contacts are also based on the GK contact model with a negligible phase angle variation.

7. Summary and conclusions

In the adopted asphalt mixture DEM model, the aggregates and mastic are directly represented. Compared with other DEM approaches, that implement a finer discretization of the asphalt material's internal structure, the adopted methodology has a reduced computational cost. An innovative 3D generalised Kelvin contact model is presented within the DEM framework to represent the viscoelastic behaviour of asphalt materials. When compared with the usually adopted Burgers contact model, the proposed GK contact model allows a better rheological representation for a higher range of frequencies and temperatures. Two different GK contact model formulations are presented and validated for mastics and asphalt mixtures in uniaxial tension-compression sinusoidal tests for a given frequency range. The main findings of the study presented here can be summarised:

• The proposed GK contact models were validated for mastics and asphalt mixtures under uniaxial loading for a given frequency range. The predicted numerical responses have a good agreement with the observed experimental response for both mastics and asphalt mixtures. For the two asphalt materials, when compared with a Burgers contact model, the GK contact model yielded better results for both rheological properties, dynamic modulus and phase angle. For the mastics, the average errors for $\left|E^{\ast }\right|$ and ϕ for the GK model (Burgers model) are 2.40% (3.46%) and 3.64% (4.17%), respectively. For the asphalts mixtures, the average errors for $\left|E^{\ast }\right|$ for the GK model (Burgers model) are on average 7.00% (7.92%). Due to the predicted nearly constant phase angle, the average errors for this property adopting the GK model (Burgers model) are 32.02% (33.29%).

• The proposed generalised contact models predicted, as expected, the same rheological response. However, the GK${}_1$ contact model has lower computational requirements (14% faster). Nevertheless, the GK${}_2$ contact model has the advantage of being easier to be adopted in the study of the non-linear behaviour in fracture tests due to its incremental numerical formulation.

• The influence of the contact parameters in the GK contact model was evaluated. The obtained numerical results clearly indicate that ${\kappa }_m$ is an influential parameter when assessing both the dynamic modulus and the phase angle. The results presented also indicate that a fine tune of the predicted phase angle numerical response can be obtained by adjusting the δ multiplier associated with ${\eta }_m$, as its influence on the dynamic modulus is negligible. The set of all viscosities of the contact model is also shown to have a significant influence on the predicted phase angle.

• The presented results show that the aggregate-to-mastic contacts in the asphalt mixture particle assembly should adopt a viscoelastic contact model. When compared with an elastic model, the GK contact model proposal resulted in a better agreement with the known experimental behaviour.

• The predicted phase angle of asphalt mixtures had a reduced variation with a changing frequency for the GK and Burgers contact models. This is possibly due to the experimental results for mastics, which are predominantly constant over the frequency range and the fact that the mastic-to-mastic contacts are calibrated with these experimental values.

The proposed 3D DEM framework that is applied in the modelling of asphalt mixtures is shown to be a valid approach. Following studies are planned on a larger experimental database to assess if it is possible to predict a higher variation in phase angle by implementing a more flexible calibration approach where a different coefficient of adjustment is adopted for each contact parameter.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is part of the research activity carried out at Civil Engineering Research and Innovation for Sustainability (CERIS) and has been funded by Fundação para a Ciência e a Tecnologia (FCT) in the framework of [project number UIDB/04625/2020].