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:
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].
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
The pavement-unevenness excitation of the front and rear wheels has the following relationship:
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:
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:
The time-domain function
q1(
t) of front-wheel pavement unevenness is substituted into the above formula to obtain the following:
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:
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:
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
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:
The left-front-wheel pavement-unevenness time-domain function
q1(
t) is substituted into the above formula thus:
According to the correlation between input
q2(
t) and
q4(
t) of the pavement unevenness of the right-front and -rear wheels,
By substituting the right-front-wheel pavement-unevenness time-domain function
q2(
t) into the above Equation (13), we can obtain
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:
Body-pitch-motion equation:
Body roll motion equation:
Equation of vertical motion of four-wheel masses:
The relationship between suspension displacement and body-mass center displacement is as follows:
In the upper formula: 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; 1, 2, 3, 4—vertical acceleration of the front and rear wheels on the left and right sides, m/s2; b—vertical speed at body center of mass, m/s; b1, b2, b3, b4—vertical speed of the left- and right-side suspension at the body position, m/s; 1, 2, 3, 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.
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 (
) [
34]; for example
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 .
The dynamic balance differential equation, without considering the physical space problem, is as follows:
The dynamic balance differential Equation (22) can be rewritten as
where
ρ is the density, and
ωl is the angular frequency.
According to Hooke’s law, the equation of physics is as follows:
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:
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:
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:
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:
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:
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:
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
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:
The expression of stress in the transformation domain:
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
For the LM-type, there are
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:
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:
First, the coefficient of the given
j layer is eliminated, and the following layer relationship is obtained:
where the matrix [
Sijj] is expressed as
For layer
j + 1, the following relationship can be obtained:
When the interface
zj between layer
j and layer
j + 1 is a completely continuous interface, the transfer relationship between the interface z
j−1 and the interface z
j is as follows:
where the transfer matrix [S
ijj:j+1] is
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:
where the matrix of [a] is
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:
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:
Combining the above two formulas can be obtained as follows:
The new transfer matrix coefficients are as follows:
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
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
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:
The lowest layer is the elastic semi-infinite space, at which the upper interface
z =
zn of the elastic base layer is
Using the boundary condition
z0 = 0, we have
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.
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.