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

Abstract

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.

1. 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 pavement design procedure [4]. However, the direct test of relaxation modulus is typically time- and cost-consuming, as the constant-strain is difficult to control in experiments, whereas it is relatively more convenient to measure the dynamic modulus and creep compliance in laboratory with high testing accuracy. Thus, relaxation modulus in most cases is obtained from dynamic modulus or creep compliance by conversion methods.

The master curve of dynamic modulus or creep compliance needs to be determined before converting to relaxation modulus. Different mathematical functions can be used to obtain master curve, such as pure power law, generalized power law, modified power law, and Prony series [5]. The pure power law, generalized power law, and modified power law have their own benefits, such as being easy to implement, or can be better used to represent the transition zone of a master curve. However, they were not selected since they are not suitable for characterizing short-term or long-term properties [6,7]. Prony series was selected in this paper since it has been widely proved to best represent the fitting equations of master curve [8]. Prony series is a mathematical method to fit a set of equidistant data based on a series of damped complex exponentials. It is widely used in linear viscoelastic material functions, since it can significantly improve the computational efficiency attributed to its physical basis in the theory of mechanical models dealing with linear springs and dashpots [9].

The interconversion between dynamic modulus and relaxation modulus was firstly proposed by Schapery [3]. Later, Cost and Becker [10] proposed the least squares method to directly fit the Prony series to experimental data. However, during these conversions, negative coefficients cannot be avoided, which may generate oscillations in the master curve. To solve the issue, Park and Kim [2] introduced constraints to ensure that all the Prony series coefficients are positive so that a smooth master curve could be plotted. Before fitting the Prony series, modification technologies such as pre-smoothing raw data was proved to portray a smoother curve over a wide range of time [4,11]. Recently, researchers have been developing methods to generate relaxation modulus with no need to use the collocation method or the least squares method, such as semi-analytical algorithm and Levenberg-Marquardt method [12,13,14].

The relaxation modulus can be obtained from the creep compliance through numerical or analytical approaches, which have been developed by several researchers [15,16,17,18]. In regards to exact analytical interconversion, the integral solution usually can be achieved in an infinite range. However, the experimental data can only be collected in a limited range. Thus, it is highly possible to get complicated expressions in the process of solving integrals [5]. Therefore, researchers proposed different approximate analytical methods to simplify the integral computation procedure, such as quasi-elastic relationship, power-law, Prony series, and interrelationship [15,16,17,19,20]. It is generally accepted that numerical method based on Prony series is easier to implement in contrast to the exact analytical method, and it is also more accurate than the approximate analytical method [21,22].

With various fitting approaches available, it is of great interest to compare the accuracy and efficiency of those approaches for the same group of testing results to help identify more appropriate method. In previous studies, despite the overall performance of converted relaxation modulus being acceptable, uncertain difference still often occurred. The performance of relaxation modulus converted from both dynamic modulus and creep compliance based on the same specimens was rarely compared. Thus, the objective of this paper was to adopt the numerical method based on Prony series to obtain relaxation modulus from dynamic modulus and creep compliance. The Prony series representations of these viscoelastic functions were introduced first. Then, rationale of the interconversion methods was elaborated. Finally, numerical examples were presented to compare the difference of two methods.

2. Methodology

2.1. Viscoelastic Response Representations

The viscoelastic responses of asphalt mixtures can be expressed by combination of spring and dashpot [23]. The generalized Maxwell model consisting of a spring and m Maxwell elements installed 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:

$$E^{*}=E_e+{\displaystyle \sum _{i=1}^m\frac{i\omega {\rho }_i{E}_i}{i\omega {\rho }_i+1}}$$

(1)

where $E_e$ is the equilibrium modulus, $E_e=\underset{t\to \infty }{\lim }E\left(t\right)$, MPa; i is the complex number, $i=\sqrt{-1}$; $\omega$ is the angular frequency, rad/s; $E_i$ are the relaxation strengths, MPa; ${\rho }_i$ are the relaxation times, ${\rho }_i={\eta }_i/E_i$, s; ${\eta }_i$ is the viscosity of dashpot in a Maxwell element.

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:

$$D\left(t\right)=D_g+{\displaystyle \sum _{j=1}^n{D}_j\left[1-e^{-\left(t/{\tau }_j\right)}\right]}$$

(2)

where $D_g$ is the glassy compliance, 1/MPa; $D_j$ are the retardation strengths, 1/MPa; ${\tau }_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:

$$E\left(t\right)=E_e+{\displaystyle \sum _{i=1}^m{E}_i}{e}^{-\left(t/{\rho }_i\right)}$$

(3)

2.2. 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]:

$$\sigma \left(t\right)={\displaystyle {\int }_0^{t}E\left(t-\tau \right)}\frac{d\epsilon \left(t\right)}{d\tau }d\tau$$

(4)

$$\epsilon \left(t\right)={\displaystyle {\int }_0^{t}D\left(t-\tau \right)}\frac{d\sigma }{d\tau }d\tau$$

(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]:

$${\displaystyle {\int }_0^{t}E\left(t-\tau \right)}\frac{dD\left(\tau \right)}{d\tau }d\tau =1$$

(6)

Then, the relationship of relaxation modulus and creep compliance in Laplace transform domain [15] can be obtained:

$$\overline{E}\left(s\right)\overline{D}\left(s\right)=1/s^2,\tilde{E}\left(s\right)\tilde{D}\left(s\right)=1$$

(7)

where $\overline{f}\left(s\right)$ is the Laplace transformation of $f\left(t\right)$, $\overline{f}\left(s\right)={\displaystyle {\int }_0^{\infty }f\left(t\right)}\cdot e^{-st}dt$; s is the Laplace operator; $\tilde{E}\left(s\right)$ is the Laplace transform of relaxation modulus, $\tilde{E}\left(s\right)=s\overline{E}\left(s\right)$; $\tilde{D}\left(s\right)$ is the Laplace transform of creep compliance, $\tilde{D}\left(s\right)=s\overline{D}\left(s\right)$.

Since the complex function of viscoelastic materials is excited by sinusoidal loading under $\omega$ angular frequency, the relationship between it and Laplace transformed function can be found as [25]:

$$E^{*}\left(\omega \right)=\tilde{E}{\left(s\right)|}_{s\to i\omega },D^{*}\left(\omega \right)=\tilde{D}{\left(s\right)|}_{s\to i\omega }$$

(8)

where E* is the complex modulus; D* is the complex compliance.

Thus, relationship between $E^{*}$ and $D^{*}$ can be given by:

$$E^{*}\left(\omega \right)D^{*}\left(\omega \right)=1$$

(9)

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.

2.2.1. 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.

$$E^{*}=E^{\prime }+i{E}^{″}$$

(10)

In Equation (10), $E^{\prime }$ 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:

$$E^{\prime }\left(\omega \right)=E_e+{\displaystyle \sum _{i=1}^m\frac{{\omega }^2{\rho }_i^2{E}_i}{{\omega }^2{\rho }_i^2+1}},E^{″}\left(\omega \right)={\displaystyle \sum _{i=1}^m\frac{\omega {\rho }_i{E}_i}{{\omega }^2{\rho }_i^2+1}}$$

(11)

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{\prime }_{expriment}$ based on experimental data;

(3)

Make the storage modulus $E{\prime }_{predict}$ determined by Equation (11) equal to the $E{\prime }_{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.

2.2.2. 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:

$${\displaystyle {\int }_0^t{E}_e}+{\displaystyle \sum _i^m{E}_i}{e}^{-\left[\left(t-\tau \right)/{\rho }_i\right]}\left(D_g\delta \left(t\right)+{\displaystyle \sum _j^n\frac{D_j}{{\tau }_j}{e}^{-\left(t/\tau \right)}}\right)d\tau =1$$

(12)

where $\delta ( )$ is the Dirac delta function.

According to the literature [15], further simplifying the integral can derive the following equations:

$$\{\begin{array}{l}{\displaystyle \sum _{i=1}^m\left[{\displaystyle \sum _{j=1}^n\frac{{\rho }_i{D}_j}{{\rho }_i-{\tau }_j}\left(e^{-\left(t/{\rho }_i\right)}-e^{-\left(t/{\tau }_j\right)}\right)+D_g{e}^{-\left(t/{\rho }_i\right)}}\right]}{E}_i=1-\left(D_g+{\displaystyle \sum _{j=1}^n{D}_j\left(1-e^{-\left(t_k/{\tau }_j\right)}\right)}\right)E_e,{\rho }_i\ne {\tau }_j\\ {\displaystyle \sum _{i=1}^m\left[{\displaystyle \sum _{j=1}^n\frac{t{D}_j}{{\tau }_j}{e}^{-\left(t/{\rho }_i\right)}+D_g{e}^{-\left(t/{\rho }_i\right)}}\right]D_j=1-\left(D_g+{\displaystyle \sum _{j=1}^n{D}_j\left(1-e^{-\left(t_k/{\tau }_j\right)}\right)}\right)E_e,{\rho }_i={\tau }_j}\end{array}$$

(13)

where,

$$E_e\equiv \underset{t\to \infty }{\lim }E\left(t\right)=\underset{s\to 0}{\lim }\tilde{E}\left(s\right)=\frac{1}{\underset{s\to 0}{\lim }\tilde{D}\left(s\right)}=\frac{1}{\underset{t\to \infty }{\lim }D\left(t\right)}=\frac{1}{D_g+{\displaystyle \sum _{j=1}^n{D}_j}}$$

(14)

Adjusting them into linear system of equations, Equation (15) can be given as:

$$\begin{gathered} \left[A\right]\left[E\right]=\left[B\right] \\ {A}_{ki}=\{\begin{array}{l}{D}_g{e}^{-\left(t_k/{\rho }_i\right)}+{\displaystyle \sum _{j=1}^n\frac{{\rho }_i{D}_j}{{\rho }_i-{\tau }_j}}\left(e^{-\left(t_k/{\rho }_i\right)}-e^{-\left(t_k/{\tau }_j\right)}\right),{\rho }_i\ne {\tau }_j\\ D_g{e}^{-\left(t_k/{\rho }_i\right)}+{\displaystyle \sum _{j=1}^n\frac{t_k{D}_j}{{\tau }_j}\left(e^{-\left(t_k/{\rho }_i\right)}\right),{\rho }_i={\tau }_j}\end{array} \\ {B}_{kj}=1-\left(D_g+{\displaystyle \sum _{j=1}^n{D}_j\left(1-e^{-\left(t_k/{\tau }_j\right)}\right)}\right)/\left(D_g+{\displaystyle \sum _{j=1}^n{D}_j}\right) \end{gathered}$$

(15)

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\left(k=1,2,...,p\right)$ 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 ${\Vert AE-B\Vert }^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]⁻¹[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.

3. Numerical Examples

3.1. 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. Aggregate Gradation of Hot Mix Asphalt AC-13.

Sieve Size (mm)	19	12.5	9.5	4.75	2.36	1.18	0.6	0.3	0.15	0.075
Passing Rate (%)	100	93.0	74.8	50.3	30.3	24.6	19.4	15.0	10.3	6.5

.

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 με. The test was conducted under four different temperatures of 4.4 °C, 21.1 °C, 37.8 °C and 54.4 °C according to AASHTO PP61 for the preparation of plotting master curve of dynamic modulus. Also, nine different frequencies of 0.1Hz, 0.2 Hz, 0.5 Hz, 1.0 Hz, 2.0 Hz, 5.0 Hz, 10.0 Hz, 20.0 Hz and 25.0 Hz were under sinusoidal loading at each temperature. The loading time is 100 s and time intervals are 5 min at each temperature. The experimental data of dynamic modulus and phase angle was determined by the average of four LVDTs and the results were the average of three repeated experiments.

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$\mathsf{\mu }\mathsf{\epsilon }$. The results were also the average of three repeated experiments.

3.2. 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. Summary of fitting parameters for the Hot Mix Asphalt AC-13.

Fitting parameters, dynamic modulus master curve (T0 = 21.1 °C)	α	β	δ	γ
	1.3530	2.8961	0.7631	0.6983
Fitting parameters, phase angle master curve (T0 = 21.1 °C)	ξ1	ξ2	ξ3
	−0.6555	−0.0008405	−0.59608
Time-temperature shifting factor ($\log \left({\alpha }_T\right)$)	4.4 °C	21.1 °C	37.8 °C	54.4 °C
	1.8454	0	−1.3743	−2.9746

.

$$\log \left|E^{*}\right|=\delta +\frac{\alpha }{\left(1+e^{\beta +\gamma \log \left(f_r\right)}\right)}$$

(16)

$$\delta \left(f_r\right)={\xi }_1+{\xi }_2\log \left(f_r\right)-{\xi }_3\frac{\pi }{2}\alpha \gamma \frac{e^{\beta +\gamma \log \left(f_r\right)}}{{\left(1+e^{\beta +\gamma \log \left(f_r\right)}\right)}^2}$$

(17)

$$\log \left(f_r\right)=\log \left(f\right)+\log \left({\alpha }_T\right)$$

(18)

where, $f_r$ is the reduced frequency; $f$ is the loading frequency; δ is the minimum modulus; $\alpha$ is the parameter of vertical span of dynamic modulus; β, γ are the shape parameters of the material; ${\xi }_1,{\xi }_2,{\xi }_3$ are all fitting parameters of experimental data; ${\alpha }_T$ is the time-temperature shifting factor.

The determination of ${\alpha }_T$ is based on Williams-Landel-Ferry (WLF) model [28] in Equation (19).

$$\log \left({\alpha }_T\right)=\frac{-C_1\left(T-T_0\right)}{C_2+(T-T_0)}$$

(19)

where C₁ and C₂ are the WLF coefficients; T is the testing temperature; T₀ 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.

3.3. Interconversions

3.3.1. 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 ${\rho }_i$ is expected to be identified first. According to reference [3], the distribution of ${\rho }_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 ${\rho }_i=2\times {10}^{i-7}s, i=1,2,...,11$. Consequently, the Prony series coefficients converted by dynamic modulus are listed in Table 2.

3.3.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 ${\rho }_i$ of dynamic modulus, the same number of elements 11 is selected to describe retardation times ${\tau }_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\times {10}^{\left(k-1\right)/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, ${\tau }_j$ was firstly determined and it was assumed that ${\tau }_j$ and ${\rho }_i$ are equivalent as suggested [33]. In this paper, ${\tau }_j$ was selected between 2 × 10⁻⁶ to 2 × 10⁻⁴ to cover the entire range of relaxation modulus. Theoretically, the relaxation times ${\rho }_i$ are the negative reciprocal of s when $\tilde{E}\left(s\right)$ approaches to infinite [15]. The ${\rho }_i$ can also be determined by root-finding method based on $\left|\tilde{E}\left(s\right)\right|$ versus $-1/s$ figure, whereas the calculated ${\rho }_i$ was almost identical by using the two methods: assuming ${\tau }_j$ is equivalent to ${\rho }_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. Prony Series Coefficients of Relaxation Modulus Obtained from Dynamic Modulus and Creep Compliance.

No.	${\mathit{\tau }}_{\mathit{j}}\left(\mathit{s}\right)$	${\mathit{D}}_{\mathit{j}}\left(\mathbf{1}/\mathbf{MPa}\right)$	${\mathit{\rho }}_{\mathit{i}}\left(\mathit{s}\right)$	${\mathit{E}}_{\mathit{i}}\left(\mathbf{MPa}\right)$
				Dynamic to Relaxation	Creep to Relaxation
					Collocation Method	Least Squares
1	2 × 10−6	2.14 × 10−6	2 × 10−6	2888.60	1542.07	2339.76
2	2 × 10−5	9.26 × 10−6	2 × 10−5	3858.85	7967.98	7236.69
3	2 × 10−4	2.07 × 10−5	2 × 10−4	4791.64	6962.31	7381.14
4	2 × 10−3	3.91 × 10−5	2 × 10−3	4175.53	6426.10	6218.63
5	2 × 10−2	1.00 × 10−4	2 × 10−2	3510.67	4146.14	4243.43
6	2 × 10−1	2.89 × 10−4	2 × 10−1	1669.80	1712.93	1750.13
7	2 × 100	8.35 × 10−4	2 × 100	530.77	616.41	626.74
8	2 × 101	2.63 × 10−3	2 × 101	144.70	197.50	200.97
9	2 × 102	6.64 × 10−3	2 × 102	42.56	64.27	67.28
10	2 × 103	1.30 × 10−2	2 × 103	18.22	34.14	32.71
11	2 × 104	2.92 × 10−2	2 × 104	9.63	7.63	8.84
	$D_g\left(1/\mathrm{MPa}\right)$	3.33×10−5	$E_e\left(\mathrm{MPa}\right)$	21.83	18.94	18.94

Note: τ_j and ρ_i are the relaxation time for creep compliance and dynamic modulus, respectively; D_j is the Prony series for creep compliance.

.

3.4. 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. 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.

3.5. 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. Prony Series Coefficients of Relaxation Modulus Obtained from Dynamic Modulus and Creep Compliance for new data points from previous publications.

No.	${\mathit{\rho }}_{\mathit{i}}\left(\mathit{s}\right)$	${\mathit{E}}_{\mathit{i}}\left(\mathbf{MPa}\right)$
		Evotherm		Rubber		Polymer
		Dynamic to Relaxation	Creep to Relaxation	Dynamic to Relaxation	Creep to Relaxation	Dynamic to Relaxation	Creep to Relaxation
1	2 × 10−6	2887.79	1628.27	2757.41	2889.77	2888.43	2677.87
2	2 × 10−5	4120.24	7975.11	3796.04	7164.79	3842.27	6895.38
3	2 × 10−4	4511.36	7458.14	4191.18	7465.54	4325.67	7378.90
4	2 × 10−3	3532.57	4518.68	3605.38	5829.68	3549.11	5891.67
5	2 × 10−2	2838.39	3347.47	3139.01	4460.67	3207.17	3560.67
6	2 × 10−1	1558.67	1965.74	1744.32	2274.53	1663.13	1633.28
7	2 × 100	632.04	788.14	672.53	793.78	652.71	755.31
8	2 × 101	202.66	234.71	198.79	227.54	213.06	261.15
9	2 × 102	71.68	69.27	47.66	78.26	71.329	85.80
10	2 × 103	25.76	22.50	39.76	42.71	30.29	39.23
11	2 × 104	3.01	7.09	8.32	9.45	6.90	11.19
$E_e\left(\mathrm{MPa}\right)$		22.38	22.90	23.15	23.66	23.14	25.46

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.

4. 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.

5. 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.