Comparative Study on the Dynamic Response of Asphalt Pavement Structures: Analysis Using the Classic Kelvin, Maxwell, and Three-Parameter Solid Models

Abstract

This study addressed the complex problems of selecting a constitutive model to objectively characterize asphalt mixtures and accurately determine their viscoelastic properties, which are influenced by numerous variables. Inaccuracies in model or parameter determination can result in significant discrepancies between the calculated and measured results of the pavement’s structural dynamic response. To address this, the research utilized the physical engineering principles of asphalt pavement structure to perform dynamic modulus tests on three types of high-content rubberized asphalt mixtures (HCRAM) within the surface layer. The research aimed to investigate the influencing factors of the dynamic modulus and establish a comprehensive master curve. This study also critically evaluated the capabilities of three viscoelastic models—the three-parameter solid model, the classical Maxwell model, and the classical Kelvin model—in depicting the dynamic modulus of HCRAM. The findings indicated a negative correlation between the dynamic modulus of the asphalt mixture and temperature, while a positive association exists between the loading frequency and temperature, with the impact of the loading frequency diminishing as the temperature increases. Notably, the three-parameter solid model was identified as the most accurate in describing the viscoelastic properties of the HCRAM. Furthermore, the dynamic response calculations revealed that most indexes in the surface layer’s dynamic response are highest when evaluated using the three-parameter viscoelastic model, underscoring its potential to enhance the pavement performance’s predictive accuracy. This research provides valuable insights into optimizing the material performance and guiding the pavement design and maintenance strategies.

1. Introduction

One of the most often utilized road materials worldwide is asphalt pavement, with higher durability, comfort, and other advantages [1]. With the acceleration of urbanization and transportation demand growth, asphalt pavement application is still promising. In the future, with the advancement of technology and the application of new materials, many researchers will devote themselves to the study of asphalt pavements in terms of materials, design, construction, and maintenance to meet the ever-increasing transportation demand [2,3,4,5].

In terms of materials, in order to enhance the pavement performance and make the materials more ecologically and energy-efficient, and ultimately to achieve the goal of a “carbon neutral” vision [6], the application of rubberized asphalt (RA) and its mixtures have been a prominent topic of study, with numerous researchers conducting in-depth investigations on it [7,8,9]. Rubberized asphalt is widely used in stress-absorbing layers and rutting-resistant sections of asphalt pavements due to its exemplary environmental performance, rutting resistance, and durability, thus prolonging the pavement’s service life [10]. For rubberized asphalt, a rubber powder content of over 25% can become a high-content rubberized asphalt, and a large amount of research has proven the feasibility of high-content rubberized asphalt [11,12]. In recent years, it has been found that high-content rubberized asphalt (HCRA) has better performance in terms of high-temperature rutting resistance and fracture resistance compared to ordinary RA [13,14]. Also, its asphalt mixtures have good rutting and fatigue resistance [15,16,17,18]. However, few domestic and international studies have been conducted on the viscoelastic characteristics of high-content rubberized asphalt mixtures (HCRAM). Therefore, this study selected HCRAM as the research object, and their viscoelastic characteristics were examined using the dynamic modulus (DM) test.

Since asphalt mixture is a common viscoelastic material, it will display various mechanical properties at different temperatures and frequencies [19,20]. In order to investigate the viscoelastic properties, many scholars have carried out DM studies on asphalt mixtures [21,22]. Currently, the research on the DM of asphalt mixtures in various countries mainly includes test methods, influencing factors, and prediction of the DM. The DM of asphalt mixtures can be measured utilizing a variety of domestic and international test techniques [23,24]. The commonly used methods in China are the DM test method in “Standard Test Methods of Bitumen and Bituminous Mixtures in Highway Engineering” (JTG E20-2011) and the Superpave Simple Performance Tester test method (SPT) [25,26]. In addition to indoor testing, on-site estimation of the dynamic modulus is also crucial. Andreas Loizos et al. [27] conducted on-site tests on the experimental road section and non-destructive testing using a drop hammer deflectometer and ground penetrating radar to study the viscoelastic response of the road surface. In addition, regarding the investigation of factors affecting the DM, Wang, H et al. [28] investigated the impacts of migration, temperature, and loading frequency on the DM and phase angle of asphalt mixtures. Zhou et al. [29] found that the void ratio has the most significant combined effect on the DM and phase angle.

Regarding the analysis of DM forecasting, Witzack [30] developed a DM prediction model for densely mixed hot mix asphalt mixtures. Using the Witzack model, Wei et al. [31] fitted the Sigmoidal model to the DM test findings. They found that the Sigmoidal model could better respond to the viscoelastic mechanical characteristics of asphalt mixtures. Dana Daneshvar and Ali Behnood [32] created a DM prediction model utilizing the random forest algorithm and confirmed its superiority. In summary, to study the viscoelastic characteristics of HCRAM, this work measured the DM of asphalt mixtures by indoor DM tests and plotted the master curve utilizing the Sigmoidal to accurately predict the DM.

Stress state analysis and structural dynamic response calculation in asphalt pavement are essential for the design, construction, and maintenance [33]. Choosing a constitutive model that can adequately characterize the material and precisely estimate its characteristics is essential to obtaining a real and accurate stress state and structural dynamic response [34]. Researchers have proposed a number of viscoelastic constitutive models to explain the viscoelastic properties of asphalt mixes [19], such as the classical Maxwell and its generalized model, the classical Kelvin and its generalized model [35], and the three-parameter solid model [36]. However, the majority of nations still design asphalt pavements using the layered elastic system mechanics theory under vertical stresses due to computational simplicity [37]. This method needs to consider the influence of horizontal forces generated by driving loads on the pavement and accurately describes the viscoelastic behavior of road materials under repeated loads, which can easily lead to consistency between the calculated and measured results of the dynamic response of the road structure. Consequently, in this study, the viscoelastic constitutive model is employed to determine the internal mechanical response of the pavement structure to objectively describe the structural and mechanical response features of pavement.

Exploring the response of pavement structures is extremely necessary for pavement design or performance evaluation. However, this study cannot solely focus on pavement structure or material characteristics. Instead, the two should be combined to study the performance of asphalt pavement under their combined action. A helpful technique for bridging the gap between experimental and pavement research is numerical simulation [38,39]. With the application of numerical simulation methods, more and more scholars are using finite element software to study mechanical responses. Ji et al. [40] used ANSYS to simulate the pavement structure under different driving loads. They studied the connection between the tensile stress, shear stress, and pavement’s asphalt layer’s thickness. Z Dong et al. [41] established a large axle load model using Abaqus-6.14 software. They proposed that for semi-rigid mobile heavy-duty asphalt pavements, the vertical strain at the top of the roadbed and the vertical stress of the pavement must be considered in the design. In addition, the viscoelastic constitutive model of asphalt mixtures can be implemented into Abaqus-6.14 software through user-defined subroutines (UMAT) to evaluate the response of asphalt structures [42]. This study chooses to use UMAT to compute the constitutive model of HCRAM. It conducts dynamic response analysis of pavement under viscoelastic behavior.

2. Study Objectives and Scope

Based on the above background, this study aimed to characterize the dynamic response of HCRA pavement structures realistically and objectively. This study relied on the physical engineering of the pavement structure of the Jingde Expressway, considering its surface materials’ viscoelastic properties. Dynamic modulus tests and analysis of factors affecting the DM were conducted on three types of HCRAM, and a dynamic modulus master curve (DMMC) was established utilizing the Sigmoidal and studied the viscoelastic properties of high content (30%, 40%, and 50%) rubber asphalt mixtures. Further, the variability in the description of the DM of the HCRAM by the classical Maxwell model, the classical Kelvin model, and the three-parameter solid model was evaluated according to the theory of viscoelasticity of the HCRAM concerning the established master curves. Finally, the UMAT subroutine was invoked to implement the viscoelastic intrinsic model numerically in Abaqus. At the same time, the pavement structure’s dynamic response was calculated, and the dynamic response and differences of pavement structures based on three viscoelastic models were studied, objectively describing the structural and mechanical response properties of the pavement and providing the accurate basis for pavement structure design.

3. Materials and Methods

3.1. Raw Material and Mix Design

3.1.1. Asphalt

The HCRA employed in this study had 30%, 40%, and 50% rubber doping. The asphalt was provided by Hebei Provincial Communications Planning, Design and Research Institute Co., Ltd. in Shijiazhuang, China, which is consistent with the asphalt materials used in China’s Jingde Highway’s upper, middle, and lower layers, respectively. The PG grading performance is presented in

Table 1. RA PG grading properties.

Test Items	30% RA	40% RA	50% RA
PG classification	88-34	88-34	82-28
Rutting factor (88 °C)	1683.8	1421.6	1638
Rutting factor after aging (88 °C)	3429.5	2767.2	1672
Fatigue limit temperature (°C)	13	10	16
Creep strength (−24 °C)	201.9205	166.9801	343
Creep rate	0.3054	0.3147	0.335

, and the related indexes are displayed in

Table 2. Performance test results of RA.

Test Items		30% RA	40% RA	50% RA
Penetration 25 °C, 100 g, 5 s (0.1 mm)		62	66	58
Ductility 5 cm/min, 5 °C (cm)		16	18	12.8
Softening point (°C)		72.5	72	78.7
Dynamic viscosity 180 °C (Pa·s)		2.8	3.5	2.96
Elastic recovery 25 °C (%)		84	81	71
48 h softening point difference (°C)		2.0	0.5	6.4
0.1 Kpa recovery rate (%)		98.67	99.43	92.08
3.2 Kpa recovery rate (%)		90.79	92.55	86.45
Residue after TFOT (163 °C, 85 min)	−0.29	−0.3	0.3	T 0610-2011
	77	76	86.2	T 0604-2011
	11	12	10.4	T 0605-2011

[43].

3.1.2. Aggregates

The aggregate used in this study is consistent with the research of Wang et al. [44], which is limestone, and the filler is limestone mineral powder. Various characteristics indicators of aggregates were examined, and aggregates’ physical and mechanical performance indicators were finally obtained. Relevant literature can be found in [44].

3.1.3. Mix Design

The HCRAM gradation design method used in this paper mainly refers to the “Technical Specifications for Construction of Highway Asphalt Pavements” (JTG F40-2004), the ARHM-13 (The maximum nominal particle size of the asphalt rubber hot mixture is 13), ARHM-20, and ARHM-25 gradation curve presented in Figure 1. The asphalt aggregate ratio of the ARHM-13, ARHM-20, and ARHM-25 corresponding to the finalized 30%, 40%, and 50% blends of RA is 5.22%, 4.79%, and 4.36%, respectively.

3.2. Experiment Approach

3.2.1. Preparation of Test Specimens

After being mixed according to the desired gradation, the limestone aggregate was dried in an oven maintained at 110 °C until it attained a constant weight. For more than two hours, the asphalt and aggregate were placed in an oven set at 180 °C, and an indoor mixer was used to mix the HCRAM under the optimal oil-stone ratio conditions. Following demolding, a cylindrical sample measuring 150 mm ± 2 mm in height and 100 mm ± 2 mm in diameter was obtained. The samples were made using the rotary compaction method. The resulting HCRAM specimen is displayed in Figure 2.

3.2.2. Dynamic Modulus Test Methods

This study used a multifunctional material testing machine, MTS Landmark 370.25, for uniaxial compression the DM testing. The schematic diagram of the testing instrument and the uniaxial compression DM test are shown in Figure 3 and Figure 4, respectively. The selected temperatures were 15 °C, 30 °C, 45 °C, and 60 °C.

Table 3. RA PG grading properties.

Frequency	Repetitions	Frequency	Repetitions
25	200	1	20
10	200	0.5	15
5	100	0.1	15

shows the number of repeated loads and frequencies. The test was conducted according to the uniaxial compression DM test method in the “Standard Test Methods of Bitumen and Bituminous Mixtures in Highway Engineering” (JTG E20-2011).

3.3. Numerical Simulation Methods

3.3.1. Establishment of Finite Element Model

The Jingde Expressway test road in China’s pavement structure was replicated in the pavement structure model. The UMAT subroutine was used to set the parameters for the viscoelastic surface layer material. Elastic materials included the base course and soil foundation, and the material characteristics were directly determined by the internal software program. The pavement materials and structural design parameters are displayed in

Table 4. Pavement structure combination and material parameters.

Structure Layer No	Horizon	Material Type	Model	Density (kg/m3)	Compressive Modulus (MPa)	Tensile Modulus (MPa)	Pressure Poisson’s Ratio
1	Upper layer	ARHM-13	Viscoelastic model	0.04	2400	-	0.25
2	Middle layer	ARHM-20	Viscoelastic model	0.08	2400	-	0.25
3	Lower layer	ARHM-25	Viscoelastic model	0.1	2400	-	0.25
4	Upper base course	Inorganic binding material stabilizer	Linear elastic model	0.18	2300	10,500	0.25
5	Medium base course	Inorganic binding material stabilizer	Linear elastic model	0.18	2300	10,500	0.25
6	Lower base	Inorganic binding material stabilizer	Linear elastic model	0.18	2300	8500	0.25
7		Soil foundation	Linear elastic model	-	1800	280	0.40

. The model dimensions were 9 m in length, 2.75 m in width, and 5 m in depth. The construction basis of the model was as follows: due to the symmetrical pavement structure, considering computer performance and computational efficiency, only half of the pavement model was established in ABAQUS-6.14 software. A load band was drawn in advance to affect the exact condition of the pavement, with a length of 6 m and a width of 0.23 m. The load band is presented in Figure 5.

Considering computer performance, the grid division method in the x-direction was as follows: the grid width of the moving load band was set to 2.3 cm. The load was carried to the pavement’s center, and the grid width was set to 0.2 cm. The remaining grids were divided according to the principle of uniform distribution. The grid division method in the y-direction ensured that the grid shape approximates a square with a fixed grid length of 2.3 cm. The grid division method in the z-direction was as follows: near the upper part of the model, for example by the asphalt surface layer, the grid was encrypted, and the grid height was set to 1 cm. The height of the grid at the grassroots level was set to 3 cm. The grid of the soil base layer gradually became sparse according to the principle of equal distribution.

The boundary conditions of the finite element model directly affected the convergence of the model. The boundary conditions of the finite element model in this study were as follows: on the X = 0 m and X = 6 m planes, the displacement in the x-direction was limited, that is, U₁ = 0 on this vertical plane; on the Y = 0 m plane, the displacement in the y-direction was limited, i.e., U₂ = 0 on this vertical plane; on a plane with Y = 2.75 m, the displacement was set as a symmetrical boundary; at the bottom of the model, on a plane of Z = 5 m, the displacement was set as an infinite mesh.

Due to the tire width of the loading device being 0.23 m, the moving load was set to a rectangle with a length of 0.23 m and a width of 0.23 m. The method of applying moving loads was referred to the literature [45].

3.3.2. Numerical Implementation of Viscoelastic Constitutive Model

ABAQUS provides standard finite element analysis programs and has good openness, which can generate non-standard analysis programs by providing user subroutine interfaces to meet user needs. Through the user material subroutine, users can define any supplementary material model. Not only can any material constants be read as data, but ABAQUS also provides storage functionality for any state variables related to the solution at each material count point for application in these subroutines.

This work used UMAT to numerically apply the viscoelastic constitutive model of HCRAM, defining the asphalt pavement surface material as a viscoelastic material. The UMAT subroutine can be found in Appendix A (

Table A1. The UMAT subroutine and the input parameters (the “*” symbol is equivalent to “×”).

SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD,

1 RPL, DDSDDT, DRPLDE, DRPLDT,

2 STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED, CMNAME,

3 NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT,

4 CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, KSTEP, KINC)

C

INCLUDE ‘ABA_PARAM.INC’

C

CHARACTER*80 CMNAME

DIMENSION STRESS(NTENS), STATEV(NSTATV),

1 DDSDDE(NTENS, NTENS),

2 DDSDDT(NTENS), DRPLDE(NTENS),

3 STRAN(NTENS), DSTRAN(NTENS), TIME(2), PREDEF(1), DPRED(1),

4 PROPS(NPROPS), COORDS(3), DROT(3,3), DFGRD0(3,3), DFGRD1(3,3)

DIMENSION DSTRES(6), D(3,3)

C

C EVALUATE NEW STRESS TENSOR

C

EV = 0

DEV = 0

DO K1 = 1, NDI

EV = EV + STRAN(K1)

DEV = DEV + DSTRAN(K1)

END DO

C

TERM1 = PROPS(2)*(0.5*DTIME*PROPS(1) + PROPS(3))/(1 + PROPS(4))/

1 (0.5*DTIME*(PROPS(1) + PROPS(2)) + PROPS(3))

TERM2 = 3*PROPS(4)*PROPS(2)*(0.5*DTIME*PROPS(1) + PROPS(3))/

1 (1 + PROPS(4))/(1−2*PROPS(4))/(0.5*DTIME*(PROPS(1) + PROPS(2))

2 + PROPS(3))

TERM3 = DTIME*PROPS(1)*PROPS(2)/(1 + PROPS(4))/(0.5*DTIME*

1 (PROPS(1) + PROPS(2)) + PROPS(3))

TERM4 = 3*DTIME*PROPS(4)*PROPS(1)*PROPS(2)/(1 + PROPS(4))/

1 (1−2*PROPS(4))/(0.5*DTIME*(PROPS(1) + PROPS(2)) + PROPS(3))

TERM5 = DTIME*(PROPS(1) + PROPS(2))/(0.5*DTIME*(PROPS(1)+

1 PROPS(2)) + PROPS(3))

C

C

C

DO K1 = 1, NDI

DSTRES(K1) = TERM1*DSTRAN(K1) + TERM2*DEV/3 + TERM3*

1     STRAN(K1) + TERM4*EV/3−TERM5*STRESS(K1)

STRESS(K1) = STRESS(K1) + DSTRES(K1)

END DO

C

C

C

I1 = NDI

DO K1 = 1, NSHR

I1 = I1 + 1

DSTRES(I1) = TERM1/2*DSTRAN(I1) + TERM3/2*STRAN(I1)−

1     TERM5*STRESS(I1)

STRESS(I1) = STRESS(I1) + DSTRES(I1)

END DO

C

C

C

DO K1 = 1, NTENS

DO K2 = 1, NTENS

DDSDDE(K2,K1) = 0

END DO

END DO

C

DO K1 = 1, NDI

DDSDDE(K1,K1) = TERM1 + TERM2/3

END DO

C

DO K1 = 2, NDI

N2 = K1−1

DO K2 = 1, N2

DDSDDE(K2,K1) = TERM2/3

DDSDDE(K1,K2) = TERM2/3

END DO

END DO

C

I1 = NDI

DO K1 = 1, NSHR

I1 = I1 + 1

DDSDDE(I1,I1) = TERM1/2

END DO

RETURN

END

).

3.3.3. Validation Method for Numerical Simulation Effectiveness

The methods for verifying the effectiveness of numerical simulations can be divided into two types: direct verification and indirect verification. Direct verification refers to comparing the measured data from laboratories or actual engineering with the model calculated by numerical simulations. Indirect verification involves comparing the results of theoretical calculations or analytical solutions with numerical simulation results. This study adopts the second method, which uses the calculation results of Bisar3 software to verify the effectiveness of the numerical simulation method.

Because Bisar3 software is based on an elastic layered system calculation method, it is necessary to degrade the constitutive model in Abaqus-6.14 software from a viscoelastic model to an elastic model and then calculate the same road surface model to compare the calculation results of the two.

4. Results and Discussion

4.1. Analysis of DM Test Results for Asphalt Mixtures

4.1.1. Analysis of Factors Affecting the DM

(1)

The impact of loading frequency

Figure 6 displays the trend of the DM of the HCRAM with frequency. From Figure 6, the DM of the HCRAM was observed to increase with frequency and to have a tendency of fast growth in the early stage and slow growth in the later stage. Especially when the loading rate increased from 0.1 Hz to 1 Hz, the modulus of the HCRAM increased roughly linearly. However, the growth rate was prolonged when the frequency increased from 10 Hz to 25 Hz. It can be inferred that the HCRAM’s DM eventually reached a certain limit value and fails to grow indefinitely with an increase in loading rate.

The aforementioned phenomenon can be elucidated from the standpoint of the viscoelastic properties of asphalt mixtures [46], mainly because when the loading frequency is low, the HCRAM mainly exhibit viscous characteristics. When the frequency is high, the asphalt mixture mainly exhibits elastic properties, which means that when subjected to high-frequency loads, the asphalt mixture is no longer an apparent viscoelastic material. At this point, the impact of the loading rate on the DM is no longer as significant as that at low frequencies. The explanation for the hysteresis of asphalt mixture deformation is that when the asphalt mixture is subjected to external loads [47], the stress–strain changes generated inside the specimen will exhibit a certain degree of hysteresis. If the frequency is reduced, the viscosity characteristics of the asphalt mixture are apparent, and the subsequent phenomenon will be more pronounced; that is, the deformation is more minor, and the DM is larger. Conversely, if the loading frequency is reduced, the deformation is more significant, and the DM is smaller. To some extent, this explains that asphalt is more prone to deformation and rutting under low-speed traffic loads, and the pavement is more prone to cracking under high-speed traffic loads [48].

Figure 7 shows the variation trend of the phase angle of each asphalt mixture as loading frequency. This indicates that the phase angle variation trend with frequency differs with temperature. As the loading frequency increases at test temperatures of 15 °C and 30 °C, the phase angle progressively decreases. This trend can indicate that when the temperature is low, the viscous properties of the asphalt mixture decrease, and the elastic properties increase with the rise of loading frequency. At both 45 °C and 60 °C test temperatures, the phase angle of the asphalt mixture rises proportionally with the increase in frequency. This indicates that at elevated temperatures, the viscosity of the asphalt mixture intensifies while its elasticity diminishes with the rise in loading frequency.

Figure 8 shows the variation trend of the storage modulus of the HCRAM as frequency. According to Figure 8a, it is evident that the trend of the storage modulus changing with the loading frequency is similar to that of the DM, with a faster rate of change in the early stage and a slower trend in the later stage. The storage modulus characterizes the elasticity of the HCRAM, which can store and release energy when subjected to dynamic loads. At 25 Hz, the maximum DM is measured. This implies that the energy contained in the asphalt mixture is unable to be completely dissipated when subjected to high frequencies and low temperatures. The elastic properties exhibited by the asphalt mixture are apparent, allowing the stored modulus to reach its maximum value under these conditions. On the contrary, under high temperature and low-frequency conditions, the elastic performance of the material is weak, resulting in a decrease in the DM.

Figure 9 shows the variation trend of the loss modulus as frequency. This demonstrates that the trends of the DM, storage modulus, and loss modulus with frequency are all the same. The loss modulus characterizes the HCRAM’s viscous properties. Under certain temperature conditions and high loading frequencies, the viscosity and elasticity of asphalt mixtures increase, increasing the energy loss. This is manifested as a rise in the loss modulus of asphalt mixtures under dynamic loads.

(2)

The influence of temperature

Figure 10 shows the variation trend of the DM and other experimental data of the ARHM-13 with the temperature. The DM trend of ARHM-13 with a temperature fluctuation is presented in Figure 10a. In Figure 10a, it can be observed that, under various frequencies, as the temperature rises, the HCRAM’s DM drastically drops. For instance, the DM of the ARHM-13 is 66%, 89%, and 94% lower at 30 °C, 45 °C, and 60 °C when the frequency is 25 Hz compared to when it is at 15 °C. As the temperature rises, the DM of the asphalt mixture decreases, and its rutting resistance decreases. This phenomenon mostly occurs because, at high temperatures, the asphalt in the asphalt mixture mainly exhibits viscoelastic properties. As the temperature rises, the viscosity of the asphalt reduces, and its rutting resistance performance decreases. The pavement is prone to permanent shear deformation and rutting when subjected to external forces.

Figure 10b shows the trend of the phase angle with the temperature. As the figure illustrates, the variation of the phase angle with the loading frequency is more complicated. The phase angle change mainly reflects the viscoelasticity of the polymer. At low temperatures, by increasing the frequency, the phase angle change first rises and then falls, while at high temperatures, it shows an upward trend when the frequency is increased. When the loading frequency is 0.1 Hz, 0.5 Hz, and 1 Hz, the phase angle peaks at 30 °C; when the loading frequency is 5 Hz, 10 Hz, and 25 Hz, the phase angle reaches their peak at 45 °C. The primary cause for this changing pattern is that, under low temperatures and high-frequency loads, the asphalt binder exhibits significant elastic properties, with solid performance and increased viscosity with rising temperatures. Since the characteristics of the asphalt binder have a greater influence on the asphalt mixture’s performance, the phase angle rises as the temperature rises. In theory, as a composite material, the phase angle of the asphalt mixture is equal to the phase angle of the matrix material. The asphalt mixture’s phase angle is unaffected by aggregates. However, under high temperature and low-frequency action, the asphalt binder becomes soft, and the embedding force between the mineral skeleton increases, which exceeds the influence of the asphalt bonding force. As minerals are elastic materials, their phase angle is zero, so the phase angle will decrease and further stiffen [49].

Figure 10c,d demonstrate the temperature-dependent trends of the storage and loss modulus of ARHM-13. It is not difficult to find from the figures that their trends are similar to those of the DM, decreasing with increasing temperatures. This phenomenon is because the storage and loss modulus are related to the molecular chain movement of asphalt. The frictional resistance between the asphalt molecules is greatly affected by temperature, and temperature changes cause changes in asphalt performance.

4.1.2. Establishment of the Dynamic Modulus Master Curve

The DM was fitted and analyzed in this work using a Sigmoidal model [50], whose formula is as follows:

$$\lg \left(\left|E^{\ast }\right|\right)=\delta +\frac{\alpha }{1+e^{\beta +\gamma ·\lg \left({\omega }_r\right)}}$$

(1)

where $\left|E^{\ast }\right|$ is the DM; $\delta$ is the logarithm of the minimum dynamic modulus; $\beta$, $\gamma$ are the parameter that describes the shape of the sigmoidal function; $\delta +\alpha$ is the logarithm of the maximum dynamic modulus; ${\omega }_r$ is frequency reduction.

The DM obtained at different temperatures can be converted from the loading frequency $\omega$, which is the logarithm of the minimum dynamic modulus temperatures to the corresponding reduced frequency ${\omega }_r$ at the reference temperature through a shift factor $\lg {\alpha }_T$. As a result, the DM value at the converted frequency can finally be obtained. The relationship between $\lg {\alpha }_T$, $\omega$, and ${\omega }_r$ is shown as follows.

$${\omega }_r={\alpha }_T\cdot \omega$$

(2)

$$\lg {\omega }_r=\lg {\alpha }_T+\lg \omega$$

(3)

The DMMC is established using a sigmoidal model based on the DM test results at different temperatures and frequencies, with $\left|E^{\ast }\right|$ as the vertical axis and ${\omega }_r$ as the horizontal axis. For ease of comparison, 30 °C was selected as the reference temperature. The DM values at the other three temperatures were translated using Origin-2018 software to obtain the shift factor, and then the fitted sigmoidal function was used to obtain the DMMC. Figure 11 presented the fitting results of the DMMC for each HCRAM, as presented in

Table 5. DMMC of various HCRAM sigmoidal model.

Types	Temperature (°C)	Shift Factor	Dynamic Modulus Main Curve Model	R2
ARHM-13	15	2.62	$\lg \left(\left|E^{*}\right|\right)=1.9298+\frac{2.9669}{1+e^{0.2725-0.2551\lg \left({\omega }_{\mathrm{r}}\right)}}$	0.9990
	30	0
	45	−2.67
	60	−4.37
ARHM-20	15	2.81	$\lg \left(\left|E^{*}\right|\right)=2.0630+\frac{2.8337}{1+e^{0.1589-0.2449\lg \left({\omega }_{\mathrm{r}}\right)}}$	0.9943
	30	0
	45	−1.85
	60	−4.31
ARHM-25	15	2.88	$\lg \left(\left|E^{*}\right|\right)=2.2018+\frac{2.6949}{1+e^{0.1321-0.2701\lg \left({\omega }_{\mathrm{r}}\right)}}$	0.9971
	30	0
	45	−1.48
	60	−3.93

. The R² values are all higher than 0.99, indicating a very good fitting effect. That is, the model is highly consistent with the experimental results.

As demonstrated by Figure 11, the DMMC slope is comparatively small at both lower and higher frequencies., indicating that the DM tends to stabilize at this time. The change trend is not apparent and will tend to have a limited value. The shift factor is a function of temperature, and its amplitude can indicate the degree of temperature dependence of asphalt mixtures. Table 5 demonstrates that, starting from the gradation of the mixture and the rubber content, in terms of temperature dependence, ARHM-13 has 30% rubber content, <ARHM-20 has 40% rubber content, and <ARHM-25 has 50% rubber content. This result indicates that, as the rubber content and gradation increase, the temperature dependence of asphalt mixtures also increases.

Figure 12 shows the DMMC of three types of HCRAM at different temperatures. It depicted that there is no substantial disparity in the DMMC of the three types of HCRAM, especially in the high-frequency and low-frequency regions, where the three master curves tend to overlap. The rate at which the DM changes from low frequency to high frequency steadily slows down as the temperature rises, suggesting that the asphalt mixture is less impacted by frequency as the temperature rises. At constant temperature, the master curve’s position rises with increasing asphalt mixture gradation, suggesting that the ARHM-25 asphalt mixture has a higher DM.

4.2. Establishment of Intrinsic Relationships and Parameter Fitting for Viscoelastic Models

4.2.1. Establishment of Constitutive Relations for Viscoelastic Models

(1)

Classic Kelvin model

The Kelvin model comprises a spring and a sticky pot in parallel. The model structure diagram is displayed in Figure 13. In this connection method, the strain of the two components is equal, and the total stress is equal to the sum of the stresses of the spring and the adhesive pot. The constitutive equation is represented by Equation (4).

$$\sigma ={\sigma }_1+{\sigma }_2=E\epsilon +\eta {\epsilon }^{\prime }$$

(4)

By applying Fourier transform to the Kelvin constitutive model, Equation (5) can be obtained:

$$\sigma \left(\omega \right)=E\epsilon \left(\omega \right)+\eta \left(i\omega \right)\epsilon \left(\omega \right)$$

(5)

where: $i$ is an imaginary singular number; $\omega$ is the frequency.

The complex modulus of the Kelvin model can be obtained by performing a Fourier transform on Equation (5) and then using the definition of complex modulus, as shown in Equation (6):

$$Y^{\ast }\left(i\omega \right)=\frac{\sigma \left(\omega \right)}{\epsilon \left(\omega \right)}=E+\eta \left(i\omega \right)=Y_1\left(\omega \right)+i{Y}_2\left(\omega \right)$$

(6)

From this, the expression for the dynamic modulus can be obtained as shown in Equation (7):

$$\left|Y^{\ast }\left(i\omega \right)\right|=\sqrt{E^2+{\eta }^2{\omega }^2}$$

(7)

(2)

Classic Maxwell model

The Maxwell model comprises a spring and a sticky pot connected in series. The model structure is shown in Figure 14. Assuming that under stress, the strains generated by the spring and the adhesive pot are ${\epsilon }_1$ and ${\epsilon }_2$, respectively, the total deformation of the component is the sum of their deformations. The constitutive equation is represented by Equation (8).

$$\epsilon ={\epsilon }^{\prime }{}_1+{\epsilon }^{\prime }{}_2=\frac{\dot{\sigma }}{E}+\frac{\sigma }{\eta }$$

(8)

Using the same method to process the Classic Maxwell model yields the expression for dynamic modules as shown in Equation (9):

$$\left|Y^{\ast }\left(i\omega \right)\right|=\sqrt{\frac{{\left(E{\eta }^2{\omega }^2\right)}^2+{\left(E^2\eta \omega \right)}^2}{{\left(E^2+{\eta }^2{\omega }^2\right)}^2}}$$

(9)

(3)

Three-parameter solid model

A Kelvin model coupled in series with a spring element makes up the three-parameter solid model. The model structure is presented in Figure 15. The constitutive equation is represented by Equation (10).

$$\epsilon ={\epsilon }_1+{\epsilon }_2$$

(10)

For the Kelvin model:

$$\frac{\sigma }{E_1}={\epsilon }_1+\frac{{\eta }_1}{E_1}{\dot{\epsilon }}_1$$

(11)

For elastic components:

$$\frac{\sigma }{E_2}={\epsilon }_2$$

(12)

After sorting out the above equations and taking their derivatives, Equation (13) can be obtained:

$$\sigma +\frac{{\eta }_1}{E_1+E_2}\dot{\sigma }=\frac{E_1{E}_2}{E_1+E_2}\epsilon +\frac{{\eta }_1{E}_2}{E_1+E_2}\dot{\epsilon }$$

(13)

Using the same method to process the three-parameter solid model yields the expression for dynamic modules as shown in Equation (14):

$$\left|Y^{\ast }\left(i\omega \right)\right|=\sqrt{\frac{{\left(\frac{E_1{E}_2}{E_1+E_2}+E_2{\left(\frac{{\eta }_1}{E_1+E_2}\right)}^2{\omega }^2\right)}^2+{\left(\left(E_2-\frac{E_1{E}_2}{E_1+E_2}\right)\cdot \frac{{\eta }_1\omega }{E_1+E_2}\right)}^2}{{\left(1+{\left(\frac{{\eta }_1}{E_1+E_2}\right)}^2{\omega }^2\right)}^2}}$$

(14)

(4)

Jacobian Matrix and UMAT

For linear elastic models, there is a relationship between the elements in the stiffness matrix and the elastic modulus as follows:

$$\left\{\begin{array}{l}{c}_{11}=c_{22}=c_{33}=\frac{1-\mu }{\left(1-2\mu \right)\left(1+\mu \right)}E\\ c_{12}=c_{13}=c_{23}=\frac{\mu }{\left(1-2\mu \right)\left(1+\mu \right)}E\\ c_{44}=c_{55}=c_{66}=\frac{1}{2\left(1+\mu \right)}E\end{array}\right.$$

(15)

Assuming that the Poisson’s ratio of the classical viscoelastic model is a constant value, the relationship between the stiffness matrix after the Laplace transform and the corresponding values of the elastic modulus is similar to that of the linear elastic model, with:

$$\left\{\begin{array}{l}{\overline{c}}_{11}\left(s\right)={\overline{c}}_{22}\left(s\right)={\overline{c}}_{33}\left(s\right)=\frac{1-\mu }{\left(1-2\mu \right)\left(1+\mu \right)}\overline{E}\left(s\right)\\ {\overline{c}}_{12}\left(s\right)={\overline{c}}_{13}\left(s\right)={\overline{c}}_{23}\left(s\right)=\frac{\mu }{\left(1-2\mu \right)\left(1+\mu \right)}\overline{E}\left(s\right)\\ {\overline{c}}_{44}\left(s\right)={\overline{c}}_{55}\left(s\right)={\overline{c}}_{66}\left(s\right)=\frac{1}{2\left(1+\mu \right)}\overline{E}\left(s\right)\end{array}\right.$$

(16)

For common viscoelastic models, the corresponding values of their elastic modulus are presented in

Table 6. Corresponding values of elastic modulus of classical viscoelastic models.

	Classic Kelvin Model	Classic Maxwell Model	Three Parameter Solid Model
$\overline{E}\left(s\right)$	$E+\eta s$	$\frac{E\eta s}{E+\eta s}$	$\frac{E_2\left(E_1+\eta s\right)}{E_1+E_2+\eta s}$

.

Based on the above assumptions, the constitutive relationship of the three-parameter solid model can be inferred as follows:

$$\left\{\begin{array}{l}\left(1+\frac{\eta }{E_1+E_2}d\right){\sigma }_{ii}=\frac{E_2\left(E_1+\eta d\right)}{\left(1+\mu \right)\left(E_1+E_2\right)}{\epsilon }_{ii}+\frac{3\mu E_2\left(E_1+\eta d\right)}{\left(1-2\mu \right)\left(1+\mu \right)\left(E_1+E_2\right)}{\epsilon }_h\\ \left(1+\frac{\eta }{E_1+E_2}d\right){\sigma }_{ij}=\frac{E_2\left(E_1+\eta d\right)}{2\left(1+\mu \right)\left(E_1+E_2\right)}{\gamma }_{ij}\end{array}\right.$$

(17)

$$\left\{\begin{array}{l}\mathsf{\Delta }{\sigma }_{ii}=A\mathsf{\Delta }{\epsilon }_{ii}+B\mathsf{\Delta }{\epsilon }_h+C{\epsilon }_{ii}+D{\epsilon }_h-E{\sigma }_{ii}\\ \mathsf{\Delta }{\sigma }_{ij}=\frac{1}{2}A\mathsf{\Delta }{\gamma }_{ij}+\frac{1}{2}C{\gamma }_{ij}-E{\sigma }_{ij}\end{array}\right.$$

(18)

where $\left\{\begin{array}{l}A=\frac{E_2\left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]},B=\frac{3\mu E_2\left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left(1-2\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]}\\ C=\frac{2E_1{E}_2\mathsf{\Delta }t}{\left(1+\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]},D=\frac{6\mathsf{\Delta }t\mu E_1{E}_2}{\left(1+\mu \right)\left(1-2\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]}\\ E=\frac{2\left(E_1+E_2\right)\mathsf{\Delta }t}{\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta }\end{array}\right.$.

Its Jacobian matrix is as follows:

$$\left\{\begin{array}{l}{c}_{11}=\frac{\left(1-\mu \right)E_2\left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left(1-2\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]}\\ c_{12}=\frac{E_2\left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left(1-2\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]}\\ c_{44}=\frac{E_2\left(E_1\mathsf{\Delta }t+2\eta \right)}{2\left(1+\mu \right)\left[\left(E_1+E_2\right)\mathsf{\Delta }t+2\eta \right]}\end{array}\right.$$

(19)

When $E_2\to \infty$, the three-parameter solid model degenerates to the Kelvin model, and the relationship between the stress increment and the strain increment and the Jacobian matrix are as follows:

$$\left\{\begin{array}{l}\mathsf{\Delta }{\sigma }_{ii}=A\mathsf{\Delta }{\epsilon }_{ii}+B\mathsf{\Delta }{\epsilon }_h+C{\epsilon }_{ii}+D{\epsilon }_h-E{\sigma }_{ii}\\ \mathsf{\Delta }{\sigma }_{ij}=\frac{1}{2}A\mathsf{\Delta }{\gamma }_{ij}+\frac{1}{2}C{\gamma }_{ij}-E{\sigma }_{ij}\end{array}\right.$$

(20)

where $\left\{\begin{array}{l}A=\frac{E_1\mathsf{\Delta }t+2\eta }{\left(1+\mu \right)\mathsf{\Delta }t},B=\frac{3\mu \left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left(1-2\mu \right)\mathsf{\Delta }t}\\ C=\frac{2E_1}{1+\mu },D=\frac{6\mu E_1}{\left(1+\mu \right)\left(1-2\mu \right)}\\ E=2\end{array}\right.$.

Its Jacobian matrix is:

$$\left\{\begin{array}{l}{c}_{11}=\frac{\left(1-\mu \right)\left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left(1-2\mu \right)\mathsf{\Delta }t}\\ c_{12}=\frac{\mu \left(E_1\mathsf{\Delta }t+2\eta \right)}{\left(1+\mu \right)\left(1-2\mu \right)\mathsf{\Delta }t}\\ c_{44}=\frac{E_1\mathsf{\Delta }t+2\eta }{2\left(1+\mu \right)\mathsf{\Delta }t}\end{array}\right.$$

(21)

When $E_1=0$, the three-parameter solid model degenerates to the Maxwell model, and the relationship between the stress increment and the strain increment and the Jacobian matrix are as follows:

$$\left\{\begin{array}{l}\mathsf{\Delta }{\sigma }_{ii}=A\mathsf{\Delta }{\epsilon }_{ii}+B\mathsf{\Delta }{\epsilon }_h-C{\sigma }_{ii}\\ \mathsf{\Delta }{\sigma }_{ij}=\frac{A}{2}\mathsf{\Delta }{\gamma }_{ij}-C{\sigma }_{ij}\end{array}\right.$$

(22)

where $\left\{\begin{array}{l}A=\frac{2E_2\eta }{\left(1+\mu \right)\left(E_2\mathsf{\Delta }t+2\eta \right)},B=\frac{6\mu E_2\eta }{\left(1+\mu \right)\left(1-2\mu \right)\left(E_2\mathsf{\Delta }t+2\eta \right)}\\ C=\frac{2E_2\mathsf{\Delta }t}{E_2\mathsf{\Delta }t+2\eta }\end{array}\right.$.

Its Jacobian matrix is:

$$\left\{\begin{array}{l}{c}_{11}=\frac{2\left(1-\mu \right)E_2\eta }{\left(1+\mu \right)\left(1-2\mu \right)\left(E_2\mathsf{\Delta }t+2\eta \right)}\\ c_{12}=\frac{2\mu E_2\eta }{\left(1+\mu \right)\left(1-2\mu \right)\left(E_2\mathsf{\Delta }t+2\eta \right)}\\ c_{44}=\frac{E_2\eta }{\left(1+\mu \right)\left(E_2\mathsf{\Delta }t+2\eta \right)}\end{array}\right.$$

(23)

According to the above content, the code of the UMAT subroutine is shown in the Appendix A (Table A1).

4.2.2. Evaluation of Several Classical Viscoelastic Models for the DM Description of Asphalt Mixtures

This work used the DM calculation formula in Section 4.2.1 to fit the DM calculated from the DMMC of the three HCRAM established in Section 4.1.2. Figure 16 shows the fitting curve of the DM of ARHM-13 at diverse temperatures, and the viscoelastic parameters obtained by fitting each asphalt mixture are listed in

Table 7. Viscoelastic model parameters of various HCRAM.

Mixture	Temperature (°C)	Fitting Parameters of Viscoelastic Models
		Classic Kelvin Model			Classic Maxwell Model			Three-Parameter Solid Model
		E	η	R2	E	η	R2	E1	E2	η	R2
ARHM-13	15	4929.89	2.69	0.46	12,602.79	957.69	0.85	3115.43	14,552.65	43.08	0.93
	30	1683.62	1.04	0.54	4604.64	299.60	0.82	1031.70	5526.33	13.03	0.94
	45	589.37	0.32	0.55	1461.25	108.33	0.75	433.68	1755.65	4.24	0.94
	60	340.99	0.16	0.53	736.44	67.92	0.63	305.63	867.06	2.45	0.93
ARHM-20	15	6659.82	3.21	0.43	15,638.63	1395.49	0.82	4740.86	17,698.89	61.16	0.91
	30	2345.27	1.31	0.51	5999.82	440.16	0.81	1553.55	7061.87	18.81	0.93
	45	1169.86	0.64	0.53	2906.24	218.61	0.78	835.69	3451.10	8.86	0.94
	60	522.61	0.24	0.52	1123.76	105.48	0.64	466.49	1314.37	3.86	0.93
ARHM-25	15	8771.20	4.18	0.41	20,531.81	1861.68	0.83	6191.55	23,092.45	83.16	0.91
	30	2933.94	1.70	0.51	7734.54	537.87	0.82	1856.92	9148.62	23.51	0.94
	45	1634.33	0.95	0.54	4249.58	295.72	0.79	1095.27	5088.36	12.33	0.94
	60	693.00	0.33	0.54	1557.14	134.24	0.68	586.32	1843.44	4.97	0.94

.

Figure 16 illustrates how the DM determined by the three-parameter solid model is the most similar to the test data of the three models, followed by Maxwell model, while the Kelvin model has the most significant difference. This result is also verified in Table 7. It demonstrates that the viscoelastic parameters and correlation coefficients were obtained by fitting the three models in different HCRAM. In contrast, the three-parameter solid model has the largest correlation coefficient, so the fitting effect of the model is better, which aligns most closely with the outcomes of the DM test. The constitutive equation of the three-parameter solid model can comprehensively describe the instantaneous elasticity, creep, stress relaxation, and other behaviors of viscoelastic materials, showing the characteristics of solids. Therefore, using the three-parameter solid model to describe the viscoelastic mechanical property of the HCRAM is most appropriate.

4.3. Numerical Simulation Validation and Result Analysis

Using Bisar3 software and Abaqus-6.14 software to calculate the mechanical response of the pavement structure in Table 4, we determined that the load area is the unit load, and the load center, half of the load circle radius, and the load circle circumference are taken as the mechanical response analysis points, referred to as points 1, 2, and 3, respectively. Figure 17 shows a comparison of calculated mechanical response values. From Figure 17, we see that the calculation results of Abaqus are smaller than those of Bisar. However, the differences in longitudinal stress calculated by the two software are minimal, indicating that the difference between the analytical and numerical solutions is insignificant. Therefore, the numerical simulation method used in this article is practical and applicable.

4.4. Dynamic Response Analysis of Pavement Structures and Its Differential Analysis

This study mainly analyzed the structural dynamic mechanical response of the pavement of Jingde Expressway under the conditions of driving speed of 10.8 km/h, vertical load of 0.7 MPa, horizontal load of 0.3 MPa, and temperature of 15 °C. The mechanical response indexes include equivalent stress, longitudinal stress, deflection, and longitudinal strain, and the variation laws along the depth direction, longitudinal along the surface, and transverse to the pavement were studied. As illustrated in Figure 18, several distinctive aspects of mechanical response analysis were identified here to aid in the study. Based on the mechanical response findings determined by the three-parameter solid model, this chapter first examined the dynamic mechanical response of the Jingde Expressway pavement structure and then analyzed the differences of the results calculated by different viscoelastic models. To facilitate the difference analysis, only feature point A was selected for comparative analysis.

Figure 19 shows the cloud map results of three viscoelastic models. Figure 19 shows that the mechanical indicators of the pavement structure calculated based on the classical Kelvin model, classical Maxwell model, and three-parameter solid model have roughly similar changes in different directions, but there are still slight differences. In addition, the differences between the calculation results are also relatively significant. Therefore, in the following section, the differences between the sizes of the mechanical indicators will be emphasized.

4.4.1. Equivalent Stress Dynamic Response

The pavement equivalent stress is an essential parameter in pavement structure design and evaluation. Accurately calculating the pavement equivalent stress is essential to ensure the safety and durability of pavement structure [51]. Therefore, the equivalent stress dynamic response was selected for analysis, and its calculation method is as follows [52]:

$${\sigma }_e=\sqrt{\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right]/2}$$

(24)

where ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$ is the first, second, and third principal stresses.

Figure 20 shows the equivalent stress curve based on the three-parameter solid model. Figure 20a shows the variation curve of equivalent stress along the depth direction. It manifests that the base layer differs somewhat, whereas the four characteristic points on the surface layer differ significantly. The maximum value of equivalent stress at each characteristic point is in the order of A > B > D > C from large to small, and the calculated values are 0.59 MPa, 0.52 MPa, 0.17 MPa, and 0.31 MPa, respectively. Therefore, for the pavement structure, the equivalent stress at the wheel load center is highest, and the equivalent stress at the wheel gap center is the lowest, meaning that damage is more likely to happen at the wheel load center of the pavement. Figure 20b and Figure 20c, respectively, show the variation curves of equivalent stress along the longitudinal and transverse directions of the pavement. It is evident that each layer’s maximum equivalent stress ratio is arranged from large to tiny as follows: base layer bottom, surface layer bottom, and surface layer. This phenomenon occurs because of the asphalt surface layer being in direct contact with the traffic load and being subject to the impact, vertical, and horizontal forces of that load. Furthermore, the base layer is composed of a stiff material, while the surface layer is composed a viscoelastic material. In contrast, the overall stiffness of the surface layer is more minor, so the equivalent stress of the surface layer is the largest.

Figure 21 compares the equivalent stress time history curves based on different models. In Figure 21, the abscissa is the time axis, and 2 s is when the load completes a movement on the belt. Under the action of moving load, the equivalent stress of different viscoelastic models is the same with time, and the maximum value appears when t = 1 s. Based on Figure 22, it can be observed that the equivalent stress calculated by the three viscoelastic models at point A of the surface layer from large to small are the three-parameter solid model, classical Kelvin model, classical Maxwell model, and the equivalent stress response values are 0.57 MPa, 0.56 MPa, and 0.55 MPa, respectively. The reason for this phenomenon is that, among the viscoelastic model fitting parameters calculated by using the DMMC, the modulus of the classical Maxwell model is the largest, the stronger with regard to the stiffness, the stronger with regard to the resistance to the load, and the smaller with regard to the equivalent stress.

4.4.2. Longitudinal Stress Dynamic Response

Figure 23 shows the longitudinal stress curve based on the three-parameter solid model. According to Figure 23a, the equivalent stress is the same as that the maximum longitudinal stress of each characteristic point in the pavement structure is it arranged from large to small as A > B > D > C, and the maximum calculated longitudinal stress of point A is −0.39 MPa, which is 1.63 times, 2.2 times, and 2.0 times of points B, C, and D, respectively. Further, it is discovered that the pavement is in a state of compression, and the longitudinal stress at points A and B generally decreases along the pavement depth. In contrast, the longitudinal stress at points C and D first increases and then decreases. However, most of the base is in tension, and the longitudinal stress increases along the depth of the pavement. Figure 23b and Figure 23c show the variation curves of longitudinal stress along the pavement’s longitudinal and transverse directions, respectively. The figure illustrates that the equivalent stress is the same; the highest longitudinal stress appears in the surface layer and is located at the longitudinal distance of 173 mm from the center of the wheel load, and the response value is −0.65 MPa. It indicates that the most unfavorable point of the pavement is not necessarily located at the characteristic point specified in the “Specifications for Design of Highway Asphalt Pavement”. Simultaneously, it also proves that it is more reasonable to calculate the most unfavorable point using the global search method in China’s “Design Standard of Durable Asphalt Pavement for Highway with Structural Life Increment”.

Figure 24 is a comparison diagram of the longitudinal stress time history curves based on the different models. From an analysis of Figure 24, it is apparent that the time history curves of longitudinal stress calculated by the three models are the same. In the process of the load movement, the longitudinal stresses calculated by the three models show positive and negative alternating changes, that is, tension-compression alternating phenomenon. This phenomenon is because, within 0 s–0.875 s, most of the surface layer is in a tension state and has also changed to a compression state. It is preliminarily considered that the load has just begun to move and has not yet stabilized, so it fluctuates wildly. Within 0.875 s–1 s, that is, before the arrival of the wheel, the asphalt concrete is under pressure. As the tire approaches, the longitudinal compressive stress gradually increases and reaches its maximum when the tire arrives. As the tire moves away from the tire, the longitudinal compressive stress gradually decreases and then transits to the tensile stress.

Figure 25 depicts the comparison of longitudinal stresses of the different viscoelastic models. The figure clearly illustrates that the maximum longitudinal stress at point A of the surface layer calculated by the three viscoelastic models is ranked from large to small as follows: Maxwell, Kelvin, and three-parameter solid model, and the response sizes are −0.47 MPa, −0.44 MPa, and −0.39 MPa, respectively. Combined with Figure 23a, the longitudinal stress based on the classical Maxwell model is the fastest with the depth change. This is because this model can better reflect the instantaneous elastic deformation and relaxation phenomenon than other models, and its properties are closer to liquid.

4.4.3. Dynamic Response of Deflection

Deflection is an essential index for the pavement design and construction. This paper studied the dynamic response of the deflection (vertical displacement). Figure 26 demonstrates the deflection curve based on the three-parameter solid model. As shown from Figure 26, the wheel load center of the surface layer is the location of the maximum deflection value, which is 0.1 mm., indicating that this position is more prone to pavement settlement or deformation than the other three characteristic points.

Figure 27 is a comparison chart of the bending time history curves based on the different models. In Figure 26, the peak value of the curve appears at t = 1 s; that is, the load is just at point A. In addition, we can see that in 1 s–2 s, the load progressively shifts away from point A, the slowest rate of change is found in the deflection response based on the Maxwell model, and the three-parameter solid model has the fastest rate of change. In other words, the residual value of the deflection response calculated by the classical Maxwell model is the largest after the movement of a load. Figure 28 is a comparison diagram of the deflection based on the different models. It is evident from the diagram that the different viscoelastic models obtain the maximum deflection response of the different layers, and the three-parameter solid model calculates the maximum deflection of the surface layer. The base course bottom and surface layer bottom maximum values are computed by the classical Maxwell model. The surface layer, surface layer bottom, and base layer bottom have maximum deflection values of 0.1 mm, 0.0799 mm, and 0.072 mm, respectively.

4.4.4. Longitudinal Strain Dynamic Response

Strain refers to the degree of deformation of pavement under external load. The longitudinal strain of asphalt pavement is a crucial design index. Figure 29 shows the longitudinal strain curve based on the three-parameter solid model. From Figure 29a, it is demonstrated that the greatest longitudinal strain at the four characteristic points is arranged from large to small as C > D > A > B, and the response values are −41.89 με, −40.43 με, −23.49 με, −28.14 με. This phenomenon is because the characteristic point C is the wheel gap center of the road, which is squeezed by the wheel, resulting in the superposition effect of two-wheel actions [45], so the longitudinal strain of point C is considerable. Furthermore, it is apparent that the asphalt pavement is subject to the longitudinal compressive strain, with a progressive reduction in the longitudinal strain as depth increases. The longitudinal strain progressively grows as the depth deepens and reaches its greatest value at the base’s bottom. Figure 29b,c demonstrate that along the pavement’s longitudinal direction, the surface layer’s longitudinal strain has a peak value of 134.6 με at 210 mm from the wheel load center. Along the longitudinal and transverse directions of the pavement, the variation trend of longitudinal strain is the same. With the direction away from the load, the response value gradually tends to zero, and more speed is placed on the longitudinal strain attenuation rate of each structural layer compared to on the transverse strain attenuation rate.

Figure 30 is a comparison diagram of the longitudinal strain time history curves based on different models. Combined with Figure 24, the longitudinal strain exhibits a similar temporal change pattern as the longitudinal stress, and there will be tension compression alternation during the load movement. The difference is that the peak value of the curve does not appear at t = 1 s. This is because the damping effect is considered in the analysis of moving load, and the existence of damping affects the dynamic time history response, resulting in the hysteresis of the strain response. A comparison diagram of the longitudinal strain based on the several models is shown in Figure 31. From Figure 30 and Figure 31, it is evident that the longitudinal strain at point A of the surface layer based on the three-parameter solid model is the largest, which is 1.007 and 2.223 times that of the Maxwell and the Kelvin model, respectively. Moreover, its attenuation rate in the asphalt pavement is the slowest. The reason for this is that the asphalt pavement surface is defined as viscoelastic material in the finite element calculation, and the properties of the three-parameter solid model are closer to solids. Compared with other models, the strain change rate is slower when the depth is increased.

5. Conclusions

This study selected the pavement structure and materials of the Jingde Expressway in China as the research object, explored the viscoelastic characteristics of the HCRAM, and studied the mechanical response and differences of the asphalt pavement structures utilizing three viscoelastic models of the HCRAM. The following are the primary conclusions:

(1)

The relationship between temperature, frequency, and the DM was obtained through DM testing of the HCRAM. Specifically, the DM had a negative correlation with temperature and a positive correlation with loading frequency. The relationship between the storage and loss modulus and temperature and frequency had a similar pattern to that of the DM. The phase angle exhibited complicated variations with the frequency and temperature. At low temperatures, the phase angle decreased as the loading frequency increased. At high temperatures, the phase angle of the HCRAM increased as the loading frequency increased. This reveals that under the action of high-frequency and low-temperature loads, the asphalt mixtures exhibited significant elastic characteristics, and their viscosity rose proportionally with the rise in temperature.

(2)

Utilizing the DM test results, a DMMC was established based on the Sigmoidal model. The master curve demonstrated that the asphalt mixture’s reliance on temperature intensified as the rubber fraction and gradation increased. The rate at which the DM changed from low frequency to high frequency steadily slowed down as temperature rose, indicating that the impact of the loading frequency on the asphalt mixture decreased as the temperature increased.

(3)

This study found that the three-parameter solid model most accurately described the viscoelastic dynamic mechanical behavior of HCRAM. From the results of mechanical response calculations, we found that the equivalent stress, longitudinal stress, and bending response of characteristic point A were all the largest, with the response values of 0.59 MPa, −0.39 MPa, and 0.1 mm, respectively. Under the comparison of three viscoelastic models, most of the response indexes in the surface layer calculations were all the largest in the three-parameter viscoelastic model.

(4)

This study systematically investigated the dynamic mechanical response of road surfaces under different viscoelastic models. However, there has yet to be research on the fatigue and high and low-temperature performance of high-content rubber asphalt mixtures, nor has the influence of temperature, load, and speed on the dynamic mechanical response of pavement structures been considered. In addition, it is necessary to compare and analyze the numerical simulation results with the accelerated loading test data. This will be carried out in subsequent work, and further consideration of the above research content is needed.