Comparison of Relaxation Modulus Converted from Frequency- and Time-Dependent Viscoelastic Functions through Numerical Methods

Due to the difficulty of obtaining relaxation modulus directly from experiments, many interconversion methods from other viscoelastic functions to relaxation modulus were developed in previous years. The objectives of this paper were to analyze the difference of relaxation modulus converted from dynamic modulus and creep compliance and explore its potential causes. The selected methods were the numerical interconversions based on Prony series representation. For the dynamic to relaxation conversion, the time spectrum was determined by the collocation method. Meanwhile, for the creep to relaxation conversion, both the collocation method and least squares method were adopted to perform the Laplace transform. The results show that these two methods do not present a significant difference in estimating relaxation modulus. Their difference mostly exists in the transient reduced time region. Calculating the average of two methods is suggested to avoid great deviation of single experiment. To predict viscoelastic responses from creep compliance, the collocation method yields comparable results to the least squares method. Thus, simply-calculated collocation method may be preferable in practice. Further, the master curve pattern is sensitive to the Prony series coefficients. The difference in transient reduced time region may be attributed to the indeterminate Prony series coefficients.


Introduction
Asphalt mixtures exhibit viscoelastic behaviors over a wide range of time and temperatures. The responses characterizing the time-and frequency-dependent viscoelastic properties of asphalt mixtures include dynamic modulus, creep compliance and relaxation modulus. These parameters can not only be used to describe the linear viscoelastic properties, but also can indicate the nonlinear viscoelastic properties and failure properties of materials. In fact, these three properties can be interconverted between time-and frequency-domains because linear viscoelastic material properties shared an equivalent mathematical theory, including differential and integral equations [1][2][3].
Dynamic modulus is mathematically defined as the maximum dynamic stress divided by the peak recoverable axial strain, and is a key input for the prediction of field rut depth or fatigue crack amount in the current AASHTOWare pavement ME design program. Creep compliance describes the phenomenon that strain decreases with time under constant stress. In contrast, relaxation describes the phenomenon that stress decreases with time under constant strain. Relaxation modulus is required to calculate thermal stress within asphalt pavement structure in finite element analysis and in the in parallel in Figure 1a can be used to depict the total stress of asphalt mixtures. Differentiating its strain based on elastic theory [15], the Prony series representation of complex modulus (E*) can be determined as follows: where E e is the equilibrium modulus, E e = lim t→∞ E(t), MPa; i is the complex number, i = √ −1; ω is the angular frequency, rad/s; E i are the relaxation strengths, MPa; ρ i are the relaxation times, ρ i = η i /E i , s; η i is the viscosity of dashpot in a Maxwell element. Differentiating its strain based on elastic theory [15], the Prony series representation of complex modulus (E*) can be determined as follows:   (2) where g D is the glassy compliance, 1/MPa; j D are the retardation strengths, 1/MPa;  j are the retardation times, s. In contrast, relaxation modulus is the ratio of stress and strain. It is also derived from the generalized Maxwell model and its Prony series representation is given by:

Interconversion Relationship between Time-Domain and Frequency-Domain Functions
Under uniaxial, nonaging and isothermal conditions, the linear viscoelastic behavior of asphalt mixtures can be depicted by Boltzmann superposition integral [1]: (5) where σ, ε denote stress and strain, respectively; τ Means the integral variable of time. Substituting Equation (5) to Equation (4) can derive the constitutive model of relaxation modulus and creep compliance in time region [24]: Creep compliance is the ratio of strain and stress under certain temperature and stress. It can be depicted by the generalized Kelvin model consisting of a spring and a dashpot and n Kelvin elements connected in series in Figure 1b. The Prony series representation of creep compliance (D(t)) can be expressed by: where D g is the glassy compliance, 1/MPa; D j are the retardation strengths, 1/MPa; τ j are the retardation times, s. In contrast, relaxation modulus is the ratio of stress and strain. It is also derived from the generalized Maxwell model and its Prony series representation is given by:

Interconversion Relationship between Time-Domain and Frequency-Domain Functions
Under uniaxial, nonaging and isothermal conditions, the linear viscoelastic behavior of asphalt mixtures can be depicted by Boltzmann superposition integral [1]: where σ, ε denote stress and strain, respectively; τ Means the integral variable of time. Substituting Equation (5) to Equation (4) can derive the constitutive model of relaxation modulus and creep compliance in time region [24]: Then, the relationship of relaxation modulus and creep compliance in Laplace transform domain [15] can be obtained: dt; s is the Laplace operator; E(s) is the Laplace transform of relaxation modulus, E(s) = sE(s); D(s) is the Laplace transform of creep compliance, D(s) = sD(s).
Since the complex function of viscoelastic materials is excited by sinusoidal loading under ω angular frequency, the relationship between it and Laplace transformed function can be found as [25]: where E* is the complex modulus; D* is the complex compliance. Thus, relationship between E * and D * can be given by: As mentioned above, the time-domain and frequency-domain interrelationship can be derived based on linear viscoelastic theory. Thus, in the following section, relaxation modulus converted by dynamic modulus and creep compliance will be conducted according to this principle.

Interconversion from Dynamic Modulus to Relaxation Modulus
The dynamic modulus and relaxation modulus sheared same Prony series coefficients, thus, the Prony series coefficients of relaxation modulus can be directly acquired by fitting master curve of dynamic modulus. Direct fitting of dynamic modulus is extremely difficult since the equation is the sum of storage modulus and loss modulus which cannot be directly expressed. Fitting dynamic modulus includes fitting the two components of dynamic modulus: Storage modulus and loss modulus. The Prony series by fitting either storage modulus or loss modulus can be used as the Prony series for dynamic modulus since dynamic modulus, storage modulus and loss modulus shared the same Prony series as shown in Equations (1) and (11). However, as shown in Equation (11), the fitting process of storage modulus can give the result of Ee which is a required parameter to calculate relaxation modulus in Equation (3), whereas the fitting process of loss modulus cannot. Thus, the fitting of storage modulus was applied to obtain the Prony series of dynamic modulus. Here, the storage modulus in viscoelastic materials measures the stored energy, representing the elastic portion, and the loss modulus measures the energy dissipated as heat, representing the viscous portion. Equation (10) shows the relationship between complex modulus, storage modulus, and loss modulus.
In Equation (10), E is the storage modulus; E is the loss modulus.
According to the representation of complex modulus in Equation (1), the storage modulus and loss modulus can be found as: Therefore, the procedure of interconversion from dynamic modulus to relaxation modulus is available: (1) Pre-smooth and fit the experimental data of dynamic modulus and phase angle; (2) Calculate the storage modulus E expriment based on experimental data; (3) Make the storage modulus E predict determined by Equation (11) equal to the E expriment so as to obtain Prony series coefficients E e and E i ; (4) Substitute E e and E i into the relaxation modulus representation to plot its master curve over a wide time region.

Interconversion from Creep Compliance to Relaxation Modulus
As noted above, Equation (6) reveals the relationship between creep compliance and relaxation modulus in Laplace transform domain. Accordingly, relaxation modulus can be converted by creep compliance based on this formula in the following.
Substituting the Prony series representations of relaxation modulus and creep compliance into Equation (6), the integral can be expressed by: where δ( ) is the Dirac delta function.
According to the literature [15], further simplifying the integral can derive the following equations: where, Adjusting them into linear system of equations, Equation (15) can be given as: In Equation (15), matrix [E] denotes the collection of E i for relaxation modulus, which can be derived from matrix operation between matrix [A] and matrix [B]. Time spectrum t k (k = 1, 2, ..., p) denotes the upper limit of integral in Equation (12) and can be identified through the collocation method or the least squares method. In the collocation method, the number of t k must be equal to E i , whereas the number of t k could be more than E i in the least squares method. In that case, the linear system may be overdetermined equations, and the exact solution cannot be reached. Thus, the optimal solution for it may be effective by calculating in MATLAB to minimize AE − B 2 .
Consequently, the procedure of the converted relaxation modulus from creep compliance can be given as: (1) Pre-smooth the experimental data of creep compliance and fit its Prony series coefficients D g and D j ; (2) Calculate the Prony series of relaxation modulus based on [E] = [A] −1 [B] by converting D g and D j to E e and E i ; (3) Substitute E e and E i into the relaxation modulus representation to plot the master curve of it over a wide time region.

Experiments Design
In the experiment, selected specimen is hot mix asphalt AC-13. The performance grade of asphalt binder is PG70-22. The aggregate gradation of tested HMA can be seen in Table 1. The specimen was compacted by Superpave Gyratory Compactor (SGC). The cylindrical AC sample was 150 mm in diameter and 170 mm in height. Then, these samples were core-drilled and sawed to finished specimens with 100 mm in diameter and 150 mm in height. Three duplicated specimens were prepared for repeated experiments and the air voids were controlled between 4% ± 0.5%. Having prepared the specimens, they were heated up to experimental temperature and preserved for one hour, then the dynamic modulus and creep compliance were measured.
The dynamic modulus and creep compliance were obtained by uniaxial compression tests which were conducted on Material Testing System (MTS). The vertical deformation was measured by four linear variable differential transducers (LVDT) at 90 • radial intervals with a 90 mm gauge length. In the dynamic modulus test, the maximum axial and horizontal strain was controlled in 100 µε. After dynamic modulus test, the specimens were required to stand for one hour waiting for deformation recovery. Then, the creep compliance test can be conducted. The specimens were tested under −20 • C, −10 • C, 0 • C, 10 • C, 20 • C and 30 • C. And 100 s at each temperature. The horizontal strain was controlled in 100 µε. The results were also the average of three repeated experiments.

Experimental Data Processing
As designed in the experiments, the micro-strain of specimen under different frequency is around 100 µε. Moreover, the measured dynamic modulus and phase angle is constant under changeable loading. Thus, the experiments can be regarded as non-destructive ones that can maintain the specimen within linear viscoelastic range [26].
According to the time-temperature superposition principle, dynamic modulus and phase angle under different temperatures and frequencies can be translated to an integrated curve, called master curve [27]. Therefore, to obtain a dynamic modulus master curve, the curves under different temperatures are required to be horizontally shifted to one curve at a given temperature [28]. In this case, the long-term mechanical behavior of the asphalt material can be predicted without having to perform experiments for a long duration of testing time. In addition, the mechanical behavior of the material is at very short loading times (or very high frequencies) which are difficult to obtain from the experiments can be determined using the master curve [29].
Sigmoid functions as Equation (16) can be employed to pre-smooth the discrete experimental data of dynamic modulus. Similarly, Equation (17) can be used to pre-smooth the discrete experimental data of phase angle [30]. The function of calculating the reduced frequency can be obtained through Equation (18). The reduced frequency is defined as the computed frequency at the reference temperature, equivalent to the actual loading frequency at the test temperature. All of these fitted parameters can be determined by fitting the functions to experimental data and the summary of the fitted results are shown in Table 2.
where, f r is the reduced frequency; f is the loading frequency; δ is the minimum modulus; α is the parameter of vertical span of dynamic modulus; β, γ are the shape parameters of the material; ξ 1 , ξ 2 , ξ 3 are all fitting parameters of experimental data; α T is the time-temperature shifting factor. The determination of α T is based on Williams-Landel-Ferry (WLF) model [28] in Equation (19).
where C 1 and C 2 are the WLF coefficients; T is the testing temperature; T 0 is the reference temperature, which indicates that temperature at which the master curve is constructed. Thus, after the processing of experimental data, dynamic modulus and phase angle master curves at 21.1 • C Reference Temperature of asphalt mixtures can be obtained as presented in Figure 2.
Similar to the process of fitting master curves of dynamic modulus and phase angle, the creep compliance master curve can be found as in Figure 3. Considering the requirement of comparison between dynamic modulus and creep compliance interconversions, the reference temperature is required to be the same. Similar to the process of fitting master curves of dynamic modulus and phase angle, the creep compliance master curve can be found as in Figure 3. Considering the requirement of comparison between dynamic modulus and creep compliance interconversions, the reference temperature is required to be the same.

Dynamic Modulus to Relaxation Modulus
Regarding to the procedures of interconversion from dynamic modulus to relaxation modulus, storage modulus is needed to be calculated based on Equation (12). In the equation, relaxation times ρ i is expected to be identified first. According to reference [3], the distribution of ρ i could be estimated through the collocation method. Moreover, in general terms, the integrated curves could be smoother with smaller relaxation time interval. The logarithmic time interval lower than 0.5 is recommended [2]. However, the curve oscillation has already reached an acceptable range when the logarithmic time interval is less than 1 [15]. Therefore, the logarithmic time interval 1 will be employed in the following computation. Another important issue is the number of Maxwell elements, more elements yields to higher accuracy in general, but leading to more complexity at the same time [31]. Usually, 10 elements are enough for the transformation [32]. Therefore, to span enough time range of the data, 11 elements will be taken, meaning that ρ i = 2 × 10 i−7 s, i = 1, 2, ..., 11. Consequently, the Prony series coefficients converted by dynamic modulus are listed in Table 2.

Creep Compliance to Relaxation Modulus
Before converting creep compliance to relaxation modulus, the retardation times of creep compliance need to be determined by the collocation method. Compared with the relaxation times ρ i of dynamic modulus, the same number of elements 11 is selected to describe retardation times τ j of creep compliance.
Time constants t k are also required which can be determined by the collocation method or the least squares method [15]. The number of time constants t k was also selected as 11 to keep consistency with the number of retardation times. However, in the least square method, t k = 2 × 10 (k−1)/2−6 s, k = 1, 2, ..., 21. Therefore, in the collocation method, matrix A is 11 × 11 and matrix B is 11 × 1, whereas in the least square method, matrix A is 20 × 11 and matrix B is 20 × 1, both matrix calculations were completed in software MATLAB. As can be seen, the matrix utilized in the least square method is more complicated and required more computational time in contrast to the collocation method.
Regarding the interconversion from creep compliance to relaxation modulus, τ j was firstly determined and it was assumed that τ j and ρ i are equivalent as suggested [33]. In this paper, τ j was selected between 2 × 10 −6 to 2 × 10 −4 to cover the entire range of relaxation modulus. Theoretically, the relaxation times ρ i are the negative reciprocal of s when E(s) approaches to infinite [15]. The ρ i can also be determined by root-finding method based on E(s) versus −1/s figure, whereas the calculated ρ i was almost identical by using the two methods: assuming τ j is equivalent to ρ i , and root-finding method [4,34]. The Prony series coefficients of creep compliance and relaxation modulus are finally calculated in software MATLAB and results are listed in Table 3.

Discussion of Results
Based on the above procedures, the results of the two numerical interconversion methods can be analyzed. Before comparing the results of relaxation modulus converted from dynamic modulus and creep compliance, it needs to determine which method is better to represent the creep to relaxation conversion results. The collocation method was adopted in the transferring dynamic modulus to relaxation modulus, whereas both collocation method and least squares method were adopted in the transferring creep compliance to relaxation conversion. Figure 4 compares the relaxation modulus obtained by the collocation method and squares method. The relaxation modulus from collocation method was compared with relaxation modulus from least squares method based on statistical t-test at a significance level of 0.05. Results indicate that there was no statistical difference between relaxation modulus based on the two methods. In other words, the relaxation modulus converted from creep compliance through the collocation method and least squares method is almost the same.  Assuming the collocation method is more reliable, the relaxation modulus converted through it can be selected to represent the creep to relaxation conversion results. Combined with the dynamic to relaxation conversion results, these two curves can be compared in Figure 5 in semi-log scale. Figure 5a is plotted by using semi-log scale in which a clear pattern of relaxation can be observed: the relaxation modulus dramatically reduced within a short period and then tends to plateau. Assuming the collocation method is more reliable, the relaxation modulus converted through it can be selected to represent the creep to relaxation conversion results. Combined with the dynamic to relaxation conversion results, these two curves can be compared in Figure 5 in semi-log scale.

Validation of the Methodology
To validate the relaxation modulus conversion methods mentioned above, data points from another research project are adopted in this paper for validation purposes. The test specimen was gyratory compacted using loose mixture from field asphalt plant. Three types of additives were adopted to modify neat binders, namely Evotherm, rubber, and polymer. Both dynamic modulus and creep compliance of specimen was tested using the same test apparatus, test temperatures and frequencies as specified in Section 3.1. Relaxation modulus was converted by using both dynamic modulus and creep compliance. Table 4 presents the Prony series coefficients of relaxation modulus of the three mixture types. Figures 6a-c show the master curve and percentage difference caused by dynamic modulus and creep compliance, respectively. As can be seen, the two master curves of relaxation modulus obtained by dynamic modulus and creep compliance are close to each by visual observation. The difference ranges are 1% to 28% for the specimens using Evotherm modified binders, 2% to  5a is plotted by using semi-log scale in which a clear pattern of relaxation can be observed: the relaxation modulus dramatically reduced within a short period and then tends to plateau. Figure 5b is plotted by using log-log scale in which the relationship between relaxation modulus and reduced time at a wide range is available. The relaxation modulus difference based on dynamic modulus and creep compliance was also calculated. As seen, such a difference was relatively small at the short-or long-term region, whereas the difference greatly increased in the middle part of reduced time. The difference curve also presents to be fluctuated. Meanwhile, it is noted that the relaxation modulus appears a three-phase feature in the log-log scale and presents approximately linearity in the transition region. Again, statistical t-test at a significance level of 0.05 was conducted and results indicate that there was no statistical difference between relaxation modulus based on the two methods.
Since the methods adopted are numerical methods, the accuracy is directly related to the number of mechanical model elements. In the context of 11 elements, there are still obvious difference despite the results are similar overall. It is also proved that the difference mostly lies in the lower reduced time range. Due to the same specimens were tested for these two methods, the specimen to specimen variation may not be the determinant factor. A potential reason may be the error of fitting Prony series coefficients in lower reduced time region. The master curve pattern is sensitive to the Prony series coefficients. However, the experimental data is discrete and hard to get continuous and accurate processed data. Pre-smoothing the raw data can improve the performance of obtaining more consistent Prony series coefficients but also may generate some deviation in the Prony series coefficients fitting. Therefore, the difference curve may present to be fluctuated and does not have an obvious rule.

Validation of the Methodology
To validate the relaxation modulus conversion methods mentioned above, data points from another research project are adopted in this paper for validation purposes. The test specimen was gyratory compacted using loose mixture from field asphalt plant. Three types of additives were adopted to modify neat binders, namely Evotherm, rubber, and polymer. Both dynamic modulus and creep compliance of specimen was tested using the same test apparatus, test temperatures and frequencies as specified in Section 3.1. Relaxation modulus was converted by using both dynamic modulus and creep compliance. Table 4 presents the Prony series coefficients of relaxation modulus of the three mixture types. Figure 6a-c show the master curve and percentage difference caused by dynamic modulus and creep compliance, respectively. As can be seen, the two master curves of relaxation modulus obtained by dynamic modulus and creep compliance are close to each by visual observation. The difference ranges are 1% to 28% for the specimens using Evotherm modified binders, 2% to 35% for the specimens using rubber modified binders, and 7% to 32% for the specimens using polymer modified binders, respectively. The validation results confirm that the methods used in this paper are robust.

Conclusions
This study compared the relaxation modulus converted form dynamic modulus and creep compliance based on the same HMA specimens. To save testing costs and obtain more accurate viscoelastic parameters of asphalt mixture at a wider range, some conclusions can be illustrated:

Conclusions
This study compared the relaxation modulus converted form dynamic modulus and creep compliance based on the same HMA specimens. To save testing costs and obtain more accurate viscoelastic parameters of asphalt mixture at a wider range, some conclusions can be illustrated:

•
For the transition of creep compliance to relaxation modulus, using the collocation method and least squares method is quite related. However, the collocation method is not as complex as the least squares method. Thus, the collocation method may be preferable in practice.

•
The master curves of relaxation modulus converted from dynamic modulus and creep compliance do not present statistical difference at both short-and long-term time domain, and the major difference is observed at the transient time region. • Potential reasons for the difference in the transient reduced time region may be attributed to the error in Prony series coefficients. The master curve pattern is sensitive to the Prony series coefficients. However, in the pre-smoothing and fitting process, the approximation is common, leading to loss of information. The most vulnerable time range of relaxation modulus may be the transient reduced time range.

Future Studies
The limitation of this paper is that no direct experimental data of relaxation modulus can be acquired to validate the efficiency of these methods. Thus, it cannot be identified which method is more accurate based on the relaxation modulus tests in the premise of same model elements. In future research, the converted relaxation modulus may be integrated into asphalt pavement cracking prediction model to explore its reliability to be the inputs.

Conflicts of Interest:
The authors declare no conflicts of interest.