Next Article in Journal
The One-Dimensional Flow Pressure Loss Correction Model Based on the Particle Flow through Concrete Bend
Previous Article in Journal
Prediction of Leakage Flow Rate and Blow-Down in Brush Seals via 2D CFD Simulation with Porosity Correction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamic Response Analysis of Asphalt Pavement under Pavement-Unevenness Excitation

1
School of Water Conservancy and Transportation, Zhengzhou University, Zhengzhou 450001, China
2
State Key Laboratory of Tunnel Boring Machine and Intelligent Operations, Zhengzhou University, Zhengzhou 450001, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2024, 14(19), 8822; https://doi.org/10.3390/app14198822
Submission received: 29 August 2024 / Revised: 22 September 2024 / Accepted: 26 September 2024 / Published: 30 September 2024
(This article belongs to the Special Issue Advanced Structural Health Monitoring in Civil Engineering)

Abstract

:
This paper investigates and analyzes the dynamic response of asphalt pavement under pavement-unevenness excitation based on an orthogonal vector function system and efficient DVP (dual variable and position) method. Firstly, starting from the pavement unevenness of the vehicle excitation source, the pavement-unevenness excitation is established by using the filtered white-noise method, and the random load of the vehicle model is obtained by simulation. Then, based on the basic governing equation of the road-surface problem under the random load, the analytical solution of the road-surface mechanical response is obtained by using the orthogonal vector function system and DVP method. The effects of pavement-unevenness grade, vehicle speed, vehicle load, interlayer contact condition, and transverse isotropy on the mechanical response of the road surface are analyzed via the analytical results. The results show that DVP can effectively solve the dynamic response of pavements under the excitation of pavement unevenness; in addition, it can also be applied to certain situations, such as transverse isotropy of materials and interface conditions. The results show that the pavement unevenness does not affect the average stress and strain of each layer but has a significant effect on the peak value and dispersion degree. An increase in vehicle speed causes a peak in strain and a larger coefficient of variation. Poor bonding between interfaces can lead to increased stress and strain at the bottom of the surface layer.

1. Introduction

Pavement unevenness acts as a medium between the vehicle and the pavement, and the vehicle generates random loads under its excitation, which increases the mechanical response within the pavement and accelerates the deterioration of the pavement serviceability. Therefore, the response of pavement under the excitation of pavement unevenness is a concern of road engineers. In the study of this kind of problem, scholars use the linear filtering white-noise method, the harmonic superposition method (or trigonometric series method), AR and ARMA methods, and the Poisson method to establish a pavement-unevenness model. Vehicles are simplified into quarter-vehicle models, half-vehicle models, and whole-vehicle models, and the road surface is assumed to be a beam or plate model on an elastic foundation [1,2,3,4,5].
Initially, scholars regarded the vehicle–road interaction as a weakly coupled problem which could be solved by sequential decoupling. Because the vertical deformation of the road surface is usually much smaller than the displacement of vehicle suspension and tire, and the vehicle–road coupling has little influence on vehicle vibration, scholars usually assume the road surface as a rigid body to solve the vehicle random load and then apply the random load to the road-surface model to solve the mechanical response of the road surface [6]. Cebon et al. [7] used the quarter-vehicle model to study the impact of heavy vehicles on the road surface and optimized the suspension of heavy vehicles based on the road-surface results. Kim et al. [8] used the triple Fourier transform and mode superposition method to study the dynamic problem of infinite plates on the Winkler foundation under the action of vehicle dynamic load, and the results showed that the increase of pavement unevenness and vehicle wheelbase would lead to the increase of maximum deflection and stress. Uzzal et al. [9] conducted an in-depth study on the dynamic response of Eulaberly beams on the Pasternak foundation under the influence of moving loads. Huang et al. [10] conducted an in-depth study on the dynamic response of Winkler foundation thin plate under the influence of moving load and observed that the foundation stiffness, vehicle speed, and frequency would all affect the dynamic response of pavement. Lyu et al. [11] considered the vehicle load caused by pavement unevenness. The road surface was regarded as a Kirchhoff sheet, and the roadbed was made of elastic material with saturated holes. Zhang et al. [12] explored the finite integral transform approach for analytical thermal buckling solutions of composite thin plates with all edges rotationally restrained. Li et al. [13,14] simplified the vehicle as a four-degree-of-freedom, half-vehicle model and the asphalt pavement as an infinite double-layer Kirchhoff sheet on the Kelvin foundation. They solved the analytical solution of the dynamic response of the pavement through Fourier integral transformation. The effects of pavement unevenness, vehicle load, suspension stiffness, suspension damping, tire stiffness, and tire damping on the road’s dynamic response were also discussed.
In recent years, scholars have begun to think that vehicle–road interaction is a coupled system, and it is necessary to establish a coupling model to solve it. Lu et al. [15] proposed a vehicle–road coupling model based on a quarter-vehicle model and pavement unevenness. Elnashar et al. [16] used the Eulbernoulli beam pavement model to study a quarter-vehicle model and found that the coupling effect had more impact on the road-surface deflection than on the wheel displacement. Ding et al. [17] studied the vibration problem of a vehicle–pavement coupled system of Timoshenko beams on a nonlinear foundation. Snehasagar et al. [18] discussed the dynamic response of vehicle–pavement systems based on the Eulabenouli beam model on elastic foundations, taking into account the viscoelastic properties of asphalt materials. Yang et al. [19] and Li et al. [20] further improved the vehicle–road coupling model, simplified the pavement unevenness into a sinusoidal road surface and the vehicle into a seven-degree-of-freedom vehicle model, and studied the thin-plate problem on linear elastic and viscoelastic foundation successively. Krishnanunn et al. [21] proposed an iterative solution method for the vehicle–road coupling system based on a quarter vehicle and Timoshenko beam, which can quickly solve the vehicle–road coupling system.
In the above sequential decoupling or coupling analysis of vehicle–pavement interaction, the main pavement model used is a beam or plate on an elastic foundation, which is particularly suitable for studying the response of cement pavement. For asphalt pavement, the common pavement model regards it as a multilayer elastic system, which can better reflect the actual structure of asphalt pavement. Although there are many analytical methods for the mechanical response of multilayer elastic asphalt pavement under dynamic load, there are relatively few analytical solutions considering vehicle loads. Most scholars mainly use the finite-element method for simulation analysis to solve the mechanical response of asphalt pavement under vehicle loads. Kropac et al. [22] proposed a classification method for road unevenness based on power spectral density, considering the stochastic longitudinal road irregularity classification scheme, the road waviness and vehicle response, and both the vehicle vibration response and the wide-road waviness interval (road elevation, PSD slope). Based on the one-quarter vehicle vibration model, Li et al. [23] used Simulink software (R2015a) to simulate the maneuverable load of vehicles on the road surface and then established a finite-element analysis model of semi-rigid asphalt pavement to obtain the mechanical response of the road surface under maneuverable load. Xing et al. [24] used ANSYS (R18.0) to build a finite-element model of asphalt pavement. Based on MATLAB (R2018a), the dynamic response of the pavement was obtained when the vehicle load was applied to the model, and the influence of vehicle speed and load amplitude on the pavement was discussed. Zhao et al. [25] established a vehicle model to obtain vehicle dynamic load and then established a three-dimensional finite-element model to evaluate the road response under harmonic excitation. Hettiarachchi et al. [26] summarized smoothness specifications in the United States and various countries around the world and provided an overview of the comparison of smoothness measurement. The smoothness measuring devices, indices, limits, and payment schemes vary from country to country. The International Roughness Index is the most used smoothness index in the US and around the world. Inertial profilers are used to measure the pavement profile for smoothness. Wang et al. [27] used the vehicle model to solve the tire load caused by different roughness of the road surface and used the three-dimensional finite-element (FE) model to simulate the pavement response under dynamic load with non-uniform tire contact stress. The results showed that the position of rougher pavement would cause greater pavement response and accelerate pavement failure. Liu et al. [28] developed a finite-element model of tire–pavement interaction to study the effect of pavement unevenness on the mechanical response of asphalt pavement. The results show that the pavement irregularity significantly affects the mechanical response of the asphalt pavement, thereby affecting the service life of the asphalt pavement. Although the finite-element method can consider complex mechanical problems and explain the actual vehicle–pavement interaction mechanism, the finite-element method has the disadvantages of large computational cost, long calculation time, and non-convergence of results, and the analytical method is still a favored and efficient method.
Recently, based on the exact-integration method and the dual-vector technique, Zhang et al. [29] obtained the analytical solution of the isotropic displacement of each layer of material under a rectangular load. In order to avoid the influence of the exponential term of e, the fine-integration method divides each layer into several fine sub-layers in the calculation, which greatly increases the computational cost and leads to low efficiency of the algorithm. But inspired by the unique properties of fine integration, Liu et al. [30,31,32] proposed a new method, the DVP (dual variable and position) method, which can efficiently deal with multilayer elastic pavement with different pavement structures. To the best of the authors’ knowledge, the analytical response of a multilayer pavement system under the excitation of pavement unevenness based on the DVP method has not been reported, which inspired the work of this paper. This paper will use this method to solve the analytical solution of the mechanical response of the pavement under the excitation of unevenness. Firstly, starting from the pavement unevenness of the vehicle’s excitation source, the filtered white-noise method is used to establish the unevenness excitation, and the random dynamic load of the vehicle model is simulated and solved. Then, based on the basic control equation of the pavement problem under random dynamic load, the analytical solution of the pavement mechanical response is obtained using the orthogonal vector function system and the DVP method. The analytical results are used to analyze the influence of various factors on the pavement mechanical response, providing a certain technical reference for the calculation of road mechanical response.

2. Solution of Vehicle Load under Pavement-Unevenness Excitation

2.1. Establishment of Pavement-Unevenness Model

The establishment of a random pavement-unevenness excitation model is the basis for the study of random road loads exerted by vehicles. In this paper, the filtered white-noise method is used to build a time-domain model of four-wheel pavement unevenness to provide pavement excitation for the solution of maneuvering loads under the seven-degree-of-freedom vehicle model.
In this paper, the power-spectral-density function is used to characterize the classification and statistical characteristics of pavement unevenness. The fitting expression of pavement-unevenness spectral density Gq (n) within the spatial frequency n1 < n < n2 is as follows:
G q ( n ) = G q ( n 0 ) ( n n 0 ) w
where n—spatial frequency; n0—reference spatial frequency, usually n0 = 0.1 m−1; Gq(n0)—pavement-unevenness coefficient; w—frequency index, where in general w = 2. According to the power spectral density of the road surface, the uneven condition of the road surface is divided into different grades, and the quality of the road surface gradually decreases with Gq(n0) increasing [33].
q ˙ 1 ( t ) = 2 π v n q q 1 ( t ) + 2 π n 0 G q ( n 0 ) v ω ( t )
In this paper, the filtered white-noise method is used to generate random pavement-unevenness excitation. Assuming that the left-front wheel, right-front wheel, left-rear wheel, and right-rear wheel are stimulated by road q1(t), q2(t), q3(t), and q4(t) in turn, Formula (2) represents the time-domain model of the uneven pavement of the left-front wheel.
In Formula (2), qn(t) is the pavement-unevenness excitation of the wheel, nq is the cutoff space frequency under the pavement space, taking nq = 0.011 m−1; ω(t)—white-noise signal; v—vehicle speed in m/s.
Assume that the front- and rear-wheel road excitation are exactly the same but with a lag of some time. The lag time td is the ratio of the wheelbase L of the front and rear wheels to the speed v; that is
t d = L v
The pavement-unevenness excitation of the front and rear wheels has the following relationship:
q 3 ( t ) = q 1 ( t ) ( t t d )
Due to the time-lag correlation between the front and rear wheels, the Pade approximation algorithm can be adopted. The transfer function G31 between the input of pavement unevenness of the front and rear wheels is approximately as follows:
G 31 = q 3 ( s ) q 1 ( s ) = 1 t d 2 s + t d 2 12 s 2 1 + t d 2 s + t d 2 12 s 2
The Laplace inverse transformation is performed on the above Equation (5), and the time-domain relationship between the front-wheel road input and the rear-wheel road input is obtained as follows:
q ˙ 3 ( t ) = 2 t d q 3 ( t ) + 2 t d q 1 ( t ) q ˙ 1 ( t )
The time-domain function q1(t) of front-wheel pavement unevenness is substituted into the above formula to obtain the following:
q ˙ 3 ( t ) = ( 2 t d + 2 π n q v ) q 1 ( t ) 2 t d q 3 ( t ) 2 π n 0 G q ( n 0 ) v ω ( t )
Assuming that the self-spectral density of the left and right wheels of the vehicle is the same and the phase difference is 0, the coherence function of the left-front-wheel road input q1(t) and the right-front-wheel road input q2(t) is as follows:
c o h ( ω ) = s 12 ( ω ) s 1 ( ω )
where, s12(ω) is the left- and right-wheel mutual-spectral density function; s1(ω) is the self-spectral density function.
Since the frequency-response function is the ratio of cross-spectral density to self-spectral density, that is, the frequency-response function H12(ω) of the system is as follows:
H 12 ( ω ) = c o h ( ω ) = s 12 ( ω ) s 1 ( ω )
Considering the wheelbase B of the left and right wheels, the second-order Pade algorithm is adopted to obtain the transfer function G21 between the road input of the left and right wheels
G 21 = q 2 ( s ) q 1 ( s ) = 1 B 2 v s + B 2 12 v 2 s 2 1 + B 2 v s + B 2 12 v 2 s 2
Laplace inverse transformation is performed on the above formula, and the time-domain relationship between left-wheel road input q1(t) and right-wheel road input q2(t) is obtained as follows:
q ˙ 2 ( t ) = 2 v B q 2 ( t ) + 2 v B q 1 ( t ) q ˙ 1 ( t )
The left-front-wheel pavement-unevenness time-domain function q1(t) is substituted into the above formula thus:
q ˙ 2 ( t ) = ( 2 v B + 2 π n q v ) q 1 ( t ) 2 v B q 2 ( t ) 2 π n 0 G q ( n 0 ) v ω ( t )
According to the correlation between input q2(t) and q4(t) of the pavement unevenness of the right-front and -rear wheels,
q ˙ 4 ( t ) = 2 t d q 4 ( t ) + 2 t d q 2 ( t ) q ˙ 2 ( t )
By substituting the right-front-wheel pavement-unevenness time-domain function q2(t) into the above Equation (13), we can obtain
q ˙ 4 ( t ) = ( 2 v B + 2 π n q v ) q 1 ( t ) + ( 2 v B + 2 t d ) q 2 ( t ) 2 t d q 4 ( t ) + 2 π n 0 G q ( n 0 ) v ω ( t )
The differential equation of the four-wheel pavement-unevenness excitation model is deduced above, and the simulation model of four-wheel pavement excitation can be established according to the differential equation.

2.2. The Maneuverable Load of the Whole-Vehicle Model Is Solved

In this paper, a seven-degree-of-freedom vehicle model is adopted, and the four-wheel pavement unevenness is used as an incentive input to the vehicle to solve the random load exerted by the wheel on the road surface. The vehicle-dynamics model is shown in Figure 1. In the figure, mb is the body mass, and m1, m2, m3, and m4 are the mass of the four wheels, respectively. kt1, kt2, kt3, and kt4 are the elastic stiffness of the four tires. ks1, ks2, ks3, and ks4 are the elastic stiffness of the suspension on the wheel, respectively. cs1, cs2, cs3, and cs4 are the damping of the suspension on the wheel respectively. z1, z2, z3, and z4 are the vertical displacements of the four tires, respectively. zb1, zb2, zb3, and zb4 are the vertical displacement of the suspension on the tire, that is, the displacement at the four corners of the body. zb is the vertical displacement at the center of mass of the body; q1, q2, q3, and q4 are the excitation of pavement unevenness at four wheels, respectively. a is the distance between the front car and the center of mass of the vehicle; b is the distance between the rear wheel and the center of mass of the vehicle; Bf is the distance between the left and right-front wheels; Br is the distance between the left- and right-rear wheels; θb is the pitch displacement at the center of mass of the car body; φ is the lateral angle displacement at the body center of mass; v stands for direction of travel.
The dynamic equations of the system include the vertical-motion equation of the body’s center of mass, the motion equation of the body’s pitch and roll, and the vertical-motion equation of four-wheel masses.
Vertical-motion equation of body centroid:
m b z ¨ b + c 1 ( z ˙ b 1 z ˙ 1 ) + k s 1 ( z b 1 z 1 ) + c s 2 ( z ˙ b 2 z ˙ 2 ) + k s 2 ( z b 2 z 2 ) + c s 3 ( z ˙ b 3 z ˙ 3 ) + k s 3 ( z b 3 z 3 ) + c s 4 ( z ˙ b 4 z ˙ 4 ) + k s 4 ( z b 4 z 4 ) = 0
Body-pitch-motion equation:
I p θ ¨ b [ c s 1 ( z ˙ b 1 z ˙ 1 ) + k s 1 ( z b 1 z 1 ) ] a [ c s 2 ( z ˙ b 2 z ˙ 2 ) + k s 2 ( z b 2 z 2 ) ] a + [ c s 3 ( z ˙ b 3 z ˙ 3 ) + k s 3 ( z b 3 z 3 ) ] b + [ c s 4 ( z ˙ b 4 z ˙ 4 ) + k s 4 ( z b 4 z 4 ) ] b = 0
Body roll motion equation:
I r φ ¨ + [ c s 1 ( z ˙ b 1 z ˙ 1 ) + k s 1 ( z b 1 z 1 ) ] B f 2 [ c s 2 ( z ˙ b 2 z ˙ 2 ) + k s 2 ( z b 2 z 2 ) ] B f 2 + [ c s 3 ( z ˙ b 3 z ˙ 3 ) + k s 3 ( z b 3 z 3 ) ] B r 2 [ c s 4 ( z ˙ b 4 z ˙ 4 ) + k s 4 ( z b 4 z 4 ) ] B r 2 = 0
Equation of vertical motion of four-wheel masses:
m 1 z ¨ 1 c s 1 ( z ˙ b 1 z ˙ 1 ) k s 1 ( z b 1 z 1 ) + k t 1 ( z 1 q 1 ) = 0 m 2 z ¨ 2 c s 2 ( z ˙ b 2 z ˙ 2 ) k s 2 ( z b 2 z 2 ) + k t 2 ( z 2 q 2 ) = 0 m 3 z ¨ 3 c s 3 ( z ˙ b 3 z ˙ 3 ) k s 3 ( z b 3 z 3 ) + k t 3 ( z 3 q 3 ) = 0 m 4 z ¨ 4 c s 4 ( z ˙ b 4 z ˙ 4 ) k s 4 ( z b 4 z 4 ) + k t 4 ( z 4 q 4 ) = 0
The relationship between suspension displacement and body-mass center displacement is as follows:
z b 1 = z b a θ b B f 2 φ z b 2 = z b b θ b + B f 2 φ z b 3 = z b + b θ b + B r 2 φ z b 4 = z b + a θ b B r 2 φ
In the upper formula: z ¨ b—vertical acceleration at the center of mass of the body, m/s2; θ ¨ b, φ ¨ —pitch-and-roll acceleration at body center of mass, rad/s2; z ¨ 1, z ¨ 2, z ¨ 3, z ¨ 4—vertical acceleration of the front and rear wheels on the left and right sides, m/s2; z ˙ b—vertical speed at body center of mass, m/s; z ˙ b1, z ˙ b2, z ˙ b3, z ˙ b4—vertical speed of the left- and right-side suspension at the body position, m/s; z ˙ 1, z ˙ 2, z ˙ 3, z ˙ 4—vertical speed of front and rear wheels on left and right sides, m/s; Ip, Ir—moment of inertia of body pitch and roll, kg·m2.
By establishing a seven-degree-of-freedom vehicle simulation model in MATLAB/Simulink (R2021a), the four-wheel pavement-unevenness excitation established above is used as input, the pavement-unevenness grades are different grades, the vehicle parameters are Granada real vehicle model, the vehicle speed is 10 m/s, and the simulation time is 0–20 s. The random loads of the vehicles at different positions of the wheels under grade B and grade C were obtained by simulation.
F 1 = [ b 2 ( a + b ) m b + m 1 ] g + k t 1 ( z 1 q 1 ) = Re ( l = 0 L P l e i ( ω l t + θ l ) ) F 2 = [ b 2 ( a + b ) m b + m 2 ] g + k t 2 ( z 2 q 2 ) = Re ( l = 0 L P l e i ( ω l t + θ l ) ) F 3 = [ a 2 ( a + b ) m b + m 3 ] g + k t 3 ( z 3 q 3 ) = Re ( l = 0 L P l e i ( ω l t + θ l ) ) F 4 = [ a 2 ( a + b ) m b + m 4 ] g + k t 4 ( z 4 q 4 ) = Re ( l = 0 L P l e i ( ω l t + θ l ) )
In the above formula, F is the dynamic load of each wheel; Pl, ωl, and θl are amplitude, angular frequency, and phase of harmonic load, respectively. i is an imaginary unit; when l is 0, ω0 = 0 Hz, θ0 = 0 rad, representing the dead-load part of the vehicle.

3. Solution of Multilayer Pavement under Pavement-Unevenness Excitation

3.1. Pavement Model Description

As shown in Figure 2, the column coordinate system (r, θ, z) is established on the surface of the layered body, and the action on the asphalt surface depends on the maneuver load. In this paper, the load of the left-front wheel of the vehicle model is taken as an example, and the effect of a single-wheel load is superimposed by the superposition method to simulate the effect of the vehicle model on the road surface. In addition, the center of the loading circle is located at the origin of coordinates, the structure layer is numbered from 1 to n from top to bottom, and the coordinate z = z0 = 0 represents the surface layer of the structure. Take the j layer as an example; the upper-interface coordinate of the j layer is z = zj−1, the coordinate of the lower interface is z = zj, and its thickness is hj = zj−zj−1. The contact conditions between the layers can be completely continuous, completely smooth, and semi-bonded, and the bottom is an elastic half-space, where each layer is composed of isotropic or translaterally isotropic materials.

3.2. Fundamental Governing Equation

Since the maneuvering loads acting on the surface of the elastic half-space are multiple harmonic loads with angular frequency ωl, then in the cylindrical coordinate system, all variables can be expressed as functions containing factors ( l = 0 L P l e i ( ω l t + θ l ) ) [34]; for example
u x ( x , y , z , t ) = u x ( x , y , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) u y ( x , y , z , t ) = u y ( x , y , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) u z ( x , y , z , t ) = u z ( x , y , z , ω l ) l = 0 L P l e i ( ω l t + θ l )
Therefore, we only need to ask for the displacement u(x, y, z, ωl) and the stress–strain, which is independent of time and depends only on the angular frequency ωl; unless otherwise stated, the following derived components are proportional to the factor l = 0 L P l e i ( ω l t + θ l ) .
The dynamic balance differential equation, without considering the physical space problem, is as follows:
σ r r r + σ r θ r θ + σ r z z + σ r r σ θ θ r = ρ 2 u r t 2 σ r θ r + σ θ θ r θ + σ θ z z + 2 σ r θ r = ρ 2 u θ t 2 σ r z r + σ θ z r θ + σ z z z + σ r z r = ρ 2 u z t 2
Due to
2 u r ( r , θ , z , t ) t 2 = 2 u r ( r , θ , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) t 2 = ω l 2 u r ( r , θ , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) 2 u θ ( r , θ , z , t ) t 2 = 2 u θ ( r , θ , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) t 2 = ω l 2 u θ ( r , θ , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) 2 u z ( r , θ , z , t ) t 2 = 2 u z ( r , θ , z , ω l ) l = 0 L P l e i ( ω l t + θ l ) t 2 = ω l 2 u z ( r , θ , z , ω l ) l = 0 L P l e i ( ω l t + θ l )
The dynamic balance differential Equation (22) can be rewritten as
σ r r r + σ r θ r θ + σ r z z + σ r r σ θ θ r + ρ ω l 2 u r = 0 σ r θ r + σ θ θ r θ + σ θ z z + 2 σ r θ r + ρ ω l 2 u θ = 0 σ r z r + σ θ z r θ + σ z z z + σ r z r + ρ ω l 2 u z = 0
where ρ is the density, and ωl is the angular frequency.
According to Hooke’s law, the equation of physics is as follows:
σ r r = c 11 ε r r + c 12 ε θ θ + c 13 ε z z    σ θ z = 2 c 44 ε θ z σ y y = c 12 ε r r + c 11 ε θ θ + c 13 ε z z    σ r z = 2 c 44 ε r z σ z z = c 13 ε r r + c 13 ε θ θ + c 33 ε z z    σ r θ = 2 c 66 ε r θ
Among them
c 66 = ( c 11 c 12 ) / 2
In Formula (25), the parameters c11, c12, c13, c33, and c44 are the five elastic constants of the medium. Their size is derived from the following equation:
c 11 = c 33 = E ( 1 v ) / [ ( 1 + v ) ( 1 2 v ) ] c 12 = c 13 = E v / [ ( 1 + v ) ( 1 2 v ) ] c 44 = E / [ 2 ( 1 + v ) ]
where E is the elastic modulus of the material in MPa, and v is Poisson’s ratio.
The final basic governing equation is the relation between strain displacement, and the geometric equation is as follows:
ε r r = u r r ,   ε θ θ = u θ r θ + u r r ,   ε z z = u z z ε θ z = 1 2 ( u θ z + u z r θ ) ε r z = 1 2 ( u z r + u r 2 z ) ε r θ = 1 2 ( u r r θ + u θ r u θ r )

3.3. Boundary and Interlayer Bonding Conditions

When a circular harmonic load is applied on the surface of a multilayer elastic body (z = z0 = 0), the amplitude of which is Pl, the boundary conditions at the surface z = 0 can be expressed as follows:
σ z z | z = 0 = { q r R 0 r > R σ r z | z = 0 = σ θ z | z = 0 = 0 r 0
where R is the radius of the circular load and q = Pl/πR2.
When the layer is not completely continuous at the interface z = zj, the interlayer interface conditions are simplified by using the Goldman model, and the transfer relation expression at z = zj is as follows:
σ r z ( r , θ , z j ) = σ r z ( r , θ , z j + ) σ θ z ( r , θ , z j ) = σ θ z ( r , θ , z j + ) σ z z ( r , θ , z j ) = σ z z ( r , θ , z j + ) u z ( r , θ , z j ) u z ( r , θ , z j + ) = ( 1 / k r ) σ r z ( r , θ , z j ) u θ ( r , θ , z j ) u θ ( r , θ , z j + ) = ( 1 / k r ) σ θ z ( r , θ , z j ) u z ( r , θ , z j ) u z ( r , θ , z j + ) = ( 1 / k z ) σ z z ( r , θ , z j )
In the above formula, kz and kr are the interface association coefficients in the vertical and horizontal directions, respectively. When the interface association coefficients in the vertical and horizontal directions are zero, it means that the interface is completely separated from each other; when they approach infinity, it means that the interface is completely continuous.
In the practical application of the sliding interface model, only the horizontal interface bonding coefficient is considered; that is, the vertical interface bonding coefficient kz is infinite, and the interlayer contact is characterized by changing the horizontal interface bonding coefficient kr.

3.4. Analytic Solutions Based on Dual Vectors

Because the basic governing equation of the pavement problem is a partial differential equation, it is necessary to convert the response component into an ordinary differential equation expressed by an expansion coefficient by virtue of the basic property of the vector function system. In the cylindrical coordinate system, the definition and properties of the orthogonal vector function system are as follows:
L ( r , θ ; λ , m ) = e z S ( r , θ ; λ , m ) M ( r , θ ; λ , m ) = ( e r r + e θ r θ ) S ( r , θ ; λ , m ) N ( r , θ ; λ , m ) = ( e r r θ e θ r ) S ( r , θ ; λ , m )
In the formula, er, eθ, and ez are, respectively, the unit vectors of the r, θ, and z axes in the cylindrical coordinate system, and S is the scalar function, which is defined as follows:
S ( r , θ ; λ , m ) = 1 2 π J m ( λ r ) e i m θ ; m = 0 , ± 1 , ± 2 ,
where Jm(λr) represents the Bessel function of order m, m = 0 corresponds to the axisymmetric deformation, the variables λ and m are the transformation variables of r and θ under the vector function system, and i is the imaginary unit.
Therefore, by importing the displacement, stress, and other functions in the basic governing equation into the orthogonal vector function system, we can get
u ( r , θ , z ) u r e r + u θ e θ + u z e z = m 0 + [ U L ( z ) L ( r , θ ) + U M ( z ) M ( r , θ ) + U N ( z ) N ( r , θ ) ] λ d λ t ( r , θ , z ) σ r e r + σ θ e θ + σ z e z = m 0 + [ T L ( z ) L ( r , θ ) + T M ( z ) M ( r , θ ) + T N ( z ) N ( r , θ ) ] λ d λ
In the formula, UL, UM, UN, TL, TM, and TN represent the vector-function coefficients corresponding to the displacement and stress, respectively.
According to the expansion coefficient in the vector-function system, we can also get the expressions of the stress, displacement, and stress components in the physical domain.
Expression for displacement in the conversion domain:
u r = m 0 + ( U M S r + U N S r θ ) λ d λ u θ = m 0 + ( U M S r θ U N S r ) λ d λ u z = m 0 + U L S λ d λ
The expression of stress in the transformation domain:
σ r z = m 0 + ( T M S r + T N S r θ ) λ d λ σ θ z = m 0 + ( T M S r θ T N S r ) λ d λ σ z z = m 0 + T L S λ d λ
By substituting the displacement expression in the transformation domain into the relationship between strain and stress, the following ordinary differential equations for the expansion coefficients of each layer can be obtained.For the N-type, there are
d d z [ U N T N ] = [ 0 1 c 11 λ 2 c 66 ρ ω l 2 0 ] [ U N T N ]
For the LM-type, there are
d d z [ U L U M T L T M ] = [ 0 λ 2 c 13 c 33 1 c 33 0 1 0 0 1 c 44 ρ ω l 2 0 0 λ 2 0 λ 2 ( c 11 c 13 2 c 33 ) ρ ω l 2 c 13 c 33 0 ] [ U L U M T L T M ]
By means of an orthogonal vector function system, the partial differential equation of the basic governing equation of the layered pavement problem is converted to solve ordinary differential equations of vector function systems UL, UM, UN, TL, TM, and TN. Additionally, the system of ordinary differential equations is divided into two independent parts: LM-type ordinary differential equations composed of the expansion coefficients UL, UM, TM, and TL, and N-type ordinary differential equations composed of the expansion coefficients UN and TN. The LM-type ordinary differential equation is related to vertical normal stress, and the N-type ordinary differential equation is related to shear deformation. Because the load with the maneuver is vertical load, this paper only needs to solve the LM-type ordinary differential equations.

3.5. Interlayer Transfer Relationship Based on the DVP Method

In this section, the DVP method [30,31,32] is used to deduce the transfer relationship of a multilayer elastic system. Take any layer j in the layered elastic system as an example; its upper-interface coordinate is z = zj−1, and the lower-interface coordinate is z = zj. For layer j, the general solution of the LM equation in any layer can be expressed as follows:
[ U ( z ) T ( z ) ] = [ E 11 E 12 E 21 E 22 ] [ < e s 12 ( z z j ) > 0 0 < e s 34 ( z z j ) > ] [ C + C ]
In the above formula, si is the eigenvalue of the coefficient matrix in the system of ordinary differential equations of LM-type, Real (s1) ≥ Real (s2) ≥ Real (s3) ≥ Real (s4) among the eigenvalues, and [Eij] is the constant matrix composed of the eigenvectors corresponding to the eigenvalues. [C+] and [C] represent matrices of undetermined coefficients determined by the interface and boundary conditions, and the diagonal matrix is defined as follows:
< e s 12 ( z z j ) > = d i a g [ e s 1 z , e s 2 z ] < e s 34 ( z z j ) > = d i a g [ e s 3 z , e s 4 z ]
First, the coefficient of the given j layer is eliminated, and the following layer relationship is obtained:
[ U ( z j 1 ) T ( z j ) ] = [ S 11 j S 12 j S 21 j S 22 j ] [ U ( z j ) T ( z j 1 ) ]
where the matrix [Sijj] is expressed as
[ S 11 j S 12 j S 21 j S 22 j ] = [ E 11 < e s 12 h j > E 12 E 21 E 22 < e s 34 h j > ] [ E 11 E 12 < e s 34 h j > E 21 < e s 12 h j > E 22 ] 1
For layer j + 1, the following relationship can be obtained:
[ U ( z j ) T ( z j + 1 ) ] = [ S 11 j + 1 S 12 j + 1 S 21 j + 1 S 22 j + 1 ] [ U ( z j + 1 ) T ( z j ) ]
When the interface zj between layer j and layer j + 1 is a completely continuous interface, the transfer relationship between the interface zj−1 and the interface zj is as follows:
[ U ( z j 1 ) T ( z j + 1 ) ] = [ S 11 j : j + 1 S 12 j : j + 1 S 21 j : j + 1 S 22 j : j + 1 ] [ U ( z j + 1 ) T ( z j 1 ) ]
where the transfer matrix [Sijj:j+1] is
S 11 j : j + 1 = S 11 j [ I S 12 j + 1 S 21 j ] 1 S 11 j + 1 S 12 j : j + 1 = S 11 j [ I S 12 j + 1 S 21 j ] 1 S 12 j + 1 S 22 j + S 12 j S 21 j : j + 1 = S 21 j + 1 + S 22 j + 1 [ I S 12 j + 1 S 21 j ] 1 S 21 j S 11 j + 1 S 22 j : j + 1 = S 22 j + 1 [ I S 12 j + 1 S 21 j ] 1 S 22 j
When the interface is not completely continuous, the relationship between the expansion coefficients on both sides of the interface zj (zj and zj+) is as follows:
[ U ( z j ) T ( z j + ) ] = [ I α 0 I ] [ U ( z j + ) T ( z j ) ] = [ α j ] [ U ( z j + ) T ( z j ) ]
where the matrix of [a] is
[ α ] = [ α z 0 0 α r ]
That is, when the interface zj is a discontinuous interface, the transfer relation between layer j−1 and layer j can be divided into two parts. The transfer relation of the first part is transferred from layer j − 1 to the upper interface of the interface of layer j, and the coordinate of the interface is z = zj. Since the transfer relation of this part does not involve the discontinuous interface, the following relation can be obtained:
[ U ( z j 1 ) T ( z j ) ] = [ S 11 j S 12 j S 21 j S 22 j ] [ U ( z j ) T ( z j 1 ) ]
In the second part, stress and displacement are transferred on both sides of the non-completely continuous interface zj, and the transfer relationship can be obtained as follows:
[ U ( z j ) T ( z j + ) ] = [ I α 0 I ] [ U ( z j + ) T ( z j ) ]
Combining the above two formulas can be obtained as follows:
[ U ( z j 1 ) T ( z j + ) ] = [ S 11 n e w j S 12 n e w j S 21 n e w j S 22 n e w j ] [ U ( z j + ) T ( z j 1 ) ]
The new transfer matrix coefficients are as follows:
S 11 n e w j = S 11 j + S 11 j α [ I S 21 j α ] 1 S 21 j S 12 n e w j = S 12 j + S 11 j α [ I S 21 j α ] 1 S 21 j S 21 n e w j = [ I S 21 j α ] 1 S 21 j S 22 n e w j = [ I S 21 j α ] 1 S 22 j
Then, by using the known boundary conditions, the undetermined vector coefficients corresponding to any position inside the layered body are obtained. The coefficient propagates from the surface to the last layer of the interface where z = zn; we have
[ U ( z 0 ) T ( z n ) ] = [ S 11 j : n S 12 j : n S 21 j : n S 22 j : n ] [ U ( z n ) T ( z o ) ]
The boundary conditions of the surface layer and the expression of the corresponding expansion coefficient in the conversion domain were introduced in detail before, and the underlying boundary conditions will be introduced in detail below. For the underlying boundary conditions, there are two cases, as follows:
Using the boundary condition Equation (29), the vector coefficients in the surface can be expressed as
T ( z 0 ) = [ 2 π q R J 1 ( λ r ) λ , 0 ] T
When the lowest layer is a rigid base, the displacement of the bottom is 0, then the boundary condition at the bottom is U(zn) = 0, and the solution at any field point can be obtained by combining the Formula (51) as follows:
[ U ( z 0 ) T ( z n ) ] = [ S 11 1 : n S 12 1 : n S 21 1 : n S 22 1 : n ] [ 0 T ( z 0 ) ]
The lowest layer is the elastic semi-infinite space, at which the upper interface z = zn of the elastic base layer is
T ( z n ) = [ E 22 E 12 1 ] U ( z n )
Using the boundary condition z0 = 0, we have
[ U ( z 0 ) E 22 E 12 1 U ( z n ) ] = [ S 11 1 : n S 12 1 : n S 21 1 : n S 22 1 : n ] [ U ( z n ) T ( z 0 ) ]
U ( z 0 ) = { S 11 1 : n ( E 22 E 12 1 S 21 1 : n ) 1 S 22 1 : n + S 12 1 : n } T ( z 0 )
Equation (56) is the solution to the surface displacement of the above matrix.
In summary, the pavement-response solution under single-wheel vehicle random vibration load has been solved. Then, according to the vehicle random vibration load Formula (20) obtained in Section 2, the final results for the dynamic-response analysis of asphalt pavement under pavement-unevenness excitation can be obtained by superposition summation.

4. Numerical Results

4.1. Example Verification

Since the mechanical response of asphalt pavement under pavement-unevenness excitation is solved from two independent systems, the accuracy of the results in this paper will be verified from the vehicle system and the road-surface system.
In order to verify the correctness of the results of the vehicle load with maneuvering, the random load of the left-front wheel of the vehicle-simulation model established in the literature [34] was selected as a reference in this section, and the dynamic load coefficient was used as an analysis index to evaluate the vehicle load with maneuvering. The expression of the dynamic load coefficient DLC is as follows:
D L C = 1 + σ F G
In the above formula, F is the vehicle load with mobility, G is the vehicle static load, and σF is the root mean square of (FG); i.e.,
σ F = 1 n i = 1 n ( F G ) i 2
According to the vehicle parameters in the literature, the dynamic load coefficient of the vehicle driving at 16.7 m/s on the A–D grade road surface was calculated; the dynamic load coefficient results were compared with those calculated in the literature [34] and the simulation model, and the calculation results are shown in Figure 3.
As can be seen from Figure 3, the vehicle dynamic load coefficient solved by the simulation model established in this paper can achieve good consistency with the simulation results in the literature, and there is a small error, but the error of both is within the allowable range. Therefore, the correctness of the simulated vehicle loading method adopted in this chapter is verified.
In order to verify the correctness of the pavement response under dynamic load in this paper, firstly, the load is degraded into static load (ωl approaches 0), and then the three-layer pavement structure layer is calculated. Each layer of the pavement is composed of isotropic materials, the contact condition between the layers is completely continuous, and the bottom soil foundation is an elastic half-space. The pavement parameters are selected from the reference [30]. The calculated results are compared with BISAR, a pavement structure design software.
It can be clearly seen from Figure 4 that the results of displacement uz, stress σzz, and BISAR in this paper are in good agreement with slight errors, but the errors of both are within the allowable range, which verifies the correctness of the method in this paper.
After comparing the proposed method with the calculated results of BISAR, the pavement structure is compared with the analytical results of Liu Heng et al. [30] to calculate the displacements of multiple horizontal positions on the road surface, and the results are shown in Figure 5.
As can be seen from Figure 5, the results obtained in this paper are basically consistent with those obtained by Liu Heng et al. The error is less than 4%, and with the increase in horizontal distance, the error becomes smaller, which verifies the correctness of the proposed method.
The second example is selected below to further illustrate the correctness of the calculated results of this method. The random load of the vehicle is degraded into a single harmonic load, and the pavement is composed of four layers. The basic parameters of the pavement structure reference [35] are shown in Table 1. The applied load on the surface of the asphalt pavement is 0.7 MPa circular load, the applied radius R = 0.15 m, the frequency range is 5 to 65 Hz, and the increment is 5 Hz. The DVP method was used to calculate the surface displacements of four positions (x = 0.203 m, 0.457 m, 0.914 m, and 1.524 m) (without considering the time factor eiωt), and compared with the literature [35], the results are shown in Figure 6.
Figure 6 shows the displacement distribution of pavement structure at different frequencies. It can be seen from Figure 6 that the displacement results calculated by the method in this paper are basically consistent with the results in the literature, and for the displacement uz, its amplitude shows a trend of decreasing gradually with the increase of frequency; while when the loading point is far from the loading center, the displacement uz decreases with the increase of distance, which is consistent with the results in the literature, and further verifies the correctness of the method in this paper.

4.2. Comparison of Responses at Different Wheels of the Vehicle Model

The vehicle parameters are based on the real-vehicle model, as shown in Table 2. The vehicle speed is set at 20 m/s, the road level is set at Grade B, and the road-structure parameters are shown in Table 3. The superposition method is used to superposition the effects of a single wheel load. As shown in Figure 7, the response at different wheels is shown in Figure 8 and Figure 9 below.
Figure 8 shows the variation and distribution of uz and σxx on the surface of the four tires of the vehicle model. It can be observed from the figure that the mean and peak values of displacement uz and horizontal stress σxx are larger at the left- and right-front wheels than at the left- and right-rear wheels. This is because the body center of mass of the complete-vehicle model is closer to the front suspension, resulting in a greater dead load on the front wheels. At the same time, the dispersion of the left- and right-front wheels is greater, resulting in high displacement and high horizontal stress at the left- and right-front wheels. Due to the symmetry of the vehicle model, the displacement changes at the left-front wheel and the right-front wheel are similar, and their mean value and dispersion degree are basically the same.
Figure 9 reflects the variation and distribution of horizontal strain εxx and vertical strain εzz at the bottom of the surface layer of four tires of the vehicle model. It can be seen from the figure that the variation law of strain at the bottom of the surface layer is similar to that of surface displacement. The mean value and dispersion degree of horizontal strain and vertical strain at the front wheel are larger than those at the rear wheel, while the mean value and peak value at the left- and right-front wheel are basically the same, indicating that high horizontal strain and high vertical strain are more likely to occur at the bottom of the surface layer at the front wheel.
Therefore, in the following research, the action center at the left-front wheel of the vehicle model and the point below it will be selected as the research object, and the influence of pavement-unevenness grade, speed, vehicle load, interlayer contact condition and transverse isotropy on the road surface at the bottom of the surface layer will be analyzed.

4.3. Influence of Pavement-Unevenness Grade

Pavement unevenness is the cause of mechanical load, and it is divided into several grades according to the degree of pavement unevenness. The mechanical response results of the bottom layer of different pavement grades are shown in Figure 10.
Figure 10 reflects the variation and distribution of σzz, σx, εxx, and εzz at the bottom of the uneven pavement of grades A, B, C, and D. As can be seen from the figure, the peaks of σzz, σxx, εxx, and εzz at the bottom of the top layer of grade A pavement are 0.0381 MPa, 0.0115 MPa, 0.0098, and 0.0225, respectively, which increase by 20.3% compared with the level pavement surface. At the same time, these peaks increase sharply as the pavement-unevenness grade decreases, and when the road grade decreases from grade A to grade D, the εxx peaks increase by 12.3%, 51.5%, and 103.7%, respectively. σzz, σxx, εxx, and εzz distribution in the range of 25–75% and 1–99% will increase with the decline of the unevenness grade of the road surface, thus increasing its dispersion, and the probability of high stress and high strain at the bottom of the surface layer will increase. In pavement design, the horizontal strain at the bottom of the asphalt layer is usually used as the design index of pavement fatigue life, and the increase of σzz, σxx, εxx, and εzz will undoubtedly accelerate the fatigue damage of pavement.

4.4. Impact of Vehicle Speed

Speed is an important factor affecting the maneuverable load. The increase of speed will aggravate the vibration of the vehicle and greatly increase the maneuverable load of the wheel on the road surface, thus affecting the mechanical response of the road surface. Below, we explore the change rule of the mechanical response of pavement structure when the speed on B-grade pavement increases from 10 m/s to 40 m/s, and the calculation results are shown in Figure 11.
Figure 11 reflects the variation and distribution of σxx, σzz, εxx, and εzz at the bottom of the layer below different speeds. As can be seen from the figure, the increase in speed will increase the peak value of stress and strain at the bottom of the surface layer and the range of distribution between 1% and 99%, and the degree of dispersion also gradually increases, with the peak value of σxx increasing from 0.0119 MPa to 0.0145 Mpa. The peak value of σzz increases from 0.0395 Mpa to 0.0479 Mpa, εxx peak increases from 0.0101 to 0.0123, and the peak value of εzz increases from 0.0233 to 0.0283, which increases by 29.1%. The increase in speed will lead to high stress at the bottom of the surface layer, which will undoubtedly accelerate the fatigue damage of the road surface. Therefore, the road should be strictly controlled during the use of vehicles speeding to avoid premature damage to the road.

4.5. Impact of Vehicle Load

Vehicle overload is one of the important causes of road-surface damage. The influence of vehicle load on road-surface mechanical response is discussed below. The change in the vehicle’s body mass ratio mb is used to characterize the load change, with the default value of body mass mb representing full load, 0.5 times of body mass mb representing empty load, and 2.0 times of body mass mb representing overload. The mechanical response results at the bottom of the surface layer are shown in Figure 12.
Figure 12 shows the variation and distribution of σxx, σzz, εxx, and εzz at the bottom of layers with different loads. The increase in vehicle load will increase the mean value and peak value of stress and strain at the bottom of the surface layer. From no-load to overload, the peak value of σxx increases from 0.00871 MPa to 0.0217 MPa. The peak value of σzz increases from 0.0288 MPa to 0.0719 MPa, εxx peak increases from 0.00737 to 0.0184, and the peak value of εzz increases by 1.5 times from 0.0171 to 0.0425. Overload results in high stress and strain at the bottom of the surface, which affects the service life of the pavement.

4.6. Influence of Interlayer Contact Conditions

This section mainly studies the effect of different interlayer interface coefficients kr on the mechanical response of road surfaces under dynamic load. The random load generated when the vehicle speed is 20 m/s driving on the B-grade pavement is selected, the value range of kr is 107 Pa/m3 to 1012 Pa/m3, k0 = 106 Pa/m3, and kr * = kr/k0 are defined. The relevant parameters of the pavement are shown in Table 3. The discontinuous interface is located between the surface layer and the base layer, and the calculation results are shown in Figure 13.
Figure 13 shows the variation and distribution of σzz, σxx, εxx, and εzz at the bottom of the layer below the discontinuous interface, respectively. The mean and extreme values of σzz, σxx, εxx, and εzz at the bottom all decrease with the increase of the interlayer interface coefficient. As can be seen from Figure 13a,b, when the interlayer interface coefficient increases from 107 Pa/m3 to 1012 Pa/m3, the peak value of σzz decreases from 0.0466 MPa to 0.0431 MPa, decreasing by 7.5%. However, the peak value of σxx decreases from 0.0666 MPa to 0.0132 MPa, decreasing by 80.2%, and gradually approaches a fixed value at the interface coefficient of 1011 Pa/m3. The horizontal stress decreases much more than that in the vertical direction. At the same time, the dispersion of horizontal stress is much more affected by the interlayer interface coefficient than the vertical one, which indicates that the interlayer interface coefficient has a more significant effect on horizontal stress. As can be seen from Figure 13c,d, with the increase of interlayer interface coefficient, the peak value of εxx decreases from 0.0303 to 0.0111, decreasing by 63.4%, and the peak value of εzz decreases from 0.0433 to 0.0255, decreasing by 41.1%. The distribution of horizontal strain and vertical strain at the bottom of the surface layer in the range of 25–75% and 1–99% gradually decreases with the increase of interlayer interface coefficient, the degree of dispersion decreases, and the probability of high horizontal strain decreases, indicating that good interface connection conditions can slow down the fatigue damage of the pavement and delay the service life.

4.7. Effects of Transverse Isotropy

In this section, six kinds of Young’s modulus Ev in different vertical directions are selected, and the Young’s modulus Eh in horizontal direction remains unchanged (as shown in Figure 14). In addition, parameter ζ = Ev/Eh is defined, and parameter ζ equal to 1 represents isotropy, while parameter ζ is not 1 represents transverse isotropy. A speed of 20 m/s and a pavement grade of B are selected to calculate the mechanical response of the pavement with different ζ values. The results are shown in Figure 15.
Figure 15 shows the variation of σzz, σxx, εxx, and εzz with the modulus ratio of ζ at the bottom of the surface layer. It can be seen from the figure that there are significant differences between the response at the bottom of the surface layer and the change of displacement stress on the surface. In Figure 15a, the mean value and peak value of σzz increase with the increase of modulus ratio. Contrary to the change law of surface displacement, when the modulus ratio ζ increases from 0.5 to 3, its peak value increases from 0.0364 MPa to 0.0487 MPa, an increase of 33.8%. In Figure 15b, the mean value and peak value of σxx show a tendency to decrease sharply with the modulus ratio. At the same time, the degree of dispersion also decreases sharply, and the probability of high horizontal stress at a high modulus ratio is greatly reduced. It can be seen from Figure 15c,d that both the mean and peak values of εxx and εzz decrease with the increase of the modulus compared with ζ. The peak value of εxx decreases from 0.0158 to 0.00817, decreasing by 48.3%, while the peak value of εzz decreases from 0.0271 to 0.0251, decreasing by only 7.8%. The reduction of horizontal strain is much greater than that of vertical strain. At the same time, the horizontal strain εxx distribution in the range of 25%–75% and 1%–99% will decrease with the increase of the modulus ratio, and the degree of dispersion gradually becomes smaller, reducing the probability of high horizontal strain at the bottom of the surface layer, which indicates that too low of a modulus ratio will increase the probability of fatigue failure of the pavement. Therefore, suitable base and bottom base materials should be selected in pavement design to avoid this situation.

5. Conclusions

In this paper, the DVP method is used to solve the mechanical response of a multilayer elastic system under the action of vehicle loads. Firstly, starting from the pavement unevenness of the vehicle excitation source, the filtering white-noise method is used to establish the pavement-unevenness excitation, and the random load of different vehicle models is obtained by simulation. Then, based on the basic governing equation of the road-surface problem under the random load, the analytical solution of the road-surface mechanical response is obtained by using the orthogonal vector function system and the DVP method. Based on the analytical results, the influences of pavement-unevenness grade, vehicle speed, vehicle load, interlayer contact condition, and transverse isotropy on-road mechanical response were analyzed. This method has strong applicability and high efficiency and can effectively deal with the mechanical response of road-surface structure with discontinuous interface and transverse isotropy under vehicle maneuvering load, and the main conclusions are drawn as follows:
(1) Pavement unevenness does not affect the mean values of stress and strain, but it has a significant effect on the peak values and the degree of dispersion. When the pavement-unevenness level is low, it increases the peak values and coefficients of variation of stress and strain, resulting in a higher probability of high stress and strain within the pavement.
(2) The increase in vehicle speed causes the peak values and coefficients of variation of the strains to become larger, leading to frequent high strains in the pavement; in addition, the increase in vehicle loads increases the mean and peak values of the stress strains.
(3) The poor bonding between interfaces will lead to the increase of stress and strain at the bottom of the surface layer, and the increase of stress and strain in the horizontal direction is much larger than that in the vertical direction. When the base material is transverse-isotropic, the high modulus ratio will slightly increase the vertical stress and greatly reduce the horizontal stress and horizontal strain.

Author Contributions

Conceptualization: H.L. and Y.C. proposed an orthogonal vector function system and an efficient DVP method to analyze the dynamic response of asphalt pavement under pavement-unevenness excitation. Methodology: H.L. presented detailed methods and experimental procedures using numerical simulation; Software: A.W. contributed to the article by conducting numerical simulations of the experiments using relevant software; Validation: Y.C. and X.L. conducted experimental verification of the method according to the design steps; Writing—review and editing: X.L. organized and edited the initial draft of the article. All authors have read and agreed to the published version of the manuscript.

Funding

The financial support provided by the National Natural Science Foundation (Grant No. 52108424) and the Science and Technology Tackling Project of Henan Province, China (Grant No. 232102240025).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors upon request.

Conflicts of Interest

The authors are all from the School of Water Conservancy and Transportation, Zhengzhou University. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflicts of interest. The authors declare no conflicts of interest.

Abbreviations

nSpatial frequency
Gq (n0)Pavement-unevenness coefficient
qn (t)Pavement-unevenness excitation of the wheel
ω (t)White-noise function
θ ¨ b, φ ¨ Pitch-and-roll acceleration at body center of mass
s (ω)Density function
H12 (ω)Frequency response function
BWheelbase
z ˙ bVertical speed at body center of mass
z ¨ bVertical acceleration at the center of mass of the body
Ip, IrMoment of inertia of body pitch and roll
FDynamic load function of each wheel
uiDisplacement components in i-direction
εijStrain tensor
σijStress tensor
cijElastic stiffness coefficients
hiThickness of layer i
kr, kzInterface moduli in horizontal and vertical directions
RRadius of surface load circle
PlAmplitude of load
ωlAngular frequency of load
θlPhase of load
ρMaterial density
er, eθ, and ezUnit vectors of the r, θ, and z axes in the cylindrical coordinate system

References

  1. Hardy, M.S.A.; Cebon, D. Response of continuous pavements to moving dynamic loads. J. Eng. Mech. 1993, 119, 1762–1780. [Google Scholar] [CrossRef]
  2. Liu, C.; McCullough, B.F.; Oey, H.S. Response of rigid pavements due to vehicle-road interaction. J. Transp. Eng. 2000, 126, 237–242. [Google Scholar] [CrossRef]
  3. Ding, H.; Chen, L.Q.; Yang, S.P. Convergence of Galerkin truncation for dynamic response of finite beams on nonlinear foundations under a moving load. J. Sound Vib. 2012, 331, 2426–2442. [Google Scholar] [CrossRef]
  4. Chen, H.Y.; Ding, H.; Li, S.H.; Chen, L.Q. The scheme to determine the convergence term of the Galerkin method for dynamic analysis of sandwich plates on nonlinear foundations. Acta Mech. Solida Sin. 2021, 34, 1–11. [Google Scholar] [CrossRef]
  5. Chen, H.Y.; Ding, H.; Li, S.H.; Chen, L.Q. Convergent term of the Galerkin truncation for dynamic response of sandwich beams on nonlinear foundations. J. Sound Vib. 2020, 483, 115514. [Google Scholar] [CrossRef]
  6. Ma, X.; Quan, W.; Dong, Z.; Dong, Y.; Si, C. Dynamic response analysis of vehicle and asphalt pavement coupled system with the excitation of road surface unevenness. Appl. Math. Model. 2022, 104, 421–438. [Google Scholar] [CrossRef]
  7. Cebon, D. Interaction between Heavy Vehicles and Roads; SAE Technical Paper; SAE: Warrendale, PA, USA, 1993. [Google Scholar]
  8. Kim, S.M.; McCullough, B.F. Dynamic response of plate on viscous Winkler foundation to moving loads of varying amplitude. Eng. Struct. 2003, 25, 1179–1188. [Google Scholar] [CrossRef]
  9. Uzzal, R.U.A.; Bhat, R.B.; Ahmed, W. Dynamic response of a beam subjected to moving load and moving mass supported by Pasternak foundation. Shock. Vib. 2012, 19, 201–216. [Google Scholar] [CrossRef]
  10. Huang, M.H.; Thambiratnam, D.P. Dynamic response of plates on elastic foundation to moving loads. J. Eng. Mech. 2002, 128, 1016–1022. [Google Scholar] [CrossRef]
  11. Lyu, Z.; Qian, J.; Shi, Z.; Gao, Q. Dynamic responses of layered poroelastic ground under moving traffic loads considering effects of pavement roughness. Soil Dyn. Earthq. Eng. 2020, 130, 105996. [Google Scholar] [CrossRef]
  12. Zhang, J.; Zhao, Q.; Wang, X.; Ullah, S.; Zhao, D.; Civalek, O.; Xue, C.; Qi, W. New exact analytical thermal buckling solutions of composite thin plates with all edges rotationally-restrained. Mech. Adv. Mater. Struct. 2023, 31, 3518–3530. [Google Scholar] [CrossRef]
  13. Li, H.; Yang, S.; Li, S. Influence of vehicle parameters on the dynamics of pavement structure due to vehicle-road interaction. International Conference on Measuring Technoloy and Mechatronics Automation. IEEE 2009, 3, 522–525. [Google Scholar]
  14. Li, H.; Yang, S.P.; Li, S.H. Dynamical analysis of an asphalt pavement due to vehicle-road interaction. J. Vib. Shock. 2009, 28, 86–89. [Google Scholar]
  15. Lu, Z.; Yao, H. Effects of the dynamic vehicle-road interaction on the pavement vibration due to road traffic. J. Vibroengineering 2013, 15, 1291–1301. [Google Scholar]
  16. Elnashar, G.; Bhat, R.B.; Sedaghati, R. Modeling and dynamic analysis of a vehicle-flexible pavement coupled system subjected to road surface excitation. J. Mech. Sci. Technol. 2019, 33, 3115–3125. [Google Scholar] [CrossRef]
  17. Ding, H.; Yang, Y.; Chen, L.Q. Vibration of vehicle–pavement coupled system based on a Timoshenko beam on a nonlinear foundation. J. Sound Vib. 2014, 333, 6623–6636. [Google Scholar] [CrossRef]
  18. Snehasagar, G.; Krishnanunni, C.G.; Rao, B.N. Dynamics of vehicle–pavement system based on a viscoelastic Euler–Bernoulli beam model. Int. J. Pavement Eng. 2020, 21, 1669–1682. [Google Scholar] [CrossRef]
  19. Yang, S.; Li, S.; Lu, Y. Investigation on dynamical interaction between a heavy vehicle and road pavement. Veh. Syst. Dyn. 2010, 48, 923–944. [Google Scholar] [CrossRef]
  20. Li, S.; Yang, S.; Chen, L. A nonlinear vehicle-road coupled model for dynamics research. J. Comput. Nonlinear Dyn. 2013, 8, 021001. [Google Scholar] [CrossRef]
  21. Krishnanunni, C.G.; Rao, B.N. Decoupled technique for dynamic response of vehicle-pavement systems. Eng. Struct. 2019, 191, 264–279. [Google Scholar] [CrossRef]
  22. Kropac, O.; Mucka, P. Classification scheme for random longitudinal road unevenness considering road waviness and vehicle response. J. Shock. Vib. 2009, 16, 935858. [Google Scholar]
  23. Li, J.H.; He, J.; Li, X.H. Dynamic response of pavement under vehicle random load and moving constant load. J. Chang. Univ. 2015, 35, 38–45. [Google Scholar]
  24. Xing, Y.; Liu, Z. Dynamic response analysis of semi-rigid pavement under random load of vehicle. Dev. Build. Technol. 2018, 45, 80–82. [Google Scholar]
  25. Zhao, J.; Wang, H. Dynamic pavement response analysis under moving truck loads with random amplitudes. J. Transp. Eng. Part B Pavements 2020, 146, 04020020. [Google Scholar] [CrossRef]
  26. Hettiarachchi, C.; Yuan, J.; Amirkhanian, S.; Xiao, F. Measurement of pavement unevenness and evaluation through the IRI parameter—An overview. J. Meas. 2023, 206, 112284. [Google Scholar] [CrossRef]
  27. Wang, H.; Zhao, J.; Hu, X.; Zhang, X. Flexible pavement response analysis under dynamic loading at different vehicle speeds and pavement surface roughness conditions. J. Transp. Eng. Part B Pavements 2020, 146, 04020040. [Google Scholar] [CrossRef]
  28. Liu, D.; Li, G.; Chen, H.; Kang, Z.; Jiang, R. Dynamic response of pavement based on random dynamic load of vehicle. Trans. Chin. Soc. Agric. Mach. 2011, 24, 28–33. [Google Scholar]
  29. Zhang, P.; Liu, J.; Lin, G.; Wang, W. Elastic displacement fields of multi-layered transversely isotropic materials under rectangular loads. Eur. J. Environ. Civ. Eng. 2018, 22, 1060–1088. [Google Scholar] [CrossRef]
  30. Liu, H.; Pan, E.; Cai, Y. General surface loading over layered transversely isotropic pavements with imperfect interfaces. Adv. Eng. Softw. 2018, 115, 268–282. [Google Scholar] [CrossRef]
  31. Liu, H.; Pan, E. Time-harmonic loading over transversely isotropic and layered elastic half-spaces with imperfect interfaces. Soil Dyn. Earthq. Eng. 2018, 107, 35–47. [Google Scholar] [CrossRef]
  32. Liu, H.; Cai, Y.; Zhong, Y.; Pan, B. Thermo-hydro-mechanical response of a multi-layered pavement with imperfect interface based on dual variable and position method. Appl. Math. Model. 2021, 99, 704–729. [Google Scholar] [CrossRef]
  33. Si, C.; Cao, H.; Chen, E.; You, Z.; Tian, R.; Zhang, R.; Gao, J. Dynamic Response Analysis of Rutting Resistance Performance of High Modulus Asphalt Concrete Pavement. Appl. Sci. 2018, 8, 2701. [Google Scholar] [CrossRef]
  34. Shi, X.H.; Jiang, X.; Zhao, J. Response characteristics of a seven-degree-of-freedom vehicle under four-wheel random road excitation. Technol. Eng. 2018, 18, 71–78. [Google Scholar]
  35. Pan, E.; Lin, C.P.; Zhou, J. Fundamental solution of general time-harmonic loading over a transversely isotropic, elastic and layered half-space: An efficient and accurate approach. Eng. Anal. Bound. Elem. 2021, 132, 309–332. [Google Scholar] [CrossRef]
  36. Jiang, P. Simulation Analysis and Parameter Optimization Design of Automobile Suspension System; Zhejiang University: Hangzhou, China, 2006. [Google Scholar]
Figure 1. Vehicle model.
Figure 1. Vehicle model.
Applsci 14 08822 g001
Figure 2. Schematic diagram of N-layer elastic pavement structure.
Figure 2. Schematic diagram of N-layer elastic pavement structure.
Applsci 14 08822 g002
Figure 3. Comparison of the simulation model in this paper with simulation results in the literature.
Figure 3. Comparison of the simulation model in this paper with simulation results in the literature.
Applsci 14 08822 g003
Figure 4. Comparison of results between this paper and BISAR. (a) surface uz; (b) surface σzz.
Figure 4. Comparison of results between this paper and BISAR. (a) surface uz; (b) surface σzz.
Applsci 14 08822 g004
Figure 5. Comparison between the calculated results in this paper and those in the reference [30].
Figure 5. Comparison between the calculated results in this paper and those in the reference [30].
Applsci 14 08822 g005
Figure 6. Comparison between the calculated results in this paper and those in the reference [35].
Figure 6. Comparison between the calculated results in this paper and those in the reference [35].
Applsci 14 08822 g006
Figure 7. Schematic diagram of wheel-track action of the vehicle model.
Figure 7. Schematic diagram of wheel-track action of the vehicle model.
Applsci 14 08822 g007
Figure 8. Surface mechanical response of surface layers at different wheels of the vehicle model. (a) Surface uz; (b) surface σxx.
Figure 8. Surface mechanical response of surface layers at different wheels of the vehicle model. (a) Surface uz; (b) surface σxx.
Applsci 14 08822 g008
Figure 9. Mechanical response at the bottom of the surface layer at different wheels of the vehicle model. (a) The bottom of the surface layer εzz; (b) the bottom of the surface layer εxx.
Figure 9. Mechanical response at the bottom of the surface layer at different wheels of the vehicle model. (a) The bottom of the surface layer εzz; (b) the bottom of the surface layer εxx.
Applsci 14 08822 g009
Figure 10. Influence of pavement-unevenness grade on the mechanical response of surface bottom. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz. The pavement roughness coefficient Gq (n0) was selected according to different pavement roughness grades [33] A, B, C, and D, as shown in Table 4.
Figure 10. Influence of pavement-unevenness grade on the mechanical response of surface bottom. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz. The pavement roughness coefficient Gq (n0) was selected according to different pavement roughness grades [33] A, B, C, and D, as shown in Table 4.
Applsci 14 08822 g010
Figure 11. Influence of speed on strain at bottom of surface layer. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Figure 11. Influence of speed on strain at bottom of surface layer. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Applsci 14 08822 g011
Figure 12. Influence of vehicle load on mechanical response at bottom of surface layer. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Figure 12. Influence of vehicle load on mechanical response at bottom of surface layer. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Applsci 14 08822 g012
Figure 13. Influence of interface coefficient on mechanical response at bottom of surface layer. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Figure 13. Influence of interface coefficient on mechanical response at bottom of surface layer. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Applsci 14 08822 g013
Figure 14. Schematic diagram of isotropic base and base in transverse view.
Figure 14. Schematic diagram of isotropic base and base in transverse view.
Applsci 14 08822 g014
Figure 15. Influence of transverse isotropy on mechanical response of surface bottom. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Figure 15. Influence of transverse isotropy on mechanical response of surface bottom. (a) Bottom of the surface layer σzz; (b) bottom of the surface layer σxx; (c) bottom of the surface layer εxx; (d) bottom of the surface layer εzz.
Applsci 14 08822 g015
Table 1. Pavement structure parameters.
Table 1. Pavement structure parameters.
Pavement Structural LayerEv (MPa)Eh (MPa)vvvhThickness (mm)ρ (kg/m3)
Asphalt surface14007000.350.351502400
Basic level5001500.40.43602200
Subbase course200700.40.41802000
Soil matrix60600.450.451600
Table 2. Vehicle model parameters [36].
Table 2. Vehicle model parameters [36].
Real-Vehicle ParameterNumerical ValueReal-Vehicle ParameterNumerical Value
mb/kg1380cs1, cs2/kN·s·m−11.5
m1, m2/kg40.5cs3, cs4/kN·s·m−11.5
m4, m3/kg45.4ks1, ks2/kN·m−117
Ip/kg·m22444ks3, ks4/kN·m−122
Ir/kg·m2380Bf, Br/m1.48
kt1, kt2/kN·m−1192a/m1.25
kt3, kt4/kN·m−1192b/m1.51
Table 3. Pavement structure parameters.
Table 3. Pavement structure parameters.
Pavement StructureLayer Thickness (m)Modulus of Elasticity (MPa)Poisson’s RatioDensity ρ (kg/m3)
Surface course0.1520000.32.1 × 103
Basic level0.3510000.32.1 × 103
Subbase course0.25000.352.0 × 103
Soil matrix1000.41.9 × 103
Table 4. The pavement roughness coefficient Gq (n0) for different roughness grades.
Table 4. The pavement roughness coefficient Gq (n0) for different roughness grades.
Roughness GradeLower LimitGq (n0) (10−6m3)
Average Value
Upper Limit
A81636
B3264128
C128256512
D51210242048
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Liu, H.; Liu, X.; Wei, A.; Cai, Y. Dynamic Response Analysis of Asphalt Pavement under Pavement-Unevenness Excitation. Appl. Sci. 2024, 14, 8822. https://doi.org/10.3390/app14198822

AMA Style

Liu H, Liu X, Wei A, Cai Y. Dynamic Response Analysis of Asphalt Pavement under Pavement-Unevenness Excitation. Applied Sciences. 2024; 14(19):8822. https://doi.org/10.3390/app14198822

Chicago/Turabian Style

Liu, Heng, Xiaoge Liu, Ankang Wei, and Yingchun Cai. 2024. "Dynamic Response Analysis of Asphalt Pavement under Pavement-Unevenness Excitation" Applied Sciences 14, no. 19: 8822. https://doi.org/10.3390/app14198822

APA Style

Liu, H., Liu, X., Wei, A., & Cai, Y. (2024). Dynamic Response Analysis of Asphalt Pavement under Pavement-Unevenness Excitation. Applied Sciences, 14(19), 8822. https://doi.org/10.3390/app14198822

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