1. Introduction
Asphalt mixtures are considered as heterogeneous pavement materials typically composed of asphalt, aggregate, admixtures, and air voids. Furthermore, asphalt mixtures are the primary materials used worldwide for constructing pavement layers because of their remarkable advantages, such as short construction periods, long service lives, recyclability, ease of maintenance, and comfortable driving. Although aggregates, which are a type of rigid material, account for approximately 90% of all components within the asphalt mixtures, such mixtures exhibit complex rheological properties. Asphalt mixtures exhibit both viscous (fluid) and elastic (solid) behavior [
1,
2,
3], and their mechanical properties are highly dependent on temperature, loading stress level, loading time, and loading speed [
4,
5]. The theoretical and experimental studies on the mechanical behaviors of asphalt mixtures provide a crucial basis for designing pavement layers [
6,
7]. Creep behavior is amongst the mechanical characteristics of asphalt mixtures [
8], and it is strongly related to the rutting formation of the pavement structure [
3,
9]. Therefore, it is important to establish a constitutive model that can accurately predict the creep behavior of asphalt mixtures.
The improvement of constitutive models is gaining increasing interest in the field of rheology in order to better describe the mechanical behaviors of complex rheological materials. Furthermore, several rheological models of asphalt mixtures were proposed in the past [
10,
11,
12,
13,
14,
15,
16,
17,
18,
19]. These models can be divided into integral constitutive and differential constitutive models [
15,
20]. The integral viscoelastic constitutive model, which considers the loading history and memory of materials, is developed based on methods such as the Boltzmann superposition principle or the Leaderman modified superposition principle. Compared with multiple integral models, single integral models have fewer material parameters and simpler forms, which makes them easier to use in practical engineering. Schapery’s model, which is widely applied to characterize the viscoelastic properties of asphalt mixtures, is a single integral constitutive model [
21]. Im et al. [
22] described the viscoelastic and viscoplastic behavior of mixtures by combining Schapery’s viscoelastic model and Perzyna viscoplasticity with a generalized Drucker–Prager yield surface.
Classical integer differential constitutive models, such as the Maxwell model, Kelvin model, Burgers model, and generalized models, comprise springs and dashpots in series or in parallel. The mathematical expressions of the differential constitutive models can correspond to the mechanical models that have advantages of intuitionism and simplicity. The classical integer derivative constitutive models are widely used to analyze the creep behavior of asphalt mixtures. The Burgers model provides a better description of mechanical behavior than the Maxwell and Kelvin–Voigt models [
23]. However, Cheng et al. [
24] reported that, compared to the Burgers model, the generalized Kelvin (GK) model and the generalized Maxwell (GM) model have a better performance in describing the creep performance of asphalt mixtures. Although the classical constitutive models can describe the viscoelastic properties of asphalt mixtures, they yield a large error in describing the initial stage of the creep performance [
13]. With the development of viscoelasticity theory, the fractional derivative operator that can describe the historical memory dependence and spatially wide relevance of a material meets the research requirements for examining the complex viscoelastic, non-Newtonian fluid and the porous media mechanics [
25]. Furthermore, the fractional calculus theory achieves great success in describing the rheological behavior of asphalt materials [
13,
14,
15,
26,
27]. Lagos-Varas et al. [
13] established a viscoelastic model using fractional-order derivatives, and unlike the Burgers model, the fractional order model can predict the elastic jump observed at the beginning of the creep modulus. Celauro et al. [
28] demonstrated that the fractional-order Burgers model can accurately characterize the creep/recovery response of asphalt mixtures. In addition, the fractional derivative model has the advantage of exhibiting a simpler expression, utilizing fewer parameters, and providing more accurate results. Furthermore analysis revealed that a fractional derivative model comprising two elements of fractional derivatives in series with a Maxwell element can be used to model the rheological behavior of asphalt binders. This includes the dynamic viscoelastic behavior, static creep, and relaxation characteristics [
27].
The entire creep process of asphalt mixtures is characterized into three stages: primary (decelerated stage), steady (stationary stage), and tertiary (accelerated stage). Furthermore, different deformation mechanisms occur throughout the creep process of asphalt mixtures at various stresses. At a lower stress level, the asphalt mixture is gradually compacted and the aggregates are gradually stacked to form a stable skeleton structure, which results in a “consolidation effect” [
19]. This deformation mechanism is known as the creep hardening mechanism [
18], and the asphalt mortar resists deformation together with the skeleton structure during the procedure. However, at higher stress levels, the stress in the asphalt mixture exceeds the frictional resistance between the aggregates and the bonding between the aggregate and asphalt mortar, which results in the fragmentation of the asphalt mixture and the destruction of the skeleton. This consequently accelerates the viscoelastic flow rate and creep failure of the asphalt mixture in the tertiary stage. Here, the deformation mechanism is called the damage softening mechanism [
18]. Linear models (e.g., the ordinary Burgers model, generalized Kelvin model, and fractional-order Burgers model) cannot describe the tertiary stage of creep.
One potential approach is to add a plastic element into the classic differential models to construct viscoelastoplastic models to accurately describe the three-stage creep behavior of asphalt mixtures at high stresses. The Nishihara model, which adds a plastic element to the Burgers model, was developed to describe the tertiary stage of creep and was used in the mechanical response analysis of bituminous materials [
19]. A nonlinear viscoelastic-plastic creep model for asphalt mixtures was developed using variable-order fractional calculus, where time-varying viscoplastic elements were used to describe the tertiary stage of creep performance [
15]. In addition, damage occurs in the asphalt mixture during the creep process. Therefore, another feasible approach is coupling a continuum damage evolution law with a linear viscoelastic model. Sun et al. [
29] developed a damage constitutive model on the basis of the generalized Burgers model and the Perzyna viscoplastic flow theory, combined with a modified Rabotnov damage theory, which can better reflect the three stages of deformation in creep testing and the hardening and softening effects in the constant strain rate compression testing. Zhang et al. [
18] suggested that the characteristics of the three-stage permanent deformation were attributed to a competition between the damage softening and strain hardening effect, and they proposed a viscoelastoplastic damage mechanics model wherein the damage and hardening variables were introduced to modify the Burger’s model for describing the deformation of asphalt mixtures.
Most existing creep damage models are formed by the integral order difference constitutive model coupled directly with the damage factor. These models have the disadvantage of many material parameters and overly complex expressions, which are not conducive to their generalization in practical engineering. Furthermore, these models directly couple the strain caused by damage (called the damage strain) to the viscoelastic strain of the asphalt mixture, which is not conducive to studying the strain development mechanism and damage characteristics during the creep process. Therefore, this study aims to propose a creep damage model based on the fractional derivative operator that considers both the creep hardening mechanism and damage softening mechanism to characterize and analyze the creep process of asphalt mixtures at various stress levels and temperatures. Simultaneously, a unified damage evolution model is constructed for different stress levels, and the development of the strain and damage characteristics of the asphalt mixture throughout the creep process is analyzed. Finally, a method for the statistical quantification of the damage evolution of the asphalt mixture is proposed.
The paper is structured in the following way: a nonlinear viscoelastic creep damage model adopting the fractional rheology theory is developed in Part 2. The compression creep tests are conducted on AC-13 asphalt mixtures in Part 3. In Part 4, the experimental results are described and the proposed model is experimentally validated. Meanwhile, the model parameters and the creep damage evolution of the asphalt mixtures are analyzed in Part 4. Finally, several conclusions are drawn in Part 5.
4. Test Results and Analysis
4.1. Test Results
Three tests were performed at each stress level, and the results were averaged. The averaged creep strain vs. time curves for the applied stress levels at different temperatures are shown in
Figure 5.
Figure 5 shows that the creep deformation of the asphalt mixtures is stress- and temperature-dependent. Here, the creep curve at 30 °C and 0.8 MPa is used as an example illustrate the three stages of creep and the time at which each stage begins. The curve of three stages of creep and the curve of creep rate are shown in
Figure 6a,b, respectively. The creep rate curve was obtained by deriving the creep curve versus time. In
Figure 6,
t1 and
t2 represent the start times of the steady stage and tertiary stage, respectively. Additionally, during this time, the creep rate remains almost constant and the value reaches a minimum. Therefore, the change in creep rate can be used to determine which creep stage the asphalt mixture is in. Combining
Figure 5 and
Figure 6, at all temperatures, the creep process under larger stress levels exhibited all three stages of creep, whereas under smaller stress levels it only exhibited decelerated and stationary stages (i.e., consolidation effects). During the primary stage, the strain rate of the asphalt mixture gradually decreased with increasing in loading time; this is referred to as creep hardening. Creep hardening is different from strain hardening, which refers to the phenomenon wherein the material enters plastic deformation and the strength increases with an increase in deformation. The prerequisite for strain hardening is that material stresses should reach the yield strength and produce irreversible plastic deformation, i.e., the hardening is attributed to the strengthening of the yield strength. However, even under lower stress, the asphalt mixture exhibited varying degrees of decay creep; thus, creep hardening in the decelerating stage was not attributed to strain hardening, but to the compaction effect of the asphalt mixture under compressive stress. In this stage, the asphalt mixture was gradually compacted under the compressive stress, the voids were reduced and the strength gradually increased. In the accelerated stage, a rapid increase in deformation in addition to an increase in the strain rate with time was observed because of the presence and evolution of the damage, and this eventually led to a creep failure in the asphalt mixture. The strain rate at the accelerated creep stage was related to the stress and temperature: the higher the stress, the greater the strain rate; the higher the temperature, the greater the strain rate.
Thus, based on the testing results of the creep performance, this study considers that there are two mechanisms in the entire creep process based on the creep testing results: creep hardening and creep damage deterioration. When the creep hardening mechanism is dominant, the strain rate of the asphalt mixture gradually decreases with an increase in loading time, and the creep decays. When the damage deterioration mechanism is dominant, the strain rate increases with an increase in loading time, and the creep exhibits an accelerated state. When the two mechanisms occur in proximity, the strain rate remains constant and asphalt mixture is in a stable state. The point at which the steady stage enters the accelerated stage is called the flow time. The flow time depends on the level of stress applied; the higher the stress level, the earlier the damage deterioration mechanism occurs, and thus the flow time is shorter.
4.2. Model Verification and Parameter Determination
We verify whether the proposed creep damage model can accurately describe and reproduce all creep behavior presented in the tests. Equation (16) is determined by a curve fitting to the creep test data using the Levenberg–Marquardt optimization algorithm in 1stop. The parameters of the proposed model are listed in
Table 4; the testing data of creep tests under various conditions and the simulation results of the proposed model are illustrated in
Figure 7. The parameters of the correlation coefficients R
2 are greater than 0.998, which indicates that the fractional derivative creep damage model can accurately characterize the three-stage creep process of asphalt mixtures under various conditions.
Table 4 indicates that, at lower temperatures (10 °C), when the asphalt mixtures are subjected to smaller stress levels (0.8–2.0 MPa), the deformation factor
K1, which describes the creep hardening mechanism, remains approximately constant and its value is much greater than at the other temperatures. This indicates that the deformation factor
K1 is independent of the stress levels at lower stress levels (0.8–2.0 MPa) and the asphalt mixture has greater strength at lower temperatures (10 °C). This is most likely due to the fact that asphalt mortars have a greater strength and exhibit elastic behavior at lower temperatures. Additionally, in the smaller stress range, the asphalt mixture has a smaller strain, and the compaction within the asphalt mixture is not sufficient to change the skeletal structural strength of the mixture, so the deformation factor
K1 remains constant. As only consolidation effects are exhibited within the asphalt mix at these small stresses, these smaller stress ranges (0.8–2.0 MPa at 10 °C and <0.2 MPa at 30 °C) are referred to as creep consolidation stress (CCT) ranges. The creep consolidation stress range is temperature-dependent. At higher temperatures (30 °C) the asphalt mixtures exhibit a creep consolidation only at a smaller stress level (0.2 MPa); while, at lower temperatures (10 °C), the asphalt mixtures exhibit creep consolidation behavior over a wide stress range (0.8–2.0 MPa).
In the creep consolidation state, the fractional order r increases as the stress increases and the parameters (α, K2), describing the damage mechanism, are close to zero. It follows from Equation (16) that when the deformation factor K2 tends to be zero, the damage must be sufficiently small to enable the proposed constitutive model to describe the creep properties. The strain induced by the damage is too small relative to the strain of the hardening mechanism. Therefore, the Abel model can be applied to characterize the entire creep process of an asphalt mixture under smaller stress levels. This once again shows that the Abel model can describe the consolidation effect of the asphalt mixes. Under other creep conditions, both the fractional order r and deformation factor K1 increase as the stress increases for the hardening mechanism. The fractional order increases because the flow rate of the asphalt mortar in the asphalt mixture increases as the stress increases, which causes the asphalt mixture to gradually exhibit a viscous behavior. The increase in the deformation factor indicates that the resistance to deformation increases with the increasing stress. As the stress increases in the hardening mechanism, the faster the asphalt mixture is compacted, the faster the increase in the skeletal strength between the aggregates, and the smaller the deformation. For the softening mechanism, the fractional order α increases with the increasing stress, whereas the deformation factor K2 decreases with the increasing stress. The fractional order increases because the flow rate of the asphalt mortar increases as the stress increases. The decrease in the deformation factor is attributed to the presence of damage in the softening mechanism; the greater the damage to the inter-aggregate asphalt mortar as the stress increases, and the greater the deformation. For the same stress level, asphalt mixtures exhibit different viscoelastic properties at different temperatures. For example, r increases and K1 decreases as the temperature increases for the stress of 0.8 MPa. This is because, at higher temperatures, the asphalt mortar exhibits a stronger viscosity, while at lower temperatures the asphalt exhibits a stronger elasticity. When the model is applied to describe the three stages of creep, i.e., when the hardening and deterioration mechanisms are present, the fractional order of the damage softening mechanism is larger compared to the hardening mechanisms at higher temperatures (30 °C and 50 °C). This is due to the presence of damage, which makes the asphalt mixture exhibit stronger viscosity. The deformation factor of the damage softening mechanism is also larger compared to the hardening mechanism. This is because the asphalt mixture is compacted in the primary stage, which makes its resistance to deformation increase when it enters the softening mechanism. However, at 10 °C, the parameter α describing the damage mechanism is small, indicating that the damage of the asphalt mixture at this temperature tends to occur more closely to the elastic damage.
4.3. Analysis of the Relationship between Model Parameters and Stress Levels
As can be seen from
Table 4, several parameters in the fractional derivative creep damage model are related to both loading stress levels and temperatures.
Figure 8,
Figure 9 and
Figure 10 illustrate the distribution of the parameters of various stress levels at different temperatures.
Table 4 indicates that the material parameter
v is essentially equal for different stresses at the same temperature; therefore, parameter
v at a given temperature is considered the average of the parameter
v at different stresses, as indicated in
Figure 8f,
Figure 9f and
Figure 10f. When the asphalt mixtures were subjected to linear viscoelastic stress levels, the parameters
K1 remained approximately constant, and parameters
α and
K2 were close to zero. Thus, the parameters
K1,
α, and
K2 were considered constants in the range of the linear viscoelastic stresses. Under other stress levels, fractional orders
r and
α showed a good first-order exponential decay distribution regarding the stress levels, and the deformation factors,
K1 and
K2, showed a good linear distribution with the stress levels, as indicated in
Figure 8,
Figure 9 and
Figure 10. The linear viscoelastic stress range could also be obtained using
K1, as shown in
Figure 9b and
Figure 10b. The intersection point of the two straight lines, which described the relationship between
K1 and the applied stress, was the point of demarcation between the linear viscoelastic range and the nonlinear viscoelastic range. When the stress level is less than the intersection point, the asphalt mixture is in the linear viscoelastic range and vice versa in the nonlinear viscoelastic range. The creep failure times at different stresses were fitted using Equation (12) to obtain the material parameters that describe the damage at different temperatures; the results are shown in
Figure 8e,
Figure 9e and
Figure 10e.
Based on the relationship between the model parameters and stresses given above, the creep curves could be predicted for any stress level at the three temperatures. Generally, the standard load was taken as 0.7 MPa when designing the pavement structure. Therefore, the creep behavior of the asphalt mixture was predicted and compared with the experimental data at 30 °C and 0.7 MPa, as shown in
Figure 11. It was revealed that the predicted values were in good agreement with the experimental results.
4.4. Analysis of Damage Evolution
It is believed that the fractional derivative creep damage model provides a good description of the three-stage creep process in asphalt mixtures. Once the parameters of the proposed constitutive model are identified, the proposed damage evolution model can be applied to calculate the damage value
D. The damage variation curves with a loading time at various stresses and temperatures are shown in
Figure 12.
As indicated in
Figure 12, at a given temperature, under higher stresses, the damage value D increases in an approximately linear fashion during the steady stage and rapidly during the tertiary stage. In contrast, at lower stresses, the damage value
D remains at virtually zero throughout the creep process. Regardless of the applied stresses, the asphalt mixture exhibited negligible damage during the primary stage [
16,
34]. This is because, in the initial stages of loading, the asphalt mixture is compacted, which makes the mixture structurally hardened, at which point no damage occurs to the asphalt mixture. On a macroscopic scale, asphalt mixtures show a decrease in the void fraction and an increase in the resistance to deformation; this indicates creep hardening. Over time, the asphalt mixture structure continues to stabilize, and the flow of the asphalt mortar tends to slow down. However, the asphalt mortar and coarse aggregate interface are gradually damaged and develop into microcracks. Macroscopically, the asphalt mixture exhibits a stable void ratio and constant strain rate. As the asphalt mortar flow deformation continues to accumulate, microcracks coalesce to form macro cracks, the mineral skeleton of the asphalt mixture begins to destabilize, and the asphalt mixture enters an accelerated creep phase. Macroscopically, it is characterized by cracks in the asphalt mixture and a rapid reduction in load-bearing capacity until the mixture breaks down in the creep phase; in general, damage to the asphalt mixture is initiated in the steady stage and developed in the tertiary stage, which is consistent with the works reported by Zeng et al. [
16] and Al-rub et al. [
34]. The experimental data and model damage analysis show the feasability of the damage evolution proposed in this study.
Equation (11) considers only the time dependence of the damage evolution, assuming that the damage only develops over time; it does not consider the relationship between the damage evolution and the relevant mechanical quantities such as stress or strain. Equations (11) and (16) show that each moment t corresponds to a damage value and a creep strain value. Therefore, through creep tests, the relationship curve between the damage value and strain can be established for different stress levels. At different stress levels, the damage values corresponding to the strain values at each moment of the decelerated stage are small and essentially close to zero. The strain caused by the damage (damage strain,
εd) is obtained by subtracting the strain in the decelerated stage from the total creep strain in the tests. According to Equation (16),
εd can be expressed as:
Based on the time dependence of the damage evolution and damage values obtained from the creep experiments, the relationship curves between the damage value
D and
εd for different stress levels at different temperatures is established in this study, as shown in
Figure 13.
Figure 13 shows that the relationship curves between the creep
D and the
εd of the asphalt mixture at different stress levels overlap. Therefore, a unified damage evolution model between
D and
εd can be developed, where the damage parameters in the damage evolution relationship are only material- and temperature-dependent, independent of the magnitude of the stress levels. Katsuki et al. [
4] used the Weibull distribution function to describe the variation in the damage with the strain, and verified the feasibility of the model using the experimental computed tomography images from Tashman at al. [
35]. Therefore, in this study, the Weibull distribution function is used to describe the evolution of the damage with
εd, and the model expression is shown in Equation (18). By adopting the proposed damage model (Equation (18)) to fit the experimental data under different levels of stress for different temperatures in
Figure 13, the damage model curve shown in
Figure 13 is obtained. The parameters are listed in
Table 5. The parameters of the correlation coefficient R
2 are greater than 0.992, which shows that the Weibull distribution function can be used to describe the damage evolution and demonstrate the feasibility of establishing a unified damage evolution model for different stress levels. Statistically speaking, the microelement strength of randomly distributed microdefects within the material obeys the Weibull distribution under
εd:
where
p and
q are temperature-dependent material parameters, respectively.
The nonlinear damage constitutive model proposed in this study has the advantages of fewer parameters and a clear physical meaning of parameters relative to the creep model proposed by Zeng et al. [
16]. Relative to the model proposed by Zhang et al. [
19], there is no need to assume that there is critical stress. Additionally, relative to other damage models [
16,
17,
19], this study decouples the damage strain from the creep strain and establishes a unified damage evolution model for different stress levels, whose parameters are only related to material properties and temperature, and are independent of the magnitude of the applied stress.