Application of the Viscoelastic Continuum Damage Theory to Study the Fatigue Performance of Asphalt Mixtures—A Literature Review

Abstract

A viscoelastic continuum damage (VECD) model, which accounts for the effects of rate-dependent damage growth, has been successfully applied to quantify the reduction in the material integrity as a function of damage accumulation (characteristic curve) of asphalt mixtures. This paper presents an overview of the fundamentals of the VECD model and its simplifications (S-VECD) applied to the damage characterization of asphalt mixtures. It also presents a laboratory study carried out to evaluate the effects of the addition of reclaimed asphalt pavements (RAP) and a new binder (PG 64-22 or PG 58-16) on the fatigue performance of fine aggregate matrices (FAMs), in which the S-VECD theory was applied to analyze the results. The addition of RAP increased the stiffness and reduced the relaxation rate, resulting in FAMs that were stiffer and more susceptible to damage at high strain levels. The FAMs’ fatigue factors (FF_FAM) indicated that the increase in RAP from 20% to 40% decreased the fatigue life of the mixtures. A strict control of the mixture variables is required, since the intrinsic heterogeneity of asphalt mixtures can lead to different characteristic curves for the same material.

1. Introduction

An accurate prediction of the fatigue behavior of an asphalt mixture has been the goal of many studies focused on improving the asphalt mixture’s design and the flexible pavement’s performance [1,2,3,4,5,6,7,8,9]. Phenomenological models, which relate the stress or strain in the specimen with the number of cycles to failure, are simple tools to determine the fatigue behavior of an asphalt mixture in the laboratory [1,10,11,12,13]. However, this approach does not account for the complexity of the fatigue phenomenon [3]. More recently, mechanistic approaches have been employed rather than the phenomenological ones, as mechanistic models account for how damage evolves throughout the fatigue life at different loading and environmental conditions, leading to a better estimation of the fatigue behavior of the asphalt mixture [2,3,5,7].

Studies on the response of asphalt mixtures to fatigue cracking are divided into two categories: (1) the full asphalt mixture, which comprised asphalt binder, coarse aggregate, fine aggregate, and mineral filler; and (2) the fine aggregate matrix (FAM), composed of fine aggregate, mineral filler, and asphalt binder [14]. Based on the premise that the changes in the material microstructure are the beginning of the fatigue process [4,14,15], many researchers have studied the fatigue behavior of asphalt mixtures using the FAM approach [5,6,14,16,17,18]. The fine aggregate matrix represents an intermediate scale between the asphalt mastic and the full asphalt mixture, presenting an internal structure that is not affected by the coarse aggregate particles [19]. Studies at the FAM scale might be capable of providing a more realistic characterization of the fatigue response of full asphalt mixtures than the one provided by tests performed on the mastic [14]. Another advantage of using FAMs in fatigue testing is that the reduced size of the specimens requires a smaller amount of material, as compared with the amount needed to produce specimens of full asphalt mixtures, which also reduces the laboratory work needed in the preparation of specimens [20].

Asphalt concrete is a material with high viscoelasticity imparted by the binder matrix [21]. The data obtained in tests performed with both FAMs and full asphalt mixtures can be analyzed by means of the continuum mechanics theory. According to the viscoelastic continuum damage (VECD) theory, the time dependency of viscoelastic materials can be eliminated by means of correspondence principles, which transform physical variables (stress, strain, and stiffness) in pseudo variables (pseudo stress, pseudo strain, and pseudo stiffness) [15,22,23]. According to this theory, a damaged body presenting internal microcracks is assumed to be an undamaged body with a reduced pseudo stiffness (C), with the microcracks uniformly distributed within the body [22,23]. The material damage evolution is described as the function C(S), in which a reduction in pseudo stiffness (C) is related to a material internal state variable (S) [15,21,22,23].

This review paper aims to provide an overview of the fundamentals of the viscoelastic continuum damage theory applied to the damage characterization of asphalt mixtures at the scales of both full asphalt mixture and fine aggregate matrix and its simplifications, along with an overview of the mechanistic fatigue life prediction model based on continuum mechanics. In addition to this review, this paper also presents examples of laboratory tests performed to evaluate the effect of reclaimed asphalt pavements and recycling agents on the fatigue response of FAMs. Although the theory presented in this paper assumes that each material presents a single characteristic curve that is independent of the loading conditions, laboratory tests sometimes result in different characteristic curves for specimens of the same FAM. Material heterogeneity is supposed to be the reason for distinct characteristic curves for specimens of the same FAM. In an attempt to reduce the material variability, the method adopted in this study employed the average of the linear-viscoelastic properties of at least three specimens to build a unique characteristic curve of the material, following the procedure employed to perform damage tests on full asphalt mixture specimens. Figure 1 provides an overview of the topics covered in this paper.

2. Theory of Viscoelastic Continuum Damage

The work potential theory, which is based on the methods of thermodynamics of irreversible process, was developed by Schapery [15,16,17,18,19,20,21,22,23] in order to describe the mechanical behavior of elastic materials with growing damage. The theory characterizes the material using macroscale observations, quantifying the changes of the material microstructure by means of internal state variables (S). The elastic model was extended to describe the mechanical behavior of viscoelastic media, by means of elastic–viscoelastic correspondence principles that eliminate the time dependence of the viscoelastic material.

Park et al. [21], Lee [24], and Lee and Kim [25] applied Schapery’s theory to develop a constitutive model that describes the damage evolution process of asphalt mixtures for different materials, as well as loading and environmental conditions. This constitutive model was simplified by Lee et al. [7] to build a practical fatigue prediction model for specimens of asphalt mixtures under uniaxial cyclic loading. Such a model was later adapted by Kim and Little [5] for specimens of sand asphalt mixtures under torsional shear without rest periods. This improved model is considered capable of providing a reasonable representation of the fatigue life of asphalt mastics and fine aggregate matrices [5].

2.1. The Work Potential Theory

According to the thermodynamic theory, a generalized notation of the behavior of an elastic body that changes its structure is expressed by means of relationships between generalized forces, Q_j, and independent generalized displacements, q_j, as shown in Equation (1), where δq_j is the virtual displacement and δW’ is the virtual work. For different physical situations, q_j can represent strain, displacement, or rotation; and Q_j can be stress, force, or moment. For any process of interest, the existence of a strain energy function, W = W(q_j, S_m), was assumed, where S_m (m = 1, 2, 3, M) refers to the increase in the value of the internal state variable, S. The relationship between the work done on a body during the process in which damage occurs and the strain energy function is expressed by Equation (2), where f_m is the thermodynamic force (Equation (3)).

$${\mathsf{\delta }\mathrm{W}}^{\prime }{=\mathrm{Q}}_{\mathrm{j}}{\mathsf{\delta }\mathrm{q}}_{\mathrm{j}}$$

(1)

$$\mathrm{dW}=\frac{\partial \mathrm{W}}{\partial {\mathrm{q}}_{\mathrm{j}}}{\mathrm{dq}}_{\mathrm{j}}+\frac{\partial \mathrm{W}}{\partial {\mathrm{S}}_{\mathrm{m}}}{\mathrm{dS}}_{\mathrm{m}}{=\mathrm{Q}}_{\mathrm{j}}{\mathrm{dq}}_{\mathrm{j}}-{\mathrm{f}}_{\mathrm{m}}{\mathrm{dS}}_{\mathrm{m}}$$

(2)

$${\mathrm{f}}_{\mathrm{m}}=-\frac{\partial \mathrm{W}}{\partial {\mathrm{S}}_{\mathrm{m}}}$$

(3)

By integrating Equation (2) during a time interval t₁–t₂, and assuming that a state function W_S = W_S(S_m) exists such that the thermodynamic force, f_m, is given by Equation (4), the work to vary the internal state of the material from a state 1 to a state 2 is expressed by Equation (5). The variable $\dot{{\mathrm{S}}_{\mathrm{m}}}$ is the damage evolution rate. By solving the integral in Equation (5), the work is given by Equation (6). Assuming that the time t₁ = 0, the total work from t = 0 to the current time t₂ is given by Equation (7).

$${\mathrm{f}}_{\mathrm{m}}=\frac{\partial {\mathrm{W}}_{\mathrm{S}}}{\partial {\mathrm{S}}_{\mathrm{m}}}\mathrm{when}\dot{{\mathrm{S}}_{\mathrm{m}}}\ne 0$$

(4)

$${\mathsf{\Delta }\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}}^{\left(2\right)}-{\mathrm{W}}^{\left(1\right)}+{{\displaystyle \int }}_1^2{\mathrm{f}}_{\mathrm{m}}{\mathrm{dS}}_{\mathrm{m}}$$

(5)

$${\mathsf{\Delta }\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}}^{\left(2\right)}-{\mathrm{W}}^{\left(1\right) }{+\mathrm{W}}_{\mathrm{S}}{}^{\left(2\right)}-{\mathrm{W}}_{\mathrm{S}}{}^{\left(1\right)}$$

(6)

$${\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}+\mathrm{W}}_{\mathrm{S}}$$

(7)

The total work input to the elastic body by the forces Q_j is given by Equation (8), where j = 1, 2, …, J, and W_T is the total work done on the body while considering that S_m is variable in time. The elements of Schapery’s theory, expressed in terms of stress–strain relationships, can be represented as follows in Equations (9) and (10), where σ_ij is the stress tensor, ε_ij is the strain tensor, and S_m is the internal state variable. By considering Equations (3) and (4), the damage evolution law for the elastic media is represented by Equation (11), in which W_S = W_S(S_m) is the dissipated energy due to damage growth. The right-hand side of the law represents the required force for damage growth, while the left-hand side of the damage evolution law represents the available thermodynamic force to produce damage growth.

$${\mathrm{W}}_{\mathrm{T}}={{\displaystyle \int }}^{}{\mathrm{Q}}_{\mathrm{j}}{\mathrm{dq}}_{\mathrm{j}}$$

(8)

$${\mathrm{W}=\mathrm{W}(\mathsf{\epsilon }}_{\mathrm{ij}}{,\mathrm{S}}_{\mathrm{m}})$$

(9)

$${\mathrm{W}}_{\mathrm{T}}={{\displaystyle \int }}^{}{\mathsf{\sigma }}_{\mathrm{ij}}{\mathrm{d}\mathsf{\epsilon }}_{\mathrm{ij}}$$

(10)

$$-\frac{\partial \mathrm{W}}{\partial {\mathrm{S}}_{\mathrm{m}}}=\frac{\partial {\mathrm{W}}_{\mathrm{S}}}{\partial {\mathrm{S}}_{\mathrm{m}}}$$

(11)

2.2. Elastic–Viscoelastic Correspondence Principle

As demonstrated by Schapery [22], by using correspondence principles and a new damage evolution law, the elastic continuum damage theory can be extended to describe damage evolution in viscoelastic materials. For elastic materials, the stress–strain relationship is expressed by Hooke’s law (Equation (12)), where σ is the stress, E is the elasticity modulus, and ε is the strain. For viscoelastic materials, the time dependency has to be considered, and the stress is expressed by a convolution integral (Equation (13)), where τ is an increment in the value of time t, G(t) is the material relaxation modulus, and ε is the time-dependent strain in a viscoelastic material.

$$\mathsf{\sigma }=\mathrm{E}\mathsf{\epsilon }$$

(12)

$$\mathsf{\sigma }={{\displaystyle \int }}_0^{\mathrm{t}}\mathrm{G}(\mathrm{t}-\mathsf{\tau })\frac{\mathrm{d}\mathsf{\epsilon }}{\mathrm{d}\mathsf{\tau }}\mathrm{d}\mathsf{\tau }$$

(13)

The artifice proposed to eliminate the time-dependent effects was to transform the stress–strain relationships of the viscoelastic media into a pseudo domain that corresponds to a hypothetical elastic material, suggesting that the constitutive equation for viscoelastic media (Equation (14)) is identical to that for elastic media (Equation (12)). However, the stress and the strain are not necessarily physical quantities, but pseudo variables: pseudo stress (σ^R) and pseudo strain (ε^R). According to the second correspondence principle (CP II), σ^R = σ, where σ is the time-dependent stress applied to a viscoelastic material; and the pseudo strain is given by Equation (15), where ε is the time-dependent strain in a viscoelastic material, G(t) is the linear-viscoelastic relaxation modulus of the material, and E^R is the modulus of hypothetical elastic material [22,23].

$${\mathsf{\sigma }=\mathrm{E}}_{\mathrm{R}}{\mathsf{\epsilon }}^{\mathrm{R}}$$

(14)

$${\mathsf{\epsilon }}^{\mathrm{R}}=\frac{1}{{\mathrm{E}}_{\mathrm{R}}}{{\displaystyle \int }}_0^{\mathrm{t}}\mathrm{G}(\mathrm{t}-\mathsf{\tau })\frac{\mathrm{d}\mathsf{\epsilon }}{\mathrm{d}\mathsf{\tau }}\mathrm{d}\mathsf{\tau }$$

(15)

The same equations for elastic materials are used to solve viscoelastic cases by considering elastic–viscoelastic correspondence principles. Equation (9) is then expressed by Equation (16), the pseudo strain energy density (W^R) function, where the physical strain, ε, is substituted by the pseudo strain, ε^R, and S_m is the internal state variable. The stress-pseudo strain relationship is given by Equation (17), where $\mathsf{\sigma }$ is the stress, W^R is the pseudo strain energy density, and ε^R is the pseudo strain [23]. For most viscoelastic media, the available force for damage growth and the resistance against the growth are rate-dependent. For this reason, the damage evolution law for elastic materials (Equation (11)) cannot just be transformed into a damage evolution law for viscoelastic materials by the use of correspondence principles without further modification. The new damage evolution law for viscoelastic materials is given by Equation (18), where $\dot{{\mathrm{S}}_{\mathrm{m}}}$ is the damage evolution rate, W^R is the pseudo strain energy density, S_m is the internal state variable, and α_m is a material-dependent constant [26].

$${\mathrm{W}}^{\mathrm{R}}{=\mathrm{W}}^{\mathrm{R}}{(\mathsf{\epsilon }}^{\mathrm{R}}{,\mathrm{S}}_{\mathrm{m}})$$

(16)

$$\mathsf{\sigma }=\frac{\partial {\mathrm{W}}^{\mathrm{R}}}{\partial {\mathsf{\epsilon }}^{\mathrm{R}}}$$

(17)

$$\dot{{\mathrm{S}}_{\mathrm{m}}}={\left(-\frac{\partial {\mathrm{W}}^{\mathrm{R}}}{\partial {\mathrm{S}}_{\mathrm{m}}}\right)}^{{\mathsf{\alpha }}_{\mathrm{m}}}$$

(18)

2.3. Viscoelastic Continuum Damage Model Applied to Asphalt Mixtures

Kim and Little [27] successfully applied Schapery’s [15,22] theory on the development of a nonlinear-viscoelastic constitutive equation that should represent the damage growth in asphalt concrete. The equation was developed from tests on mixtures of asphalt concrete composed of crushed granite fines and asphalt binder, and it is presented by Equation (19), where σ = stresses in a body; ε^R = strains in a body (Equation (15)); ${\mathsf{\epsilon }}_{\mathrm{L}}^{\mathrm{R}}$ = maximum pseudo strain in the past history; and S_p = damage parameter [28]. The form of the damage parameter based on pseudo strain is shown as follows in Equation (20), where p = (1 + N)k; N = the exponent of the power law between stress and strain, σ~|ε^R|^N; k = 2(1 + 1/m); m = the exponent of the power law between creep compliance and time; and D(t) = D₁t^m. When repetitive loading is applied, numerical integration can be used to obtain S_p, as shown in Equation (21), assuming that dε^R/dt is constant within the range of the experimental data points. A uniaxial tensile testing was employed to generate all data for the construction of the constitutive equation, which satisfactorily predicted effects due to multilevel loading, the sequence of multilevel loading, and various durations of rest periods. From the tests, it was proved that the history dependence of asphalt concrete with negligible damage growth can successfully be eliminated by the CP II, σ^R = σ.

$${\mathsf{\sigma }=\mathsf{\sigma }(\mathsf{\epsilon }}^{\mathrm{R}}{,\mathsf{\epsilon }}_{\mathrm{L}}^{\mathrm{R}}{,\mathrm{S}}_{\mathrm{p}})$$

(19)

$${\mathrm{S}}_{\mathrm{p}}={\left(\underset{0}{\overset{\mathrm{t}}{{\displaystyle \int }}}{\left|{\mathsf{\epsilon }}^{\mathrm{R}}\right|}^{\mathrm{p}}\mathrm{dt}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathrm{p}$}\right.}$$

(20)

$${{[\mathrm{S}}_{\mathrm{p}}\left(\mathrm{t}\right)]}^{\mathrm{p}}={{\displaystyle \int }}_0^{{\mathrm{t}}_{\mathrm{j}-1}}{\left|{\mathsf{\epsilon }}^{\mathrm{R}}\right|}^{\mathrm{p}}\mathrm{dt}-{{\displaystyle \int }}_{{\mathrm{t}}_{\mathrm{j}-1}}^{\mathrm{t}}{\left|{\mathsf{\epsilon }}^{\mathrm{R}}\right|}^{\mathrm{p}}\mathrm{dt}$$

(21)

Park et al. [21] also proposed a uniaxial viscoelastic damage model to characterize the behavior of asphalt mixtures with time-dependent damage growth, subjected to different strain rates under uniaxial stress. The pseudo strain energy density (W^R) function is given by Equation (22), where C is a function of the damage parameter S, and ε^R is the pseudo strain. For linear-viscoelastic behavior and fixed damage, the stress, σ, can be written as in Equation (23). The damage evolution law (Equation (18)) is reduced to the single equation for S (Equation (24)), where S is the damage parameter, α is a material-dependent constant, and the over dot denotes a time derivative.

$${\mathrm{W}}^{\mathrm{R}}=\frac{1}{2}\mathrm{C}\left(\mathrm{S}\right){\left({\mathsf{\epsilon }}^{\mathrm{R}}\right)}^2$$

(22)

$$\mathsf{\sigma }\equiv \frac{\partial {\mathrm{W}}^{\mathrm{R}}}{\partial {\mathsf{\epsilon }}^{\mathrm{R}}}{=\mathrm{C}\left(\mathrm{S}\right)\mathsf{\epsilon }}^{\mathrm{R}}$$

(23)

$$\dot{{\mathrm{S}}_{\mathrm{m}}}={\left(-\frac{\partial {\mathrm{W}}^{\mathrm{R}}}{\partial \mathrm{S}}\right)}^{\mathsf{\alpha }}$$

(24)

By substituting Equation (22) into Equation (24), Park et al. [21] obtained a relationship between S and a new damage parameter, S∗, given by Equation (25), where S∗ is a function of the strain history (Equation (26)). The variable k is a free constant considered in their study, such that the maximum values of S and S∗ are numerically equal. Thus, Equation (23) was replaced by Equation (27). In Equations (25)–(27), α is a material-dependent constant, ε^R is the pseudo strain, and C is a function of the damage parameter S.

$${\mathrm{S}}^{\ast }=\mathrm{k}{\left[{{\displaystyle \int }}_0^{\mathrm{S}}\frac{\mathrm{dS}}{{\left(-0.5\mathrm{dC}/\mathrm{dS}\right)}^{\mathsf{\alpha }}}\right]}^{1/\left(2\mathsf{\alpha }\right)}$$

(25)

$${\mathrm{S}}^{\ast }\equiv \mathrm{k}{\left[{{\displaystyle \int }}_0^{\mathrm{t}}{|\mathsf{\epsilon }}^{\mathrm{R}}{|}^{2\mathsf{\alpha }}\mathrm{dt}\right]}^{1/\left(2\mathsf{\alpha }\right)}$$

(26)

$${\mathsf{\sigma }=\mathrm{C}(\mathrm{S}}^{\ast }{)\mathsf{\epsilon }}^{\mathrm{R}}$$

(27)

The function C(S) and the α value must be determined. However, by fitting Equation (23) to an experimental stress–pseudo strain curve, it is possible to obtain a C modulus that depends on ε^R and strain rate, but not on S. In order to find the dependence of the C modulus on S, the use of the damage evolution law (Equation (18)) is required. However, the evolution law itself requires prior knowledge of C(S), making this procedure inefficient to find C and its dependence on S. The method proposed by Park et al. [21] to overcome this problem was to determine a transformed damage variable, Ŝ (Equation (28)), that may be obtained from the numerical scheme presented in Equation (29), where ${\mathsf{\epsilon }}_{\mathrm{i}}^{\mathrm{R}}$(i = 1, 2, 3, N) denotes pseudo strain levels, C’ ≡ dC/dŜ, and Ŝ(0) = 0. This method allows one to obtain the function C(Ŝ) from experimental stress–pseudo strain curves, and then the function C(S) can be obtained from Equation (28), by replacing the transformed damage variable Ŝ for the original damage variable, S, where α is a material-dependent constant.

$$\mathsf{Ŝ}\mathsf{\equiv }\frac{1}{(1+1/\alpha )}{\mathrm{S}}^{(1+1/\alpha )}$$

(28)

$$\mathsf{Ŝ}\left({\mathsf{\epsilon }}_{\mathrm{i}+1}^{\mathrm{R}}\right)=\mathsf{Ŝ}\left({\mathsf{\epsilon }}_{\mathrm{i}}^{\mathrm{R}}\right)+\frac{\mathrm{C}\left({\mathsf{\epsilon }}_{\mathrm{i}+1}^{\mathrm{R}}\right)-\mathrm{C}\left({\mathsf{\epsilon }}_{\mathrm{i}}^{\mathrm{R}}\right)}{\Delta \mathsf{Ŝ}}$$

(29)

The constant α, related to the material’s creep or relaxation properties, also must be determined. According to Schapery [15], depending on the characteristics of the failure zone on a crack tip, α is expressed by Equation (30), if the material’s fracture energy and failure stress are constant, or α is expressed by Equation (31), if the fracture process zone size and the material’s fracture energy are constant, where the material relaxation rate, m, is given by Equation (32), where D(t) = uniaxial creep compliance, E(t) = uniaxial relaxation modulus, and t = time. Park et al. [21] employed the relationship presented in Equation (30) by assuming that the damage in the specimen is closely related to the growth of microcracks. The constant α was obtained by a method of successive approximations until the model better fit the experimental observations.

$$\mathsf{\alpha }=(1+\frac{1}{\mathrm{m}})$$

(30)

$$\mathsf{\alpha }=\frac{1}{\mathrm{m}}$$

(31)

$$\mathrm{m}=\frac{\mathrm{logD}\left(\mathrm{t}\right)}{\log \left(\mathrm{t}\right)}\mathrm{or}\mathrm{m}=\frac{-\mathrm{logE}\left(\mathrm{t}\right)}{\log \left(\mathrm{t}\right)}$$

(32)

Lee and Kim [24,25] also proposed a solution to the damage evolution law by studying the mechanical behavior of asphalt concrete. They conducted uniaxial tensile cyclic loading tests with different loading amplitudes, under controlled-strain and controlled-stress modes of loading. The researchers observed a decrease in the slope of each σ–ε^R cycle as the number of loading increased, and they found it necessary to define the secant pseudo stiffness, S^R, to represent this change in the slope of stress–pseudo-strain loops (Equation (33)), where ${\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}$ is the peak pseudo strain in each stress–pseudo-strain cycle and σ_m is the stress that corresponds to ${\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}$. To minimize the sample-to-sample variability, the pseudo stiffness was divided by the initial pseudo stiffness, I, resulting in the normalized pseudo stiffness, C, represented by Equation (34). By considering Equation (33) and Equation (34), the constitutive equation for viscoelastic materials with growing damage is expressed by Equation (35). The normalized pseudo stiffness, C(S_m), is a function of the internal state variables, S_m, and represents the microstructural changes of the body. The researchers assumed an internal state variable, S₁, to determine the change in pseudo stiffness due to growing damage, and the work function (W^R) for viscoelastic materials is given by Equation (36), where C₁(S₁) is a function that represents S^R.

$${\mathrm{S}}^{\mathrm{R}}=\frac{{\mathsf{\sigma }}_{\mathrm{m}}}{{\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}}$$

(33)

$$\mathrm{C}=\frac{{\mathrm{S}}^{\mathrm{R}}}{\mathrm{I}}$$

(34)

$${\mathsf{\sigma }=\mathrm{IC}(\mathrm{S}}_{\mathrm{m}}{)\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}$$

(35)

$${\mathrm{W}}^{\mathrm{R}}=\frac{\mathrm{I}}{2}{\mathrm{C}}_1\left({\mathrm{S}}_1\right){\left({\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}\right)}^2$$

(36)

Despite that the material function C₁ can be determined by using experimental data and the damage evolution law (Equation (18)), this procedure is not convenient to find C₁ and its dependence on S₁, because the evolution law requires prior knowledge of C₁(S₁). The method presented to overcome this problem was to use a chain rule (Equation (37)) to eliminate the S on the right-hand side of the evolution law, and by means of mathematical substitutions (Equation (38)) the numerical approximation given by Equation (39) was obtained, where ε^R is the pseudo strain, t is the time, I is a factor to normalize the pseudo stiffness, and α is a material-dependent constant (Equations (30) and (31)). The function C₁(S₁) can be obtained by cross-plotting the C values obtained from Equation (35) against the S values obtained from Equation (39), and by performing a regression on the data (Equation (40)), where C₁₀, C₁₁, and C₁₂ are regression coefficients; or an exponential function (Equation (41)), where a and b are the calibration constants.

$$\frac{\mathrm{dC}}{\mathrm{dS}}=\frac{\mathrm{dC}}{\mathrm{dt}}\frac{\mathrm{dt}}{\mathrm{dS}}$$

(37)

$$\frac{\mathrm{dS}}{\mathrm{dt}}={\left[-\frac{1}{2}\frac{\mathrm{dC}}{\mathrm{dS}}{\left({\mathsf{\epsilon }}^{\mathrm{R}}\right)}^2\right]}^{\raisebox{1ex}{$\mathsf{\alpha }$}\!\left/ \!\raisebox{-1ex}{$\left(1+\mathsf{\alpha }\right)$}\right.}$$

(38)

$$\mathrm{S}\equiv {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left[\frac{\mathrm{I}}{2}{\left({\mathsf{\epsilon }}^{\mathrm{R}}\right)}^2\left({\mathrm{C}}_{\mathrm{i}-1}-{\mathrm{C}}_{\mathrm{i}}\right)\right]}^{\mathsf{\alpha }/(1+\mathsf{\alpha })}{\left({\mathrm{t}}_{\mathrm{i}}-{\mathrm{t}}_{\mathrm{i}-1}\right)}^{1/(1+\mathsf{\alpha })}$$

(39)

$${\mathrm{C}}_1\left({\mathrm{S}}_1\right){=\mathrm{C}}_{10}-{\mathrm{C}}_{11}{{(\mathrm{S}}_1)}^{{\mathrm{C}}_{12}}$$

(40)

$${\mathrm{C}}_1\left({\mathrm{S}}_1\right){=\exp (\mathrm{aS}}_1{}^{\mathrm{b}})$$

(41)

For the material parameter α, Lee and Kim [29] observed that Equation (30) was adequate for the controlled strain mode of loading, while Equation (31) was a better assumption for the controlled stress mode of loading. These observations suggested that the material’s fracture energy and failure stress are constant under the controlled strain mode, whereas the material’s fracture energy and the fracture process zone size are constant under the controlled stress mode.

2.4. Pseudo-Strain Calculation

The constant cyclic shear strain, ε(t), as a function of time, can be simply represented as an analytical harmonic sinusoidal function (Equation (42)), where ε₀ is the shear strain amplitude, ω is the angular velocity, θ is a regression constant, and H(t) is the Heaviside step function. Substituting Equation (42) within the definition of pseudo strain (Equation (15)) and assuming that E_R is 1, the pseudo strain at the current time can be analytically represented by Equation (43), where ϕ is the phase angle and |G*| is the linear-viscoelastic dynamic shear modulus. As can be seen, the pseudo strain at any time can be predicted with a well-defined strain history as a function of time and two material properties: dynamic modulus and phase angle. In the study of fatigue, only the peak pseudo strain within each cycle is used, and the pseudo strain reaches the peak pseudo strain in each cycle when the sine function in Equation (43)) becomes 1, as represented by Equation (44)), where ${\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}$ is the peak pseudo strain at each cycle [8,30].

$$\mathsf{\epsilon }\left(\mathrm{t}\right)=\mathsf{\epsilon }0\sin (\mathsf{\omega }\mathrm{t}+\mathsf{\theta })\mathrm{H}\left(\mathrm{t}\right)$$

(42)

$${\mathsf{\epsilon }}^{\mathrm{R}}\left(\mathrm{t}\right)=\mathsf{\epsilon }0\left|{\mathrm{G}}^{\ast }\right|\sin (\mathsf{\omega }\mathrm{t}+\mathsf{\theta }+\varphi )$$

(43)

$${\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}\left(\mathrm{t}\right)=\mathsf{\epsilon }0\left|{\mathrm{G}}^{\ast }\right|$$

(44)

2.5. Simplified Viscoelastic Continuum Damage Model

Different researchers have worked on developing a simplified viscoelastic continuum damage (S-VECD) model in order to characterize asphalt mixtures in a quick and easy way using cyclic fatigue test data [9,30,31,32,33,34,35]. The modification in the continuum damage equations [8,15,21,22,23,24,25,27,29] proposed by Christensen and Bonaquist [31] can be used to predict the fatigue behavior of asphalt concrete mixtures based on their volumetric composition and degree of compaction, and the rheological type of the asphalt binder used in the mixture. In a later study, Christensen and Bonaquist [32] presented two new concepts for inclusion in the continuum damage analysis of fatigue data on hot-mix asphalt (HMA) mixtures: the concept of reduced loading cycles and the concept of effective strain. Kutay et al. [30] developed and validated a simple equation for calculation of the continuum damage parameter (S) for the specific case of fixed-frequency cyclic tests. Underwood et al. [33,34] reviewed and discussed the limitations of the simplified models proposed by Kutay et al. [30], and Christensen and Bonaquist [31] suggested an improved S-VECD model.

The original application of the VECD model to cyclic data requires that the pseudo strain, pseudo stiffness, and damage are calculated and tracked for the entire loading history. Although the analysis of tests with a large number of data points is not impossible with modern computers, it depends on advanced computational schemes. Further, experimental difficulties related to data storage and electrical interference (noise and phase distortion) can lead to significant errors. In an effort to minimize these shortcomings, simplified mechanics analyses were developed by researchers at North Carolina State University (NCSU) and by the Federal Highway Administration (FHWA). In the methods developed at NCSU, the pseudo strain is calculated for the entire loading history (Equation (15)), and present the same sequence of calculations for damage during the first load path (Equation (39)). This early portion of the damage calculation is referred to as the transient calculation. After the first loading cycle, however, the sequence of calculations is referred to as cyclic calculations, and a simplified calculation method is used. The NCSU simplified methods also consider the specimen-to-specimen variation by the I factor, which is defined by the initial slope of the stress–pseudo-strain plot obtained in a stress-controlled test within the linear–viscoelastic (LVE) range. In some of the NCSU simplified models, the permanent pseudo strain (${\mathsf{\epsilon }}_{\mathrm{S}}^{\mathrm{R}}$) is also taken into account for the damage calculation [33,34]. The main difference between the NCSU S-VECD models and the one proposed by the FHWA is the simplification of the pseudo strain calculation (Equation (45)). This simplification is used as reference for the steady-state assumption, because it is considered theoretically accurate for the zero-mean-stress steady-state conditions. In addition, in the FHWA S-VECD model, the pseudo strain is calculated based on the peak-to-peak values, and assuming that the damage growth occurs during the tension and compression load [33,34].

$${\mathsf{\epsilon }}_0^{\mathrm{R}}{=\mathsf{\epsilon }}_0{\left|{\mathrm{G}}^{\ast }\right|}_{\mathrm{LVE}}$$

(45)

Characteristic curves of five asphalt mixtures built by using the NCSU S-VECD models and the FHWA S-VECD model were compared. Based on the results of this comparison, a new simplified model was developed while considering the advantages of each S-VECD model [33,34]. Taking into account the modified correspondence principle [22] and zero-mean-stress steady-state condition considered in the FHWA S-VECD, a new S-VECD model was proposed [33,34]. The relationships for pseudo-strain amplitude, pseudo-stress amplitude, and pseudo stiffness at the k^th load cycle for a stress-controlled cyclic test were simplified, as shown in Equations (46) and (47). In this FHWA S-VECD method, damage is expected to grow according to Equation (48).

$${\mathsf{\sigma }}_{\mathfrak{0}\mathfrak{,}\mathrm{k}}^{\mathrm{P}}{\mathsf{\sigma }}_{\mathfrak{0}\mathfrak{,}\mathrm{k}}$$

(46)

$${\mathrm{C}}_{\mathrm{k}}\left(\mathrm{S}\right)\equiv \frac{{\mathsf{\sigma }}_{0,\mathrm{k},\mathrm{pp}}^{\mathrm{R}}}{\mathrm{I}·{\mathsf{\epsilon }}_{0,\mathrm{k},\mathrm{pp}}^{\mathrm{R}}}=\frac{{\left|{\mathrm{G}}^{\ast }\right|}_{\mathrm{k}}}{\mathrm{I}·{\left|{\mathrm{G}}^{\ast }\right|}_{\mathrm{LVE}}}$$

(47)

$$\mathrm{dS}\equiv {\left[\frac{\mathrm{I}}{2}{\left({\mathsf{\epsilon }}_{0,\mathrm{k},\mathrm{pp}}^{\mathrm{R}}\right)}^2\left({\mathrm{C}}_{\mathrm{k}-1}-{\mathrm{C}}_{\mathrm{k}}\right)\right]}^{\mathsf{\alpha }/(1+\mathsf{\alpha })}{\left({\mathrm{t}}_{\mathrm{k}}-{\mathrm{t}}_{\mathrm{k}-1}\right)}^{1/(1+\mathsf{\alpha })}$$

(48)

The parameters ${\mathsf{\sigma }}_{0,\mathrm{k}}$, ${\mathsf{\sigma }}_{0,\mathrm{k}}^{\mathrm{R}}$, ${\mathsf{\epsilon }}_{0,\mathrm{k}}^{\mathrm{R}}$, |G*|_k, and |G*|_LVE| represent, respectively, the shear stress amplitude, the shear pseudo-stress amplitude, the shear pseudo-strain amplitude, the dynamic shear modulus at the k^th loading cycle, and the linear-viscoelastic dynamic shear modulus. The comma subscript, pp, denotes the pseudo strain computed based on peak-to-peak values. The I factor is used as a correction factor in order to normalize test results while considering the sample-to-sample variability of the initial dynamic shear modulus. Hou [9] conducted a study in order to verify the S-VECD model proposed by Underwood et al. [33] by applying it to various types of asphalt concrete mixtures under various conditions. Hou [9] observed that the model could be applied to accurately predict the fatigue life of asphalt concrete under cyclic loading at multiple temperatures and strain levels. Lee et al. [36] recommended that the power law model (Equation (40)) should be fitted to the characteristic curves obtained from the S-VECD model after filtering the data to produce damage (S) increments of 5000. Lee et al. [36] also proposed a shape factor, which is calculated as the ratio of the area above the damage characteristic curve for an individual test replicate (i.e., A_measured) to the area above the fitted model damage characteristic curve obtained by Equation (40) (i.e., A_predicted). Replicates with shape factors greater than 1.1 or lower than 0.9 are considered outliers and should not be considered for analysis.

2.6. Mechanistic Fatigue Life Prediction Model

There are two main approaches to characterize the material fatigue behavior in the laboratory: the phenomenological approach and the mechanistic approach. The phenomenological models relate the stress or strain in the specimen with the number of cycles to failure. This approach is relatively simple, but it does not account for the complexity of the fatigue phenomenon. On the other hand, the mechanistic approach accounts for how the damage evolves throughout the fatigue life at different loads and environmental conditions, leading to a better estimation of the fatigue behavior of asphalt mixtures.

The mechanistic VECD fatigue prediction model has been applied as a modern method to estimate the mechanical behavior of asphalt mixtures. Different versions of the formulations have been developed by various researchers [5,7,37,38]. In this paper, the model developed by Lee et al. [7] and Kim and Little [5] is presented, and it will be applied to compare the fatigue lives of the FAMs evaluated in this study. The model estimates the number of cycles required to degrade the material, N_f, to a certain pseudo-stiffness level, C, or to reach a certain amount of damage, S_f, at an arbitrary frequency, f, and controlled pseudo-strain amplitude, ε^R. The parameters C₁₁ and C₁₂ are obtained from Equation (40) [5,7]. Kim and Little [5] performed torsional shear cyclic tests in sand asphalt mixtures, and compared the measured fatigue lives of the materials with the values predicted from the fatigue model (Equations (49)–(51)). They concluded that the model parameters might provide a reasonable representation of the fatigue response.

$$\mathrm{A}=\mathrm{f}{\left\{\frac{1}{2}{\mathrm{C}}_{11}{\mathrm{C}}_{12}\right\}}^{\mathsf{\alpha }}{\left\{1+\mathsf{\alpha }\left(1-{\mathrm{C}}_{12}\right)\right\}}^{-1}{\mathrm{S}}_{\mathrm{f}}^{\left[1+\mathsf{\alpha }(1-{\mathrm{C}}_{12})\right]}$$

(49)

$$\mathrm{B}=2\mathsf{\alpha }$$

(50)

$${\mathrm{N}}_{\mathrm{f}}={\mathrm{A}\left[{\mathsf{\epsilon }}_{\mathrm{R}}\right]}^{-\mathrm{B}}$$

(51)

2.7. Fatigue Failure Criterion

Studies have shown that the traditional fatigue failure relationship based on constant amplitude loading [39], used in a variety of layered elastic pavement design methods, may underestimate field fatigue life by as much as 100 times [40]. This difference between laboratory and field fatigue performance may be attributed to differences between laboratory and field loading conditions, such as rest periods and varying load magnitude in the field. By applying the VECD approach, it was observed that fatigue failure was influenced by both relaxation and damage mechanisms, and due to this, the traditional failure criterion (50% reduction in the initial stiffness) was modified to the 50% reduction in the initial secant pseudo-stiffness criterion [25]. However, a 50% reduction in the initial secant pseudo-stiffness may not represent the fatigue failure of all materials [5,41,42]. Other criteria based on changes in phase angle during fatigue testing and the reduced energy ratio concept have been defended in studies by Rowe [10], Reese [11], Hopman et al. [12], and Rowe and Bouldin [13].

Research efforts have been conducted to identify a proper failure criterion based on the VECD approach and energy-based concept [5,43,44,45,46,47,48,49]. Daniel et al. [43] made a comparison between the VECD and dissipated energy (DE) approaches using uniaxial direct-tension fatigue tests. The two approaches were also compared to the traditional phenomenological approach relating the initial strain to the number of cycles to 50% reduction in initial stiffness. The number of cycles to failure showed strong correlation between the VECD and the DE failure criteria. Findings of a study conducted by Kim and Little [5] showed that three damage parameters were effective to characterize fatigue damage during torsional loading: (i) a decay in pseudo stiffness; (ii) a loss of nonlinear dynamic modulus; and (iii) a change in dissipated strain energy.

In a study by Hou [9], the pseudo stiffness at failure was plotted against test reduced frequency for multiple mixtures to provide an empirical observation of all the tested mixtures, in order to determine the failure criterion. It was found that: (i) the pseudo stiffness at failure increased with the reduced frequency, and (ii) the rate of change in the pseudo stiffness at failure as a function of reduced frequency increased. The observations led to a reduced-frequency piecewise function that was applied for failure criterion development. Tarefder et al. [44] conducted a study in order to compare the traditional stiffness [50] and the energy-based fatigue failure criterion [13] to the fatigue failure criterion based on the VECD approach [40]. It was shown that the traditional approach (50% reduction in the initial stiffness) was conservative, and that the fatigue life predicted by the VECD approach was always shorter than the one predicted by the energy ratio approach.

Sabouri and coworkers proposed a new energy-based failure criterion for the continuum damage model that should consistently predict failure of a material that reaches macrofracture. The proposed failure criterion relates the average of the dissipated pseudo energy values in cycles, denoted as G^R, and the number of cycles to fatigue failure. The criterion was able to predict the fatigue life of asphalt concrete mixtures across different temperatures and strain amplitudes [45,46], and has been applied to predict fatigue behavior of reclaimed asphalt pavement (RAP) mixtures and non-RAP mixtures, mixtures prepared with modified and unmodified binders, warm-mix asphalt mixtures, and long-term aged mixtures [51,52,53,54,55,56,57,58,59,60]. Keshavarzi and Kim [49] extended the G^R criterion concepts to determine when fracture occurs in monotonic failure tests, such as thermal stress restrained specimen tests (TSRSTs).

Mensching et al. [61] suggested that more studies are needed to determine a G^R index that is more representative of field conditions, because there was not a clear trend between RAP increase and a change in the mixture fatigue performance when the damage characteristic curves and G^R failure criterion were used. More recently, Wang and Kim [48] developed and validated a new energy-based failure criterion based on the S-VECD model (D^R). The advantages of the D^R failure criterion in comparison to the previous G^R failure criterion are: (i) the D^R can be computed for each fatigue test and then used to check the sample-to-sample variability for each test, (ii) D^R is obtained in arithmetic scale rather than in log–log scale, and (iii) the number of tests needed to obtain D^R is fewer than for the G^R failure criterion [48]. In a study by Wang et al. [62], a three-dimensional finite element program (FlexPAVETM) was used to simulate the fatigue performance of field test sections. The fatigue damage of the sections was predicted using the FlexPAVETM software by considering both G^R and D^R criteria. The D^R failure criterion was found to yield more realistic fatigue-cracking performance predictions than the G^R failure criterion [62], and has been checked for RAP and aged mixtures [63,64].

Etheridge et al. [47] studied the correlations among D^R, an index parameter referred to as apparent damage capacity (S_app) [65], and mix design factors such as nominal maximum aggregate size (NMAS), asphalt binder type, and binder content. The S_app parameter was shown to be able to predict the fatigue-cracking propensity of asphalt mixtures [47,65,66], and it was found to have a strong relationship with polymer modification and the NMAS of asphalt mixtures [47]. Based on experimental data and the Georgia Department of Transportation’s practical guidelines for specific mixtures, a study by Etheridge et al. [47] developed S_app threshold values for different traffic levels. In a study by Zhang et al. [67], the D^R parameter showed good correlation with three new performance indices from the linear amplitude sweep (LAS) test: (i) strain tolerance (ε_T), (ii) strain energy tolerance (ε_E), and (iii) average reduction in integrity to failure (I^R).

2.8. Linear Viscoelasticity

The linear-viscoelastic characterization of asphalt mixtures is carried out by measuring relaxation or creep properties. The creep and relaxation tests measure the response of the materials to a constant load or displacement over time; i.e., they describe the material properties in the time domain within the linear-viscoelastic region. In the creep test, a constant load is applied to the specimen over time at a constant temperature, and the creep compliance, D(t), is defined as the ratio of accumulated strain to the constant stress magnitude at a specific time. In the relaxation test, the specimen is submitted to a constant strain for a given period of time at a constant temperature, and the relaxation modulus, E(t), is defined as the ratio of the stress evolution to the constant strain magnitude. However, in some cases, it is not possible to obtain an accurate response of the material by carrying out a short-time test with transient excitation; e.g., static loading. To overcome this limitation, tests with steady-state sinusoidal excitation are adopted; e.g., dynamic loading.

Viscoelastic materials under dynamic loading conditions provide frequency-domain dynamic properties, such as (i) phase angle, ϕ(ω), which represents the gap between the stress and strain due to the time-dependency of viscoelastic materials; (ii) storage shear modulus, G′(ω), which represents the elastic characteristics of the material; and (iii) loss shear modulus, G″(ω), which corresponds to the viscous behavior of the material. The combined form of storage shear modulus and loss shear modulus results in Equation (52) for the phase angle, where ω is the angular frequency; in Equation (53) for the dynamic shear modulus, |G*(ω)|, where τ_max is the maximum shear stress at each cycle and γ_max is the applied cyclic shear strain amplitude; and in Equation (54) for the complex shear modulus, G*(ω), where i is equal to √(−1).

$$\mathsf{\phi }\left(\mathsf{\omega }\right){=\tan }^{-1}\left[\frac{{\mathrm{G}}^{″}\left(\mathsf{\omega }\right)}{{\mathrm{G}}^{\prime }\left(\mathsf{\omega }\right)}\right]$$

(52)

$$\left|{\mathrm{G}}^{\ast }(\mathsf{\omega })\right|=\frac{{\mathsf{\tau }}_{\max }}{{\mathsf{\gamma }}_{\max }=\sqrt{{\left[{\mathrm{G}}^{\prime }\left(\mathsf{\omega }\right)\right]}^2+{\left[{\mathrm{G}}^{″}\left(\mathsf{\omega }\right)\right]}^2}}$$

(53)

$${\mathrm{G}}^{\ast }(\mathsf{\omega })={\mathrm{G}}^{\prime }(\mathsf{\omega })+\mathrm{i}{\mathrm{G}}^{″}\left(\mathsf{\omega }\right)$$

(54)

A dynamic frequency sweep test within the linear-viscoelastic range is carried out in the dynamic shear rheometer (DSR) to define the linear-viscoelastic relaxation behavior of the material. A curve-fitting function for the linear-viscoelastic modulus and frequency is required to determine the linear-viscoelastic relaxation modulus from a test in the frequency domain.

The mechanical behavior of viscoelastic materials can be represented by an association of springs and dashpots. The spring represents the elastic behavior of the material, which obeys Hooke’s law: σ = Gε, where σ is the stress, ε is the strain, and G is the elasticity modulus. The viscous behavior of the material is represented by the dashpot. Viscous elements are mathematically represented by Newton’s law of viscosity: σ = η$\dot{\mathsf{\epsilon }}$, where η is the viscosity coefficient and $\dot{\mathsf{\epsilon }}$ is the strain rate. Different arrangements of springs and dashpots can represent different material behaviors. The Maxwell model is a two-component model composed of a series connection of a spring and a dashpot, and corresponds to a viscoelastic fluid [68]. The Kelvin–Voigt Model is also a two-component model represented by a spring and a dashpot combined in a parallel fashion, and represents a viscoelastic solid [68,69]. Burger’s model is a four-component model and consists of a series arrangement of a Maxwell and a Kelvin–Voigt element [69]. While the Maxwell and Burger’s elements can exhibit typical rheological behavior of asphalt binders under constant stress–strain tests, the Kelvin element cannot describe relaxation and should not be used to represent asphalt binders in any hot-mix asphalt models [69].

In this study, the generalized Maxwell model, commonly referred to as the Prony series, is used as a curve-fitting function for the viscoelastic materials due to its capability of describing the different stages of the behavior of viscoelastic materials [70]. In the generalized Maxwell model, the Maxwell model elements are combined in a parallel fashion as shown in Figure 2. The relationship between the tensile stress σ(t) and time is shown in Equation (55) [68].

$$\mathsf{\sigma }(\mathrm{t})=\stackrel{\cdot }{\mathsf{\epsilon }}\mathsf{\eta }(1-{\mathrm{e}}^{-\frac{\mathrm{Et}}{\mathsf{\eta }}})$$

(55)

Prony series representation of storage and loss moduli as a function of frequency (frequency domain) was presented by Christensen [68], which is given by Equations (56) and (57), where G_e is the equilibrium modulus, G_i is the elastic modulus, ρ_i is the relaxation time, ω is the angular frequency, and n is the number of elements of the Prony series needed to fit the analytical representation to the experimental data.

$${\mathrm{G}}^{\prime }\left(\mathsf{\omega }\right){=\mathrm{G}}_{\mathrm{e}}+{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{G}}_{\mathrm{i}}{\mathsf{\omega }}^2{\mathsf{\rho }}_{\mathrm{i}}^2}{{\mathsf{\omega }}^2{\mathsf{\rho }}_{\mathrm{i}}^2+1}$$

(56)

$${\mathrm{G}}^{\prime \prime }\left(\mathsf{\omega }\right)={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{G}}_{\mathrm{i}}{\mathsf{\omega }\mathsf{\rho }}_{\mathrm{i}}}{{\mathsf{\omega }}^2{\mathsf{\rho }}_{\mathrm{i}}^2+1}$$

(57)

The spring constants (G_i) are defined by means of the collocation method, a matching process between the analytical representation and the experimental data for a certain number of points. Considering the Prony series parameters found by the collocation method, the static relaxation shear modulus as a function of time (time domain) can be predicted from the dynamic shear modulus as a function of frequency (frequency domain) by Equation (58). The relaxation property (m-value) is determined as the slope of the relaxation modulus curve, in logarithm scale, and is used in the VECD approach to determine the damage evolution rate of the material. This material property can be obtained by adjusting a power law function (Equation (59)) to the relaxation curve predicted by the Prony series, where G₀ and G₁ are material constants, t is the time, and m is the slope of the relaxation curve in the time domain.

$${\mathrm{G}\left(\mathrm{t}\right)=\mathrm{G}}_{\infty }+{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{G}}_{\mathrm{i}}{\mathrm{e}}^{-\frac{\mathrm{t}}{{\mathsf{\rho }}_{\mathrm{i}}}}$$

(58)

$${\mathrm{G}\left(\mathrm{t}\right)=\mathrm{G}}_0{+\mathrm{G}}_1{\mathrm{t}}^{-\mathrm{m}}$$

(59)

The experimental characterization of the dynamic modulus of asphalt concrete mixtures follows the test protocol AASHTO TP 62 [71]. In the absence of a specific procedure for FAMs, researchers have been adopted different criteria to obtain the dynamic modulus within the linear-viscoelastic region of each material. For example, Ng [72] and Klug et al. [73] defined the linear-viscoelasticity region as the range of stresses under which the material presents a deviation of less than 10% of its initial stiffness, and a stress of 15 kPa was applied to the specimens to obtain the linear-viscoelastic properties. In turn, Karki et al. [16] subjected each specimen to an oscillatory strain amplitude of 0.008% in shear for 30 s at 10 Hz and 25 °C. Ideally, the linear-viscoelastic dynamic shear modulus, |G*|_LVE, is equal to the dynamic shear modulus value that is measured by applying a small level of strain or stress within the linear-viscoelastic range of the material [16].

3. Studies on the Application of the VECD Theory

3.1. Full Asphalt Mixture Approach

Based on previous studies developed to characterize the damage behavior of asphalt mixtures [21,24,25,27,74], Daniel and Kim [75] proposed a testing procedure for fatigue characterization of asphalt concrete specimens under monotonic (constant-crosshead rate tests to failure) and cyclic (controlled-crosshead strain amplitude cyclic fatigue testing in tension) loading, which consists of two steps: (i) to perform a frequency sweep test on at least three replicate specimens at different frequencies and temperatures, for linear-viscoelastic material characterization (phase angle and relaxation modulus); and (ii) to perform a constant strain rate test to failure on all replicate specimens at a single rate at the desired temperature, for damage characterization. Pseudo strains are calculated according to Equation (15). Normalized pseudo stiffness, C, and damage parameter, S, for all times are calculated by Equation (34) and Equation (39), respectively. The values are cross-plotted to construct the characteristic curve that describes the reduction in material integrity as damage grows in the specimen (CxS) and to determine the functional coefficients C₁ and C₂ (Equation (40)). Daniel and Kim [75] showed that a single CxS curve can be obtained for each material, regardless of the applied loading conditions (cyclic vs. monotonic, amplitude/rate, frequency). However, Lundström and Isacsson [76] indicated that it was difficult to generally predict fatigue results based on characteristic curves obtained from monotonic tests. A later study conducted by Keshavarzi and Kim [77] applied the viscoelastic continuum damage (VECD) theory to simulate asphalt concrete behavior under monotonic loading. In that study, direct tension monotonic testing that incorporated a constant crosshead displacement rate and four temperatures was used to simulate thermal cracking of asphalt concrete prepared with four reclaimed asphalt pavement (RAP) proportions. The predictions of monotonic simulation matched the measured data of the monotonic tests very well up to the point of maximum stress. More recently, Cheng et al. [66] observed that the asphalt mixture CxS curves were independent of the strain level, but affected by the loading waveform.

Daniel and Kim [75] also verified that the characteristic curve at any temperature below 20 °C can be found by utilizing the time–temperature superposition principle (t-TS) and the concept of reduced time. Later, Chehab et al. [78] demonstrated that t-TS can be extended from material’s linear-viscoelastic range to high damage levels. Findings from compression tests [79,80] complemented the findings from Chehab et al. [79] for tension tests. Underwood et al. [81] also found that the t-TS principle with growing damage was applicable to mixtures with modified asphalt binders. The analysis of the fatigue behavior using the constant-crosshead rate test method and the t-TS principle [82] has been successfully employed to evaluate mixtures prepared with RAS, mixtures prepared with RAP, warm-mix asphalt (WMA) mixtures, mixtures prepared with modified binders, and other factors such as aggregate gradation, air voids, moisture, and aging [9,56,57,58,59,60,61,62,63,64,83,84,85,86,87,88,89,90,91,92].

In the studies carried out by Kim et al. [93] and Underwood et al. [33,35], a simplified VECD model was implemented in a finite element package (FEP++) to predict the fatigue performance of asphalt mixtures tested at test road project sites. In this VECD-FEP++ approach, the viscoelastic nature of asphalt concrete (AC) mixtures with growing damage is addressed using the VECD model, whereas the finite element program (FEP++) accounts for other important characteristics, such as temperature, layer thickness, stiffness gradient, etc. Comparisons between the field fatigue performance of the test road pavements to those predicted by the VECD-FEP++ simulations showed a generally positive relationship between model predictions and field observations. For a quick fatigue assessment, Underwood et al. [33,35,38] implemented the simplified VECD (S-VECD) model in the FEP++, instead of the original VECD model. The fatigue life of the pavements predicted using the S-VECD-FEP++ was found to agree well with the measured field response (R² = 0.8473% for the no-Terpolymer scenario, and R² = 0.9932 for the with-Terpolymer scenario). The S-VECD-FEP++ model was shown to be able to capture the effects of structure, climatic region, unbound layer modulus, and asphalt mixture properties, and to distinguish between top-down and bottom-up cracking.

The flexural bending test, also known as the beam fatigue test [50], is another testing method used to characterize the damage behavior of asphalt mixtures. This test measures the fatigue life of a compacted asphalt beam subjected to repeated flexural bending. According to De Mello et al. [94], what most likely happens in the field, and which is common to flexural fatigue tests, is a stress/strain field varying through the section from maximum compression to maximum tension in opposite sides instead of a homogeneous state of stress/strain throughout the sample section in cylindrical specimens subjected to uniaxial loading. The VECD approach was applied to flexural fatigue tests in studies with different goals, for example: (i) simplifying the calculation of damage parameters in the VECD model by considering the peak-to-peak values of stress and strain [95], (ii) evaluating which factors can influence the fatigue behavior more significantly [94], (iii) comparing fatigue cracking characteristics of a fine mix and a coarse mix [92], and (iv) studying the influence of RAP content on fine aggregate matrix (FAM) mixes [96]. Zhang et al. [96] also concluded that the linear amplitude sweep (LAS) test of FAM mixes under flexural bending mode can provide acceptable data with good repeatability as an alternative test method for tests with cylindrical samples in the dynamic shear rheometer (DSR).

3.2. Fine Aggregate Matrix Approach

Microstructural discontinuities, such as air voids and microcracks, coalesce and propagate due to repeated dynamic loading from the traffic of heavy vehicles and environmental loads, giving rise to the fatigue cracking process. The fatigue cracking reduces the structural performance of the pavement, with a negative impact on its service life. This process develops under two circumstances: (i) after adhesive failure, when the crack occurs at the aggregate–mortar interface; and/or (ii) after cohesive failure, when the crack develops within the mortar. Based on such an interpretation of the cracking phenomenon, Kim et al. [8] began to study the fatigue behavior of asphalt mixtures using the fine aggregate matrix (FAM) approach, and developed a protocol to evaluate the FAM’s properties.

The FAM is the matrix phase of the asphalt concrete composed of fine aggregates, filler, binder, and air voids, and represents the intermediate scale between the asphalt mastic and the full asphalt mixture. In studies with FAM mixtures, the primary assumption is that it reproduces the internal structure of the fine portion of the aggregate gradation of a full hot-mix asphalt (HMA) mixture. Another important assumption is that the physicochemical interactions between the aggregate and binder are replicated in the FAM specimen. Studies with a FAM gained notoriety for having a good agreement between the FAM and asphalt concrete (AC) properties, which was observed for the moisture characterization [97], fatigue cracking, and permanent deformation characterization [98,99,100,101,102].

Caro et al. [97] carried out surface energy measurements and dynamic mechanical analyzer tests on specimens of four FAM mixtures, and compared the results with the ones obtained by assessing the moisture susceptibility of the corresponding full asphalt mixtures by means of the saturation aging tensile stiffness (SATS) test [103]. A good agreement was observed between the results obtained for the fine-graded asphalt mixtures and the ones obtained for the dense asphalt mixtures. In both approaches, the granite mixture treated with 2% (by weight) hydrated lime was the most resistant to moisture damage, and the mixture containing only granite was the most susceptible to damage. However, the authors reported some differences between the FAM and HMA results for the sample containing crushed granite and limestone dust.

In order to investigate the inherent fatigue cracking resistance of modified asphalt binders (PPA, SBS, PPA + SBS, and PPA + Elvaloy), Motamed et al. [98] submitted FAM samples produced with modified asphalt binders and glass beads to torsional loading (controlled strain mode–275 kPa) at 10 Hz and 16 °C. The rationale for the use of glass beads in substitution for the mineral aggregate was to simulate the same stress state to which the binder was submitted in the asphalt concrete structure. The FAMs and the asphalt concrete mixtures fatigue lives (number of cycles to achieve 50% of the initial modulus) were compared, and it was observed that the FAM mixtures presented the same rank order for fatigue life of the asphalt concrete mixtures produced with the same modified asphalt binders [98].

Gudipudi and Underwood [18] observed a good agreement for the damage characteristic curves (C vs. S) between FAM and asphalt concrete for the tests carried out at 10 and 19 °C, but the C-values at failure for the FAMs were lower as compared to those of the asphalt concrete. It was not possible to compare results from tests with FAM and asphalt concrete carried out at 25 °C once the damage curve for both materials (FAM and asphalt concrete) presented a significant variation that could be related to viscoplasticity or another mechanism [18]. Coutinho [99] found the same rank order between the fatigue resistance of AC mixtures and the fatigue resistance of FAM mixtures when the FAM mixtures were subjected to stress-controlled time-sweep tests. Im et al. [100] observed a strong correlation between the linear and nonlinear viscoelastic and viscoplastic deformation characteristics of the asphalt concrete and its corresponding FAM. Underwood and Kim [104] evaluated the effect of different compositions for the four material scales (binder, mastic, FAM, and asphalt concrete) using linear-viscoelastic properties, such as dynamic shear modulus (|G*|) and phase angle (δ). The authors concluded that the dynamic modulus and the phase angle for the FAM materials were much more similar to the full mixture data than were the mastic materials. The study showed that the materials at different scales presented differing levels of sensitivity to changes in the blending parameters. The materials at the FAM scale presented a sensitivity that was more in line with that observed for asphalt concrete mixtures under all of the tested conditions.

Palvadi et al. [17] validated the VECD theory to characterize damage in FAM specimens based on the similarity of the characteristic curves (CxS) for a given FAM for both monotonic and cyclic loading modes and different amplitudes. Palvadi et al. [17] also proposed a test procedure to investigate the healing characteristics of FAM specimens. This test procedure consists of four rest periods (5, 10, 20, and 40 min) and three levels of stiffness (20%, 30%, and 40% reduction in C). In this method, four specimens of each FAM were tested in order to apply a specific rest period in a specific specimen. Palvadi et al. [17] concluded that the healing percentage of each FAM is a material characteristic, once that the values for this parameter were similar, regardless of both the sequence of application of the rest period and the damage level. In an attempt to improve the procedure proposed by Palvadi et al. [17], Karki et al. [16,105] developed an integrated testing procedure that was capable of quantifying damage and healing characteristics using a single specimen without separating the damage and healing tests. Karki et al. [16,105] were the first to apply the simplified viscoelastic continuum damage (S-VECD) theory to characterize FAMs. The researchers highlighted that the characteristic curve (CxS) is a unique material property due to the similarity of the curves for a given material, regardless of the loading conditions (different amplitudes and frequencies) and the introduction of rest periods during the test.

Gudipudi and Underwood [18] analyzed the fundamental similarities or differences between FAM and AC scales by means of the S-VECD theory. They observed similarities in material properties between the two material scales, and the CxS curves for a particular FAM and its corresponding AC mix were very similar. However, the C-value when failure occurred was generally lower in the FAM as compared to the mix. This result suggested that the FAM can reach a greater damage accumulation before failure occurs. The use of FAM testing for material characterization and ranking of AC mixtures has a great potential, if the material fabrication protocols are accurately followed. Freire et al. [106] applied the S-VECD theory to evaluate the effect of different maximum nominal aggregate sizes (MNS) of the mineral aggregate particles on the FAM fatigue resistance to identify which one best represented the damage characteristics of the asphalt mixture. The FAM mixtures were prepared with three different MNS (4.00, 2.00, and 1.18 mm), and their fatigue characteristics were evaluated and compared to the ones of a hot-mix asphalt (HMA) mixture prepared with an MNS of 12.5 mm. The main finding was that the Wöhler curves created for the FAMs produced with mineral aggregate particles of 2.00 mm and the corresponding asphalt mixture presented similar trends. The authors pointed out that a direct comparison of the absolute results, instead of the trends observed, could not be done directly because the parameters for the HMA mixture were obtained by axial loading, whereas the FAM parameters were obtained by shear loading.

Some researchers adapted the linear amplitude sweep test (LAS) method proposed by Johnson [107] to characterize the fatigue resistance of the FAMs using the VECD approach. The investigations evaluated the effect of (i) different particle size distributions [99], (ii) different nominal maximum aggregate size [108,109], (iii) the thermal and photochemical aging [110,111], and (iv) the effect of RAP content and rejuvenating agents (RAs) [112,113]. However, Freire et al. [109] did not recommend the use of the LAS test to analyze the fatigue resistance of the FAM mixes due to the difficulty to achieve the failure, once the torque capacity of the DSR was low and unable to take the sample to failure. The authors observed that for the highest strain amplitudes of the LAS test, the equipment needed to work near its capacity due to the high stiffness of the FAM specimens.

Regarding the FAM mixes containing RAP and recycled asphalt shingle (RAS), researchers [6,73,96,101,112,113,114] concluded that the use of these materials decreased the fatigue life of the mixture due to the hard binder present in the RAS and RAP. The use of RAs (petroleum tech, green tech, and agriculture tech) in the FAM mixes containing RAS and RAP was investigated by Zhang et al. [96], Nabizadeh [101], Nabizadeh et al. [114], and Zhu et al. [6] as an alternative to increase the fatigue life of the FAMs. Nabizadeh [101] and Zhang et al. [96] concluded that the RAs resulted in softer mixtures with improved fatigue life (especially for the FAMs with high RAP contents). Zhu et al. [6] observed the same behavior in the case of the FAM with RAS mixed with another petroleum-based RA. The combination of the WMA additive with the petroleum tech rejuvenator was evaluated by Nabizadeh [101], and this combination resulted in the softest FAM compared with other rejuvenators (green tech and agriculture tech).

With the aim of investigating the fatigue cracking of the asphalt binders at the FAM scale without the physicochemical interaction with the mineral aggregate, Motamed et al. [98] used rigid particles, such as glass beads, in substitution for the mineral aggregate to produce the FAM specimens. This new technique resulted in similar fatigue cracking characteristics between the FAMs and the asphalt mixtures produced with the same asphalt binder. The researchers concluded that the glass beads could be used in substitution for the mineral aggregate when the binder properties were the main issue of interest. More recently, Li et al. [115] used a combined fatigue–healing method based on the VECD model to evaluate fatigue and self-healing properties of three rock asphalt composites. Li et al. [115] indicated that the replacement of a portion of the virgin asphalt binder by the rock asphalts enhanced the fatigue cracking resistance of the FAM mixtures, and the influence on fatigue life was dependent on both the type and concentration of the rock asphalt.

Warm fine aggregate mixtures (W-FAM) fabricated using different WMA additives were compared with an HMA (control mixture) in a study by Sadeq et al. [116]. The control mixture presented higher dissipated pseudo-strain energy (DPSE) than the W-FAM, and the fatigue life in the VECD analysis approach was not statistically significant among the control mixture and the W-FAM mixtures, indicating that the WMA mixtures had fatigue resistance comparable to the hot-mix asphalt mixtures. In studies by Sadek et al. [117] and Sharma and Swamy [118], a probabilistic analysis approach was applied to the fatigue life prediction model deduced from the VECD theory. The inherent variability of asphaltic materials exhibited in the fatigue test results led the researchers to develop a new probabilistic approach. Probabilistic approaches present the ability to account for uncertainties associated with fatigue tests, models, and model parameters. In a study conducted by Sadek et al. [119], the efficacy of using the probabilistic approach in the analysis of the viscoelastic continuum damage (VECD) and fatigue life was examined for hot and warm fine aggregate mixtures (H-FAM and W-FAM). The probabilistic VECD approach had the advantage of providing more reliable predictions of fatigue life that accounted for uncertainty in determining the model parameters, instead of the deterministic approach. The probabilistic analysis results showed that the W-FAMs had shorter fatigue lives than the one obtained by the control H-FAM mixture; however, their fatigue lives presented more consistency and less uncertainty.

4. Analysis Protocol of Tests with FAM Using the S-VECD Approach

The outputs of tests carried out using fine aggregate matrix (FAM) specimens can be analyzed by means of the theory of continuum mechanics, by performing a linear-viscoelastic characterization followed by a damage evolution characterization.

Table 1. Summary of the VECD theory applied to asphalt mixtures—fingerprint and damage tests.

Fingerprint Test			FAM1 Sample 1
1. Data obtained from the fingerprint	Dynamic shear modulus within the linear-viscoelastic region (|G*LVE|), relaxation rate (m) (Figure 3a,b)	-	|G*|LVE = 9.93 × 108 m = 0.476
2. Prony series fitted to the storage modulus values to obtain ${\mathrm{G}}_{\mathrm{e}}$, ρi, and ${\mathrm{G}}_{\mathrm{i}}$.	${{\mathrm{G}}^{\prime }\left(\mathsf{\omega }\right)=\mathrm{G}}_{\mathrm{e}}+{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}\frac{{\mathrm{G}}_{\mathrm{i}}{\mathsf{\omega }}^2{\mathsf{\rho }}_{\mathrm{i}}{}^2}{{\mathsf{\omega }}^2{\mathsf{\rho }}_{\mathrm{i}}{}^2+1}$	Equation (56)	${\mathrm{G}}_{\mathrm{e}}$, ρi, and ${\mathrm{G}}_{\mathrm{i}}$ are calculated by using Solver
3. Laplace transform to convert data from frequency to time domain	${\mathrm{G}\left(\mathrm{t}\right)=\mathrm{G}}_{\infty }+{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{G}}_{\mathrm{i}}{\mathrm{e}}^{-\mathrm{t}/{\mathsf{\rho }}_{\mathrm{i}}}$	Equation (58)	$\mathrm{G}\left(\mathrm{t}\right)$ values obtained by using the parameter of the Prony series
4. Model adjusted to data to obtain the material parameter, m	G(t) = G0+G1t−m (Figure 3b)	Equation (59)	G0 = 1 G1 = 2.65 × 108 m = 0.476
5. Equation to obtain the material parameter, α	α = (1 + 1/m) for strain-controlled tests α = 1/m for stress-controlled tests	Equation (30) Equation (31)	α = 2.10
Damage Test			FAM1 sample 1 k = 5
1. Data obtained from the damage tests	Complex modulus (|G*|), phase angle (φ), and strain (ε), at each cycle k	-	|G*|k = 1 = 7.87 × 108 |G*|k = 5 = 6.91 × 108 φ k = 5 = 49.5 ε k = 5 = 0.051%
2. Peak pseudo strain during each cycle k, ${\mathsf{\epsilon }}_{\mathrm{m}}^{\mathrm{R}}\left(\mathrm{t}\right)$	${\mathsf{\epsilon }}_0^{\mathrm{R}}{\left(\mathrm{t}\right)=\mathsf{\epsilon }}_0{|\mathrm{G}}^{\ast }{|}_{\mathrm{LVE}}$	Equation (45)	${\mathsf{\epsilon }}_0^{\mathrm{R}}{\left(\mathrm{t}\right)=\mathsf{\epsilon }}_{0\mathrm{k}=5}=5.03\times 10$5
3. Pseudo stress, σR	σR = σ	Elastic-viscoelastic correspondence principle	σR = 350 kPa
4. Initial pseudo stiffness, I	$\mathrm{I}=\frac{{|\mathrm{G}}^{\ast }{|}_{\mathrm{k}=1}}{{|\mathrm{G}}^{\ast }{|}_{\mathrm{LVE}}}$	Dynamic shear modulus at the first cycle	$\mathrm{I}=\frac{7{.87\times 10}^8}{9{.93\times 10}^8}=0.792$
5. Normalized pseudo stiffness at each cycle k, C	${\mathrm{C}}_{\mathrm{k}}\left(\mathrm{S}\right)=\frac{{|\mathrm{G}}^{\ast }{|}_{\mathrm{k}=\mathrm{n}}}{{\mathrm{I}|\mathrm{G}}^{\ast }{|}_{\mathrm{LVE}}}$	Equation (47)	${\mathrm{C}}_5=\frac{6.91}{0.792\times 7.87}=0.88$
6. Damage parameter, S	$\mathrm{dS}\equiv {\left[\frac{\mathrm{I}}{2}{\left({\mathsf{\epsilon }}_{0,\mathrm{k},\mathrm{pp}}^{\mathrm{R}}\right)}^2\left({\mathrm{C}}_{\mathrm{k}-1}-{\mathrm{C}}_{\mathrm{k}}\right)\right]}^{\mathsf{\alpha }/(1+\mathsf{\alpha })}{\left({\mathrm{t}}_{\mathrm{k}}-{\mathrm{t}}_{\mathrm{k}-1}\right)}^{1/(1+\mathsf{\alpha })}$	Equation (48)	${\mathrm{S}}_5\equiv 4{.71\times 10}^7$
7. Characteristic curve, CxS	The characteristic curves, CxS, are obtained by cross-plotting the C values against the S values at each cycle k	Figure 4a	C5 = 0.88 vs. S5 = 4.71 × 107 at cycle k = 5

presents the resulting data obtained from fingerprint and damage tests, along with a summary of the equations employed to determine the linear-viscoelastic and damage properties, and

Table 2. Summary of the VECD theory applied to asphalt mixtures—fatigue life prediction.

Fatigue Life Prediction Model
1. Power law fitted to the CxS curve to obtain C10, C11, and C12	${\mathrm{C}}_1\left({\mathrm{S}}_1\right){=\mathrm{C}}_{10}-{\mathrm{C}}_{11}{{(\mathrm{S}}_1)}^{{\mathrm{C}}_{12}}$	Equation (40)	Figure 4a
2. Calculation of the parameter A	$\mathrm{A}=\mathrm{f}{\left\{\frac{1}{2}{\mathrm{C}}_1{\mathrm{C}}_2\right\}}^{\mathsf{\alpha }}{\left\{1+\mathsf{\alpha }\left(1-{\mathrm{C}}_2\right)\right\}}^{-1}{\mathrm{S}}_{\mathrm{f}}^{\left[1+\mathsf{\alpha }(1-{\mathrm{C}}_2)\right]}$	Equation (49)	A = 7.35 × 1029
3. Calculation of the parameter B	B = 2α	Equation (50)	B = 4.35
4. Prediction fatigue life curve	Nf = A[εR]−B	Equation (51)	Figure 4b
5. Nf (strain = 0.005%)	Nf(0.005%) = 7.35 × 1029 [4.46 × 104]−4.35	Equation (51)	4.59 × 109
6. Nf (strain = 0.2%)	Nf(0.2%) = 7.35 × 1029 [1.78 × 106]−4.35	Equation (51)	500.25

presents the equations employed to build the FAM prediction fatigue models.

5. Application of the S-VECD Theory: Laboratory Tests and Discussion of Results

5.1. Experimental Method—Materials and Preparation of the FAM Specimens

The effect of the addition of different proportions of reclaimed asphalt pavements (RAP) and virgin binder on the fatigue performance of fine aggregate matrices was evaluated by performing tests according to the protocol proposed by Kim and Little [5], and by analyzing the obtained data applying the simplified viscoelastic continuum damage (S-VECD) theory, following the procedures described in Table 1 and Table 2. Fingerprint and damage tests under the stress-controlled mode of loading were carried out using cylindrical fine aggregate matrix (FAM) specimens (40 mm in height and 12 mm in diameter). The FAMs were produced with two RAP proportions (20% and 40%). The binder content was adjusted by using either a performance grade (PG) 58-16 binder or a PG 64-22 binder. One source of virgin mineral aggregate was used—a basalt rock obtained from Bandeirantes Quarry, located in the city of Sao Carlos, state of Sao Paulo, Brazil. One RAP source obtained from roads around the city of Sao Carlos was used in the study. The mineral aggregate and asphalt binder were characterized following the standard procedures from the American Society for Testing and Materials (ASTM D7928, ASTM C128, ASTM D70, ASTM D6373, ASTM D6648, and ASTM D7175) and from the American Association of State Highway and Transportation Officials (AASHTO T84, AASHTO T85, and AASHTO T209). The mineral aggregate and asphalt binder characteristics are presented in

Table 3. Material characteristics.

Basalt Rock
Quarry identification	Bandeirantes
Specific gravity of coarse aggregates (g/cm³)	2.904		AASHTO T 85
Specific gravity of fine aggregates (g/cm³)	2.999		AASHTO T 84
Specific gravity of filler (g/cm³)	2.769		ASTM D7928
Absorption (%)	0.6		ASTM C128
RAP Material
Quarry location	São Carlos/SP
Maximum specific gravity (g/cm³)	2.596		AASHTO T209
Asphalt Binders
Performance grade (PG)	PG 58-16	PG 64-22	ASTM D6373
Specific gravity (g/dm³)	1.015	1.004	ASTM D70
Continuous grade—virgin (°C)	61.07	66.84	ASTM D7175
Continuous grade—short-term aged (°C)	65.52	66.94	ASTM D7175
Continuous grade—long-term aged (S [60]1)	−20.7	−26.9	ASTM D6648
Continuous grade—long-term aged (m [60]2)	−20.2	−27.4	ASTM D6648

¹ Creep stiffness at 60 s (MPa); ² slope at 60 s.

. Aggregate particles passing sieve #10 (2.00 mm) were used. The FAM aggregate gradation met the specification of the Brazilian Department of Terrestrial Infrastructure (DNIT 031/2004-ES). A dense-graded aggregate gradation represented by the middle of the C band was chosen, as it is a typical mineral composition used for road construction in Brazil. The RAP binder content was 8.3%. The binder content to be used to cover the virgin fine aggregate particles was estimated by means of the specific surface method [120], and the result was 9.0%. The difference between the actual binder content and the binder content calculated by means of the specific surface method was compensated by adding a complementary proportion of virgin asphalt binder (PG 58-16 or PG 64-22).

The results presented here refer to the following FAMs: (i) FAM1: 20% of RAP and binder PG 64-22; (ii) FAM2: 40% of RAP and binder PG 64-22; (iii) FAM3: 20% of RAP and binder PG 58-16; and (iv) FAM4: 40% of RAP and binder PG 58-16. The mixtures were compacted in a specially designed mold with internal dimensions of 40 mm in length and 12.8 mm in diameter in order to produce the cylindrical FAM specimens.

5.2. Experimental Method: Fingerprint Test

Fingerprint tests were carried out in an Anton Paar MCR 302 dynamic shear rheometer (DSR). The linear-viscoelastic properties of the FAMs were obtained by performing a fingerprint test at a stress of 15 kPa. This stress was selected within the linear-viscoelastic (LVE) range in order to avoid deformations larger than 100 µstrain at the end of the test. The frequencies were the following: 30, 26, 22, 18, 14, 10, 6, 4, 2, 1, 0.5, 0.2, 0.1, 0.05, and 0.01 (Hz). The software used to control the DSR was the Rheoplus/32 V3.62, and the data were exported to a Microsoft Excel spreadsheet. A Prony series (Equation (55)) was fitted to the storage modulus G′(ω) values obtained in the fingerprint test by using the Solver tool in Microsoft Excel in order to obtain the equation parameters. The static relaxation shear modulus as a function of time, G(t), was predicted by applying a Laplace transform (Equation (57)), which employed the Prony series parameters previously obtained. For the sake of exemplification, Figure 3a shows the results of G′(ω) versus frequency (ω) obtained for the specimen 1-FAM1, and Figure 3b shows a log–log chart of the G(t) values, where the slope of the log–log G(t) curve (Equation (32)) was obtained by adjusting a power law function (Equation (58)). The material damage evolution rate, α, was calculated according to Equation (31), as the damage tests were conducted under controlled stress [29]. The linear-viscoelastic complex modulus (|G*|_LVE) was the averaged |G*| values measured at 1 Hz in the fingerprint.

5.3. Experimental Method: Damage Test

Damage tests were carried out in the Anton Paar MCR 302 dynamic shear rheometer (DSR). Oscillatory tests under controlled stress at a frequency of 1 Hz and a temperature of 25 °C were performed in order to evaluate the damage evolution process of the FAM specimens. The damage tests were performed by using the same sample tested in the fingerprint test, except for the cases in which the specimen presented deformations larger than 100 µstrain at the last frequency of the test. In these cases, those specimens were discarded and new specimens were tested. The applied stress adopted to damage the specimens was 350 kPa. The data were exported from Rheoplus/32 V3.62 software to a Microsoft Excel spreadsheet. The pseudo stiffness (C) and the accumulated damage (S) values were calculated by means of the equations of the S-VECD theory (Equations (46)–(48)). The characteristic curves, CxS, were built according to the model proposed by Lee and Kim [25]. The material fatigue life, N_f, which is the number of axle load repetitions capable of leading the material to failure, was predicted by employing the mechanistic fatigue life prediction model (Equation (51)) developed by Lee et al. [7] and Kim and Little [5]. The fatigue models were adjusted for a 50% reduction in the material’s pseudo stiffness. Figure 4a,b shows an example of the CxS curve and the fatigue curve obtained for a single specimen. The detailed description of the procedure devised to deal with replicates of the same FAM and to build the CxS and fatigue curves is presented in the next section. Figure 4a presents the CxS curve built from laboratory data and the power law fitted to the data for the specimen 1-FAM1. Figure 4b presents the fatigue curve for the specimen 1-FAM1.

5.4. Analysis and Discussion of Results

In general, distinct specimens of the same FAM can present different values for the linear-viscoelastic and damage properties. The heterogeneity in the material properties for specimens of the same material may result in different characteristic curves. In order to overcome this issue and to accurately predict the damage behavior of the mixtures, the average properties of at least three specimens was considered in the model. By doing so, it was possible to obtain similar curves for distinct specimens of the same FAM. After averaging the properties and plotting the material CxS curves, an average CxS curve that best represented the material damage behaviour was built for each FAM.

Table 4. Linear-viscoelastic properties and fatigue model parameters of the FAMs.

Material	Sample	Air Voids	Linear-Viscoelastic Properties						CxS Parameters
			m	M (av)	cv (%)	|G*|LVE (kPa)	|G*|LVE (av) (kPa)	cv (%)	Shape Factor	DR
FAM1	s1	5.15	0.476	0.461	3.2	9.93 × 108	8.92 × 108	11.3	1.05	0.382
	s2	5.28	0.490		6.2	7.53 × 108		−15.6	0.92	0.373
	s3	4.96	0.422		−8.6	9.46 × 108		6.1	0.97	0.350
	s4	4.99	0.467		1.2	8.82 × 108		−1.1	0.94	0.345
	s5	4.99	0.452		−2.0	8.85 × 108		−0.8	1.08	0.394
FAM2	s1	5.19	0.342	0.375	−8.7	1.57 × 109	1.70 × 109	−7.6	0.8	-
	s2	5.05	0.386		2.9	1.61 × 109		−5.3	0.99	0.332
	s3	5.09	0.380		1.3	1.96 × 108		15.3	1.03	0.332
	s4	5.05	0.392		4.5	1.66 × 109		−2.4	0.98	0.325
FAM3	s1	5.16	0.463	0.460	0.7	5.04 × 108	5.95 × 108	−15.3	1.04	0.395
	s2	4.67	0.465		1.1	5.90 × 108		−0.8	0.96	0.415
	s3	5.47	0.448		−2.5	6.61 × 108		11.1	0.99	0.406
	s4	5.09	0.463		0.7	6.24 × 108		4.9	0.91	0.360
FAM4	s1	4.57	0.398	0.406	−1.8	1.54 × 109	1.57 × 109	−2.1	1.03	0.455
	s2	4.68	0.402		−1.0	1.63 × 109		3.6	1.00	0.448
	s3	4.89	0.417		2.7	1.55 × 109		−1.5	0.97	0.423

presents the air voids and the linear-viscoelastic properties (linear-viscoelastic complex modulus, |G*|_LVE, and the material relaxation rate, m) for each FAM specimen. Table 4 also shows the average linear-viscoelastic properties of each FAM and the coefficient of variation of each specimen. The specimen shape factor and the D^R criterion were tools that were able to check the specimen-to-specimen variability, and they are also presented in Table 4. As recommended by Lee et al. [7], the specimen-to-specimen variability can be quantified and minimized by means of a shape factor. By following such an idea, a shape factor was defined in this study as the ratio of the area above a particular specimen’s CxS curve to the area above the average CxS curve of the FAM. A shape factor close to 1 indicated a valid specimen; i.e., the results of a particular specimen were very close to the FAM average curve. In the study by Lee at al. [7], replicates with shape factors greater than 1.1 or lower than 0.9 were considered outliers. By following this premise, curves with shape factors of 1.0 ± 0.1 were accepted and used to build the characteristic curve of the FAMs. In the calculation of the shape factors, a decay until a pseudo stiffness equal to 0.3 was considered. The shape factors of all specimens are presented in Table 4, where it is possible to observe shape factors between 0.92 and 1.08, except for FAM2 s1, which presented a shape factor of 0.8. The other tool applied in this study to check the specimen-to-specimen variability was the D^R criterion, which is defined as the average reduction in pseudo stiffness up to failure [48]. Wang and Kim [48] found a variation of ±0.04 in the D^R value of specimens of the same mixture. The authors also carried out sensitivity studies of the pavement performance analysis using the S-VECD model with the D^R failure criterion, and they found that with this variation of ±0.04 in the D^R value, the predicted fatigue damage for the pavements did not differ significantly. The variation between the specimens of the FAMs evaluated in this study was ±0.05.

Regarding the effect of the RAP content, the linear-viscoelastic properties in Table 4 show that the FAMs containing 40% of RAP (2 and 4) presented higher |G*|_LVE values as compared to the mixtures containing 20% of RAP (1 and 3). For the same RAP proportion (20% or 40%), the highest complex modulus was observed for the FAMs containing the binder PG 64-22. Regarding the relaxation rate, m, the increase in the RAP content from 20 to 40% led to a decrease in this rate, with FAMs 1 and 3 (20% of RAP) presenting higher m values as compared to FAMs 2 and 4 (40% of RAP). Among the FAMs prepared with 20% of RAP, FAM1 (prepared with binder PG 64-22) and FAM3 (prepared with binder PG 58-16) presented equivalent m values (FAM3 presented a slight 0.2% reduction as compared to the FAM1 rate). Among the FAMs prepared with 40% of RAP, FAM4 (prepared with binder PG 58-16) presented an increase of 8.5% in the m value as compared to FAM2 (prepared with binder PG 64-22).

Figure 5 illustrates the CxS curves for each specimen built by using its own properties, while Figure 6 shows the CxS curves for each specimen built by using the average linear-viscoelastic properties. When comparing Figure 5 to Figure 6, it is possible to observe that the CxS curves tended to converge when the average properties were employed. However, this procedure did not yet result in a single CxS curve for each material. For such a reason, an average CxS curve for each FAM was built, which represented the material damage behavior, by using the Solver tool in Microsoft Excel. These curves are illustrated in Figure 6. It can be observed that FAM2 s1 is not presented in Figure 6, and this was due to its shape factor result of 0.8, which was lower than the limit of 0.9. For such a reason, the specimen FAM2 s1 was not included in the specimen set used to build the average characteristic curve of FAM2.

The FAM fatigue models (Figure 7b) were adjusted to a damage level (S) corresponding to a reduction of 50% in the material pseudo stiffness (C). These S values were obtained from the final (average) CxS curves (Figure 7a). By means of the fatigue models, it was possible to compare the effect of the different proportions of RAP (20% and 40%) and the effect of the addition of the asphalt binders (PG 64-22 and PG 58-16). The fatigue curves (Figure 7b) indicate that the materials behaved differently at different strain levels. For this reason, two strain levels were adopted for this analysis (0.005% and 0.20%). The parameters A and B of the fatigue models are presented in

Table 5. Parameters A and B of the fatigue models and FAM fatigue factors.

FAM	A	B	Nf (0.005%)	Rank Order	Nf (0.20%)	Rank Order	FFFAM	Rank Order
FAM1	7.35 × 1029	4.35	4.59 × 109	4	500.25	2	2.32	3
FAM2	2.68 × 1036	5.35	1.17 × 1010	3	31.80	4	2.30	4
FAM3	5.25 × 1029	4.35	1.77 × 1010	1	1885.85	1	2.49	1
FAM4	2.16 × 1034	4.93	1.52 × 1010	2	189.63	3	2.39	2

, along with the fatigue lives for each strain level and the mixture ordering (with the first positions occupied by the materials with the longest fatigue lives). Another procedure that can be applied to rank the FAM fatigue lives is the mixture fatigue factor (MFF) proposed by Nascimento et al. [121]. The MFF is determined by calculating the area below the curve “number of cycles vs strain” between two specific strains. Table 5 also presents the FAM fatigue factors (FF_FAM) and the rank order of the FAMs according to this factor, by assuming that the higher the factor, the higher the fatigue performance of the material.

FAMs 1, 2, 3, and 4 were compared in order to evaluate the effect of the addition of different proportions of RAP (20% and 40%) to the mixture. FAMs 2 and 4, containing 40% RAP, present higher stiffness as compared to FAMs 1 and 3, which contained 20% RAP, as well as lower relaxation rates, higher damage evolution rates, and higher parameters A and B of the fatigue model (Table 4 and Table 5). By comparing FAMs 1 and 2, it was possible to evaluate the effect of the increase in the RAP content from 20% to 40%, when a binder PG 64-22 was added to adjust the binder content of the RAP material. FAM2 (40% RAP) presented a higher stiffness and higher damage evolution rate than FAM2 (20% RAP). At low strains, FAM2 presented an increase of 2.6 times in its fatigue life when compared to FAM1; however, at high strains, the fatigue life of FAM2 was about 6.4% of the fatigue life of FAM1.

By comparing FAMs 3 and 4, it was possible to evaluate the effect of the increase in the RAP content from 20% to 40% when a binder PG 58-16 was added to adjust the binder content of the RAP material. FAM4 (40% RAP) presented a higher stiffness and higher damage evolution rate than FAM3 (20% RAP). At low strains, FAM4 presented a decrease of 14.5% in its fatigue life when compared to FAM3, and at high strains, the fatigue life of FAM4 was about 10% of the fatigue life of FAM3. In summary, the increase in the RAP content from 20% to 40% increased the fatigue life at low strains only if the binder content was corrected with a binder PG 64-22 and reduced it at high strains.

FAMs 1, 2, 3, and 4 were compared in order to evaluate the effect of adding the binder PG 64-22 or the binder PG 58-16 to adjust the binder content of the RAP material. When the FAMs containing 20% of RAP (FAMs 1 and 3) were compared, the results indicated a stiffness 33% lower for FAM3 (RAP binder content adjusted with the binder PG 58-16), and a similar damage evolution rate for both FAMs 1 and 3. At low strains, FAM3 presented a fatigue life 3.9 times longer than FAM1 (RAP binder content corrected with the binder 64-22), and at high strains, the fatigue life of FAM3 was 3.77 times longer than that of FAM1. When the FAMs containing 40% of RAP (FAMs 2 and 4) were compared, it is possible to observe that the stiffness and the damage evolution rate of FAM2 (RAP binder content adjusted with the binder PG 64-22) were about 8% higher than the values for FAM4 (RAP binder content adjusted with the binder PG 58-16). At low strains, FAM2 presented a fatigue life that was 17% of the one for FAM4, and at high strains, FAM4 has a fatigue life that was 1.3 times longer than the one expected for FAM2. In summary, the use of the binder PG 58-16 to adjust the FAM binder content increased the fatigue life at both low and high strains.

The rank order at low strains showed that FAM3, prepared with 20% RAP and the binder content of the RAP material adjusted with the binder PG 58-16, presented the longest fatigue life; and FAM1, prepared with 20% RAP and the binder content of the RAP material adjusted with the binder PG 64-22, presented the shortest fatigue life. FAM4, prepared with 40% RAP and the binder content of the RAP material adjusted with a binder PG 58-16, occupied the second position in the rank order, followed by FAM2, prepared with 40% RAP and the binder content of the RAP material adjusted with a binder PG 64-22, occupying the third position. The rank order at high strains showed FAM3 occupying the first position, followed by FAMs 1, 4, and 2. Concerning the fatigue factors (FFFAM), the rank order showed that FAMs 3 and 4, prepared with 20% and 40% of RAP, respectively, and the binder content of the RAP material adjusted with the binder PG 58-16, occupied the first and second positions, respectively. FAMs 1 and 2, prepared with 20% and 40% of RAP, respectively, and the binder content of the RAP material adjusted with the binder PG 64-22, occupied the third and fourth positions, respectively.

Concerning the D^R values presented in Table 4, it can be observed that FAM4, prepared with 40% RAP and the binder content of the RAP material adjusted with the binder PG 58-16, presented the higher D^R values, followed by FAM3, prepared with 20% RAP and the binder content of the RAP material adjusted with the binder PG 58-16. It can also be observed that the lower values of the D^R parameter were the ones of FAM2, which was prepared with 40% RAP and the binder content of the RAP material adjusted with the binder PG 64-22. The trend observed by Wang and Kim (2019) was that the D^R value decreased as the RAP content increased, which was the same trend observed in this study for the two mixtures prepared with the binder PG 64-22. However, for the two mixtures prepared with the binder PG 58-16, the D^R values increased with the increase in RAP. Wang and Kim (2019) emphasized that the D^R value alone cannot be used to compare the fatigue performance of different asphalt mixtures.

The overall findings indicated that: (i) the increase in RAP from 20% to 40% increased the fatigue life of the mixture prepared with the binder PG 64-22, decreased the fatigue life of the mixture prepared with the binder PG 58-16 at low strains, and decreased the fatigue life of the mixtures at high strains; (ii) the addition of the binder PG 58-16 instead of the binder PG 64-22 resulted in an increase in the mixture fatigue life at both low and high strain levels; and (iii) the rank order of the FAM fatigue factors indicated that the increase in RAP from 20% to 40% decreased the fatigue life of the mixtures. Figure 8 depicts the experimental method followed in this study.

6. Conclusions

This paper presented a comprehensive literature review of the fundamentals of the continuum damage approach employed in the characterization of the fatigue resistance of asphalt mixtures. It also presented results of a laboratory study carried out to evaluate to effects of the addition of different proportions of reclaimed asphalt pavements (RAP) and two asphalt binders of distinct performance grades on the fatigue response of fine aggregate matrix (FAM) mixtures. In this study, the simplified VECD (S-VECD) theory was applied to analyze the resulting data from the damage tests, along with a detailed description of the calculations as an illustrative example of how to handle the equations of the S-VECD theory for this sort of analysis.

• The findings of the experimental study with RAP and binders PG 64-22 and PG 58-16 indicated that the FAMs containing 40% of RAP (2 and 4) presented higher |G*|_LVE values and higher damage evolution ratios as compared to the FAMs containing 20% of RAP (1 and 3). Out of the FAMs prepared with 20% of RAP (FAM1 and FAM3), the highest |G*|_LVE was observed for the FAM containing binder PG 64-22 (FAM1), and the damage evolution ratios were the same for both FAMs, which was an expected result, once the presence of the softest binder (PG 58-16) was supposed to reduce the stiffness of the FAMs (3 and 4).
• Regarding the prediction of the fatigue lives of the materials evaluated in the experimental study, the addition of RAP increased the parameter A of the fatigue model (related to the initial stiffness of the material and how the stiffness changed with the evolution of the damage) and the parameter B (related to the damage evolution rate)—the resulting fatigue lives of the FAMs prepared with 20% RAP were longer than the ones obtained for the FAMs prepared with 40% of RAP. The fatigue performance was directly related to the specimen stiffness: the higher the stiffness, the higher its susceptibility to damage and the lower the relaxation rates (which resulted in higher damage accumulation rates). The best solution to adjust the binder content of FAMs produced with 20% and 40% of RAP was the use of the binder PG 58-16. The FAM tests combined with the S-VECD theory as a tool to analyze the results was a practical approach, and is widely used to evaluate all sorts of variables of an asphalt mixture. However, some variables, such as low temperatures and/or high percentages of RAP, turn the mixtures into overly stiff materials, and the tests can be unpractical due to the limits of the rheometer torque. Equipment with a higher torque capability could accelerate the test duration.
• The improvement of computational simulations of the test protocols is an important subject for future works, and could contribute to a better understanding of the mechanisms and variables involved in the fatigue process, and could also help overcome the rheometer limitations.
• Comparisons between fractures mechanics and continuum mechanics results could also be an interesting topic to improve the VECD model in order to account for the different types of damage: adhesive or cohesive.
• Regarding materials science and development of advanced/new materials, the FAM approach combined with the S-VECD approach offered several new possibilities in terms of material performance evaluation and material development. Some examples can be mentioned concerning the fatigue performance: (i) the evaluation of the impact of higher RAP contents added to new AC mixtures; (ii) the evaluation of the impact of recycling agents at different contents, including petroleum-based materials, vegetable-based oils, and recycled oils; (iii) the assessment of the aging impact on fatigue; (iv) the assessment of moisture damage on fatigue resistance; (v) the assessment of new asphalt modifiers, including hybrid modification using virgin and recycled materials; and (vi) the evaluation of the effect of distinct aggregate types and aggregate gradations, among others. Several doubts related to these subjects can be countered by carrying out tests at the FAM scale and using the S-VECD approach. However, one must keep in mind that such a development also depends on a larger number of experiments on the correlation between the fatigue performance at the two scales (FAM and full asphalt mixtures). Such experiments are essential for the development and popularization of these very promising techniques.