of the mechanical behavior of asphalt mixtures using fractional rheology to model their viscoelasticity

of aggregate mineral particles attached with an asphalt

Because most materials exhibit both fluid (viscous) and solid (elastic) behavior, multivariable non-linear mathematical and mechanical models are based on Hooke's law of elasticity and Newton's law of fluid viscosity [7].As a consequence of these laws and to better describe the behavior of complex viscoelastic materials, a number of rheological models have been proposed in the past, such as Maxwell, Kelvin-Voigt, standard linear solids, Burgers, and generalized Maxwell or Kelvin-Voigt, among others [8][9][10].
According to Chang et al. [11] the Burgers model describes the mechanical properties of asphalt mixtures better than the Maxwell and Kelvin-Voigt models.Saltan et al. [12] noted that the Burgers model is appropriate for determining the viscous component of stiffness.In contrast, Ma et al. [13] reported that it is difficult to control the variability of experimental tests, due to the complicated interference of environmental variables.Jaczewski et al. [9] stated that no model commonly used in the literature is ideal for predicting the deformations of asphalt mixtures, since each model has advantages and disadvantages.In addition, Celauro et al. [14] mentioned that the loading-unloading cycle cannot be described with the same set of model's parameters, since the prediction results are inaccurate in the complete cycle.
One of the most important phenomena produced in asphalt deformation is the creep and recovery of the loading-unloading cycle [7].According to Alrashydah et al. [15] the creep process is a significant mechanical property for pavement design.Consequently, creep allows one to predict the mechanical response by means of linear regression models.Li et al. [10] showed that the variation of loads and temperatures has also significant effects on the viscoelasticity of the mixtures.Ma et al. [16] indicate that the creep deformation is proportional to the increase in the content of air voids in the mixture.Xie et al. [17] determined that in the attenuation creep stage of an asphaltic mastic the deformation rate decreases gradually while in the steady creep stage the creep rate is relatively stable.In addition, the greater the load, the greater the deformation rate and the creep deformation are in every stage.On the other hand, Yao et al. [18] applied a repeated creep-recovery method to modified binders, demonstrating that these properties are relevant for evaluating delayed elasticity.From the same method, Liu et al. [19] have classified the deformation capacity of modified binders for different traffic levels.Han et al. [20] studied the elasticity of different binders by means of an extended recovery during 10 cycles of loading and unloading.They concluded that creep deformation and subsequent recovery are important for the design of an asphalt mixture.
Several experimental investigations have been presented with the use of models based on differential equations including fractional derivatives.Fractional calculus is a fast developing branch of mathematical physics that deals with integrals and derivatives of traditional functions of any arbitrary real or complex order.This mathematical approach is able to reveal important physical meanings for the viscoelasticity of asphalt mixtures and other materials [21].Within this framework, Gao et al. [22] verified that modified models that provide fractional functions, such as Kelvin and Maxwell, possess fewer rheological parameters and they agree well with experimental measurements.The fractional Burgers model proposed by Oeser et al. [23] requires four parameters to describe creep strain in an asphalt mixture.It was reported that this model exhibits better agreement with experimental tests than previous classical models.Although the fractional calculation improves the physical representation of the deformations, it does not completely characterize an asphalt mixture.The main practical limitation of the aforementioned models is that they do not consider a detailed characterization for each element that makes up a mixture and thereby they do not examine the asphalt binder and the aggregate separately.
In the following sections of this paper, a new approach related to the Burgers model [24] is presented to describe the time-dependent behavior of the mixture.The proposed rheological model of asphalt viscoelasticity is introduced in Section 2, where the asphalt binder and the aggregate contributions are considered separately.Equations are obtained using derivatives of fractional order and they were implemented into computer codes using MATLAB©.Creep, recovery, and relaxation phenomena in an asphalt mixture are analyzed using the new model and compared with the Burgers model.Results of the experimental validation of the model are presented in Section 3. Finally, the main conclusions are summarized in Section 4.

A rheological model of asphalt viscoelasticity
Currently, the viscoelasticity of asphalt mixtures and binder is characterized using the Burgers model, also called the "four elements model," that describes a viscoelastic material [25].This model was not specifically developed to explain the deformation produced in an asphalt mixture, since the model does not detail differences between constituent elements.Therefore, this study proposes a mechanical model of asphalt viscoelasticity that assumes an aggregate particle enclosed by an asphalt material, i.e., an elastic element surrounded by an elastic-viscous set.
Note that an asphalt mixture has a centralized mechanical model in itself that reveals the behavior of its deformations and details the independent character of the rheological properties of both the aggregate and asphalt binder, producing compression when subjected to mechanical tests [26,27].
Fig. 1 shows this relationship from the aggregate characterization, which is represented by the elasticity constant ξ 2 and the elasticity and viscosity constants representative of the asphalt binder, ξ 1 and η, respectively.

P o s t -p r i n t
The fractional differential equation that governs the previous model was obtained by analyzing two strains in series [28], where the first strain was developed with the classic Maxwell model [29]: where ε(t) is the strain, σ(t) is the stress, and D t α and D t γ are defined as the fractional derivatives with respect to time t and fractional exponents α and γ, respectively.The subscript of ε is the index of the corresponding elements in the model shown in Fig. 1.
The second strain differs from the Kelvin-Voigt model [30], since it has two purely elastic springs and a viscous damper connected in parallel to consider in more detail the elasticities of an asphalt mixture: where β is also a fractional exponent.
It is important to emphasize that the fractional exponents must have the same ranges, between 0 and 1, to satisfy the classical equations.The complete system is connected in series, so that the total strain is the sum of the strains in each component [31]: After substituting the equations and transforming the derivatives in an algebraic form, the differential equation of the proposed asphalt mixture viscoelastic model is (see Appendix A): Applying the Laplace transform to Eq. ( 4) and assuming zero initial conditions yields the following stress-strain constitutive equation: where .

Creep phenomenon
The creep phenomenon appears when the material is subjected to a constant stress at a certain time, given by a single step loading function σ = σ 0 H(t), where H(t) is the Heaviside unit step function [32].Doing this, and returning to the previous differential Eq. (5), gives [33]: (5) P o s t -p r i n t As the parameters α and γ must be equivalent for the first deformations, these are cancelled when starting the first elastic deformation.Furthermore, the creep modulus F(t) can be obtained by inverse transforming Eq. ( 6) [34][35][36]: where Γ(.) is the Gamma function.In Eq. ( 7), the third term of the creep modulus equation is called the infinite Mittag-Leffler series, which depends on parameter β, and describes the transformation process of the binder.This can reach a value of 1, since it acts at the same time as elasticity [37].
The change of curvature, or the inflection value, is known as the relaxation time and corresponds to the complete system (aggregate-binder).Parameter α appears in the second term of Eq. ( 7), which is interesting since it represents the last state of the creep phenomenon, indicating that if this parameter reaches a value of 1, it would represent a completely viscous material.However, this assumption is far from reality, due to the dependence of the aggregate on the asphalt mixture (note that this parameter depends exclusively on the capacity of the binder, and not on the complete mixture) [38].Therefore, it is possible to avoid setting boundary conditions, assuming α, β = 1 as values of the fractional exponents and, consequently, obtaining the classical equation from the proposed model: We can compare Eq. ( 8) with the Burgers model, where the classic creep equation is given by the following mathematical expression [39]: In Eq. ( 9), the first term of the equation corresponding to the creep of the material (from the point of view of the Burgers model) represents the elasticity of a single material, which is correct if the aim is to study the creep phenomenon of a binder or of any viscoelastic material separately.However, if an asphalt mixture is tested, the Burgers model states that the first strain is caused by the aggregate-binder.This could be incorrect, because when a mixture is subjected to a load the first deformation is always suffered exclusively by the binder, since this a softer material.
Based on the above, the proposed viscoelastic model in asphalt mixtures adequately describes the first strain of the asphaltic material, which corresponds exclusively to the binder strain.Thus, to continue analyzing the other terms of the proposed model, numerical examples are given in the following subsections to compare the effect caused by the fractional parameters α and β.

Fractional parameter β
Parameter β represents a change in the mixture state, physically describing the transformation from an elastic to a viscous state [40].Its mathematical domain is limited in its lower part by the value 0 until it reaches a maximum value of 1.The above indicates that mechanically, the model presents two springs and a damper connected in parallel in its curvature.This demonstrate a certain incompatibility between the proposed model and the Burgers model, deduced from classical mechanical reasoning based on two hypotheses: (1) The proposed model is characterized by a change in the viscoelastic curvature, considering the second elastic jump of the entire asphalt mixture; (2) The space and time dependent model requires fractional calculus.
Fig. 2 parametrically shows how the fractional exponent acquires different curves by varying the range of β, both in the Burgers model and the model proposed in this study.In this example, elasticity for the binder, for the aggregate, and viscosity for the binder were used.These values have merely been chosen in this example to illustrate more clearly the differences between the Burgers model and the proposed approach. ) P o s t -p r i n t Fig. 2a (the Burgers model) shows that the strain reaches deformations of 0.02 mm, while with the proposed model (Fig. 2b) the maximum deformation value does not exceed 0.007 mm.This is because the model considers the aggregate and the binder separately.The Burgers model shows a point of inflection when its deformation is 0.01 mm, which corresponds to the elastic jump if parameter β is 0. In contrast, the new model reaches a value of 0.005 mm, since it adds the elasticity of the aggregate.However, it is important to note that each value of β differs for the various binders, temperatures, and load frequencies.
Therefore, the largest difference between both models is based on the relaxation time (η/M), since it has a sum of elasticities that affect a strain (with slope values close to the proposed viscosity).The model proposed in this study presents a curve with a parallel trend, typical of an elastic material, thus representing the participation of the aggregate in the asphalt mixture.

Fractional parameter
From the mechanical point of view, parameter α corresponds to the last damper of the model (see Fig. 1), both in the Burgers model and the proposed model [41].This damper represents the last phase of the creep phenomenon in which the binder reaches full viscosity, meaning that its final state will never be viscous due to its dependence on the aggregate.
Therefore, its relationship depends exclusively on the viscous constant of the binder, which does not mean that at this point the aggregate does not contribute to the total strain.Each binder will present specific slopes and α values depending on the temperature and load frequency.
Thus, it is interesting to know the dependence on time and temperature of these parameters to appropriately characterize the model.
To show the dependence of the creep modulus with parameter α, a fixed value of β = 1 was set and the creep of the material was plotted, as shown in Fig. 3.This value of β was set to simulate the binder undergoing a viscoelastic transformation following the spontaneous elastic jump, which is typical of spring-damper mechanical work.The results obtained confirm that parameter α acts when the material begins to deform with a constant slope.The elastic jump of the binder, with a deformation value of 0.02 mm for an α = 1, gives rise to a final strain of 0.2247 mm (at 10 s).Additionally, if α = 0, its instantaneous strain become constant due to the transformation of the fractional exponent (see Fig. 3).The above is due to the fact that if α = 0, this represents a spring, giving the sum of two instantaneous jumps with initial values of 0.04 mm, followed by growth until reaching the maximum value of 0.0467 mm.
The application of the load is directly proportional to exponent α.Thus, parameters α and β can be determined by comparing the model with the results obtained from experimental tests.Once these values are obtained, other physical phenomena may be predicted and analyzed (such as the dynamic module and/or relaxation test).The importance of this model lays in that we may design an asphalt mixture with a modified binder, so we can build an optimized and more resistant pavement layer.

Recovery phenomenon
The recovery phenomenon begins when the initial stress applied to the material in the creep test is released [41].This time-dependent recovery is a characteristic of each material and depends on both the type of load that is applied and the temperature variation to which the material sample is subjected [42].
In an experiment, when the applied load exceeded the elastic limit of the material, the sample became plastic, modifying the physical properties.To describe this process, it is necessary to consider the fractional differential Eq. ( 4), eliminating the concept of initial stress .However, this is where a controversy arises, since it is not easy to understand how fractional equations work when the initial conditions are not zero.In addition, applying the Laplace transform is complicated, since there is currently no physical meaning for these initial conditions.Therefore, the right-hand side of Eq. ( 4) is zero.After applying the Laplace transform, Eq. ( 4) becomes: where α and β adopt values between m − 1 and m, where m is the nearest integer to the value of α or β.Analyzing Eq. ( 10) using first the Maxwell model and then the modified Kelvin-Voigt model gives the following equations, respectively [43]: Note that the fractional parameters can reach maximum values of 1, so m − 1 = 0. Therefore, the summations in Eq. ( 11) are canceled and only the initial strain conditions remain.This initial strain is defined as for the system connected in parallel.Now, we can express the recovery phase of the asphalt viscoelastic model by adding the elastic jump from the Maxwell model, denoted as , which gives: P o s t -p r i n t Fig. 4 shows the Mittag-Leffler function for the recovery of asphalt mixtures (the summation in Eq. ( 12)).When the value of β tends to 0, the model becomes a pure spring that indicates a constant recovery time, which is not typical of the asphaltic material.

Relaxation phenomenon
The relaxation phenomenon originates when the material is subjected to constant strain [45].To do this, a test was carried out in which the samples were deformed during a finite interval of time to analyze when the effort began to vary over the course of time.To describe this phenomenon, a constant strain was used.At any time t, deformation began with a finite jump, so that the differential equation given in Eq. ( 4) can be Laplace transformed [46].
To describe the relaxation, the deformation process was carried out by adding the relaxations of the model in series, i.e., when the material began to be governed by the Maxwell model and the modified Kelvin-Voigt model [47].Then, in the Laplace domain we get: Finally, Eq. ( 13) can be inverse transformed term by term to obtain the time-dependent relaxation modulus R(t) [36,41] for the asphalt mixture viscoelastic model:

P o s t -p r i n t
The model given by Eq. ( 14) differs from the typical function of the rheological relaxation of viscoelastic models due to the contribution of the fractional derivatives.This is the case of the second term of Eq. ( 14), which comes from the modified Kelvin-Voigt model, whose fractional parameter is β and validates the following relationships: -If β reaches a value of 1, it is an undefined case, since the Gamma function is indeterminate at zero.Its solution is based on the application of a limit when the parameter tends to that number, resulting in the classical relaxation equation whose first parameter is the sum of elasticities of the aggregate-binder system connected in parallel, plus the Mittag-Leffler function: -If β reaches much lower values (close to 0), β contributes considerably to the result of the final relaxation function, explaining the viscosity capacity of the binder to dissipate energy.Additionally, note that each value within the established range contributes to the final relaxation modulus.
The proposed model is exclusively linked to the capacity of the binder to release stress.Fig. 5 shows the stress given by the second term in Eq. ( 14) as a function of both time and β.
The figure shows that as β goes from 0 to 1, there is a release of stress over time that starts at a value of η = 50 and corresponds to the viscous constant value.When β reaches 1, stress release is zero, which is represented by a straight line with the value of 0 over time, thus confirming the classical model.For values of β < 1, there is residual stress release.This means that knowing the value of β associated with each binder is important for real applications.
In the third term of Eq. ( 14), another important point is the presence of the Mittag-Leffler function that becomes the classical exponential of the Maxwell model when α = 1.This is clearly observed in Fig. 6, which shows the variation of the third term of Eq. ( 14) as a function of both α and time. ( Fig. 5 Contribution of the second term of Eq. ( 14) to the relaxation modulus as a function of β.

P o s t -p r i n t
The classic Maxwell model explains just a part of the relaxation phenomenon.In the third term of Eq. ( 14), as the value of α tends to 0, the stress release becomes smaller and when α reaches 0, the binder becomes theoretically purely elastic although, in practice, this does not occur.In addition, when α reaches values close to 0.2, it undergoes a sudden change in relaxation, reaching the linear limit with a value corresponding to half the elastic capacity of the binder.

Experimental verification
In order to verify the proposed model, several creep-recovery laboratory tests were carried out at various temperatures on two different bituminous binders.12 Marshall samples measuring 101.6 mm in diameter and 65 mm in height were manufactured and compacted to 75 blows per side.Semi-dense AC16S gradation was used, with the variation of a conventional asphalt binder with penetration grade B50/70 and a binder modified with crumb rubber PMB45/80-65 that were used for manufacturing the samples.The softening points of the conventional and modified binder were 51.6 °C and 72.3 °C, respectively.Ophitic nature coarse aggregate and limestone nature filler were also employed.
Uniaxial static creep and recovery tests were carried out on all the samples.A constant Heaviside load of 5 kN was applied in a set period to determine the rheological properties of the materials [13,16].Tests were conducted for a temperature range between −10 °C and 35 °C (with a scale of 15 °C).
Creep-recovery tests on both binders measured at different temperatures are shown in Fig. 7. Experimental results (solid line) are compared with those obtained by using the proposed model with fitted parameters (dashed line).The theoretical model agrees well with the experimental data at all temperatures for both binders, both during the creep phase and during recovery.In addition, the model accurately represents the increasing stiffness with decreasing temperature and the lower accumulated deformations in the modified binder compared to the conventional one.The model can show that as the testing temperature increases, a greater viscoelastic deformation occur, each time generating a higher elastic jump and transforming the nonlinear creep.
Fig. 6 Contribution of the third term of Eq. ( 14) to the relaxation modulus as a function of α.

P o s t -p r i n t
The results of the fractional simulation at different temperatures are shown in Table 1.The creep phenomenon develops the instantaneous elastic jump represented by the first spring of the model (ξ 1 ).Unlike the present approach, classical models cannot represent this elastic jump with real values since it is not associated with the binder (represented by the first spring).The mixture with the conventional binder increases the elastic jump with the temperature, which causes the decreasing of the elastic modulus, transforming its behaviour to the viscoelastic range.The mixture with the PMB 45/80-65 binder keeps the elastic jump relatively constant generating a more stable value of elasticity with the variation of temperature in the linear viscoelasticity.P o s t -p r i n t After the elastic jump of the binder, nonlinear viscoelastic curvature is developed with a second elastic jump, which considers the elasticity of the complete mixture given by the sum of the elasticities of binder and aggregate, M = ξ 1 + ξ 2 [see Eq. ( 5)].The result predicted by the model makes it possible to derive the Young's modulus for the aggregate, whose value is 635.0MPa.This result explains why the aggregate is not deformable in a mixture, delivering volumetric and load dissipation properties capable of increasing the second elastic jump and delaying the viscosity.Unlike the Burgers model, the proposed one maintains the properties of elasticity (ξ 1 ) and dynamic viscosity (η) of the asphalt binder constant for the parallel and series deformations, originating a physical meaning to evaluate the deformations in the whole test cycle.In addition to the elastic jump, a time-dependent deformation related to parameter β of the model is produced, which represents the fluctuation between the elastic and viscous state.The value of β in the mixture with modified binder increases more than that in the mixture with conventional binder, generating more open viscoelastic curves and allowing a recoverable creep.
Finally, the parameter α corresponds to the last fractional damper and it decreases with the temperature due to the viscoelastic transformation previously arisen with β.When the recoverable deformation is rapid (close to zero) the slope is greater.However, if the transition is slow or with greater development of recoverable viscoelasticity (higher β values) the result is a flatter final slope.The value of α is characteristic of the non-recoverable creep and does not present large differences between mixtures either with conventional or modified binders.
The behaviour of the asphalt mixture would be represented as a Newtonian fluid (value of α = 1) when working with classical calculation, which is not the real rheological situation since an asphalt mixture never behaves that way.In the Burgers model, this is solved by adding an extra parameter of dynamic viscosity despite of losing the actual physical representation of the set of elements.Therefore, the fractional calculation helps understand the viscoelastic deformations and clarify the real physics when designing a pavement.
The dynamic viscosity η depends on the temperature and is related only to asphalt binders.The proposed model fits the strain curves with the same value of η for the whole test cycle.For low temperatures, this value is high due to the solid consistency with high stiffness that the material attains.For higher temperatures (35 °C) its physical state is more viscous, and the values of η are lower.The proposed model satisfactorily represents the lower dynamic viscosity of the conventional binder compared with the modified binder, which produces a lower elastic deformation in the asphalt mixture.
The second part of the test shows the recovery of the asphalt mixture, which fits the experimental data with the rheological values obtained in the creep process.The parameter represents the elastic jump of the recovery phenomenon.An increase in the temperature results in an increase of , which is reflected in the final values of the accumulated plastic deformation.At high temperatures, the jump of elastic recovery in the mixtures is higher.It is important to notice that at a fixed temperature, the rheological values of the binder (ξ 1 and η), the aggregate (ξ 2 ) and the parameter β obtained in the creep process were the same as those in the recovery phase.This fact reveals that, unlike the existing models, the loading-unloading cycle can be described with the same set of model's parameters.The value of β increases with temperature causing more open curves that decay faster to the starting point.The mixture with PMB 45/80-65 binder has lower accumulated plastic deformations compared to the mixture with conventional binder, which is related to the increase of β by the modified binder.This behaviour allows one to predict the phenomenon of recoverable creep of mixtures.

Conclusions
From the results of this study, the following conclusions can be made: -The new model of asphalt mixture viscoelasticity proposed here represents in a more precise way the characteristics of the materials that form the mixture.In this study, we proposed a model of asphalt viscoelasticity that involves a mechanical arrangement that characterizes an asphalt mixture in a more clear and detailed way than classic viscoelastic models, which do not explicitly consider the materials that make up the mixture.Unlike the classical models, the new model is characterized by two springs and a damper connected in parallel that represents an aggregate particle enclosed by an elastic-viscous binder element.The new model better represents practical cases of asphalt mixtures used in the construction of pavement layers.
-An elastic jump at the beginning of the creep modulus is represented by the new model that is not shown by the Burgers model.This elastic jump depends on the characteristics of the binder.After this elastic jump, the binder transforms its state from elastic to viscoelastic and the complete mixture acts.
-The recovery phenomenon depends on the parameter β of the new model.Subjecting the material to small loads causes a gradual loss in recovery for several reasons, such as the loading-unloading cycle of a test sample and high loads that may affect final conditions.
-The use of fractional derivatives in the new model results in new terms in the equation for the relaxation modulus.These new terms are useful to explain better the total energy release caused by the binder (associated with the value of β) and the whole mixture (associated with the value of α).
-Rheological parameters in the model ξ 1 , ξ 2 , and η allow representing the properties of each one of the materials in the asphalt mixture, correlating the proposed model with the deformation behavior that occurred in the laboratory.In addition, fractional parameters α and β determine the nonlinear viscoelasticity of asphalt mixtures, detailing the ranges of recoverable and non-recoverable creep for different binders and working temperatures.
-Further work is planned to compare the theory presented in this paper with experimental tests that will include different asphalt mixtures and loading frequencies.

Declarations of interest
P o s t -p r i n t

Fig. 1
Fig. 1 Schematic diagram of the proposed model of asphalt mixture.

Fig. 2
Fig.2reveals differences between the results obtained for both models.When parameter β reaches values close to 0, the binder resists certain loads and its transition to the viscous state is zero.From the mechanical point of view, this behavior corresponds to a spring, instead of a damper, resulting in a model comprised of three springs connected in parallel.In contrast, as parameter β reaches values close to 1, the binder begins behave like a viscoelastic material with exponential growth, resulting in a model with two springs and a damper connected in parallel.

Fig. 2
Fig. 2 Results of strain as a function of both time and : a) Burgers model; b) proposed model.

Fig. 3
Fig. 3 Results of strain as a function of both time and for the proposed model with β = 1.α

Fig. 4
Fig. 4 Mittag-Leffler function for different values of showing the recovery phenomenon of asphalt mixtures.

Table 1
Creep-recovery of asphalt mixtures for a load of 5 kN.