Next Article in Journal
Effect of Twinning on Angle-Resolved Photoemission Spectroscopy Analysis of Ni49.7Mn29.1Ga21.2(100) Heusler Alloy
Next Article in Special Issue
Influence of Basalt Fibers on the Crack Resistance of Asphalt Mixtures and Mechanism Analysis
Previous Article in Journal
Suspended Metasurface for Broadband High-Efficiency Vortex Beam Generation
Previous Article in Special Issue
Fatigue Models for Airfield Concrete Pavement: Literature Review and Discussion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Fast-Acquiring High-Quality Prony Series Parameters of Asphalt Concrete through Viscoelastic Continuous Spectral Models

1
City Institute, Dalian University of Technology, Dalian 116600, China
2
School of Transportation and Logistics, Dalian University of Technology, Dalian 116024, China
*
Author to whom correspondence should be addressed.
Materials 2022, 15(3), 716; https://doi.org/10.3390/ma15030716
Submission received: 6 December 2021 / Revised: 14 January 2022 / Accepted: 17 January 2022 / Published: 18 January 2022
(This article belongs to the Special Issue Experimental Testing and Constitutive Modelling of Pavement Materials)

Abstract

:
Prony series representations have been extensively applied to characterizing the time-domain linear viscoelastic (LVE) material functions for asphalt concrete. However, existing methods that can generate high-quality Prony series parameters (i.e., discrete spectra) mostly involve complicated programming algorithms, which poses a challenge for quick access of Prony series parameters. Also, very limited research has been devoted to establishing methods for simultaneously determining both retardation and relaxation spectra. To resolve these issues, this study presented a practical approach to fast acquiring high-quality Prony series parameters for both relaxation modulus and creep compliance of asphalt concrete by using the complex modulus test data. The approach adopts the analytical representations of the continuous relaxation and retardation spectra from the Havriliak-Negami (HN) and 2S2P1D complex modulus models to directly determine the discrete spectra, and the elastic constants, Ee and Dg, for both LVE modulus and compliance functions are further calculated by fitting the corresponding generalized Maxwell model representations to smoothed data from the storage modulus representations of the HN and 2S2P1D complex modulus models. In this way, all the procedures in the proposed method can be easily implemented in Microsoft Excel. The results showed that the HN and 2S2P1D models yielded slightly different continuous spectral patterns at shorter relaxation times and longer retardation times. However, at the region covered by the test data, the continuous spectra of the two complex modulus models were very close to each other. Thus, the two models can generate comparable Prony series parameters within the time or frequency range covered by the test data. Considering that the quality of the resulting Prony series parameters are closely related to the master curve models used for presmoothing, the HN and 2S2P1D models were compared with the conventional Sigmoidal model. Additionally, the Black diagram was recommended for examining the quality of the complex modulus test data before constructing the master curves.

1. Introduction

Asphalt concrete, which has been paved on most roadways in the world, is a typical particulate composite with a viscoelastic matrix. In engineering applications, it is commonly regarded as a linear viscoelastic (LVE) material [1,2,3]. As such, many mechanical tests based on the LVE theory, like the static relaxation and creep tests and the dynamic complex modulus test, can be used for characterizing its LVE behavior. Theoretically, the properties from these tests such as the relaxation modulus, creep compliance and complex modulus are equivalent [4,5,6]; however, for a practical purpose, the uniaxial compressive complex modulus test has been widely accepted as a standard LVE material characterization test. After the complex modulus test data is obtained, it is usually required to extract the LVE information from the test data through mathematical models.
In the LVE theory, the generalized Maxwell model and generalized Kelvin model appear to be the most commonly used models for describing both time- and frequency-domain material functions, and they have been implemented into many commercial numerical simulation programs, e.g., ABAQUS, ANSYS and COMSOL Multiphysics [7,8,9,10,11]. This can be primarily attributed to their high computational efficiency and wide applicability [3,12,13]. The two models are composed of linear springs and dashpots linked in different configurations, and mathematically yield the so-called Prony series expressions for the relaxation modulus and creep compliance in time domain [3,12,14]. The Prony series expressions are not only very convenient to be converted analytically into the frequency-domain complex modulus and compliance, but also can considerably facilitate the computation of the convolution integrals for the LVE constitutive equations due to the presence of decaying exponential terms. Therefore, accurate and efficient identification of the Prony series parameters (i.e., discrete relaxation and retardation spectra) is crucial to the subsequent performance analysis and prediction of asphalt pavement or mixtures.
To date, researchers have proposed various methods for determining the Prony series parameters. Several representative approaches that apply directly to raw data in the time or frequency domain have been widely used for LVE materials, e.g., the collocation method by Schapery [15], the multidata method by Cost and Becker [16], and the windowing method by Tschoegl and Emri [17]. Nonetheless, these classic schemes would encounter difficulties when utilized for asphalt concrete. Two major issues, namely negative spectrum strengths and local spectrum oscillations, occur frequently due to the narrowband nature of the Prony series terms and significant scatters in test data. To address these problems, presmoothing techniques have been introduced by using broadband functions, like the Sigmoidal model [2], power-law series [3], Huet-Sayegh model [18], Havriliak-Negami (HN) model [14] and 2S2P1D model [19,20]. The use of these broadband functions not only improves the quality of the test data, but facilitates the data shift in accordance with the time-temperature superposition principle (TTSP) during the construction of master curves. In view of the equivalence of the LVE material functions, some interconversion algorithms [13,21,22] were also presented to calculate the Prony series parameters of the relaxation functions from the retardation functions, or vice versa.
On the other hand, the continuous spectrum-based methods attract increasing attention from the asphalt paving research community in that they are able to eliminate negative spectrum strengths and excessive parameters. Levenberg [23] developed a continuous relaxation spectrum model for asphalt concrete by using a lognormal distribution function. However, this model is symmetrical on the logarithmic timescale and thus may not be appropriate for all mixtures. Zhao et al. [2] established a confining pressure dependent continuous relaxation spectrum by considering the relationship between the relaxation spectrum and storage modulus. Luo et al. [24] and Lv et al. [25] respectively deduced continuous relaxation spectra from a modified power law-based relaxation modulus model and a generalized Sigmoidal model-based storage modulus model. Nevertheless, these works were all concentrated on the relaxation spectrum, and thus may be inconvenient for those who need fast solutions for both retardation and relaxation functions. Aiming at this issue, Sun et al. [26] presented a numerical approach to determining a continuous spectrum from the other. Bhattacharjee et al. [27] and Zhang et al. [28] calculated the two continuous spectra from storage modulus and storage compliance separately based on the Sigmoidal function and the generalized Sigmoidal function; however, due to the inconsistency of the model parameters of the storage modulus and storage compliance, the LVE relationship cannot be strictly satisfied.
Although there have been so many methods developed for determining the Prony series parameters as mentioned above, most of them involve complicated programming algorithms, which poses a challenge for quick access of Prony series parameters. Furthermore, very limited research has been devoted to establishing approaches for determining both retardation and relaxation spectra at the same time. To deal with these problems, this study gave a practical approach by adopting analytical representations of the continuous relaxation and retardation spectra from two complex-valued models, and all the procedures in the method can be easily implemented in Microsoft Excel.

2. Materials and Complex Modulus Test

Two dense-graded asphalt mixtures, denoted as Mix-13.2 and Mix-9.5 herein, were prepared for the complex modulus testing. Mix-13.2 had a nominal maximum aggregate size (NMAS) of 13.2 mm. The coarse aggregates and asphalt binder used were limestone and PG 58-22 unmodified asphalt, respectively. Mix-9.5 had a NMAS of 9.5 mm. The coarse aggregates and asphalt binder used were granite and PG 64-22 neat asphalt, respectively. The asphalt contents of Mix-13.2 and Mix-9.5 were 3.9% and 5.7%, respectively. Figure 1 presents the aggregate gradations of the two asphalt mixtures.
The complex modulus tests were performed on all specimens of the two mixtures in accordance with the standard testing method AASHTO T342 [29]. The two mixtures were first compacted using the Superpave Gyratory Compactor and then trimmed into the final cylindrical specimens (150 mm in height and 100 mm in diameter) containing an air void content of 4 ± 1%. For each mixture, three replicate specimens were fabricated.
The complex modulus testing was conducted on a universal testing machine (UTM). A stress-controlled compressive mode was employed for all the complex modulus tests. For Mix-13.2, five testing temperatures (−10, 5, 20, 35 and 50 °C) and seven loading frequencies (0.1, 0.5, 1, 5, 10, 20 and 25 Hz) were adopted, and for Mix-9.5, five testing temperatures (−16, 4, 24, 40 and 50 °C) and six loading frequencies (0.1, 0.5, 1, 5, 10 and 25 Hz) were adopted. During testing, the strain was kept between 50~150 με and the accumulated strain was controlled below 1500 με to ensure LVE measurements. By means of the obtained stress and strain data, two quantities, i.e., the dynamic modulus |E*| and phase angle φ, can be calculated as follows:
E = σ 0 ε 0
φ = Δ t t p × 360
where σ0 and ε0 are the amplitudes of the axial stress and strain; Δt is the time lag of the strain curve behind the stress curve; tp is the loading period. The dynamic modulus |E*|, which is the absolute value of the complex modulus E*, characterizes the resistance to deformation of a viscoelastic material, whereas the phase angle φ characterizes the extent to which the viscoelastic material behaves like a viscous liquid (φ = 90°) or an elastic solid (φ = 0°).

3. Methodology

3.1. Viscoelastic Master Curve Models

In this study, two models, i.e., the HN model and the 2S2P1D model, were employed to build the complex modulus master curves of the asphalt mixtures. The HN model has five parameters and is represented by [30]:
E ω = E g + E e E g 1 + i ω τ 0 α β
where Eg is the glassy modulus; Ee is the equilibrium modulus; ω = 2πf is the angular frequency; f is the frequency; i = 1 ; α, β and τ0 are model parameters and they respectively control the width, asymmetry and horizontal position of the relaxation spectrum.
The real part E′ and imaginary part E″ of the complex modulus E* have the following relationship with the dynamic modulus |E*| and phase angle φ:
E = E + i E = E cos φ + i E sin φ
where E′ is the storage modulus; E″ is the loss modulus. Thus, the dynamic modulus |E*| and phase angle φ can be calculated using the storage modulus E′ and the loss modulus E″, as follows
E = E 2 + E 2   and   φ = arctan E E
From the HN model, the representations of the storage modulus and loss modulus can be analytically separated out according to De Moivre’s formula, as the following [30]:
E = E g + E e E g cos β ψ 1 + 2 ω α τ 0 α cos α π / 2 + ω 2 α τ 0 2 α β / 2
E = E g E e sin β ψ 1 + 2 ω α τ 0 α cos α π / 2 + ω 2 α τ 0 2 α β / 2
ψ = arctan ω α τ 0 α sin α π / 2 1 + ω α τ 0 α cos α π / 2
Obviously, with the analytical expressions of the storage and loss moduli, those for the dynamic modulus and phase angle are also available according to Equation (5).
The 2S2P1D model, composed of two spring elements, two parabolic elements and a dashpot element, possesses seven parameters and has the following mathematical form [19]:
E ω = E e + E g E e 1 + α i ω τ 0 k + i ω τ 0 h + i ω β τ 0 1
where α, k and h (0 < k < h < 1) are the parameters of the two parabolic elements; τ0 is a temperature-dependent parameter; β is the parameter associated with the Newtonian viscosity of the dashpot element, η = (EgEe) βτ0.
According to De Moivre’s formula, the storage and loss moduli representations of the 2S2P1D model can also be derived, as follows:
E = E e + E g E e 1 + A 1 + A 2 + B 2
E = E e E g B 1 + A 2 + B 2
A = α ω τ 0 k cos k π / 2 + ω τ 0 h cos h π / 2
B = α ω τ 0 k sin k π / 2 ω τ 0 h sin h π / 2 ω β τ 0 1
Since the storage and loss moduli of the HN and 2S2P1D models can all be derived analytically from the corresponding complex-valued models, they accurately meet the Kronig–Kramers relation that correlates the real and imaginary parts of the response to a harmonic load to each other theoretically [31].
Besides, in the viscoelastic theory, when the representation of the complex modulus E* is known, the complex compliance D* can be analytically obtained by taking the inverse of the complex modulus, as follows:
D = D i D = 1 / E
D = E E 2 + E 2
D = E E 2 + E 2
where D′ is the storage compliance; D″ is the loss compliance.

3.2. Continuous Relaxation and Retardation Spectra

For a LVE material, the modulus functions in the time and frequency domains can be uniformly expressed using the continuous relaxation spectrum H(ρ) through integral forms [31]:
E ω = E e + H ρ ω 2 ρ 2 1 + ω 2 ρ 2 d ln ρ
E ω = H ρ ω ρ 1 + ω 2 ρ 2 d ln ρ
E t = E e + H ρ e t / ρ d ln ρ
where ρ is the relaxation time; E(t) is the relaxation modulus; t is the loading time.
Similarly, the compliance functions of a LVE material in the time and frequency domains can be uniformly expressed using the continuous retardation spectrum L(τ) through integral forms [31]:
D ω = D g + L τ 1 1 + ω 2 τ 2 d ln τ
D ω = L τ ω τ 1 + ω 2 τ 2 d ln τ
D t = D g + L τ 1 e t / τ d ln τ
where τ is the retardation time; D(t) is the creep compliance.
The continuous relaxation spectrum H(ρ) and continuous retardation spectrum L(τ) essentially contain identical time- and frequency-dependent material information; thus, they are equivalent of characterizing the LVE behavior of a material. As can be seen from Equations (17)–(22), once the continuous relaxation and retardation spectra are determined, both modulus and compliance functions in the time and frequency domain can be attained. In accordance with the LVE theory, the continuous spectra have the following relationships with the complex modulus E* [31]:
H ρ = ± π 1 Im E i ω i ω ρ 1 e ± i π = ± π 1 Im E ρ 1 e ± i π
L τ = π 1 Im D i ω i ω τ 1 e ± i π = π 1 Im D τ 1 e ± i π = π 1 Im E τ 1 e ± i π 1
where Im represents the operation of retaining the imaginary part of a complex-valued function.
For the HN model, Havriliak and Negami [30] presented the analytical expression of H(ρ) through Equation (23), as follows:
H ρ = E g E e ρ / τ 0 α β sin β ϕ π 1 + ρ / τ 0 2 α + 2 ρ / τ 0 α cos α π β / 2
ϕ = arctan sin ( α π ) ρ / τ 0 α + cos α π
Analogously, Sun et al. [26] derived the close-form solution for L(τ) of the HN model through Equation (24), as follows:
L τ = Ω E g E e sin ( β ϕ ) π E e E g cos β ϕ + E g Ω 2 + E e E g sin β ϕ 2
Ω = 1 + τ 0 / τ 2 α + 2 τ 0 / τ α cos α π β / 2
For the 2S2P1D model, Alavi et al. [32] obtained the analytical representation of H(ρ) through Equation (23), as follows:
H ρ = E g E e X π X 2 + Y 2
X = 1 + α τ 0 k ρ k cos k π + τ 0 h ρ h cos h π β 1 τ 0 1 ρ
Y = α τ 0 k ρ k sin k π + τ 0 h ρ h sin h π
Sun et al. [33] successfully deduced the analytical expression for L(τ) of the 2S2P1D, as follows:
L ( τ ) = Y ( E g E e ) ( X 2 + Y 2 ) π E e ( X 2 + Y 2 ) + X ( E g E e ) 2 + Y ( E g E e ) 2
It is noted that ρ in X and Y of Equation (32) should be replaced by τ.
Evidently, for the HN and 2S2P1D models, the corresponding relaxation modulus E(t) and creep compliance D(t) in the time domain can be readily calculated with the continuous relaxation and retardation spectra according to Equations (19) and (22). Further, in terms of the Boltzmann superposition integrals, the constitutive relationships for the LVE material can be determined [31].

3.3. Construction of Master Curves

Asphalt concrete is a typical thermorheologically simple material in the LVE region; therefore, the master curves for various LVE material functions in both frequency and time domains can be constructed in accordance with the time–temperature superposition principle (TTSP). During this process, viscoelastic test data measured at different temperatures is shifted horizontally along the frequency or time axis on the logarithmic scale, thus generating a smooth master curve at a given reference temperature Tr. By means of the constructed master curve, the LVE behavior over a wider range of loading time or frequency than that offered by the test instrument can be predicted. The reduced angular frequency ωr and reduced time tr for the shifted test data are represented by:
ω r = ω × α T
t r = t α T
where αT is the time-temperature shift factor. The time–temperature shift factors can be represented using a function of temperature, e.g., the Williams–Landel–Ferry (WLF) or the Arrhenius equation, or in a non-functional form. To avoid the effect of the functional expression of αT, the non-functional method was adopted for constructing the master curve of the complex modulus in the present study.
The parameters of the complex modulus model and time-temperature shift factors were determined simultaneously through a nonlinear optimization process. To fully extract the LVE information, both dynamic modulus and phase angle test data were taken into account, and the target error function to minimize was as the following:
F = 1 N i = 1 N E m , i E c , i E m , i 2 + 1 N i = 1 N φ m , i φ c , i φ m , i 2
where N is the number of the dynamic modulus or phase angle data points; E m , i and φ m , i are the measured values for the dynamic modulus and phase angle, respectively; E c , i and φ c , i are the calculated values for the dynamic modulus and phase angle from the master curve model used, respectively. The optimization operation can be easily completed using the Solver in Microsoft Excel. Before this, initial values for both master curve model parameters and shift factors should be given. The reference temperatures for Mix-13.2 and Mix-9.5 were set to 20 and 24 °C, respectively.

3.4. Determination of Prony Series Parameters

As stated above, once the parameters of the complex-valued models, like HN and 2S2P1D models, are known, the corresponding H(ρ) and L(τ) can be automatically determined due to the existence of their analytical expressions with the same parameters as the original complex modulus models. Although all the modulus and compliance functions can further be straightforward calculated with H(ρ) and L(τ) through Equations (17)–(22), the integral forms based on the continuous spectra are actually inconvenient to implement in numerical simulation techniques, e.g., the finite element method. Instead, the Prony series expressions on the basis of discrete spectra have been extensively utilized due to their advantage at computation efficiency.
The relaxation modulus expression derived from the generalized Maxwell model and the creep compliance expression derived from the generalized Kelvin model are two typical Prony series representations. For the generalized Maxwell model, the modulus functions can be formulated by [31]:
E t = E e + j = 1 n E j e t / ρ j
E ω = E e + j = 1 n E j ω 2 ρ j 2 1 + ω 2 ρ j 2
E ω = j = 1 n E j ω ρ j 1 + ω 2 ρ j 2
where Ej is the modulus of the spring or the relaxation strength; ρj = ηj/Ej is the discrete relaxation time; ηj is the viscosity of the dashpot; the set of Prony series parameters [ρj, Ej] is called the discrete relaxation spectrum.
For the generalized Kelvin model, the compliance functions can be represented by [31]:
D t = D g + j = 1 n D j 1 e t / τ j
D ( ω ) = D g + j = 1 n D j 1 1 + ω 2 τ j 2
D ( ω ) = j = 1 n D j ω τ j 1 + ω 2 τ j 2
where Dj is the compliance of the spring or the retardation strength; τj = λjDj is the discrete retardation time; λj is the viscosity of the dashpot; the set of Prony series parameters [τj, Dj] is called the discrete retardation spectrum.
In fact, when the discrete relaxation and retardation spectra become infinitely dense, they evolve into the so-called continuous relaxation and retardation spectra. As such, Equations (36)–(41) can be interpreted as discretizations of Equations (17)–(22). For the storage modulus E′(ω), H(ρ)dlnρ in Equation (17) represents the contribution of the model to the modulus function in the interval of lnρ and lnρ+dlnρ, which leads to the following derivation:
E ω = E e + H ρ ω 2 ρ 2 1 + ω 2 ρ 2 d ln ρ E e + j = 1 n H ρ j × Δ ln ρ j ω 2 ρ j 2 1 + ω 2 ρ j 2 = E e + j = 1 n E j ω 2 ρ j 2 1 + ω 2 ρ j 2
Likewise, for the storage compliance D′(ω), the integral form based on the continuous retardation spectrum and the series expression based on the discrete retardation spectrum have the following relationship:
D ω = D g + L τ 1 1 + ω 2 τ 2 d ln τ D g + j = 1 n L τ j × Δ ln τ j 1 1 + ω 2 τ j 2 = D g + j = 1 n D j 1 1 + ω 2 τ j 2
During the determination of the Prony series parameters, the discrete time constants, ρj and τj, are commonly preselected. Specifically, they are set to values with equal intervals on the logarithmic scale according to Equation (44):
ρ i = τ i = b × 10 d + i / M
where b and d are specified according to the logarithmic time range covered by the shifted test data, and generally b = 1; M is the number of the discrete times assumed in each decade on the logarithmic scale.
It can be observed that with the discrete time constants (ρj and τj) known, the Prony series coefficients, namely the relaxation and retardation strengths (Ej and Dj), can be quickly and easily calculated using the following equations:
E i = H ρ i × Δ ln ρ i
D i = L τ i × Δ ln τ i
Δ ln ρ i = Δ ln τ i = 1 M ln 10
Finally, the remaining two elastic constants Ee and Dg can be determined by fitting Equations (37) and (40) to the corresponding real part expression of the original complex modulus model, like Equation (6) or (10), over the range covered by the shifted test data through the Excel Solver. In such a manner, all the Prony series parameters can be fast acquired. Actually, Ee and Dg have an analytical relationship, as follows:
D g = 1 E g = 1 E e + j = 1 n E j
Therefore, once the elastic constant Ee along with the discrete relaxation strengths Ej is available, Dg can be obtained accordingly.

4. Results and Discussion

4.1. Examination of Test Data Quality of Asphalt Concrete

Before constructing the master curves, it is crucial to examine the quality of the complex modulus measurements. In the present study, the Black diagram [34] was employed to conduct this manipulation, in which the dynamic modulus |E*| is plotted against the phase angle φ in a single plane. Since for a thermorheologically simple material, all the components of the complex modulus are the functions of the reduced angular frequency, any two of them can form a unique curve in a complex plane. In the Black diagram, the angular frequency axis can be treated as an additional axis perpendicular to the complex plane in accordance with the right-hand rule. Thus, the testing temperatures would have no effect on the analysis of the overlapping behavior of the test data during the construction of master curves in the Black diagram. A smoother Black curve generally represents a higher quality of the test data. In such a manner, the Black diagram allows an effective and efficient detection of inconsistency with thermorheological simplicity.
Figure 2 shows the resulting Black diagrams for the two asphalt mixtures. It can be observed that in both diagrams, the complex modulus test data obtained at different temperatures basically formed unique curves, indicating the compliance with thermorheological simplicity under the test conditions as well as the applicability of the TTSP. In addition, the test results at lower temperatures exhibited better overlapping behavior, whereas those at higher temperatures showed slightly higher dispersion. This is mainly because nonlinear behaviors (e.g., the viscoplastic deformation) of asphalt concrete occur more easily at higher temperatures, which impact the measurement of LVE responses of the material to a certain degree.

4.2. Analysis of Results from the Developed Method

Figure 3 and Figure 4 present the master curves of the dynamic modulus and phase angle respectively developed from the HN and 2S2P1D models. As observed, both models fitted to the test data of the two mixtures very well. Table 1 and Table 2 list the resulting model parameters and fitting errors. For Mix-13.2, the two complex modulus models contributed to very close fitting errors, whereas for Mix-9.5, the 2S2P1D model yielded slightly lower fitting error than that from the HN model. This may be because the 2S2P1D model has more parameters and thus higher flexibility. Besides, the time-temperature shift factors calculated using the HN and 2S2P1D methods were found very close as well, as shown in Figure 5.
As mentioned previously, with the obtained complex modulus model parameters, both continuous relaxation and retardation spectra can be analytically developed [see Equations (25)–(32)]. Figure 6 shows the continuous spectra of the two asphalt mixtures. It can be observed that for both mixtures, the HN and 2S2P1D models exhibited slightly different continuous spectral patterns, particularly at shorter relaxation times and longer retardation times. However, at the time range of 10−8 to 104 s, which is approximately corresponding to the angular frequency range of 10−4 to 108 rad/s, that is, the region mostly covered by the test data (Figure 3 and Figure 4), the continuous spectra for the HN and 2S2P1D models were very close to each other.
Based on the continuous spectra developed, the corresponding discrete spectra can be fast determined using Equations (44)–(47). Although the relaxation and retardation times of the discrete spectra can be preset at any time regions of interest with any widths, it is a common practice that they are selected at regions covered by test data [3,12,35]. In this way, the numbers of the Prony series parameters can be reduced reasonably without losing significant computation accuracy. Consequently, the range of the discrete spectra was selected at 10−8 to 104 s in this study.
Figure 7 and Figure 8 display the calculated discrete relaxation and retardation spectra for the two mixtures. Three densities of the discrete spectrum lines, namely, M = 1, 2 and 3, were considered. As can be seen, the resulting discrete spectra from both HN and 2S2P1D models were very smooth without any local oscillations. Also, since the discrete spectrum strengths were all calculated from the corresponding positive continuous spectra, no negative strength values were produced.
To establish the Prony series representations for the relaxation modulus E(t) and creep compliance D(t) in Equations (36) and (39), the elastic constants Ee and Dg need to be further determined. To this end, Equation (37) for the storage modulus E′ was fitted to Equations (6) and (10) separately for the real parts of the HN and 2S2P1D models. Before fitting, smoothed data points were generated from Equations (6) and (10), equally spaced on the logarithmic scale within the region covered by the test data. With Ee determined, Dg can be fast obtained by Equation (48).
Figure 9 gives the developed master curves of the storage modulus for Mix-9.5 using the discrete relaxation spectrum, i.e., using the generalized Maxwell model, with M = 1 from the HN and 2S2P1D models, respectively. Obviously, both methods yielded satisfactory results over the region where the spectrum lines were selected. Similar observations were made for Mix-13.2. It should be mentioned that, traditionally, one spectrum line per decade (M = 1) is extensively accepted for generating the Prony series representation. The higher density of the spectrum lines would generate higher accuracy for fitting but would produce more Prony series parameters. Thus, in the following sections, only the results for M = 1 are presented.
Figure 10 gives the master curves of the relaxation modulus and creep compliance for Mix-9.5 in the Prony series forms from the HN and 2S2P1D models. To verify the quality of the calculated Prony series parameters, the corresponding curves developed through the continuous spectra are also presented. Figure 11 gives the relative errors between the master curves from the Prony series forms and continuous spectra for Mix-9.5. It should be mentioned that since the spectrum lines were selected only at the time range covered by the test data, only the relative errors at 10−8 to 104 s were calculated. To achieve the infinite integrals in Equations (19) and (22), an integral interval of 10−40 to 10+40 s was employed to approximately represent the infinite one through the trapezoidal rule, in which 100 increments per decade were equidistantly selected on the logarithmic time scale. It can be seen that the curves from the Prony series parameters were in good agreement with those from the continuous spectra for both HN and 2S2P1D models over the region where the spectrum lines were selected, thus demonstrating the effectiveness of the proposed method in this study. Equally desirable results were also found for Mix-13.2.

4.3. Comparison to the Conventional Sigmoidal Model Method

Considering that the quality of the resulting Prony series parameters are dependent on the master curve models used for presmoothing, the results obtained were compared with those from the Sigmoidal model, which has been adopted by MEPDG [36]. The Sigmoidal model with four parameters can be expressed by:
lg E = a 1 + a 2 1 + e a 3 + a 4 lg ω
where a1 is the on the minimum logarithmic value of the dynamic modulus; a2 is the difference of the maximum and minimum logarithmic values of the dynamic modulus; a3 and a4 are model parameters governing the curve shape.
Unlike the HN and 2S2P1D models, the Sigmoidal function is a real-valued model for the dynamic modulus, and thus does not have an accurate analytical model for the corresponding phase angle. To deal with this issue, Rowe [37] developed a representation for the phase angle using an approximate Kronig–Kramers relation [38], as follows:
φ 90 × d lg E d lg ω = 90 a 2 a 4 e a 3 + a 4 lg ω 1 + e a 3 + a 4 lg ω 2
Table 3 shows the calculated Sigmoidal model parameters and fitting errors for the two mixtures. It can be observed that both HN and 2S2P1D models generated lower fitting errors than the Sigmoidal model, indicating their higher applicability to the complex modulus test data. To gain an in-depth insight into their advantages, the master curves of the dynamic modulus and phase angle for the three models were plotted in Figure 12 and Figure 13. It can be found that for the dynamic modulus, the curves from both HN and 2S2P1D models are non-centrosymmetric, that is, they offer asymmetric inflection points, whereas the Sigmoidal model is centrosymmetric on the log-log scale. As a result, the HN and 2S2P1D models exhibit higher flexibility than the Sigmoidal model in modeling the dynamic modulus of asphalt concrete. Additionally, the phase angle master curves of the HN and 2S2P1D models are non-axisymmetric, while that for the Sigmoidal is axisymmetric. Evidently, non-axisymmetric curves are more suitable for simulating the phase angle test data of asphalt concrete.
As a representation for the dynamic modulus, the Sigmoidal model does not have corresponding close-form solutions for H(ρ) or L(τ). Thus, the Prony series parameters for the relaxation modulus and creep compliance cannot be analytically yielded. To obtain the Prony series parameters, only the numerical approach can be used, in which the storage modulus representation from the generalized Maxwell model in Equation (37) is directly fitted to smoothed data produced from Equations (49) and (50). Similarly, the Prony series for the creep compliance also needs to be numerically computed. In this regard, the complex-valued models adopted in this study, like the HN and 2S2P1D models, have the prominent advantage over real-valued models.
Figure 14 displays the Black diagrams plotted using the Sigmoidal model for the two asphalt mixtures. For a comparison purpose, the generalized Maxwell model developed using the storage modulus data generated from the Sigmoidal model is also shown. To guarantee a good consistency of the generalized Maxwell model to the smoothed storage modulus data, the recursive fitting method developed by Sun et al. [14] was utilized. It can be clearly seen from Figure 14 that the curves from the two models diverge around the peaks of the phase angle, which indicates a noncompliance of the Sigmoidal model method with the LVE theory. This is ascribed to the use of approximate Kronig–Kramers relation. In this connection, the HN and 2S2P1D models employed in the presented approach can accurately satisfy the Kronig–Kramers relation due to the presence of the analytical representations of both the real and imaginary parts of the complex modulus.

5. Summary and Conclusions

This study presented a practical approach to fast acquiring high-quality Prony series parameters for both relaxation modulus and creep compliance of asphalt concrete based on the complex modulus test data. The approach can directly determine Prony series parameters through the analytical representations of the continuous relaxation and retardation spectra from the HN and 2S2P1D complex modulus models. With the model parameters determined in constructing dynamic modulus and phase angle master curves, the Prony series parameters can be immediately obtained with required accuracy. The elastic constants, Ee and Dg, for both LVE modulus and compliance functions can further be readily calculated through the smoothed data from the storage modulus representations of the HN and 2S2P1D complex modulus models. To offer an in-depth interpretation for the approach, the performance of the HN, 2S2P1D and conventional Sigmoidal models in fitting the complex modulus master curves were compared. Based on the results and analysis from this study, main conclusions can be drawn as follows:
(1)
The HN and 2S2P1D models yielded slightly different continuous spectral patterns at shorter relaxation times and longer retardation times. However, at the region covered by the test data, the continuous spectra of the two complex modulus models were very close to each other. Thus, the two models can generate comparable Prony series parameters within the time or frequency range covered by test data.
(2)
By means of the positive analytical expressions of the continuous spectra, local spectrum oscillations and undesirable negative spectrum strengths were successfully eliminated, thus generating high-quality Prony series parameters.
(3)
The HN and 2S2P1D models provide non-centrosymmetric curve patterns for the dynamic modulus master curves on the log-log scale and non-axisymmetric curve patterns for the phase angle master curves on the logarithmic angular frequency scale. Therefore, they performed better than the traditional Sigmoidal model in fitting to the complex modulus test data.
(4)
The Black diagram is recommended for examining the quality of the complex modulus test data before constructing the master curves, because it can effectively avoid the effect of testing temperatures.
(5)
The analytical expressions of the storage and loss moduli for both HN and 2S2P1D models accurately meet the Kronig–Kramers relation, and therefore the master curves constructed are consistent with the LVE theory.
(6)
All the procedures in the proposed method can be easily achieved even only by Microsoft Excel, successfully avoiding sophisticated expertise for programming in implementation process. Thus, the proposed method furnishes a practical way to fast acquiring high-quality Prony series parameters.
Further studies are required to develop predictive models of the complex modulus master curve of asphalt concrete based on the HN and 2S2P1D models by using statistical relationships between the model parameters and the constituent and volumetric properties of asphalt concrete, and the work is ongoing.

Author Contributions

Conceptualization, Y.Z. and Y.S.; methodology, Y.Z. and Y.S.; writing—original draft preparation, Y.Z.; writing—review and editing, Y.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Doctoral Start-up Foundation of Liaoning Province (2019-BS-048) and Fundamental Research Funds for the Central Universities [DUT19RC(4)023].

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing is not applicable to this article.

Acknowledgments

This study was sponsored by Doctoral Start-up Foundation of Liaoning Province (2019-BS-048) and Fundamental Research Funds for the Central Universities [DUT19RC(4)023]. The supports are gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hernandez-Fernandez, N.; Ossa-Lopez, A.; Harvey, J.T. Viscoelastic characterisation of asphalt concrete under different loading conditions. Int. J. Pavement Eng. 2021, 1–14. [Google Scholar] [CrossRef]
  2. Zhao, Y.; Liu, H.; Liu, W. Characterization of linear viscoelastic properties of asphalt concrete subjected to confining pressure. Mech. Time-Depend. Mater. 2012, 17, 449–463. [Google Scholar] [CrossRef]
  3. Park, S.W.; Kim, Y.R. Fitting Prony-series viscoelastic models with power-law presmoothing. J. Mater. Civ. Eng. 2001, 13, 26–32. [Google Scholar] [CrossRef]
  4. Nguyen, H.T.T.; Do, T.-T.; Tran, V.-T.; Phan, T.-N.; Pham, T.-A.; Nguyen, M.L. Determination of creep compliance of asphalt mixtures at intermediate and high temperature using creep-recovery test. Road Mater. Pavement Des. 2021, 22 (Suppl. 1), S514–S535. [Google Scholar] [CrossRef]
  5. Nguyen, H.T.T.; Nguyen, D.-L.; Tran, V.-T.; Nguyen, M.-L. Finite element implementation of Huet-Sayegh and 2S2P1D models for analysis of asphalt pavement structures in time domain. Road Mater. Pavement Des. 2020, 23, 22–46. [Google Scholar] [CrossRef]
  6. Li, L.; Wu, C.; Cheng, Y.; Ai, Y.; Li, H.; Tan, X. Comparative analysis of viscoelastic properties of open graded friction course under dynamic and static loads. Polymers 2021, 13, 1250. [Google Scholar] [CrossRef]
  7. Hristov, J.Y. Linear viscoelastic responses: The Prony decomposition naturally leads into the Caputo-Fabrizio fractional operator. Front. Phys. 2018, 6, 135. [Google Scholar] [CrossRef]
  8. Xu, Q.; Engquist, B.; Solaimanian, M.; Yan, K. A new nonlinear viscoelastic model and mathematical solution of solids for improving prediction accuracy. Sci. Rep. 2020, 10, 2202. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Gu, L.; Zhang, W.; Ma, T.; Qiu, X.; Xu, J. Numerical simulation of viscoelastic behavior of asphalt mixture using fractional constitutive model. J. Eng. Mech. 2021, 147, 04021027. [Google Scholar] [CrossRef]
  10. Londono, J.G.; Shen, R.; Waisman, H. Temperature-dependent viscoelastic model for asphalt–concrete implemented within a novel nonlocal damage framework. J. Eng. Mech. 2020, 146, 04019119. [Google Scholar] [CrossRef]
  11. Yu, D.; Yu, X.; Gu, Y. Establishment of linkages between empirical and mechanical models for asphalt mixtures through relaxation spectra determination. Constr. Build. Mater. 2020, 242, 118095. [Google Scholar] [CrossRef]
  12. Zhao, Y.; Ni, Y.; Zeng, W. A consistent approach for characterising asphalt concrete based on generalised Maxwell or Kelvin model. Road Mater. Pavement Des. 2014, 15, 674–690. [Google Scholar] [CrossRef]
  13. Park, S.W.; Schapery, R.A. Methods of interconversion between linear viscoelastic material functions. Part I—A numerical method based on Prony series. Int. J. Solids Struct. 1999, 36, 1653–1675. [Google Scholar] [CrossRef]
  14. Sun, Y.; Huang, B.; Chen, J. A unified procedure for rapidly determining asphalt concrete discrete relaxation and retardation spectra. Constr. Build. Mater. 2015, 93, 35–48. [Google Scholar] [CrossRef]
  15. Schapery, R.A. A Simple Collocation Method for Fitting Viscoelastic Models to Experimental Data. GALCIT SM 61-32A; California Institute of Technology: Pasadena, CA, USA, 1961. [Google Scholar]
  16. Cost, T.L.; Becker, E.B. A multidata method of approximate Laplace transform inversion. Int. J. Numer. Meth. Eng. 1970, 2, 207–219. [Google Scholar] [CrossRef]
  17. Tschoegl, N.W.; Emri, I. Generating line spectra from experimental responses. Part II: Storage and loss functions. Rheol. Acta 1993, 32, 322–327. [Google Scholar] [CrossRef]
  18. Xu, Q.; Solaimanian, M. Modelling linear viscoelastic properties of asphalt concrete by the Huet–Sayegh model. Int. J. Pavement Eng. 2009, 10, 401–422. [Google Scholar] [CrossRef]
  19. Olard, F.; di Benedetto, H. General “2S2P1D” model and relation between the linear viscoelastic behaviours of bituminous binders and mixes. Road Mater. Pavement Des. 2003, 4, 185–224. [Google Scholar]
  20. Tiouajni, S.; di Benedetto, H.; Sauzéat, C.; Pouget, S. Approximation of Linear Viscoelastic Model in the 3 Dimensional Case with Mechanical Analogues of Finite Size. Road Mater. Pavement Des. 2011, 12, 897–930. [Google Scholar]
  21. Tschoegl, N.W.; Emri, I. Generating line spectra from experimental responses. III. Interconversion between relaxation and retardation behavior. Int. J. Polym. Mater. 1992, 18, 117–127. [Google Scholar] [CrossRef]
  22. Schapery, R.A.; Park, S.W. Methods of interconversion between linear viscoelastic material functions. Part II—An approximate analytical method. Int. J. Solids Struct. 1999, 36, 1677–1699. [Google Scholar] [CrossRef]
  23. Levenberg, E. Smoothing asphalt concrete complex modulus test data. J. Mater. Civ. Eng. 2011, 23, 606–611. [Google Scholar] [CrossRef]
  24. Luo, R.; Lv, H.; Liu, H. Development of Prony series models based on continuous relaxation spectrums for relaxation moduli determined using creep tests. Constr. Build. Mater. 2018, 168, 758–770. [Google Scholar] [CrossRef]
  25. Lv, H.; Liu, H.; Tan, Y.; Sun, Z. Improved methodology for identifying Prony series coefficients based on continuous relaxation spectrum method. Mater. Struct. 2019, 52, 86. [Google Scholar] [CrossRef]
  26. Sun, Y.; Chen, J.; Huang, B. Characterization of asphalt concrete linear viscoelastic behavior utilizing Havriliak–Negami complex modulus model. Constr. Build. Mater. 2015, 99, 226–234. [Google Scholar] [CrossRef]
  27. Bhattacharjee, S.; Swamy, A.K.; Daniel, J.S. Continuous relaxation and retardation spectrum method for viscoelastic characterization of asphalt concrete. Mech. Time-Depend. Mater. 2011, 16, 287–305. [Google Scholar] [CrossRef]
  28. Zhang, F.; Wang, L.; Li, C.; Xing, Y. The discrete and continuous retardation and relaxation spectrum method for viscoelastic characterization of warm mix crumb rubber-modified asphalt mixtures. Materials 2020, 13, 3723. [Google Scholar] [CrossRef] [PubMed]
  29. AASHTO. Standard Method of Test for Determining Dynamic Modulus of Hot Mix Asphalt (HMA), AASHTO T 342-11; AASHTO: Washington, DC, USA, 2011. [Google Scholar]
  30. Havriliak, S.; Negami, S. A complex plane representation of dielectric and mechanical relaxation processes in some polymers. Polymer 1967, 8, 161–210. [Google Scholar] [CrossRef]
  31. Tschoegl, N.W. The Phenomenological Theory of Linear Viscoelastic Behavior: An Introduction; Springer: New York, NY, USA, 1989. [Google Scholar]
  32. Alavi, M.Z.; Hajj, E.Y.; Morian, N.E. Approach for quantifying the effect of binder oxidative aging on the viscoelastic properties of asphalt mixtures. Transp. Res. Rec. 2013, 2373, 109–120. [Google Scholar] [CrossRef]
  33. Sun, Y.; Huang, B.; Chen, J.; Jia, X.; Ding, Y. Characterizing rheological behavior of asphalt binder over a complete range of pavement service loading frequency and temperature. Constr. Build. Mater. 2016, 123, 661–672. [Google Scholar] [CrossRef]
  34. Airey, G.D. Use of black diagrams to identify inconsistencies in rheological data. Road Mater. Pavement Des. 2011, 3, 403–424. [Google Scholar] [CrossRef]
  35. Mun, S.; Geem, Z.W. Determination of viscoelastic and damage properties of hot mix asphalt concrete using a harmony search algorithm. Mech. Mater. 2009, 41, 339–353. [Google Scholar] [CrossRef]
  36. ARA (Applied Research Associates), ERES Consultants Division. Guide for Mechanistic-Empirical Design of New and Rehabilitated Pavement Structures, NCHRP 1-37A Final Report; ERES Consultants Division, Transportation Research Board: Washington, DC, USA, 2004. [Google Scholar]
  37. Rowe, G. Phase angle determination and interrelationships within bituminous materials. In Proceedings of the 7th International RILEM Symposium on Advanced Testing and Characterization of Bituminous Materials, Rhodes, Greece, 27–29 May 2009; pp. 43–52. [Google Scholar]
  38. Booij, H.C.; Thoone, G.P.J.M. Generalization of Kramers-Kronig transforms and some approximations of relations between viscoelastic quantities. Rheol. Acta 1982, 21, 15–24. [Google Scholar] [CrossRef]
Figure 1. Aggregate gradations of Mix-13.2 and Mix-9.5.
Figure 1. Aggregate gradations of Mix-13.2 and Mix-9.5.
Materials 15 00716 g001
Figure 2. Black diagram of the complex modulus test data set: (a) Mix-13.2; (b) Mix-9.5.
Figure 2. Black diagram of the complex modulus test data set: (a) Mix-13.2; (b) Mix-9.5.
Materials 15 00716 g002
Figure 3. Master curves of dynamic modulus and phase angle from the HN model: (a) Mix-13.2; (b) Mix-9.5.
Figure 3. Master curves of dynamic modulus and phase angle from the HN model: (a) Mix-13.2; (b) Mix-9.5.
Materials 15 00716 g003
Figure 4. Master curves of dynamic modulus and phase angle from the 2S2P1D model: (a) Mix-13.2; (b) Mix-9.5.
Figure 4. Master curves of dynamic modulus and phase angle from the 2S2P1D model: (a) Mix-13.2; (b) Mix-9.5.
Materials 15 00716 g004aMaterials 15 00716 g004b
Figure 5. Time-temperature shift factors obtained using the HN and 2S2P1D methods.
Figure 5. Time-temperature shift factors obtained using the HN and 2S2P1D methods.
Materials 15 00716 g005
Figure 6. Continuous spectra of the two asphalt mixtures: (a) Relaxation; (b) Retardation.
Figure 6. Continuous spectra of the two asphalt mixtures: (a) Relaxation; (b) Retardation.
Materials 15 00716 g006
Figure 7. Discrete spectra of Mix-13.2: (a) Relaxation; (b) Retardation
Figure 7. Discrete spectra of Mix-13.2: (a) Relaxation; (b) Retardation
Materials 15 00716 g007
Figure 8. Discrete spectra of Mix-9.5: (a) Relaxation; (b) Retardation.
Figure 8. Discrete spectra of Mix-9.5: (a) Relaxation; (b) Retardation.
Materials 15 00716 g008aMaterials 15 00716 g008b
Figure 9. Master curves of storage modulus for Mix-9.5 using the discrete relaxation spectrum with M = 1 from: (a) HN model; (b) 2S2P1D model.
Figure 9. Master curves of storage modulus for Mix-9.5 using the discrete relaxation spectrum with M = 1 from: (a) HN model; (b) 2S2P1D model.
Materials 15 00716 g009aMaterials 15 00716 g009b
Figure 10. Master curves of relaxation modulus and creep compliance for Mix-9.5 in the Prony series forms from: (a) HN model; (b) 2S2P1D model.
Figure 10. Master curves of relaxation modulus and creep compliance for Mix-9.5 in the Prony series forms from: (a) HN model; (b) 2S2P1D model.
Materials 15 00716 g010aMaterials 15 00716 g010b
Figure 11. Relative errors between the master curves from Prony series forms and continuous spectra for Mix-9.5.
Figure 11. Relative errors between the master curves from Prony series forms and continuous spectra for Mix-9.5.
Materials 15 00716 g011
Figure 12. Comparison of master curves of dynamic modulus and phase angle for Mix-13.2: (a) dynamic modulus; (b) phase angle.
Figure 12. Comparison of master curves of dynamic modulus and phase angle for Mix-13.2: (a) dynamic modulus; (b) phase angle.
Materials 15 00716 g012aMaterials 15 00716 g012b
Figure 13. Comparison of master curves of dynamic modulus and phase angle for Mix-9.5: (a) dynamic modulus; (b) phase angle.
Figure 13. Comparison of master curves of dynamic modulus and phase angle for Mix-9.5: (a) dynamic modulus; (b) phase angle.
Materials 15 00716 g013
Figure 14. Black diagrams from the Sigmoidal model for the two asphalt mixtures: (a) Mix-13.2; (b) Mix-9.5.
Figure 14. Black diagrams from the Sigmoidal model for the two asphalt mixtures: (a) Mix-13.2; (b) Mix-9.5.
Materials 15 00716 g014
Table 1. HN model parameters and fitting errors.
Table 1. HN model parameters and fitting errors.
Mix TypeTr/°C;Eg/MPaEe/MPaαβτ0/sF/%
Mix-13.22073,13292.00.3980.1930.0131.943
Mix-9.52443,70737.50.4310.1750.0112.461
Table 2. 2S2P1D model parameters and fitting errors.
Table 2. 2S2P1D model parameters and fitting errors.
Mix TypeTr/°C;Eg/MPaEe/MPaαkhβτ0/sF/%
Mix-13.22080,329126.21.8050.1040.41238,4002.485 × 10−41.961
Mix-9.52434,13249.12.3290.2050.49341,6661.558 × 10−32.204
Table 3. Sigmoidal model parameters and fitting errors.
Table 3. Sigmoidal model parameters and fitting errors.
Mix TypeTr/°C;a1a2a3a4F/%
Mix-13.2201.6473.089−0.246−0.4602.089
Mix-9.5241.1523.370−0.233−0.4592.471
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, Y.; Sun, Y. Fast-Acquiring High-Quality Prony Series Parameters of Asphalt Concrete through Viscoelastic Continuous Spectral Models. Materials 2022, 15, 716. https://doi.org/10.3390/ma15030716

AMA Style

Zhang Y, Sun Y. Fast-Acquiring High-Quality Prony Series Parameters of Asphalt Concrete through Viscoelastic Continuous Spectral Models. Materials. 2022; 15(3):716. https://doi.org/10.3390/ma15030716

Chicago/Turabian Style

Zhang, Yan, and Yiren Sun. 2022. "Fast-Acquiring High-Quality Prony Series Parameters of Asphalt Concrete through Viscoelastic Continuous Spectral Models" Materials 15, no. 3: 716. https://doi.org/10.3390/ma15030716

APA Style

Zhang, Y., & Sun, Y. (2022). Fast-Acquiring High-Quality Prony Series Parameters of Asphalt Concrete through Viscoelastic Continuous Spectral Models. Materials, 15(3), 716. https://doi.org/10.3390/ma15030716

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop