3D DEM model simulation of asphalt mastics with sunflower oil

A three-dimensional particle model, based on the asphalt mastic micro-structure representation following a discrete element model framework, was developed to investigate the influence of sunflower oil (rejuvenator) on the rheological properties of asphalt mastic. Dynamic shear rheometer tests in laboratory, for a frequency range of 0.1–20 Hz and temperatures in a range between 20 and 80 ∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\circ }$$\end{document}C, were carried out in order to assess the viscoelastic behaviour of asphalt mastics containing different oil-to-bitumen content by mass proportions (2.5–20%). Master curves were constructed for two reference temperatures (30 ∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\circ }$$\end{document}C and 50 ∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\circ }$$\end{document}C). Experimental results showed that the increase in sunflower oil content resulted in a progressive decrease in viscosity. However, the rheological behaviour of the mastic containing the highest oil amount could not be properly represented in master curves, indicating that the specimen had a different rheological behaviour when compared with the lower oil contents responses. Numerical simulations of rheometer tests were carried out with an asphalt mastic particle assembly that emulated the experimental procedure. The viscoelastic contacts within the asphalt mastic assembly were simulated with a generalized Kelvin contact model. A calibration procedure was derived based on the fitting of laboratory data. The simulations were shown to have a good agreement with laboratory values. The average errors for the dynamic shear modulus and phase angle were 3.4% and 4.0%, respectively, considering both temperatures of analysis. Finally, a considerable improvement was accomplished in comparison with the numerical response obtained with the often-used Burgers model.


Introduction
Flexible pavements, the most used road pavement system [1], result from the use of asphalt mixture in the top layers, Rui Micaelo and Nuno Azevedo have contributed equally to this work.a material that combines coarse and fine aggregates, bitumen, and air voids.During the manufacturing process, fine aggregates, filler, and bitumen mix together originating the asphalt mastic, a material that influences the dynamic modulus, cohesion, and relaxation aspects of mixtures.Regardless of being designed for a lifetime between 15 and 30 years, the asphalt pavement will face damage after years due to ultraviolet radiation, moisture, repeated traffic loads, and ageing of the binder.These conditions are responsible for stiffening the binder and reducing its relaxation capacity [2], leading to ravelling, the appearance of micro-cracks, and other pathologies, especially reflective cracking, rutting, and water damage as a consequence of the brittleness and less adhesive aged bitumen [3] that eventually affect its mechanical performance [4].
Asphalt mixture is a self-healing material.However, the in-service conditions of temperature and continuous traffic load do not promote the healing mechanism [5], which may take several days.Besides, the self-healing property reduces over time due to the ageing of the bitumen, low temperatures, large size cracks, and material type.Based on the need for reducing maintenance costs, two low-cost and effective techniques have been studied in the past years to increase the healing property and its speed in asphalt pavements: the use of induction heating and encapsulated agents in the asphalt structure [6].In the thermally induced process, alternating electromagnetic fields of variable frequency, in association with a very fast increase in temperature, are applied to the surface of the road pavement to allow metallic particles previously embedded in the mixture to soften the binder and flow into the existing micro-cracks.However, the increase in temperature is related to ageing [7].
Several studies have investigated the use of capsules containing rejuvenators [1,8,9].Some authors have used vegetable oils, such as sunflower oil, a low-viscosity material, due to low health and environmental hazards, availability, and safety while managed [10].A more sustainable solution is using waste cooking oil [11,12].It has been shown that capsules can resist asphalt mixture manufacture, temperature does not affect their mechanical strength [13,14], and they present a beneficial effect on the safety and durability of asphalt mixtures [4].In such a technique, capsules remain inactive in the pavement until a mechanical trigger, i.e., an external load, breaks the capsule's shell, leading to the release of the oil content and the reduction of the bitumen's viscosity to promote the flow of the binder to the cracks with the help of capillarity [1].Al-Mansoori et al. [8] investigated the effect of self-healing adopting different contents of capsules by total weight of the asphalt mixture.The authors were able to relate temperature and healing ratios and quantify the capsule content to obtain higher healing levels.Garcia et al. [9] verified a gradual breakage of capsules during the action of compressive forces and softening of the binder in the presence of sunflower oil during the rest periods in mixtures submitted to cyclic tests.
Up to now, the improvement in the self-healing levels with the use of capsules has been reported as positive.However, their effects on the physical and mechanical properties of bituminous materials need to be investigated.Al-Mansoori et al. [1] analysed the physical, mechanical, and self-healing properties of asphalt mixtures containing different healing amounts.It was shown that capsules did not affect the water sensitivity, particle loss, and permanent deformation of asphalt mixtures.Similarly, Norambuena-Contreras et al. [10] concluded that mechanical properties, such as stiffness modulus, tensile strength, fatigue resistance, and water sensitivity, of asphalt mixtures with and without encapsulated sunflower oil yielded similar responses to experimental tests.Additionally, the air void content and density of samples generated equivalent results.
The use of encapsulated rejuvenators in asphalt mixtures has been extensively investigated.However, the analyses have been limited to a laboratory perspective.For the past decades, rigid particle models based on the discrete element model (DEM) have been adopted in the micro-mechanical analysis of asphalt materials in the study of viscoelasticity and damage behaviour.In this study, a micro-mechanical discrete model has been implemented to describe the viscoelastic behaviour of asphalt mastics containing different amounts of sunflower oil by adopting contact models to represent the interaction between the constitutive materials and an irregular micro-structure.In the following section, some studies developed in the context of the DEM in the study of asphalt materials are mentioned.

Background
Discrete element modelling is a numerical method that repeatedly applies the Newton's second law and the forcedisplacement law to calculate and update the positions, velocities, and displacement of particles in a time-stepping algorithm.The model assumes a small time step to avoid error propagation during simulations.User-defined specifications, i.e., boundary and loading conditions and the implementation of contact models, must be developed within the numerical model.The methodology was first adopted in the analysis of rock mechanics [15,16].In the past decades, the DEM has been adopted to investigate the behaviour of asphalt materials under different conditions, such as numerical creep simulations to evaluate permanent deformation, especially at high temperatures, and simulations for describing the fracture and rheological properties of unmodified and modified binders over a wide range of loading frequencies and temperatures.
One of the first applications with the DEM for asphalt materials was carried out by Rothenburg et al. [17] with the aim of studying pavement rutting.The interactions among particles were modelled with a linear viscoelastic model.The aggregate portion accounted for a major influence on the mechanical response of asphalt mixtures.In the creep testing field, Zelelew and Papagiannakis [18] simulated the creep behaviour of asphalt concrete in a two-dimensional (2D) approach adopting a Burgers model.The viscoelastic properties of the mastic phase were captured in dynamic shear rheometer (DSR) tests.The developed model underestimated the creep compliance at an initial stage and was wellpredicted in the following step.Thus, it was suggested the implementation of generalized models for a better outcome.The constituent materials of asphalt mixtures individually influence the creep behaviour of asphalt structures.Aggregates, usually comprising approximately 85% of the total volume, have a great importance on such property.In this sense, Liu and You [19] carried out a study of the angularity, orientation angles, and spheric index on the creep stiffness of asphalt mixtures.One major conclusion was that creep stiffness was directly related to the orientation angle of aggregates.Recently, Wang et al. [20] investigated the dynamic creep mechanism on asphalt mixtures under repeated loading in two dimensions.Adopting a Burgers model, the average error was less than 10%.
Some efforts to study the fracture behaviour in asphalt materials were conducted with the use of the DEM [21][22][23].Kim and Buttlar [21] investigated the aggregate size and type, temperature, stress, and displacement properties of asphalt concrete at low temperatures.An elastic model associated with a bilinear model defined the aggregate-mastic and mastic-mastic contacts.As a result, the aggregate size and strength were defined as important properties in the analysis of the softening behaviour.Dan et al. [22] investigated the cracking mechanism of asphalt mixtures in indirect tensile tests under different loading speeds.The loading speed significantly affected the peak tensile stress, and the failure process was concentrated near the loading walls.In order to investigate the three-dimensional (3D) block cracking in asphalt pavements, Wang and Buttlar [23] implemented a viscoelastic cohesive-softening contact model with the DEM.The analysis was carried out to verify the influence of the dimension of pavement segments, relaxation capacity, and ageing.Most cracks appeared on the surface due to large temperature variations, high brittleness, and high stiffness levels.
The rheological behaviour of asphalt materials under dynamic tests has been extensively evaluated in recent years [24][25][26].For instance, Abbas et al. [24] simulated asphalt mixtures under sinusoidal loads by adopting a Burgers model.The dynamic modulus was under-predicted for unmodified specimens, while over-predicted for modified ones.With the aim of reducing the errors in the numerical analysis, Ren and Sun [25] developed a 2D generalized Maxwell model to simulate the rheological response of asphalt mixtures under dynamic modulus tests, static creep tests, and bending tests at low temperatures.The model generated a response within 8% with respect to the experimental data.In a following study, a 2D generalized Kelvin model was implemented to simulate dynamic tests [26].The numerical results were compared with the response obtained with other viscoelastic models.The generalized Kelvin model better reflected the dynamic response at higher temperatures, while the generalized Maxwell model was more suitable for lower temperatures.
The viscoelastic properties of asphalt mastics have great importance to the specifications of asphalt pavements, e.g., stiffness, fatigue resistance, rutting resistance, and bleeding.Hence, Abbas et al. [27] studied the effect of the binder stiffness and mineral filler volume fraction on the asphalt mastic stiffness.Adopting an elastic model in a 2D numerical scheme, results indicated the dynamic shear modulus of the mastic highly depended on the dynamic modulus of the binder.One common test to analyse the time-temperature properties of mastics is the dynamic shear rheometer test.Tests with frequency sweep or temperature sweep configurations are possible.Therefore, Vignali et al. [28] analysed the effect of cement and limestone contents on the rutting resistance of cold mastics over temperatures varying from 0 to 90 • C. Limestone filler contributed to the resistance against permanent deformations by increasing the sample's stiffness at high temperatures.This result was later confirmed by Zhang et al. [29] when adopting a modified Burgers model to consider volume concentration and particle-interaction reinforcements of fillers.Dondi et al. [30] adopted a Burgers model and correctly predicted the viscous response and internal force interaction of binders under DSR tests in a three-dimensional model.Recently, Zhang et al. [31] developed a 2D micro-mechanical model to simulate low-temperature splitting tests and simulate the self-healing ability in asphalt mixtures.It was verified that the increase in accumulated damage gradually decreased the healing capacity in asphalt mixtures.
Previous studies have enlightened that the analysis of permanent deformation, rheological properties, and fracture of asphalt materials are the main concern for a micromechanical perspective with the DEM.More importantly, many studies have used 2D models, although they usually underestimate interlock effects in the particle assembly.Besides, many authors have adopted simple viscoelastic contact models, e.g., the Burgers model, to describe the existing viscoelastic behaviour.However, this contact model type is not the most suitable to describe such behaviour over a wide range of loading frequencies.Moreover, the self-healing mechanism on these composites has been evaluated under experimental tests.
In this sense, a 3D DEM model was developed to investigate the viscoelastic behaviour and the effect of sunflower oil in asphalt mastics under dynamic shear rheometer tests at temperatures varying from 20 to 80 • C. A generalized Kelvin (GK) contact model was adopted to represent the viscoelastic contacts within mastic particles in the virtual sample.Additionally, this study was extended to the analysis of different oil-to-bitumen content by mass proportions (2.5-20%) in mastic samples.Finally, the numerical results obtained for the GK model were compared to the corresponding response obtained when using a Burgers model.

Generalized Kelvin contact model
The generalized Kelvin 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.In this model, the resultant displacement (u) derives from the sum of the elastic displacement of the free stiffness (u el ), the irreversible displacement of the viscosity unit (u vp ), and the sum of the delayed elastic displacement from the Kelvin elements (u i ve ).As the elements are placed in series, the contact model's force ( f ), the free stiffness element's force ( f el ), the free viscosity element's force ( f vp ), and the Kelvin chain elements' force ( f i ve ) are equal to one another [32].The developed model is constructed on an incremental formulation, in which the increment in displacement results in the increment in contact force of the model [33].The delayed elastic displacement and the irreversible displacement derive from the relationship between force and the rate of creep compliance (J ) to relate the force-displacement behaviour.Accordingly, the equation governing the displacement of the time-dependent elements is defined with the use of an integral hereditary formulation that admits an initial value followed by an integration over time: where subscript ξ assumes ve and vp to correspond to the viscoelastic and viscoplastic elements of the model, respectively, and t is an integration variable.The initial creep compliance of the respective elements is expressed in the following equation: Therefore, it is possible to determine the increment in the displacement of the viscoelastic and viscoplastic elements by deriving Eq. 1. Consequently, the displacement governing the two components at time t + t can be determined, while the displacement corresponding to the elastic portion derives from the direct relationship between the contact force and the stiffness κ m .The respective elastic, viscoelastic, and viscoplastic displacements are: where C el , C i ve , and C vp are the inverse of the respective contact parameters κ m , η i , and η m , and τ i is the retardation time defined as the relationship between viscosity and stiffness of the ith-element of the Kelvin chain: Hence, the total displacement increment of the generalized Kelvin contact model is given by: Substituting the elastic (Eq.3), viscoelastic (Eq.4), and viscoplastic (Eq.5) components into Eq.( 7) and reorganizing the terms, the displacement increment is calculated.However, the resultant equation can be simplified into a pseudo-elastic expression relating the increment in displacement to the increment in contact force, as shown: where C t corresponds to the total matrix of flexibility of the generalized Kelvin contact model, calculated in the following equation: Moreover, the matrix of flexibility corresponding to the time-dependent portion of the contact model can be defined as: Accordingly, the pseudo-elastic displacement increment ū is calculated as: where the displacement increment u 1 is estimated as: Lastly, the contact force f t+ t results from the sum between the force at time t and its increment defined in Eq. 8: 3 Experimental analysis

Materials
Laboratory tests were conducted on asphalt mastics to assess the influence of sunflower oil (rejuvenator) on the rheological properties of the mixture binder.In this study, five different oil-to-bitumen content by mass-0%; 2.5%; 5%; 10%; 20%-were used, corresponding to asphalt mastic samples M-0 (control sample), M−2.5, M-5, M-10, and M-20, respectively.The proportions of the rejuvenator were chosen following previous investigations.Norambuena-Contreras et al. [10] used a 5.23% oil-to-bitumen content by mass in bitumen, corresponding to 0.50% of rejuvenator by total mass of asphalt mixture, in the analysis of the mechanical and crack-healing properties of stone mastic asphalt.The healing process of fractured samples was improved in comparison with samples without a healing agent.The same amount of rejuvenator was adopted by Garcia-Hernández et al. [6].In fact, this quantity was proposed by Al-Mansoori et al. [1] after analysing the effect of three oil-to-bitumen content by mass in bitumen (1.1%; 2.8%; 5.5%) on the mechanical properties and self-healing ability of asphalt mixtures.
A paving grade bitumen 35/50 (density 1.03 g/cm 3 ) and limestone filler RC480 (density 2.72 g/cm 3 ) with 83% of particles passing the 0.063 mm sieve, composed of calcium carbonate (CaC03), were used to fabricate the control asphalt mastic sample (M-0).The specifications of the bitumen used were the following: needle penetration at 25 • C of 46×10 −1 mm and a softening point of 51 • C. A commercial sunflower oil with density, smoke point, and flash point corresponding to 0.92 g/cm 3 , 227 • C and 316 • C [8,14], was used.
Previous to testing, the bitumen was submitted to the Rolling Thin Film Oven Test (RTFOT) process to simulate the ageing that occurs during the fabrication and paving of asphalt pavements, and therefore, to study the effect of the vegetable oil in an aged material.Accordingly, the binder film was heated at 163 • C with a constant supply of air for 75 min.After the ageing procedure, the bitumen was conditioned at room temperature in metal tins.Then, to fabricate the mastics, the previously dried filler and the bitumen were heated at 150 • C for 60 min, and then, the two constituents were mixed with mechanical agitation at 400 rpm for 20 min.Asphalt mastics were produced with a filler volumetric concentration of 36.6%, which is near the midpoint of concentration values for dense asphalt mixtures [34].The mastics were later transferred to circular silicon moulds with dimensions corresponding to the DSR specimen testing requirement.The specimens were cooled at room temperature and de-moulded directly to the DSR plates.

Dynamic tests
During the past decades, the study of the viscoelastic behaviour of asphalt binders has progressively evolved.Different experimental tests are available to obtain the rheological properties of asphalt materials, such as the bending beam rheometer test, the direct tension tester, and the dynamic shear rheometer test.Therefore, the DSR test was adopted in this study to capture the dynamic shear modulus and phase angle of the mastic samples under sinusoidal shear loading at intermediate and high service temperatures, typically between 5 and 85 • C, and different loading frequencies.Figure 2 shows the rotational rheometer Bohlin Gemini HRnano utilized for experimental tests and the sinusoidal loading mechanism.The test is conducted in such a way that the spindle applies a torque to the testing specimen, oscillating from point A to point B, returning to point A, and subsequently, the spindle oscillates to point C, returning to the origin afterwards.This movement comprises a complete cycle and may be repeated continuously during the analysis.
The sample is sandwiched between two circular parallel plates (discs).The bottom plate remains fixed in all directions, while the upper plate sinusoidally oscillates above the testing specimen.In addition, three essential pieces of information are initially required to properly work: the thickness (gap), the diameter of the sample, and the loading definition for testing.Based on the European Standard EN 14770 [35], for moderate temperatures, it is recommended to adopt an 8-mm-diameter disc with a 2-mm gap.On the other hand, a 25-mm-diameter disc with a 1-mm gap should be adopted for elevated temperatures.In these experimental tests, an oscillatory shear was imposed on the asphalt mastic and the resulting stress was monitored.This solution is often common for asphalt materials rather than adopting an imposed stress [36].
The mastics were subjected to a sinusoidal load with a strain control mechanism of 0.50 × 10 −4 m/m.As a consequence, the strains remain small enough to assure the analysis is conducted in the linear viscoelastic zone.The dynamic shear modulus and phase angle were recorded across 0.01 to 20 Hz loading frequencies at testing temperatures varying from 50 to 20 • C with a -10 • C temperature rate and across 0.01-10 Hz at loading temperatures varying from 50 to 80 • C with a 10 • C temperature rate.For each mastic specimen, at least two replicates were prepared for each testing frequency.The mean values of the replicates were adopted as the laboratory results.
The following equations are used to determine the dynamic shear modulus and the phase angle: where |G * | represents the dynamic shear modulus and indicates the resistance to deformation when repeatedly sheared at a testing frequency, τ max is the absolute value of the peak shear stress, γ max is the absolute value of the peak shear strain, φ is the phase angle, t corresponds to the time difference between two adjacent peak moments, and T is the loading period.
The peak shear stress and the peak shear strain at the minimum and maximum points of a complete sinusoidal loading cycle are calculated according to: where T is the applied torque, r is the radius of the sample, θ is the rotation angle, and h is the height of the mastic.

Resulting laboratory master curves
The Time-Temperature Superposition Principle (TTSP) method has been used for asphalt materials in the past decades.You et al. [37] adopted the method to reduce computational time in 3D simulations for mixtures under creep tests.The principle is applied for thermorheologically simple materials, such as the asphalt mastic, for the linear viscoelastic region in dynamic simulations.It is based on the direct relationship between time and temperature, which assumes an equivalence between elevated temperatures and the reduction in loading time, as such for lower temperatures and the increase in loading time.
After defining a reference temperature (Tr), the dynamic shear measurements are shifted based on a determined value for α t (shift factor) to achieve a smooth and unique master curve.The reduced frequency ( f r ) was determined from the test frequency ( f ) and the shift factor as: In this study, the shift factors were manually attributed.Based on this representation, it is possible to assess the rheological behaviour of asphalt mastics in a wider range of loading frequencies.Constructing a master curve at a moderate temperature improves the understanding of the viscoelastic behaviour of specimens at lower loading frequencies, while a higher reference temperature generates insight at higher loading frequencies.Therefore, master curves were constructed at two temperatures: 30 • C and 50 • C.During the experimental analysis, the viscoelastic behaviour of samples was verified at a common temperature (50 • C) using both diameter discs (8 mm and 25 mm) in the DSR tests.It is worth noting that the results at this temperature were selected according to which best suited the properties of each disc [35].Thus, in the case of samples M-0 and M−2.5, the results with the 8-mm-diameter plate were adopted.On the contrary, the 25-mm-diameter plate better represented the rheological response of samples M-5, M-10, and M-20 due to the reduced dynamic shear modulus by the action of the oil, as the 25-mm plate is more suitable for a limited range of 100 kPa [35].Therefore, the experimental results obtained for the different asphalt mastics considering the two reference temperatures are shown in Figs. 3 and 4.
From Figs. 3 and 4, it is possible to verify a wider range of loading frequencies that the rheological properties may be analysed.In both cases, the control sample M-0 presented higher values of dynamic shear modulus and lower values for the phase angle.As expected, by adding the rejuvenator to the mastic, a considerable decrease in dynamic shear modulus and an increase in phase angle were evident.This behaviour was directly related to the oil content added to the binder.The rejuvenator considerably reduced the viscosity of the asphalt mastics.In the asphalt pavement, this would allow the binder to flow into the micro-cracks to extend its service life.The lower values for |G * | and higher values for φ obtained at lower frequencies (corresponding to slow-moving vehicles) for mastics containing oil indicate a reduction in the resistance to permanent deformation.Hence, the oil percentage in the sample must be appropriate.
In the case of the 30 • C (50 • C) master curves, the behaviour of mastics was extended to loading frequencies varying from 0.00001 to 178.2 Hz (0.0004-7962.1 Hz), depending on the specimen type.Comparing the |G * | response among the testing samples for a 30 • C master curve, a considerable reduction was produced with the addition of the rejuvenator along the frequency range.A similar decrease occurred for a 50 • C master curve.This is evidenced in Table 1 when comparing the dynamic shear modulus of mastic specimens at 10 Hz at 50 • C. It was verified a reduction of 57% when adding 2.5% of sunflower oil content to the binder.When doubling the amount of rejuvenator, only 8.6% of the original value endured.For samples M-10 and M-20, the dynamic shear modulus was reduced to 2.7% and 0.3%, respectively.Laboratory measurements of phase angle for asphalt mastic with a 20% sunflower oil content produced an inconsistent pattern.The result might be a consequence of the high oil content, and at this stage, the mastic may not be a thermorheologically simple material to be represented in master curves.In accordance, several studies have recommended an oil-tobitumen (aged) ratio within the range of 5.5% for a mastic to recover its initial properties [10,41], which corroborates to the 20% oil-to-bitumen ratio being an excessive amount.

Macroscopic numerical formulation
In this study, a macroscopic viscoelastic generalized Kelvin model (Fig. 5), which reflects the general behaviour of the material, was adopted to fit the experimental results.The macroscopic parameters must be converted to contact parameters for each asphalt mastic due to the impossibility of determining the contact parameters for the DEM model directly from the experimental analysis.
The stress and strain relationships for the macroscopic GK model result from the association of the elastic portion, viscoelastic portion, and viscoplastic portion constituting the model.Accordingly, the resultant strain is determined from the sum of the strain generated among the elements of the model ( el , i ve , vp ).On the contrary, the stress referred to the individual portions is equal to one another.The stress of the composing elements is: where σ el , σ i ve , and σ vp refer to the stress unit of the elastic, viscoelastic, and viscoplastic elements of the macroscopic model, respectively.The application of a dynamic stress to the model and the resultant dynamic strain are shown, respectively, in the following equations: where σ 0 is the stress value at time equals zero, ε 0 is the strain value at time equals zero, ω is the angular frequency, and t is the elapsed time.After substituting Eqs.22 and 23 into Eqs.19-21, the following expressions can be obtained: σ 0 e iωt = iωC m ε vp e i(ωt−φ) ( The viscoelastic response for the GK model due to dynamic shear stress is represented using the shear modulus (G * ), the ratio between the incremental increase in stress (Eq.22) and the corresponding increase in strain (Eq.23).It comprises a real component (G ), referred to as the storage shear modulus, and an imaginary component (G ), referred to as the loss shear modulus: After substituting Eqs.24-26 into Eq.27, it is possible to determine the shear modulus.Considering its reciprocal value, the shear compliance (J * ) is: where J is the real part of the shear compliance, referred to as the storage shear compliance, and J is the imaginary component, referred to as the storage shear compliance.The components are calculated in the following: Thus, the dynamic shear modulus, the storage shear modulus, the loss shear modulus, and the phase angle are: The fitting procedure takes into consideration a minimization function that is based on the relative difference between the predicted values of the real and imaginary components of the dynamic shear modulus and the results obtained in laboratory at each testing frequency: where G θ and G θ are the real and imaginary components of the predicted dynamic shear modulus at the θ -th frequency, G 0 θ and G 0 θ are the real and imaginary components of the measured dynamic shear modulus at the θ -th frequency, and n is the number of data points.

Resultant calibration of input parameters
The macroscopic parameters for the GK model and Burgers model were determined after fitting laboratory measurements through Eqs.29-35.The parameters were calculated for each of the two master curves (30 • C and 50 • C), considering the wide loading frequency for a proper calibration in each scenario.For the GK model, the number of elements (7-9) was selected according to the analysis of the error predictions with respect to the previous equations.The referred minimization process was adopted to define an initial guess for the generalized contact parameters.Therefore, the set of macroscopic parameters for both the GK model and Burgers model are presented in Tables 2, 3, and 4. The particularities and influence of the materials constituting each asphalt mastic are indirectly described in the respective set of parameters.
The resultant calibrations of the macroscopic parameters obtained with the viscoelastic models for all asphalt mastics are illustrated in Figs. 6 and 7.It is possible to verify that a good agreement between the resulting fitting and the experimental results was obtained with the use of the GK model at both reference temperatures.This is more evident when compared to the corresponding response generated with the use of the Burgers model.Hence, it is expected to achieve a similar response in the DEM numerical simulations.
The input parameters for the DEM simulations are determined based on the stress-strain relationships of a viscoelastic beam simulating two particles in contact.Thus, the micro-scale parameters of the GK contact model in the normal direction can be obtained from the following expressions: where L is the sum of the radius of the two contacting particles, A is the cross-sectional area, δ is the coefficient of adjustment to suit the macro-scale properties into contact parameters, which takes into account geometric relations of the particle assembly and contact properties, n indicates the normal direction, and subscript ξ assumes i for the ith ele-   ment of the generalized model or m to describe the elastic and viscoplastic elements of the model.The contact parameters for the tangent direction are determined based on the product between the values computed for the normal direction and the contact stiffness ratio α, the ratio between the stiffness of the tangent and normal directions.

Numerical specimen preparation
In this study, discrete element modelling simulations replicated the laboratory DSR tests performed for mastics.For this purpose, a micro-structure asphalt mastic was created using a random generation method.The method has been chosen for the fact that is less time-consuming, more affordable, and does not require special equipment and human resources, comparatively to other particle generation methodologies.The numerical sample preparation scheme has in prospect the use of the Laguerre-Voronoi diagrams to generate polygonalshaped particles.A similar generation model has been adopted in 2D [38].
The process of filling a determined volume with asphalt mastic particles is based on the definition of an initial particle size distribution.At this stage, input parameters, such as the maximum and minimum particles diameters, radius distribution, porosity, density of particles, and geometric characteristics, are required to create the centre of gravity  of the particles considering an imposed uniform probability distribution.Asphalt mastic particles are introduced from the largest diameter size to the smallest in the virtual assembly.In fact, particles are generated with half of their diameter's size in order to avoid excessive particle overlap that may compromise a final stability.After achieving the desired number of particles, a radius expansion procedure is performed to achieve the required particle size.
Once the initial positions have been assigned, a DEM cohesionless cycle is carried out.In this step, particles are allowed to interact and overlap with their neighbouring ones.After reaching an equilibrium state, the contact forces throughout the assembly assume zero and the positioning particle process is concluded.The polygonal-shaped particles are created following a triangulation of their centres of gravity with the use of a 3D weighted Delaunay scheme to generate the dual 3D Voronoi diagram [39].In this model, particles interact with their neighbouring ones by the Voronoi tessellation.The contact area and the corresponding contact location are defined by the area of the associated Voronoi cell facet.This technique increases the number of contacts per particle in the assembly compared to traditional mod-  els [38].However, the increase in the contact points result in a considerable increase in the time of numerical simulations.Therefore, a reasonable particle size distribution for the asphalt mastic needs to be chosen in order to limit computational costs.
The DEM simulations following laboratory DSR loading tests were performed in a three-dimensional particle assembly representing asphalt mastics M-0, M−2.5, M-5, M-10, and M-20.It is worth noting that the temperature and magnitude of dynamic modulus existing in experimental tests for defining geometric specifications are not applicable for numerical simulations, and consequently, the results obtained for both particle geometries are the same.Accordingly, simulations were carried out with the use of a single virtual sample.A cylindrical numerical sample (Fig. 8) was generated with dimensions corresponding to 8 mm in diameter and 2 mm in height.The specimen comprises 3642 randomly distributed particles sandwiched between two rigid walls (simulating the plates of the DSR equipment) with diameters ranging from 0.35 to 0.40 mm.Moreover, the assembly contains 22,516 contacts, resulting in an average number of contacts per particle of approximately 6.2.Previous studies have suggested a minimum of four contacts per particle for a 3D-particle assembly to model the internal contact structure [17,40].The numerical sample presents two contact types: interactions between asphalt mastic particles and the walls and contacts within asphalt mastic particles.With the aim of representing these contacts in the DEM model, a linear elastic contact model was assigned to describe the mastic-to-wall contacts, while the viscoelastic behaviour within asphalt mastic particles was characterized with the GK contact model and Burgers model.

Modelling results
In the first step of simulations, dynamic loading tests were carried out to calibrate the coefficient of adjustment (δ) for each mastic considering both reference temperatures.The same laboratory conditions were considered in the numerical simulations.A rotational velocity was applied to the upper wall while the bottom wall remained fixed in all directions.The analyses were performed for loading frequencies with respect to the master curves of each mastic shown in Figs. 3  and 4. The parameter δ for the different mastics was between the small range of 1.80 and 2.00.Table 5 shows the coefficient δ for the studying group.
For each loading frequency, the dynamic shear modulus and phase angle are calculated based on Eqs.14-17, taking into consideration a complete cycle for numerical calculations.The DEM results obtained adopting the GK contact model were compared to the results of the Burgers model and the laboratory data.The following equations were used to evaluate the relative errors in the measurements: where As can be seen, the GK contact model predicted reasonably well the viscoelastic behaviour of the different mastics for both temperatures.On the other hand, the Burgers model satisfactorily described the evolution of dynamic shear modulus along the loading frequencies; however, the model could not properly represent the phase angle.Other studies have suggested a similar improvement in the analysis of bituminous materials under dynamic loading with the use of generalized models when compared to the response for the Burgers model [25,26].Moreover, the relative average errors for the numerical simulations are presented in Table 6.
Similarly to the analysis of the shear dynamic modulus, the GK model yielded a better numerical response for the phase angle when compared to the Burgers model.The relative errors for such property in the group of asphalt mastics for the GK model (Burgers) were on average 4.2% (38.7%) for Tr = 30 • C, unlike relative errors of 3.9% (19.8%) for Tr = 50 • C. In regard to the maximum errors, a considerable increase was verified when adopting the Burgers model.For instance, for Tr = 30 • C, the maximum relative errors for the GK model (Burgers) for mastics M-0, M−2.5, M-5, M-10, and M-20 were on average 14.9% (98.2%), 8.2% (88.2%), 8.9% (84.9%), 16.8% (81.4%), and 14.3% (20.3%).In the case of Tr = 50 • C, a similar agreement was observed, with percentages for the GK model (Burgers) of 6.7% (89.3%), 13.2% (78.5%), 8.8% (86.5%), 17.2% (58.9%), and 25.8% (37.6%).It can be noticed that the highest average and maximum errors in the measurement of the phase angle adopting the GK model were verified in the case of asphalt mastic with the highest oil content.This was due to the inconsistent pattern of the experimental results that represented a higher complexity in the numerical modelling.Nevertheless, a better prediction was achieved for the GK model.In general, the GK model was able to capture the changing in the viscoelastic properties of the asphalt mastics with the addition of sunflower oil along the considered loading frequencies.The contact model considerably reduced the average and maximum errors during simulations when compared to the outcome generated for the Burgers model.Although the GK model has proved to better represent both rheological properties, it was verified that the Burgers model produced a good response when characterizing the dynamic shear modulus.Nevertheless, the generated responses at intermediate and high loading frequencies were very divergent for the phase angle.The computational runtimes for the GK model varied from 31 s to 59 h.The numerical simulations were carried out with an Intel i9@3.70GHz multi-core processor.The GK contact model numerical simulations were two times slower than the simulations carried out with the Burgers model, which is still reasonable taking into account the increase in accuracy that is obtained.

Conclusions
In this study, the influence of a rejuvenator (sunflower oil) on the rheological properties of asphalt mastics was investigated through laboratory tests.With the aim of understanding the micro-mechanical consequence of such impact, a three-dimensional DEM model was5 developed with the implementation of a generalized Kelvin contact model to represent the viscoelastic interactions within the virtual mastic sample in numerical simulations following experimental testing guidelines.The main results are summarized in the following: • Different oil-to-bitumen ratios by mass (2.5%, 5%, 10%, and 20%) in asphalt mastic were analysed in DSR tests and the results compared to the control sample (0%).The results indicated a considerable reduction in the viscosity of the mastics promoted with the addition of the rejuvenator.This was verified with the reduction in dynamic shear modulus and the increase in phase angle.The rate of reduction in the viscosity increased considerably faster as increasing the amount of rejuvenator.Thus, it is important to correctly specify the content of rejuvenator in the sample.• The experimental results were correctly represented with the use of the Time-Temperature Superposition Principle in two temperatures (30 • C and 50 • C).This evidenced that the specimens containing sunflower oil were rheological simple materials.Therefore, the mastics were able to be analysed in a wider range of loading frequencies.
However, the M-20 specimen could not be well represented with this methodology, indicating that the amount of rejuvenator was excessive and the mastic sample was not a rheological simple material.• A viscoelastic generalized Kelvin contact model was implemented to characterize the viscoelastic behaviour of asphalt mastics.Numerical simulations following experimental tests were carried out.The GK-predicted response was compared with the response predicted using the Burgers model.The GK model generated a better response for the dynamic shear modulus and phase angle for all asphalt mastic specimens at both reference temperatures.The average errors for the dynamic shear modulus and phase angle for the GK model (Burgers) were 3.4% (19.6%) and 4.0% (29.2%), respectively, considering the different mastics at both temperatures.• The results obtained with the Burgers model showed partial agreement with the lab-based data (shear modulus).
Therefore, the Burgers model should only be adopted for initial calibrations for further developments with the GK model.
The implemented GK contact model has shown to be a promising approach to the encapsulated rejuvenator in following the modelling of the self-healing mechanisms of damaged asphalt mixtures.

Fig. 1
Fig. 1 Generalized Kelvin contact model: a directions and b model illustration

)Fig. 2
Fig. 2 Dynamic shear rheometer: a equipment and b illustration

Fig. 8
Fig. 8 Asphalt mastic sample: a experimental and b virtual Error |G * | (%) and Error φ (%) are the percentage errors of the dynamic shear modulus and phase angle, |G * | DEM and |G * | exp are the numerical and experimental values for the dynamic shear modulus, and φ DEM and φ exp are the numerical and experimental values for the phase angle.The numerical results for the rheological response at 30 • C and 50 • C obtained for the GK model and the corresponding values for the Burgers model with respect to the experimental analysis are shown in Figs. 9 and 10, respectively.The time step adopted in the DEM simulations varied from 8.2 × 10 −8 to 1.4 × 10 −6 s/step according to the asphalt mastic sample, studying temperature, and contact model type.Due to computational efforts, numerical simulations were conducted for frequencies higher than 2 Hz.

Fig. 9
Fig. 9 Laboratory data and simulation results for GK model and Burgers model (Tr = 30 • C) for asphalt mastics: a M-0, b M−2.5, c M-5, d M-10, and e M-20

Fig. 10
Fig. 10 Laboratory data and simulation results for GK model and Burgers model (Tr = 50 • C) for asphalt mastics: a M-0, b M−2.5, c M-5, d M-10, and e M-20

Table 1
Experimental results at 10 Hz for mastics at Tr = 50 • C

Table 2
Fitted parameters for macroscopic generalized Kelvin model at Tr = 30 • C

Table 3
Fitted parameters for macroscopic generalized Kelvin model at Tr = 50 • C

Table 4
Fitted parameters for macroscopic Burgers model at Tr = 30 • C and 50 • C

Table 5
Coefficient δ for mastics

Table 6
Relative errors between DEM numerical simulations and experimental data for the GK model and Burgers model