A Creep Model of Asphalt Mixture Based on Variable Order Fractional Derivative

: In order to quantitatively describe the time ‐ varying mechanical properties of asphalt mixture during creep process, a nonlinear viscoelastoplastic creep model was proposed, by using variable ‐ order fractional calculus. The differential order of the variable ‐ order fractional element of the model is no longer constant, but a variable that changes with time, which reflects the changes of the mechanical properties of the material during the creep process. Whereas the tertiary creep phase is modeled by the viscoplastic element with time ‐ varying viscosity, which is attributed to damage evolution. The uniaxial creep compression tests of AC ‐ 13C asphalt mixture under different stress levels (0.7 MPa, 0.9 MPa, 1.1 MPa, 1.3 MPa, 1.5 MPa, 1.7 MPa) were carried out with MTS ‐ 809 testing machine at 25 °C, and the test results were analyzed by the model using Levenberg– Marquardt optimization algorithm. It is shown that creep damage occurs when the applied stress exceeds a certain critical value, and the damage incubation time depends on the applied stress level. The higher stress decreases the damage incubation time. The model is in good agreement with the experimental results, and is applicable to describe the entire creep process, which consists of primary, steady and tertiary stages. Moreover, the variation of the model parameter can describe the change of viscoelastic properties of the material during the creep process. The differential order of the variable ‐ order fractional element is constant during the primary creep stage, indicating that the creep behavior of the asphalt mixture is linear viscoelastic in small strain range. For the same stress level, the fractional order of the steady creep stage is greater than that of the primary creep stage, and it increases with the increasing stress level, which shows that the viscous behavior in the steady creep is more remarkable than that in the primary creep, and the higher the stress level, the more prominent the viscous performance exhibits.


Introduction
Asphalt mixture is a typical viscoelastic-plastic multiphase composite material, which has complex rheological properties. Asphalt pavement has been widely used in the highway construction because of the advantages of no connection joints, low noise and short construction period, and more than 95% of the current highway pavement in China is asphalt pavement. Due to the increasing traffic volume and the significant increase in axle load, the failure of asphalt pavement becomes increasingly serious. The major pavement diseases are rutting, potholes, and cracking. Creep is one of the most important mechanical properties of the asphalt mixture, and it is directly relevant to the rutting behavior. Therefore, it is of great significance to predict the creep behavior of the asphalt mixture by establishing its reasonable constitutive model.
Up to now, many constitutive models have been proposed and applied to study the dynamic viscoelastic and static viscoelastic behaviors of asphalt mixture, which can be divided into two categories in general: integral constitutive model and differential constitutive model. Considering the fact that the multiple integral constitutive models with too many parameters are rarely adopted by practical engineering, the study of the integral constitutive model mainly focuses on the single integral constitutive model. As for linear viscoelasticity, the Boltzmann superposition principle is adopted; as for nonlinear viscoelasticity, the modified Boltzmann superposition principle and the reduced time models are adopted [1,2]. The Schapery single integral constitutive model has been widely used to analyze of nonlinear viscoelastic mechanical behavior of asphalt mixtures [2][3][4][5]. In the study of the differential constitutive models, the classical models, such as the Kelvin model, the Maxwell model and the Burgers model, consist of the series and/or parallel combination of springs and dashpots, and are generally used to describe the linear viscoelasticity with the advantages of simply, intuition and clear concept. Cheng et al. [6] used the generalized Kelvin model, the generalized Maxwell model and the Burgers model to describe the creep behavior of asphalt mixture at various temperatures, and showed that the generalized Kelvin model and the generalized Maxwell model are superior to the Burgers model in describing the variation of viscoelastic properties of asphalt mixtures with loading time. In order to describe the nonlinear viscoelastic mechanical behavior, the influences of viscoplasticity and damage were taken into consideration, and the corresponding modified models were proposed on the basis of linear viscoelastic constitutive models. Zhang et al. [7] proposed a viscoelastic damage constitutive model for asphalt mixture under cyclic loading, which was constructed by coupling the Burgers model with the modified Chaboche damage model. Bai et al. [4] combined the Schaperyʹs nonlinear viscoelastic model and the modified Schwartzʹs viscoplastic model to characterize viscoelastic and viscoplastic behaviors of compressed asphalt mixture. Darabi et al. [5] suggested a coupled nonlinear viscoelastic, viscoplastic and hardening-relaxation constitutive model for asphalt mixtures, and verified it by the dynamic moduli and repeated creep-recovery tests. Moreover, the differential operators in the above mentioned classical differential constitutive models are local integer differentials. However, the fractional differential operators are nonlocal and applicable to describe the memory characteristics of viscoelastic materials [8]. It has become a research focus in the study of viscoelasticity, because of the advantages of few parameters and high accuracy [8][9][10]. Fini et al. [11] used a simple fractional viscoelastic model to evaluate the performance of asphalt binders at low service temperatures through the damping ratio and dissipated energy ratio. Celauro et al. [12] analyzed the creep and creep recovery test data and adopted the fractional Burgers model to characterize the viscoelastic behaviors of asphalt mixture. Lagos-Varas et al. [13] proposed a fractional rheological model consisting of a fractional Maxwell model and a modified fractional Kelvin model to describe the static viscoelasticity of asphalt mixture. Luo et al. [14] proposed a modified fractional Zener model to quantitatively describe the symmetric dynamic viscoelastic properties of asphalt mixture. Nevertheless, the mentioned fractional differential constitutive models have their limitations, such as no change in the fractional differential order and without considering the rheological property change in creep process.
Variable order fractional calculus developed from the basis of classic fractional calculus theory [15][16][17], has attracted much attention in recent years. Wu et al. [18] adopted the variable-order fractional Maxwell model to analyze the creep behavior of salt rock at low stress level. Tang et al. [19] further proposed a creep model derived from variable-order fractional calculus in order to describe the nonlinear creep deformation of rock.
Although some researchers have proposed viscoelastic and viscoplastic constitutive models for asphalt mixture on the basis of constant-order fractional derivative, there is no report about the analysis of the mechanical behavior of asphalt mixture based on variable-order fractional differential theory. In this paper, a new viscoelastoplastic creep damage model based on a variable-order fractional differential will be constructed to study the creep behavior of asphalt mixture at various stress levels. The model can describe the whole creep process including primary, steady and tertiary stages, it uses variable-order fractional elements to model the primary and steady creep, and uses time-varying viscoplastic element to consider the damage induced tertiary creep stage. The paper is organized as follows: the basic constitutive models of the constant order fractional viscoelastic element and the variable-order fractional viscoelastic element are briefly derived in Section 2. In order to model the damage induced tertiary creep, the dashpot element with variable viscosity is introduced in Section 3. The new viscoelastoplastic creep damage model based on variable-order fractional differential is then constructed in Section 4 and verified with experiments in Section 5. Finally, some conclusions are drawn in Section 6.

Constant Fractional Viscoelastic Element
There are various definitions of fractional calculus, and the definition of fractional calculus based on Riemann-Liouville (R-L) is a widely used one. The R-L fractional integral of function with order is defined by [15]: where 0 1, and Γ is the gamma function. The R-L fractional differential of function with order is defined as the inverse operation of the fractional integral: where n is the minimum positive integer greater than . 1 if 0 1, and the R-L fractional differential of function with order is then written as: In the classical viscoelastic constitutive modeling, the series and/or parallel combinations of springs and dashpots are usually used to describe viscoelasticity. In order to describe the viscoelastic mechanical behavior of materials more flexibly, Abel established the relationship between the stress and the fractional derivative of the strain, which is modeled by the fractional viscoelastic element as shown in Figure 1. The constitutive relation of the fractional viscoelastic element reads: where , representing the viscosity coefficient of the material, E is the material modulus, is the characteristic time. It can be seen that when →0, the fractional dashpot element reduces into a Hooke spring, and when →1, the fractional dashpot element reduces into a Newtonian dashpot.

Let
( is the Heaviside function), the creep equation for the fractional viscoelastic element can be obtained by integrating Equation (4):

Variable-Order Fractional Viscoelastic Element
Unlike constant fractional calculus, the fractional order of variable-order fractional calculus no longer remains constant, but changes with time. Let 0 1 , then the fractional calculus of function with order is given below [15]: where is a piecewise constant in time domain, and we can assume that , 1,2, … , when . Then, the constitutive equation of the fractional viscoelastic element in the time period ( ) is given by: is the viscosity coefficient corresponding to . Under creep condition, , the creep strain of the fractional viscoelastic element can be obtained from Equation (8): From Equation (1), the -th order integration of a constant C is given by: Thus, Equation (9) becomes: Based on Equation (11), the creep strain of the fractional viscoelastic element after n−1 time segments can be written as:

Dashpot Element with Variable Viscosity
Damage will occur when the material is under creep at higher stresses; the accumulated damage causes accelerated creep to failure. Therefore, the influence of damage evolution should be taken into consideration for the constitutive description of the tertiary creep stage. The mechanical property of the classical dashpot element is characterized by the constant viscosity coefficient , and its constitutive equation is expressed by . Zhou et al. [20] demonstrated that the damage accumulation and meso-cracking evolution in the tertiary creep stage would lead to the continuous decrease in the viscosity coefficient of rock, which satisfies the exponential decay law. Assuming that the creep tertiary stage starts when t , the time dependent viscosity coefficient can be expressed as: where is the initial viscosity coefficient in the tertiary creep stage, is the parameter describing the damage evolution rate. Replacing in the constant viscosity dashpot element with in Equation (13), the constitutive equation of the time-varying viscosity dashpot element as shown in Figure 2 can be expressed as: Let , then the creep strain of the time-varying viscosity dashpot element is given by:

Variable Order Fractional Viscoelastoplastic Creep Model of Asphalt Mixture
As shown in Figure 3, the creep behavior of asphalt mixture can be divided into two types according to the stress level. From the curve (a) in Figure 3, the asphalt mixture mainly exhibits linear viscoelasticity at low stress level, and the creep process only includes two stages: primary stage (0 ) and steady stage ( ). At higher stress levels, damage will occur in the creep process after a certain time ( ), and, in turn, damage evolution will accelerate the creep. The time ( ) is called the damage incubation time, which is determined by the stress level. It is well accepted that different constant applied stresses correspond to different damage incubation times, and the higher the applied stress, the shorter the damage incubation time required. It is expected, therefore, that there exist two different critical stresses at a given temperature: one being the so-called critical damage initiation stress, d , which is the minimum value of stress below which no damage occurs for any elapsed loading time, and the other, the critical fracture stress, f , at which the material fractures instantaneously, leading to a zero damage incubation time. Thus, the damage incubation time can be expressed by the model as [21]: where A and n are material constants determined from the test data fitting. As seen from curve (b) in Figure 3, when the creep stress is higher than the damage initiation stress d , besides the primary stage (0 ) and steady stage ( ), the tertiary stage will occur for During the primary stage, the strain rate gradually decreases with time; during the steady stage, the strain rate basically keeps constant, and the strain increases with a constant speed or tends to be a stable value; during the tertiary stage, because of the damage evolution, the viscosity of the material decreases with time and the strain rate increases with time, which results in a rapid increase of deformation and ultimate failure. According to the above two types of creep behavior and their characteristics, a fractional viscoelastoplastic creep model consisting of variable order fractional viscoelastic element and time-varying viscosity dashpot element will be constructed to analyze and describe the creep behavior of asphalt mixture.
The creep process was divided into primary stage, steady stage and tertiary stage according to the time history, and each stage was described by different viscoelastic elements or viscoplastic elements. Considering the instantaneous elasticity, viscoelasticity and viscoplasticity of the asphalt mixture, a variable order fractional viscoelastoplastic creep model was proposed. Figure 4 shows the schematic diagram of the model. If the creep process only consists of the primary stage and steady stage, the time-varying viscoplastic element will be removed from the model, which is then reduced to the variable order fractional viscoelastic creep model. Compared with the classic viscoelastic and or viscoplastic model, the viscoelastic element and/or viscoplastic element in the present model only play their roles in the corresponding creep stages. At low stress levels (σ ), the asphalt mixture only exhibits the primary creep and steady creep. As seen in Figure 4, at the moment of loading, the asphalt mixture produces an instantaneous elastic deformation, which can be described by the spring element in the model: In the primary stage, due to the low elastic modulus of the asphalt component in the mixture, viscous flow will occur, which will reduce the aggregate gap. In addition, the asphalt mixture will be consolidated and hardened, and its ability to resist deformation will be enhanced, so its creep rate decreases. This stage can be modeled by the first variable order fractional viscoelastic element, as marked by (1) in the dashed box, of the model shown in Figure 4: where is the intersection time of the primary stage and the steady stage. With the passage of time, the porosity of the asphalt mixture becomes smaller, the mixture components bear the load together, and the creep rate gradually turns to be stable. This behavior will be simulated by the second variable order fractional viscoelastic element, as marked by (2) in the dashed box, of the model shown in Figure 4: Therefore, the creep model of two-stage creep of asphalt mixture can be expressed as: where , means taking the smallest value from the variables a and b, 〈 〉 is the Macauley bracket, indicating 〈 〉 | | /2.
When the stress exceeds the damage initiation stress, σ d , the creep curve of the asphalt mixture includes three typical stages. Since there is no damage in the first two stages, it can still be described by Equation (20). Whereas, for the tertiary stage, the damage evolution causes the viscosity coefficient of the mixture to decrease continuously, which results in an increase in creep rate until the overall failure. Such behavior at this stage is simulated by a viscoplastic element with time-varying viscosity, which is a Bingham-like element composed of a friction element and a variable-viscosity dashpot in parallel, as marked by time-varying viscoplasticity in Figure 4. Similar to Equation (15), the creep equation of the time-varying viscoplastic element becomes: where is the start time of the tertiary creep. Therefore, by combining the creep equation for the primary stage and steady stage, the creep model of the three-stage creep of the asphalt mixture can be expressed as below:

Creep Test of Asphalt Mixture
As for the design of asphalt mixture under heavy traffic condition, Zhonghai AH-70 heavy traffic road asphalt was selected to prepare the AC-13C asphalt mixture samples. The asphalt mixture gradation was designed according to Standard Test Methods of Bitumen and Bituminous Mixtures for Highway Engineering [22], as listed in Table 1, and the asphalt-aggregate mass ratio was 4.57%. A cylinder with a diameter of 150 mm and a height of 150 mm was made by the CONTROLS 76-PV2522 asphalt mixture gyratory compactor, and then was cored and cut into a sample with a diameter of 100 mm and a height of 100 mm, as shown in Figure 5.  Uniaxial compression creep tests were carried out on the AC-13C asphalt mixture samples at 25 °C by using the MTS-809 testing machine, as shown in Figure 6. Considering the influence of the stress level on the creep behavior of the asphalt mixture, Constant stresses of 0.7 MPa, 0.9 MPa, 1.1 MPa, 1.3 MPa, 1.5 MPa and 1.7 MPa were applied, to consider the effect of stress on creep behavior. Three tests were carried out under each stress level and then averaged. Figure 7 shows the averaged creep test data at various stresses.

Results and Disscussions
In order to distinguish the linear/nonlinear characteristics of the creep behavior, the creep strain is divided by the corresponding applied stress to obtain the creep compliance curve, as shown in Figure 8. It can be seen that the creep compliance curves for stresses of 0.7 MPa, 0.9 MPa and 1.1 MPa are almost coincident, indicating that the creep compliances at these stress levels are stress independent, and the asphalt mixture exhibits linear viscoelastic creep behavior; when the stresses are higher than 1.1 MPa, the creep compliance curves depend on the stress level, and the asphalt mixture shows nonlinear characteristics. Assuming that the nonlinear creep is caused by damage evolution, the creep damage initiation stress can be taken as d 1.
1 MPa. Figure 9 shows the creep rate curve and creep acceleration curve at various stress levels, and it can be seen that the primary creep stages under all stresses are completed at 300 s. While, for the stresses of 1.  Zhu et al. [23] conducted uniaxial compression tests on AC-13C asphalt mixture, and obtained the compression strength ranging 3.70 MPa ~ 4.38 MPa. In this study, the compression strength is set to be the average value, i.e., f 4.0 MPa. As shown in Figure 10, using Equation (16)   By means of the Levenberg-Marquardt optimization algorithm in 1stopt, nonlinear fitting of the test data in the primary stage, the steady stage and the tertiary stage for various stress levels were completed by use of Equation (22), and the corresponding parameters were determined, as listed in Table 2. It can be seen that in the stress range considered, the model parameters of the primary stage creep are the same constants, regardless of the applied stress, which indicates that when the strain is small, the asphalt mixture behaves linearly viscoelastic; under the same applied stress, the fractional differential order for the steady creep stage is greater than that for the primary stage creep, i.e., . Moreover, in the steady creep stage, the fractional differential order increases with increasing stress level, which demonstrates that the viscoelastic behavior of the material is more inclined to viscous flow when the applied stress increases, indicating the shortening of the characteristic time and the increase of the creep rate. In the tertiary stage, the parameter  increases with the applied stress level, which is consistent with the fact that the creep damage rate increases with stress level. Figure 11 shows the comparison between the model and the test results. It can be seen that the model is in good agreement with the test results, which verifies that the model can effectively describe the whole primary-steady-tertiary creep process of the asphalt mixture under different stress levels.

Conclusions
Based on the fact that the asphalt mixture shows different viscoelastic characteristics in its different creep stages, a variable order fractional viscoelastoplastic creep model was proposed to model the instantaneous elasticity, variable order fractional viscoelasticity and time-varying viscoplasticity. Some conclusions can be drawn as follows: (1) When the creep stress on asphalt mixture exceeds a certain critical value, the damage will occur, and the damage incubation time depends on the applied stress level. The higher the stress, the shorter the damage incubation time. (2) The proposed model agrees well with the results of the creep tests, which verifies that the model can effectively describe the primary-steady-tertiary creep process of the asphalt mixture under different stress levels. (3) The parameters of the model have a clear physical meaning, and the change of the fractional differential order quantitatively describes the change of the viscoelastic properties of the material at different creep stages. The fractional differential order for the steady creep stage is greater than that for the primary stage creep. Moreover, the fractional differential order for the the steady stage increases with increasing stress level.