# Prediction of Viscoelastic Pavement Responses under Moving Load and Nonuniform Tire Contact Stresses Using 2.5-D Finite Element Method

Authors: Chaoyang Wu; Hao Wang; Jingnan Zhao; Xin Jiang; Qiu Yanjun; Bekhzad Yusupov
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This study developed two-and-half dimensional (2.5-D) finite element method (FEM) to predict viscoelastic pavement responses under moving loads and nonuniform tire contact stresses. The accuracy of 2.5-D FEM was validated with two analytical solutions for elastic and viscoelastic conditions. Compared to three-dimensional (3-D) FEM, the computational efficiency of the 2.5-D method was greatly improved. The effects of loading pattern and speed on pavement surface deflection and strain responses were analyzed for asphalt pavements with four different asphalt layer thicknesses. The analyzed pavement responses included surface deflections, maximum tensile strains in the asphalt layer, and maximum compressive strains on top of subgrade. The loading patterns have influence on the mechanical responses. According to the equivalent rule, the point load, rectangle type, and sinusoid-shape contact stresses were studied. It was found that the point load caused much greater pavement responses than that of the area-based loading. When the tire loading was simplified as uniform contact stress in rectangular area, the maximum tensile strains in the asphalt layer varied with the width/length ratio of contact area. Additionally, it was shown that the dynamic responses of pavement structure induced by the sinusoid-shape contact stresses and realistic nonuniform stresses were quite similar to each other in all the cases. The pavement strain responses decreased as the speed increased due to viscoelastic behavior of asphalt layer. The study results indicate that asphalt pavement responses under moving load can be calculated using the proposed 2.5-D FEM in a fast manner for mechanistic-empirical pavement design and analysis.


---

## Body

## 1. Introduction

It is important to calculate pavement responses under traffic loading for pavement design and analysis. Two primary approaches, analytical solution and finite element modeling, have been used for calculation of pavement responses. The Boussinesq solution is first derived to solve stress and strain in the half-space system under vertical force [1]. The development of Burmister theory enables the consideration of multilayer system in pavement analysis [2]. Recently, many software applications have been developed based on multilayer elastic layer theory for pavement structural analysis, such as BISAR, KENLAYER, and WESLEA [1–5]. In these early analytical methods, material linearity, isotropy, stationary loading, and no-slip interface between layers are usually assumed in order to solve the governing equations of pavement system with multiple layers.Due to viscoelastic nature of asphalt material, it is well known that mechanical responses of asphalt pavement under moving vehicular loading are affected by tire-pavement interaction and speed [6, 7]. Analytical solutions of viscoelastic pavement responses vary with different levels of complexity depending on the assumptions of pavement structure (such as finite beam, infinite plate on Winkler foundation, or multiple layers) and moving loads (constant, harmonic, or random loads) [3, 4, 8–12]. The solutions can be closed-form solutions in analytical form based on the corresponding principle or semianalytical form solved using numerical techniques.Although analytical solutions are fast and efficient to calculate pavement responses, the scope of its application and analysis accuracy are limited due to adopting many simplified assumptions (such as full bonded layer interface, linear and isotropic material properties, etc.). In addition, the analytical solutions may not exist for complex problems and cannot be customized for most cases. On the other hand, FEM has been widely used by researchers for pavement analysis, including axisymmetric FEM [13] and 3-D FEM [14–16]. Although the 3-D FEM model can be a complex and costly analysis tool, it provides the needed versatility and flexibility to accurately simulate nonlinear material behavior, complex layer interface condition, and nonuniform distribution of tire loading. Compared to analytical solutions, FEM has strong high versatility and can better adapt to different conditions of complex pavement systems [7, 17, 18]. However, due to the nonuniform contact stress distribution at tire-pavement interface, the meshes of pavement models need to be refined to capture the localized stress patterns corresponding to moving loads. Additionally, the time step needs to be divided into very small increments in order to obtain the converged solutions of dynamic problems. Therefore, the calculation process of FEM costs much longer time than analytical solutions.The semianalytical solution is an approach combined with the advantage of both accuracy and fast calculation speed. Siddharthan et al. used a continuum, finite-layer model to evaluate pavement responses under moving surface load [19]. The moving load was decomposed into multiple single harmonic pressure distributions, and viscoelastic layer properties were defined by dynamic shear modulus and internal damping ratio. The model was later used to develop the 3D-MOVE pavement analysis software. Taiebat and Carter proposed a semianalytical finite element method to investigate the laterally loaded pile problem in consolidating elastic-plastic soil [20]. Chabot et al. developed the ViscoRoute 2.0 software based on Duhamel’s semianalytical multilayer model to calculate pavement responses under moving load within rectangular or elliptic contact area [21]. Lee developed ViscoWave based on the Laplace and Hankel transform to solve viscoelastic pavement responses under impulsive loading [22]. On the other hand, Eslaminia and Guddati developed the Fourier finite element (FFE) method to solve pavement problems from elastic domain to viscoelastic domain with Prony series as representation of viscoelastic material [23]. Liu et al. applied the Fourier transformation of space along the direction perpendicular to traffic direction for calculating pavement responses under moving loads and found the agreements between the results from the semifinite element method (SFEM) and the ones obtained using ABAQUS [24, 25].The 2.5-D FEM is based on Fourier transformation of space along the traveling direction in wavenumber domain as well as the time in frequency domain at the same time. The roadway has very long length while its width is relatively small. Thus, it is possible to simplify the 3-D problem to be 2.5-D case by Fourier transformation along the traffic direction. Compared to 3-D FEM, 2.5-D FEM only needs mesh discretization in transverse section. Therefore, its computational efficiency is greatly improved. Previous works have applied 2.5-D FEM for solving the problems related to high-speed railway and soil foundation [26–30]. In particular, Yang and Hung applied the 2.5-D method to analyze environmental vibration control for high-speed railway [26]. Bian et al. discussed the critical speed of track-embankment-ground system using 2.5-D FEM [30]. Gao et al. investigated dynamic responses of rail tracks and saturated porous ground subjected to moving loads caused by high-speed trains [31]. However, these research studies focused on rail track and ground behavior with the aim to reduce vibration caused by high-speed train. Few research studies have utilized 2.5-D FEM for pavement applications.In this study, the mathematical principle of 2.5-D FEM was implemented for pavement response analysis subject to any arbitrary pattern of moving tire loading. The predicted pavement deflections and strain responses were validated with analytical solutions. The effects of tire-pavement contact stresses and speed on pavement responses at different pavement structures were investigated. This study further developed 2.5-D FEM and brought new contributions in the following two aspects. First, the traditional 2.5-D FEM used elastic materials, but in this duty linear viscoelastic material was considered for asphalt layer. This is important for calculating asphalt pavement responses under moving load. Secondly, the arbitrary shape of loads (such as nonuniform contact stress distribution) at tire-pavement interface was considered in the model.


## 2. Theory and Background

Figure1 shows the typical flow chart that illustrates the solving process implemented in 2.5-D FEM. The solving process of this method started from presenting governing equation of elastic or viscoelastic and the relationship of strain and stress in time domain of 3-D FEM. Then, the correspondence form in wavenumber domain and frequency domain was obtained. In the next step, the solutions in wavenumber domain and frequency domain were determined. Finally, the results were transformed from frequency domain to time domain.Figure 1
Flowchart of solving process in 2.5-D FEM.


### 2.1. Governing Equations

Several assumptions are made to derive the governing equations in 2.5-D FEM model. The material is assumed with isotropic properties. The dimension of pavement structure is assumed to be infinite along the traffic direction. Pavement layers are assumed to be fully bonded at the interface. The asphalt layer can be modeled as elastic material or viscoelastic material, while the other layers are modeled as elastic material. The model is based on small strain assumption.The illustration of 2.5-D FEM model converted from 3-D model is shown in Figure2. The Fourier transform and its inverse transform are predefined in equation (1) and (2):(1)Fstξx,y,z,ω=∫−∞+∞∫−∞+∞Fx,y,z,teiξxxe−iωtdxdt,(2)Fx,y,z,t=14π2∫−∞+∞∫−∞+∞Fstξx,y,z,ωe−iξxxeiωtdξxdω,where F is a general function which represents the mechanical response in pavement structure; the superscript “s” indicated the mechanical response in the wavenumber domain; the superscript “t” indicates mechanical response in the frequency domain; the superscript “st” indicates mechanical response in the wavenumber domain and the frequency domain; ξx represents wavenumber along the x axis (traffic direction); and w represents angle frequency.Figure 2
Illustration of transformation from 3-D model to 2.5-D model. (a) 3-D model. (b) 2.5-D model.


#### 2.1.1. Elastic Governing Equations

The governing equation of elastic ground is shown in equation (3) in tensor form:(3)σij=2μεij+λδijεkk,where λ and μ are Lame constants; εiji,j=1,2,3 is strain in the ground; σiji,j=1,2,3 is the stress in the ground; and δiji,j=1,2,3 is Kronecker delta function.The equilibrium equation for elastic layers is formulated by the following expressions:(4)σij,j+fi=ρu¨i,(5)εij=12ui,j+uj,i,where ρ is the density; uii=1,2,3 is displacement vector; and fi=i=1,2,3 is body force.By substituting equations (4) and (5) into equation (3), the Navier elastic equation can be obtained, as shown in the following equation:(6)μui,jj+λ+μuj,ji+fi=ρu¨i.Navier’s equation in the frequency domain and the wavenumber domain can be described as follows:(7)μcμi,jjt+λc+μcμj,jit+ω2ρμit+fit=0,where λc and μc are Lame constants can be expressed by the following equation:(8)λc=λ1+2βi,μc=1+2βi,where β represents damp ratio and i is imaginary unit.Based on small strain theory, the strain-displacement relationship can be expressed as follows:(9)ε=εxxεyyεzzγxyγyzγzx=∂∂x000∂∂y000∂∂z∂∂y∂∂x00∂∂z∂∂y∂∂z0∂∂x=uvw=Bu.By applying the predefined Fourier transformation in equation (1) to (9), the expression can be yielded in the frequency domain and the wave number domain. The strain-displacement relationship is given by(10)εst=εxxstεyystεzzstγxystγyzstγzxst=−iζx000∂∂y000∂∂z∂∂y−iζx00∂∂z∂∂y∂∂z0−iζxustvstwst=Bstust.The stress-strain relationship is given by(11)σst=Dεst.The model can be discretized into several quadrilateral elements, and each element node has three degrees of freedom, as shown in equations (12)–(18):(12)Niη,ζ=141+ηηi1+ζζi,i=1,2,3,4,(13)Niηj,ζj=1,i=j,0,i≠j,(14)Uy,z=∑j=14Njy,zdj=NU,(15)Kst−ω2MUst=Fst,(16)M=∑eρ∫−11∫−11NTNJdηdζ,(17)Kst=∑eρ∫−11∫−11B∗stNTDBstNJdηdζ,(18)Fst=∑eρ∫−11∫−11NTfJdηdζ,where J is determinant of the matrix; N is the matrix of shape function; and D is stiffness matrix.


#### 2.1.2. Viscoelastic Governing Equations

The Fourier transform is a linear transform, so this rule could be applied to all linear equations. Prony series has been widely used to capture the modulus change of linear viscoelastic materials with reduced time. For viscoelastic case, the equilibrium equation is given as(19)BTσ=f,where f is the body force vector.The strain-displacement relationship is the same as the one shown in equation (9).For stress and strain relationship, the viscoelastic model is different from elastic equation, as shown in equation (20) [32]. For stress and strain relationship, the viscoelastic model is different from elastic equation, as shown in the following equation [32]:(20)σ=∫0tDt−τdεdτdτ,where Dt−τ is stiffness history matrix.By applying Fourier transformation, the stress and strain in the frequency domain are given by(21)σt=iωDtiωεt,where Dtiω is stiffness history matrix in frequency domain.The loadP moves on the top of surface layer at a speed of v, and the surface boundary condition is expressed by(22)tr=σzzτyzτxyT=Py,t−xv.The property of material is assumed to be isotropic. By applying the energy variation principle to equation (19), equation (23) can be obtained:(23)δUω,ξx=∫Ωσstδεstdydz−∫yminymaxtrstδusty,zz=0dy=0,where U is the total potential energy of the system and δεst is virtual strain in wavenumber domain and frequency domain.The next step is to follow the general finite element method which has been performed in equation (12)–(18). However, the stiffness matrix is different from the one using elastic theory, as shown in the following equation:(24)Kst=∑eρ∫−11∫−11B∗NTDtiωBNJdηdζ,(25)Fst=∫yminymaxNy,zz=0pstdy.The Prony series as a general form of Maxwell model for linear viscoelasticity is used here, as shown in equation (25).(26)E∗ω=E∞+∑i=1NiωτiEi1+iωτi,where E∞ is long-term equilibrium relaxation modulus; E∗ω is dynamic modulus; τi and Ei are the parameters of Prony Series; τi is the reduced time; and Ei is the elastic component.The relationship between the shift factor and the temperature can be approximated by Williams–Landel–Ferry (WLF) function.(27)logaT=−C1T−T0C2+T−T0,where T0 is reference temperature (20°C in this study); aT is actual temperature corresponding to the shift factor with τ0=τ/aT; and C1 and C2 are regression parameters.


### 2.2. Arbitrary Shape of Moving Force

According to Fourier transformation, the dynamic vehicle load could be decoupled into a series of simple harmonic motions. For the sake of simplicity, only one self-oscillation frequency was considered in this case. It is assumed that the vehicle moves at speed ofv and with width of Ly along y direction (perpendicular to traffic direction) on pavement surface. The moving load of the vehicle can be expressed in the following equation:(28)fx,y,z,t=Φx−vtδy−Lyeiω0t,where Φx−vt is the load distribution along x direction and ω0 is the self-oscillation of load.The load in the frequency and wavenumber domain is shown in the following equation:(29)fstξx,y,z,ω=2πvδξx−ω−ω0vΦsξx.A typical track with five axles is employed in our case, as shown in Figure3. The track moves along the x direction at a constant speed of v and with self-oscillation frequency ω0, and the equation of moving load can be presented as shown in the following equation:(30)fx,y,z,t=p1δx−vt+p2δx−vt+L1+p3δx−vt+L1+L2+p4δx−vt+L1+L2+L3+p5δx−vt+L1+L2+L3+L4eiω0tΦx−vt,where P1,P2,P3,P4,andP5 represents the axle load of truck; L1,L2,L3,andL4 represents the distance between difference axles; and δ is the Dirac function.Figure 3
Example of typical truck with five axle loads.By employing the predefined Fourier transformation, the load in the frequency and wavenumber domain can be given in the following equation:(31)fstξx,y,z,ω=2πvδξx−ω−ω0vp1+p2∗e−iL1ξx+p3∗e−iL1+L2ξx+p4∗e−iL1+L2+L3ξx+p5∗e−iL1+L2+L3+L4ξxΦsξx.According to the definition of Dirac function, equation (32) can be obtained:(32)ξx=ω−ω0v.This means the wavenumber has specific relationship with the frequency. Then, the double dimensional Fourier transform can be simplified as one-dimensional Fourier transform, and the computational cost could be reduced significantly.Tire-pavement interaction can be modeled as different load patterns, such as point load and area-based load with different distribution patterns, as shown in Figure4.Figure 4
Illustration of different loading patterns: uniform stress in rectangular area, sinusoid-shape stress distribution, and point load. (a) Rectangle; (b) sinusoid; (c) point load.
(a)
(b)
(c)The load of distribution along traffic direction could be transformed in the wavenumber domain, as shown in equation (33)–(35).(33)ΦRx=b,x≤a,0,x>a,ΦRsξx=2∗b∗sinaξxξx,(34)ΦSx=qpeaksinπ∗x2L+π2,x≤L,0,x>L,ΦSsξx=qpeak∗sinπ/2L+ξx∗Lπ/2L+ξx+sinξx−π/2L∗Lξx−π/2L,(35)Φpsξx=2∗π∗Fv,where the subscripts r, s, and p represent rectangle shape, sinusoidal, and point load, respectively.In 3-D finite element method, all the external and body forces need be converted to equivalent nodal forces. However, in the 2.5-D method, the load was solved in frequency domain, and most of the time the distribution of regular loads had analytical forms after Fourier transformation. In this case, the analytical solution could be used directly. There are also some cases that where the distribution functions of loads along the traffic direction have no analytical solution after Fourier transformations. In this case, the method of equivalent area can be used for approximate solutions, as shown in Figure5. The loads with the elliptical or circular area are equivalent to the loads distributed in multiple rectangle-shaped areas divided from the elliptical or circular area. For exact solutions, the node forces can be calculated and treated as multiple point loads. Then, the ideal result can be acquired using equation (31). This makes the implementation of arbitrary-shape load in the 2.5-D FEM model possible.Figure 5
Illustration of ellipse and circle shapes that can be equivalent to rectangular distributed shapes.
(a)
(b)


### 2.3. Mesh and Boundary Condition

Due to the orthogonal properties of Fourier transformation, the 3-D FE model of pavement was decomposed to the 2.5-D FE model in frequency domain, while each element node in the 2.5-D FE model had movements in three dimensions. The nodes were unrestricted to move at the traffic direction. The nodes at the bottom of pavement model were fixed in this model. In order to obtain more accurate result in dynamic problem. According to previous research results, the maximum element size, minimum element size, and model size should be decided by the wavelength [33].In the dynamic problem, wave reflection will occur from the boundary, which will influence the results. The boundary condition should be selected properly. In this model, pavement structure lays on unbounded soil medium. Due to the limited size of the finite element model, it cannot handle this boundary condition directly. Several methods have been used in the literature to simulate the infinite foundation such as boundary element [34], absorbing boundary [27], thin-layer element [30], and infinite element [26, 35]. Among all these methods, the absorbing boundary showed good accuracy, and it is easy to be used in the practical application [27].Figure6 illustrates the mesh and boundary condition of the 2.5-D FEM model. In general, the minimum element size is 0.01 m by 0.0254 m at the pavement near surface, and the maximum element size is 0.5 m by 0.5 m in the subgrade. The model size was 6 m by 6 m. The left and right boundaries were absorbing boundaries, and the bottom boundary was fixed. The absorbing boundary depends on material properties and vibration frequency, as shown in the following equation:(36)ffj=kiu˙j,j=x,y,z,where ff is node force and u˙ is deformation velocity.Figure 6
General model of pavement structure and loaded area.This could be transformed by the predefined Fourier transformation in the frequency domain, as shown in equations (37)–(39).(37)ffjt=−iωkjujt,j=x,y,z,in which(38)kx=ky=ρvs,kz=ρvp,where vs is velocity of shear waves and vp is velocity of compression waves, and(39)vs=μρ,vp=λ+2μρ.The absorbing boundary can be assembled to stiffness matrix either to equation (17) (elastic material) or equation (24) (viscoelastic material).


### 2.4. Strain, Stress, Velocity, and Acceleration

According to the basic theory of the finite element method, the strain, stress, velocity, and acceleration can be calculated by the node displacements. Solving equation (15), the node displacement can be obtained. The strain, stress, velocity, and acceleration in the frequency domain can be obtained through equations (40)–(42), respectively.(40)σst=Dεst=DBstNuest,where uest is a vector concluding the node displacement of element.(41)Vst=iωμst,where Vst is the node velocity of any element.(42)ast=−ω2μst,where ast is the node acceleration of any element.


## 2.1. Governing Equations

Several assumptions are made to derive the governing equations in 2.5-D FEM model. The material is assumed with isotropic properties. The dimension of pavement structure is assumed to be infinite along the traffic direction. Pavement layers are assumed to be fully bonded at the interface. The asphalt layer can be modeled as elastic material or viscoelastic material, while the other layers are modeled as elastic material. The model is based on small strain assumption.The illustration of 2.5-D FEM model converted from 3-D model is shown in Figure2. The Fourier transform and its inverse transform are predefined in equation (1) and (2):(1)Fstξx,y,z,ω=∫−∞+∞∫−∞+∞Fx,y,z,teiξxxe−iωtdxdt,(2)Fx,y,z,t=14π2∫−∞+∞∫−∞+∞Fstξx,y,z,ωe−iξxxeiωtdξxdω,where F is a general function which represents the mechanical response in pavement structure; the superscript “s” indicated the mechanical response in the wavenumber domain; the superscript “t” indicates mechanical response in the frequency domain; the superscript “st” indicates mechanical response in the wavenumber domain and the frequency domain; ξx represents wavenumber along the x axis (traffic direction); and w represents angle frequency.Figure 2
Illustration of transformation from 3-D model to 2.5-D model. (a) 3-D model. (b) 2.5-D model.


### 2.1.1. Elastic Governing Equations

The governing equation of elastic ground is shown in equation (3) in tensor form:(3)σij=2μεij+λδijεkk,where λ and μ are Lame constants; εiji,j=1,2,3 is strain in the ground; σiji,j=1,2,3 is the stress in the ground; and δiji,j=1,2,3 is Kronecker delta function.The equilibrium equation for elastic layers is formulated by the following expressions:(4)σij,j+fi=ρu¨i,(5)εij=12ui,j+uj,i,where ρ is the density; uii=1,2,3 is displacement vector; and fi=i=1,2,3 is body force.By substituting equations (4) and (5) into equation (3), the Navier elastic equation can be obtained, as shown in the following equation:(6)μui,jj+λ+μuj,ji+fi=ρu¨i.Navier’s equation in the frequency domain and the wavenumber domain can be described as follows:(7)μcμi,jjt+λc+μcμj,jit+ω2ρμit+fit=0,where λc and μc are Lame constants can be expressed by the following equation:(8)λc=λ1+2βi,μc=1+2βi,where β represents damp ratio and i is imaginary unit.Based on small strain theory, the strain-displacement relationship can be expressed as follows:(9)ε=εxxεyyεzzγxyγyzγzx=∂∂x000∂∂y000∂∂z∂∂y∂∂x00∂∂z∂∂y∂∂z0∂∂x=uvw=Bu.By applying the predefined Fourier transformation in equation (1) to (9), the expression can be yielded in the frequency domain and the wave number domain. The strain-displacement relationship is given by(10)εst=εxxstεyystεzzstγxystγyzstγzxst=−iζx000∂∂y000∂∂z∂∂y−iζx00∂∂z∂∂y∂∂z0−iζxustvstwst=Bstust.The stress-strain relationship is given by(11)σst=Dεst.The model can be discretized into several quadrilateral elements, and each element node has three degrees of freedom, as shown in equations (12)–(18):(12)Niη,ζ=141+ηηi1+ζζi,i=1,2,3,4,(13)Niηj,ζj=1,i=j,0,i≠j,(14)Uy,z=∑j=14Njy,zdj=NU,(15)Kst−ω2MUst=Fst,(16)M=∑eρ∫−11∫−11NTNJdηdζ,(17)Kst=∑eρ∫−11∫−11B∗stNTDBstNJdηdζ,(18)Fst=∑eρ∫−11∫−11NTfJdηdζ,where J is determinant of the matrix; N is the matrix of shape function; and D is stiffness matrix.


### 2.1.2. Viscoelastic Governing Equations

The Fourier transform is a linear transform, so this rule could be applied to all linear equations. Prony series has been widely used to capture the modulus change of linear viscoelastic materials with reduced time. For viscoelastic case, the equilibrium equation is given as(19)BTσ=f,where f is the body force vector.The strain-displacement relationship is the same as the one shown in equation (9).For stress and strain relationship, the viscoelastic model is different from elastic equation, as shown in equation (20) [32]. For stress and strain relationship, the viscoelastic model is different from elastic equation, as shown in the following equation [32]:(20)σ=∫0tDt−τdεdτdτ,where Dt−τ is stiffness history matrix.By applying Fourier transformation, the stress and strain in the frequency domain are given by(21)σt=iωDtiωεt,where Dtiω is stiffness history matrix in frequency domain.The loadP moves on the top of surface layer at a speed of v, and the surface boundary condition is expressed by(22)tr=σzzτyzτxyT=Py,t−xv.The property of material is assumed to be isotropic. By applying the energy variation principle to equation (19), equation (23) can be obtained:(23)δUω,ξx=∫Ωσstδεstdydz−∫yminymaxtrstδusty,zz=0dy=0,where U is the total potential energy of the system and δεst is virtual strain in wavenumber domain and frequency domain.The next step is to follow the general finite element method which has been performed in equation (12)–(18). However, the stiffness matrix is different from the one using elastic theory, as shown in the following equation:(24)Kst=∑eρ∫−11∫−11B∗NTDtiωBNJdηdζ,(25)Fst=∫yminymaxNy,zz=0pstdy.The Prony series as a general form of Maxwell model for linear viscoelasticity is used here, as shown in equation (25).(26)E∗ω=E∞+∑i=1NiωτiEi1+iωτi,where E∞ is long-term equilibrium relaxation modulus; E∗ω is dynamic modulus; τi and Ei are the parameters of Prony Series; τi is the reduced time; and Ei is the elastic component.The relationship between the shift factor and the temperature can be approximated by Williams–Landel–Ferry (WLF) function.(27)logaT=−C1T−T0C2+T−T0,where T0 is reference temperature (20°C in this study); aT is actual temperature corresponding to the shift factor with τ0=τ/aT; and C1 and C2 are regression parameters.


## 2.1.1. Elastic Governing Equations

The governing equation of elastic ground is shown in equation (3) in tensor form:(3)σij=2μεij+λδijεkk,where λ and μ are Lame constants; εiji,j=1,2,3 is strain in the ground; σiji,j=1,2,3 is the stress in the ground; and δiji,j=1,2,3 is Kronecker delta function.The equilibrium equation for elastic layers is formulated by the following expressions:(4)σij,j+fi=ρu¨i,(5)εij=12ui,j+uj,i,where ρ is the density; uii=1,2,3 is displacement vector; and fi=i=1,2,3 is body force.By substituting equations (4) and (5) into equation (3), the Navier elastic equation can be obtained, as shown in the following equation:(6)μui,jj+λ+μuj,ji+fi=ρu¨i.Navier’s equation in the frequency domain and the wavenumber domain can be described as follows:(7)μcμi,jjt+λc+μcμj,jit+ω2ρμit+fit=0,where λc and μc are Lame constants can be expressed by the following equation:(8)λc=λ1+2βi,μc=1+2βi,where β represents damp ratio and i is imaginary unit.Based on small strain theory, the strain-displacement relationship can be expressed as follows:(9)ε=εxxεyyεzzγxyγyzγzx=∂∂x000∂∂y000∂∂z∂∂y∂∂x00∂∂z∂∂y∂∂z0∂∂x=uvw=Bu.By applying the predefined Fourier transformation in equation (1) to (9), the expression can be yielded in the frequency domain and the wave number domain. The strain-displacement relationship is given by(10)εst=εxxstεyystεzzstγxystγyzstγzxst=−iζx000∂∂y000∂∂z∂∂y−iζx00∂∂z∂∂y∂∂z0−iζxustvstwst=Bstust.The stress-strain relationship is given by(11)σst=Dεst.The model can be discretized into several quadrilateral elements, and each element node has three degrees of freedom, as shown in equations (12)–(18):(12)Niη,ζ=141+ηηi1+ζζi,i=1,2,3,4,(13)Niηj,ζj=1,i=j,0,i≠j,(14)Uy,z=∑j=14Njy,zdj=NU,(15)Kst−ω2MUst=Fst,(16)M=∑eρ∫−11∫−11NTNJdηdζ,(17)Kst=∑eρ∫−11∫−11B∗stNTDBstNJdηdζ,(18)Fst=∑eρ∫−11∫−11NTfJdηdζ,where J is determinant of the matrix; N is the matrix of shape function; and D is stiffness matrix.


## 2.1.2. Viscoelastic Governing Equations

The Fourier transform is a linear transform, so this rule could be applied to all linear equations. Prony series has been widely used to capture the modulus change of linear viscoelastic materials with reduced time. For viscoelastic case, the equilibrium equation is given as(19)BTσ=f,where f is the body force vector.The strain-displacement relationship is the same as the one shown in equation (9).For stress and strain relationship, the viscoelastic model is different from elastic equation, as shown in equation (20) [32]. For stress and strain relationship, the viscoelastic model is different from elastic equation, as shown in the following equation [32]:(20)σ=∫0tDt−τdεdτdτ,where Dt−τ is stiffness history matrix.By applying Fourier transformation, the stress and strain in the frequency domain are given by(21)σt=iωDtiωεt,where Dtiω is stiffness history matrix in frequency domain.The loadP moves on the top of surface layer at a speed of v, and the surface boundary condition is expressed by(22)tr=σzzτyzτxyT=Py,t−xv.The property of material is assumed to be isotropic. By applying the energy variation principle to equation (19), equation (23) can be obtained:(23)δUω,ξx=∫Ωσstδεstdydz−∫yminymaxtrstδusty,zz=0dy=0,where U is the total potential energy of the system and δεst is virtual strain in wavenumber domain and frequency domain.The next step is to follow the general finite element method which has been performed in equation (12)–(18). However, the stiffness matrix is different from the one using elastic theory, as shown in the following equation:(24)Kst=∑eρ∫−11∫−11B∗NTDtiωBNJdηdζ,(25)Fst=∫yminymaxNy,zz=0pstdy.The Prony series as a general form of Maxwell model for linear viscoelasticity is used here, as shown in equation (25).(26)E∗ω=E∞+∑i=1NiωτiEi1+iωτi,where E∞ is long-term equilibrium relaxation modulus; E∗ω is dynamic modulus; τi and Ei are the parameters of Prony Series; τi is the reduced time; and Ei is the elastic component.The relationship between the shift factor and the temperature can be approximated by Williams–Landel–Ferry (WLF) function.(27)logaT=−C1T−T0C2+T−T0,where T0 is reference temperature (20°C in this study); aT is actual temperature corresponding to the shift factor with τ0=τ/aT; and C1 and C2 are regression parameters.


## 2.2. Arbitrary Shape of Moving Force

According to Fourier transformation, the dynamic vehicle load could be decoupled into a series of simple harmonic motions. For the sake of simplicity, only one self-oscillation frequency was considered in this case. It is assumed that the vehicle moves at speed ofv and with width of Ly along y direction (perpendicular to traffic direction) on pavement surface. The moving load of the vehicle can be expressed in the following equation:(28)fx,y,z,t=Φx−vtδy−Lyeiω0t,where Φx−vt is the load distribution along x direction and ω0 is the self-oscillation of load.The load in the frequency and wavenumber domain is shown in the following equation:(29)fstξx,y,z,ω=2πvδξx−ω−ω0vΦsξx.A typical track with five axles is employed in our case, as shown in Figure3. The track moves along the x direction at a constant speed of v and with self-oscillation frequency ω0, and the equation of moving load can be presented as shown in the following equation:(30)fx,y,z,t=p1δx−vt+p2δx−vt+L1+p3δx−vt+L1+L2+p4δx−vt+L1+L2+L3+p5δx−vt+L1+L2+L3+L4eiω0tΦx−vt,where P1,P2,P3,P4,andP5 represents the axle load of truck; L1,L2,L3,andL4 represents the distance between difference axles; and δ is the Dirac function.Figure 3
Example of typical truck with five axle loads.By employing the predefined Fourier transformation, the load in the frequency and wavenumber domain can be given in the following equation:(31)fstξx,y,z,ω=2πvδξx−ω−ω0vp1+p2∗e−iL1ξx+p3∗e−iL1+L2ξx+p4∗e−iL1+L2+L3ξx+p5∗e−iL1+L2+L3+L4ξxΦsξx.According to the definition of Dirac function, equation (32) can be obtained:(32)ξx=ω−ω0v.This means the wavenumber has specific relationship with the frequency. Then, the double dimensional Fourier transform can be simplified as one-dimensional Fourier transform, and the computational cost could be reduced significantly.Tire-pavement interaction can be modeled as different load patterns, such as point load and area-based load with different distribution patterns, as shown in Figure4.Figure 4
Illustration of different loading patterns: uniform stress in rectangular area, sinusoid-shape stress distribution, and point load. (a) Rectangle; (b) sinusoid; (c) point load.
(a)
(b)
(c)The load of distribution along traffic direction could be transformed in the wavenumber domain, as shown in equation (33)–(35).(33)ΦRx=b,x≤a,0,x>a,ΦRsξx=2∗b∗sinaξxξx,(34)ΦSx=qpeaksinπ∗x2L+π2,x≤L,0,x>L,ΦSsξx=qpeak∗sinπ/2L+ξx∗Lπ/2L+ξx+sinξx−π/2L∗Lξx−π/2L,(35)Φpsξx=2∗π∗Fv,where the subscripts r, s, and p represent rectangle shape, sinusoidal, and point load, respectively.In 3-D finite element method, all the external and body forces need be converted to equivalent nodal forces. However, in the 2.5-D method, the load was solved in frequency domain, and most of the time the distribution of regular loads had analytical forms after Fourier transformation. In this case, the analytical solution could be used directly. There are also some cases that where the distribution functions of loads along the traffic direction have no analytical solution after Fourier transformations. In this case, the method of equivalent area can be used for approximate solutions, as shown in Figure5. The loads with the elliptical or circular area are equivalent to the loads distributed in multiple rectangle-shaped areas divided from the elliptical or circular area. For exact solutions, the node forces can be calculated and treated as multiple point loads. Then, the ideal result can be acquired using equation (31). This makes the implementation of arbitrary-shape load in the 2.5-D FEM model possible.Figure 5
Illustration of ellipse and circle shapes that can be equivalent to rectangular distributed shapes.
(a)
(b)


## 2.3. Mesh and Boundary Condition

Due to the orthogonal properties of Fourier transformation, the 3-D FE model of pavement was decomposed to the 2.5-D FE model in frequency domain, while each element node in the 2.5-D FE model had movements in three dimensions. The nodes were unrestricted to move at the traffic direction. The nodes at the bottom of pavement model were fixed in this model. In order to obtain more accurate result in dynamic problem. According to previous research results, the maximum element size, minimum element size, and model size should be decided by the wavelength [33].In the dynamic problem, wave reflection will occur from the boundary, which will influence the results. The boundary condition should be selected properly. In this model, pavement structure lays on unbounded soil medium. Due to the limited size of the finite element model, it cannot handle this boundary condition directly. Several methods have been used in the literature to simulate the infinite foundation such as boundary element [34], absorbing boundary [27], thin-layer element [30], and infinite element [26, 35]. Among all these methods, the absorbing boundary showed good accuracy, and it is easy to be used in the practical application [27].Figure6 illustrates the mesh and boundary condition of the 2.5-D FEM model. In general, the minimum element size is 0.01 m by 0.0254 m at the pavement near surface, and the maximum element size is 0.5 m by 0.5 m in the subgrade. The model size was 6 m by 6 m. The left and right boundaries were absorbing boundaries, and the bottom boundary was fixed. The absorbing boundary depends on material properties and vibration frequency, as shown in the following equation:(36)ffj=kiu˙j,j=x,y,z,where ff is node force and u˙ is deformation velocity.Figure 6
General model of pavement structure and loaded area.This could be transformed by the predefined Fourier transformation in the frequency domain, as shown in equations (37)–(39).(37)ffjt=−iωkjujt,j=x,y,z,in which(38)kx=ky=ρvs,kz=ρvp,where vs is velocity of shear waves and vp is velocity of compression waves, and(39)vs=μρ,vp=λ+2μρ.The absorbing boundary can be assembled to stiffness matrix either to equation (17) (elastic material) or equation (24) (viscoelastic material).


## 2.4. Strain, Stress, Velocity, and Acceleration

According to the basic theory of the finite element method, the strain, stress, velocity, and acceleration can be calculated by the node displacements. Solving equation (15), the node displacement can be obtained. The strain, stress, velocity, and acceleration in the frequency domain can be obtained through equations (40)–(42), respectively.(40)σst=Dεst=DBstNuest,where uest is a vector concluding the node displacement of element.(41)Vst=iωμst,where Vst is the node velocity of any element.(42)ast=−ω2μst,where ast is the node acceleration of any element.


## 3. Model Verification

### 3.1. Verification with Elastic Analytical Solution

The model validation was first conducted by comparing the deflections computed using the analytical solution and the 2.5-D FEM method. A single stratum on rigid bedrock represents half-space foundation. The applied loadP was a point load moving on pavement surface along the positive x direction. The density of ground ρ was assumed to be 2000 kg/m3. Poisson’s ratio and damping ratio were assumed to be 0.25 and 0.05, respectively. The velocity of shear wave υs was 100 m/s as the input property related to material modulus.The proposed model was validated against the analytical solution by Eason [36] to verify its capability in simulating dynamic responses in homogenous elastic foundation. To demonstrate the method, the moving load was applied with constant speed of 70 m/s on pavement surface. The data extraction point was set as (0, 0, −1 m). Figure 7(a) and 7(b) show the normalized displacement obtained from the 2.5-D method and the analytical solution in the longitudinal and vertical directions, respectively. In addition to the normalized displacement, the normalized vertical stress, which was determined by the above mentioned two methods, is shown in Figure 7(c). It was noted that the plotted displacements and stresses were normalized by 2πρvs2/P and π/P, respectively. The maximum discrepancy between two solutions was found less than 10%.Figure 7
The normalized (a) longitudinal displacements; (b) vertical displacements; and (c) vertical stress in the ground (0, 0, −1 m) under moving point load at 70 m/s.
(a)
(b)
(c)


### 3.2. Verification with Viscoelastic Analytical Solution

To further validate the developed 2.5-D FEM model, another analytical model developed by Dong and Ma [9] was employed to compare the results. In this case, the viscoelastic pavement responses under moving load were solved. A three-layer flexible pavement structure was used in the analysis. The asphalt surface layer was modeled as viscoelastic material described with Prony series. The semirigid base layer and subgrade were modeled as elastic material with the modulus of 1000 MPa and 200 MPa, respectively. Four different speeds (10 m/s, 20 m/s, 25 m/s, and 40 m/s) were considered. The dual-tire loading was simulated as nonuniform tire contact stresses under each tire rib applied on pavement surface.The strains at the bottom of surface layer in three directions at two speeds (10 m/s and 40 m/s) were compared between the results obtained by analytical solutions and 2.5-D FEM model, respectively. Figure8 shows the pulse shapes and magnitudes of vertical strains at the bottom of surface layer. The results derived by 2.5-D FEM were in found good agreements with the analytical results. The speed effect was well observed with the shorter pulse period and the smaller peak strain values at the higher speed. The horizontal strains at the bottom of the asphalt layer are also presented in Figure 8. It revealed that the transverse strain was always in tensile when the tire was approaching or leaving the point of interest, while the longitudinal strains were compressive when the tire was far from the point of interest. These trends were consistent with field measurements obtained from the embedded strain gauges [37].Figure 8
Validation of strains with analytical solutions [9] at (a) 10 m/s and (b) 40 m/s.
(a)
(b)The 2.5-D method is one kind of semianalytical solution. The computation efficiency was greatly improved by means of the Fourier transformation. In other words, the 3-D problem was reduced to several 2-D problems. As a result, the number of element nodes and elements were greatly reduced. For example, if the verification case of viscoelastic solution of pavement was conducted using commercial software such as ABAQUS [38], it required 155 minutes on a workstation (Intel(R) Xeon(R) Silver 4116 CPUx2 (24 Cores), 64 GB Memory), while the 2.5-D method could obtain the same result in 226 seconds (i7 6700HQ, 8 GB Memory) in single thread.


## 3.1. Verification with Elastic Analytical Solution

The model validation was first conducted by comparing the deflections computed using the analytical solution and the 2.5-D FEM method. A single stratum on rigid bedrock represents half-space foundation. The applied loadP was a point load moving on pavement surface along the positive x direction. The density of ground ρ was assumed to be 2000 kg/m3. Poisson’s ratio and damping ratio were assumed to be 0.25 and 0.05, respectively. The velocity of shear wave υs was 100 m/s as the input property related to material modulus.The proposed model was validated against the analytical solution by Eason [36] to verify its capability in simulating dynamic responses in homogenous elastic foundation. To demonstrate the method, the moving load was applied with constant speed of 70 m/s on pavement surface. The data extraction point was set as (0, 0, −1 m). Figure 7(a) and 7(b) show the normalized displacement obtained from the 2.5-D method and the analytical solution in the longitudinal and vertical directions, respectively. In addition to the normalized displacement, the normalized vertical stress, which was determined by the above mentioned two methods, is shown in Figure 7(c). It was noted that the plotted displacements and stresses were normalized by 2πρvs2/P and π/P, respectively. The maximum discrepancy between two solutions was found less than 10%.Figure 7
The normalized (a) longitudinal displacements; (b) vertical displacements; and (c) vertical stress in the ground (0, 0, −1 m) under moving point load at 70 m/s.
(a)
(b)
(c)


## 3.2. Verification with Viscoelastic Analytical Solution

To further validate the developed 2.5-D FEM model, another analytical model developed by Dong and Ma [9] was employed to compare the results. In this case, the viscoelastic pavement responses under moving load were solved. A three-layer flexible pavement structure was used in the analysis. The asphalt surface layer was modeled as viscoelastic material described with Prony series. The semirigid base layer and subgrade were modeled as elastic material with the modulus of 1000 MPa and 200 MPa, respectively. Four different speeds (10 m/s, 20 m/s, 25 m/s, and 40 m/s) were considered. The dual-tire loading was simulated as nonuniform tire contact stresses under each tire rib applied on pavement surface.The strains at the bottom of surface layer in three directions at two speeds (10 m/s and 40 m/s) were compared between the results obtained by analytical solutions and 2.5-D FEM model, respectively. Figure8 shows the pulse shapes and magnitudes of vertical strains at the bottom of surface layer. The results derived by 2.5-D FEM were in found good agreements with the analytical results. The speed effect was well observed with the shorter pulse period and the smaller peak strain values at the higher speed. The horizontal strains at the bottom of the asphalt layer are also presented in Figure 8. It revealed that the transverse strain was always in tensile when the tire was approaching or leaving the point of interest, while the longitudinal strains were compressive when the tire was far from the point of interest. These trends were consistent with field measurements obtained from the embedded strain gauges [37].Figure 8
Validation of strains with analytical solutions [9] at (a) 10 m/s and (b) 40 m/s.
(a)
(b)The 2.5-D method is one kind of semianalytical solution. The computation efficiency was greatly improved by means of the Fourier transformation. In other words, the 3-D problem was reduced to several 2-D problems. As a result, the number of element nodes and elements were greatly reduced. For example, if the verification case of viscoelastic solution of pavement was conducted using commercial software such as ABAQUS [38], it required 155 minutes on a workstation (Intel(R) Xeon(R) Silver 4116 CPUx2 (24 Cores), 64 GB Memory), while the 2.5-D method could obtain the same result in 226 seconds (i7 6700HQ, 8 GB Memory) in single thread.


## 4. Case Study Using 2.5-D FEM

### 4.1. Pavement Structure and Material Properties

In the sensitivity analysis, the asphalt layer thickness varied from 10.16 cm to 25.4 cm, to consider the variation of pavement structure, which represents common asphalt pavement in local and major highways in the USA. The viscoelastic properties of pavement layers were kept the same as those used in the second validation case [9]. The pavement material properties of each layer are shown in Tables 1 and 2. The base and subgrade modulus was determined based on one field section constructed at MnRoad for pavement response measurement [39]. Considering that fatigue cracking and rutting are the major failure mechanism of asphalt pavement, the pavement responses considered in the analysis are tensile strains at the bottom of asphalt layer and compressive strains on top of subgrade. The same material properties were used for all the analysis cases with different tire loading patterns and speeds.Table 1
Material parameters of each pavement layer.
Pavement structure
Modulus (MPa)
Poisson’s ratio
Density (kg/m3)
Thickness (cm)
Asphalt course
Viscoelastic
0.35
2400
10.16–25.40
Base course
108
0.35
1762
20.32
Subgrade
58.6
0.40
1762
—Table 2
Viscoelastic parameters of asphalt course in Prony Series.
Relaxation time
Prony series parameters
1000
0.000541
100
0.000902
10
0.002705
1
0.007663
0.1
0.025243
0.01
0.071131
0.001
0.148753
0.0001
0.185896
0.00001
0.457167
Instantaneous modulus
10,000 MPa


### 4.2. Effect of Tire Contact Stress Patterns

Four different tire contact stress patterns were considered in the analysis, including point load, uniform contact stress in rectangular area, sinusoid-shape stress at each tire rib, and nonuniform contact stress at each tire rib. The illustrations of four contact stress patterns are shown in Figure9. The total load on dual-tire is 35.5 kN, and the tire inflation pressure is 724 kPa. The loading time is 2.56 seconds, and Fourier transformation is conducted using 0.005-second time increment. The point load assumed concentrated forces applied on the specific nodes at each tire center. The uniform contact stress distribution assumed the contact stress equal to tire inflation pressure, and the corresponding contact area was calculated with the specific width/length ratio. The nonuniform contact stress distribution is based on the measured contact stress at each tire rib provided by tire manufacturer and the predicted results from tire-pavement interaction model [37, 40, 41]. The sinusoid-shape contact stress mimics the realistic contact stress distribution using the fitted sinusoidal function, which basically kept the nonuniform pattern of contact stresses.Figure 9
Tire contact stress patterns: (a) point load, (b) uniform stress in rectangular area, (c) sinusoidal distribution, and (d) nonuniform stresses under each rib (unit: kPa).
(a)
(b)
(c)
(d)Figure10 shows the calculated deflections at pavement surface in the center of dual tires using different loading patterns. As expected, the deflections decreased with the increase of asphalt layer thickness. In general, the point load and uniform contact stresses caused very similar deflections, which were slightly greater than those under sinusoid-shaped and nonuniform contact stresses. Additionally, the deflections induced by sinusoid-shaped and nonuniform contact stresses were close to each other.Figure 10
Comparisons of deflections at pavement surface under different types of loading patterns.For asphalt pavement structures analyzed here, the critical failure mechanisms are considered to be fatigue (bottom-up) cracking and subgrade rutting, which are caused by tensile strains in the asphalt layer and compressive strains on top of subgrade, respectively. Figure11 shows the maximum horizontal tensile strains at the bottom of asphalt layer, respectively, in the transverse and longitudinal direction. Similar with deflections, the tensile strains decreased with the increasing thickness of asphalt layer. The point load caused much greater tensile strains than the area-based loading patterns, and their differences gradually decreased when asphalt layer thickness increased. As asphalt layer thickness varies from 10.16 cm to 25.4 cm, the longitudinal tensile strains induced by nonuniform tire contact stresses were 31.3%, 16.1%, 9.2%, and 6.1% smaller as compared to the ones under point load, while the transverse tensile strains induced by nonuniform tire contact stresses were 45.2%, 26.4%, 15.3%, and 9.1% smaller. On the other hand, the sinusoid-shape contact stresses caused slightly greater tensile strains as compared to the nonuniform contact stresses for thin asphalt layer.Figure 11
Comparisons of maximum transverse and longitudinal tensile strains in asphalt layer under different loading patterns.The maximum compressive strains on top of subgrade under different loading patterns are compared in Figure12. As expected, compressive strains decreased significantly as asphalt layer thickness increased. The decreasing trend was not linear and became less significant as asphalt layer was thicker. As for the effect of loading pattern, the compressive strains induced by nonuniform contact stresses were 9.7%, 2.9%, 2.4%, and 2.1% smaller than those caused by point load as asphalt layer thickness varied from 10.16 cm to 25.14 cm. On the other hand, there were no obvious differences between the compressive strains calculated using uniform and nonuniform tire contact stresses. This indicates that the effect of localized tire contact stress diminishes at the deep pavement depth.Figure 12
Comparisons of maximum compressive strains on top of subgrade different loading patterns (negative values means compression).


### 4.3. Effect of Width/Length Ratio in Rectangular Contact Area

In the uniform contact assumption in rectangular area, the width/length ratio is an important parameter that needs to be assumed to determine the shape of contact area including the regular tire and wide tire. In the analysis, the width/length ratio was changed at different values to see its effect on pavement responses. Figure13 shows the effects of width/length ratio on maximum tensile strains at the bottom of asphalt layer, respectively, for 10.16 cm, 15.24 cm, 20.32 cm, and 25.4 cm asphalt layer.Figure 13
Effects of width/length ratio on tensile strains for (a) 10.16 cm, (b) 15.24 cm, (c) 20.32 cm, and (d) 25.4 cm asphalt layer.
(a)
(b)
(c)
(d)The results indicated that as the width/length ratio of rectangular contact area increased, the longitudinal tensile strain increased but the transverse tensile strain decreased. This was because the distribution of load was assigned more in the traffic direction as compared to the transverse direction. The effect of width/length ratio on tensile strains became less significant as asphalt layer thickness increased. In general, as the width/length ratio of rectangular area was 1.3, the tensile strains caused by the uniform contact stress were close to the ones caused by nonuniform contact stresses. This indicates that if the uniform contact stress is used as the simplified assumption of tire-pavement interaction, the width/length ratio needs to be carefully selected in order to reduce the discrepancy as compared to the realistic tire-pavement contact stress pattern.


### 4.4. Effect of Moving Speed

The range of traffic speed is 10–40 m/s (26–144 km/h), which covers the typical speed of traffic. The effect of speed on asphalt pavement response was analyzed on pavement sections with four different asphalt layer thicknesses. In the analysis, the realistic nonuniform contact stresses were applied on pavement surface. Figure14 shows pavement surface deflections at the center of dual tires at different speeds. The results show that pavement surface deflections decreased slightly as speed increased regardless of pavement thickness effect. The calculated deflections at 20 m/s, 30 m/s, and 40 m/s were 1.51%, 2.09%, and 2.1% smaller than the one at 10 m/s. This indicates that the speed effect on surface deflection is not significant.Figure 14
Pavement surface deflections at different speeds.Figure15 shows maximum tensile strains at the bottom of asphalt layer at different speeds. The effect of speed on tensile strains was found more significant as compared to the case of surface deflection. As the speed increased from 10 m/s to 20 m/s and 30 m/s, the decrease of tensile strains was significant. However, the decrease of tensile strain was minor as the speed further increased from 30 m/s to 40 m/s. The strain variations with speeds were due to the fact that dynamic modulus of asphalt layer increased with the increasing speed (loading frequency) due to its viscoelastic behavior. These prediction results are consistent with the field testing results reported in reference [42].Figure 15
Comparisons of tensile strains at the bottom of asphalt layer at different speeds.The compressive strains on top of subgrade at different speeds are plotted in Figure16. Although compressive strain decreased with the increase of speed, the decreasing trend was not as significant as the strains in the asphalt layer. For example, the compressive strain at 40 m/s was 3.8% smaller than the one at 10 m/s, while the tensile strains in the asphalt layer at 40 m/s were 15.3% smaller than the one at 10 m/s. This is because the frequency dependency of asphalt layer modulus has more significant effect on pavement responses in the asphalt layer.Figure 16
Comparisons of compressive strains on top of subgrade at different speeds.


## 4.1. Pavement Structure and Material Properties

In the sensitivity analysis, the asphalt layer thickness varied from 10.16 cm to 25.4 cm, to consider the variation of pavement structure, which represents common asphalt pavement in local and major highways in the USA. The viscoelastic properties of pavement layers were kept the same as those used in the second validation case [9]. The pavement material properties of each layer are shown in Tables 1 and 2. The base and subgrade modulus was determined based on one field section constructed at MnRoad for pavement response measurement [39]. Considering that fatigue cracking and rutting are the major failure mechanism of asphalt pavement, the pavement responses considered in the analysis are tensile strains at the bottom of asphalt layer and compressive strains on top of subgrade. The same material properties were used for all the analysis cases with different tire loading patterns and speeds.Table 1
Material parameters of each pavement layer.
Pavement structure
Modulus (MPa)
Poisson’s ratio
Density (kg/m3)
Thickness (cm)
Asphalt course
Viscoelastic
0.35
2400
10.16–25.40
Base course
108
0.35
1762
20.32
Subgrade
58.6
0.40
1762
—Table 2
Viscoelastic parameters of asphalt course in Prony Series.
Relaxation time
Prony series parameters
1000
0.000541
100
0.000902
10
0.002705
1
0.007663
0.1
0.025243
0.01
0.071131
0.001
0.148753
0.0001
0.185896
0.00001
0.457167
Instantaneous modulus
10,000 MPa


## 4.2. Effect of Tire Contact Stress Patterns

Four different tire contact stress patterns were considered in the analysis, including point load, uniform contact stress in rectangular area, sinusoid-shape stress at each tire rib, and nonuniform contact stress at each tire rib. The illustrations of four contact stress patterns are shown in Figure9. The total load on dual-tire is 35.5 kN, and the tire inflation pressure is 724 kPa. The loading time is 2.56 seconds, and Fourier transformation is conducted using 0.005-second time increment. The point load assumed concentrated forces applied on the specific nodes at each tire center. The uniform contact stress distribution assumed the contact stress equal to tire inflation pressure, and the corresponding contact area was calculated with the specific width/length ratio. The nonuniform contact stress distribution is based on the measured contact stress at each tire rib provided by tire manufacturer and the predicted results from tire-pavement interaction model [37, 40, 41]. The sinusoid-shape contact stress mimics the realistic contact stress distribution using the fitted sinusoidal function, which basically kept the nonuniform pattern of contact stresses.Figure 9
Tire contact stress patterns: (a) point load, (b) uniform stress in rectangular area, (c) sinusoidal distribution, and (d) nonuniform stresses under each rib (unit: kPa).
(a)
(b)
(c)
(d)Figure10 shows the calculated deflections at pavement surface in the center of dual tires using different loading patterns. As expected, the deflections decreased with the increase of asphalt layer thickness. In general, the point load and uniform contact stresses caused very similar deflections, which were slightly greater than those under sinusoid-shaped and nonuniform contact stresses. Additionally, the deflections induced by sinusoid-shaped and nonuniform contact stresses were close to each other.Figure 10
Comparisons of deflections at pavement surface under different types of loading patterns.For asphalt pavement structures analyzed here, the critical failure mechanisms are considered to be fatigue (bottom-up) cracking and subgrade rutting, which are caused by tensile strains in the asphalt layer and compressive strains on top of subgrade, respectively. Figure11 shows the maximum horizontal tensile strains at the bottom of asphalt layer, respectively, in the transverse and longitudinal direction. Similar with deflections, the tensile strains decreased with the increasing thickness of asphalt layer. The point load caused much greater tensile strains than the area-based loading patterns, and their differences gradually decreased when asphalt layer thickness increased. As asphalt layer thickness varies from 10.16 cm to 25.4 cm, the longitudinal tensile strains induced by nonuniform tire contact stresses were 31.3%, 16.1%, 9.2%, and 6.1% smaller as compared to the ones under point load, while the transverse tensile strains induced by nonuniform tire contact stresses were 45.2%, 26.4%, 15.3%, and 9.1% smaller. On the other hand, the sinusoid-shape contact stresses caused slightly greater tensile strains as compared to the nonuniform contact stresses for thin asphalt layer.Figure 11
Comparisons of maximum transverse and longitudinal tensile strains in asphalt layer under different loading patterns.The maximum compressive strains on top of subgrade under different loading patterns are compared in Figure12. As expected, compressive strains decreased significantly as asphalt layer thickness increased. The decreasing trend was not linear and became less significant as asphalt layer was thicker. As for the effect of loading pattern, the compressive strains induced by nonuniform contact stresses were 9.7%, 2.9%, 2.4%, and 2.1% smaller than those caused by point load as asphalt layer thickness varied from 10.16 cm to 25.14 cm. On the other hand, there were no obvious differences between the compressive strains calculated using uniform and nonuniform tire contact stresses. This indicates that the effect of localized tire contact stress diminishes at the deep pavement depth.Figure 12
Comparisons of maximum compressive strains on top of subgrade different loading patterns (negative values means compression).


## 4.3. Effect of Width/Length Ratio in Rectangular Contact Area

In the uniform contact assumption in rectangular area, the width/length ratio is an important parameter that needs to be assumed to determine the shape of contact area including the regular tire and wide tire. In the analysis, the width/length ratio was changed at different values to see its effect on pavement responses. Figure13 shows the effects of width/length ratio on maximum tensile strains at the bottom of asphalt layer, respectively, for 10.16 cm, 15.24 cm, 20.32 cm, and 25.4 cm asphalt layer.Figure 13
Effects of width/length ratio on tensile strains for (a) 10.16 cm, (b) 15.24 cm, (c) 20.32 cm, and (d) 25.4 cm asphalt layer.
(a)
(b)
(c)
(d)The results indicated that as the width/length ratio of rectangular contact area increased, the longitudinal tensile strain increased but the transverse tensile strain decreased. This was because the distribution of load was assigned more in the traffic direction as compared to the transverse direction. The effect of width/length ratio on tensile strains became less significant as asphalt layer thickness increased. In general, as the width/length ratio of rectangular area was 1.3, the tensile strains caused by the uniform contact stress were close to the ones caused by nonuniform contact stresses. This indicates that if the uniform contact stress is used as the simplified assumption of tire-pavement interaction, the width/length ratio needs to be carefully selected in order to reduce the discrepancy as compared to the realistic tire-pavement contact stress pattern.


## 4.4. Effect of Moving Speed

The range of traffic speed is 10–40 m/s (26–144 km/h), which covers the typical speed of traffic. The effect of speed on asphalt pavement response was analyzed on pavement sections with four different asphalt layer thicknesses. In the analysis, the realistic nonuniform contact stresses were applied on pavement surface. Figure14 shows pavement surface deflections at the center of dual tires at different speeds. The results show that pavement surface deflections decreased slightly as speed increased regardless of pavement thickness effect. The calculated deflections at 20 m/s, 30 m/s, and 40 m/s were 1.51%, 2.09%, and 2.1% smaller than the one at 10 m/s. This indicates that the speed effect on surface deflection is not significant.Figure 14
Pavement surface deflections at different speeds.Figure15 shows maximum tensile strains at the bottom of asphalt layer at different speeds. The effect of speed on tensile strains was found more significant as compared to the case of surface deflection. As the speed increased from 10 m/s to 20 m/s and 30 m/s, the decrease of tensile strains was significant. However, the decrease of tensile strain was minor as the speed further increased from 30 m/s to 40 m/s. The strain variations with speeds were due to the fact that dynamic modulus of asphalt layer increased with the increasing speed (loading frequency) due to its viscoelastic behavior. These prediction results are consistent with the field testing results reported in reference [42].Figure 15
Comparisons of tensile strains at the bottom of asphalt layer at different speeds.The compressive strains on top of subgrade at different speeds are plotted in Figure16. Although compressive strain decreased with the increase of speed, the decreasing trend was not as significant as the strains in the asphalt layer. For example, the compressive strain at 40 m/s was 3.8% smaller than the one at 10 m/s, while the tensile strains in the asphalt layer at 40 m/s were 15.3% smaller than the one at 10 m/s. This is because the frequency dependency of asphalt layer modulus has more significant effect on pavement responses in the asphalt layer.Figure 16
Comparisons of compressive strains on top of subgrade at different speeds.


## 5. Conclusions

This study analyzed pavement responses under moving tire loading using 2.5-D FEM, which can significantly reduce computation time than 3-D FEM but maintain the proper accuracy for pavement applications. The 2.5-D approach is suitable for the situation that the speed of moving load and the pavement structure do not change along the travel direction. The accuracy of 2.5-D FEM was verified with two analytical solutions for elastic and viscoelastic pavement responses, respectively. The effects of loading pattern and speed on pavement surface deflection and strain responses were analyzed. The analyzed pavement responses included surface deflections, maximum tensile strains in the asphalt layer, and maximum compressive strains on top of subgrade.It was found that the point loading caused much greater pavement responses than the area-based loading, and thus it should not be used in pavement analysis and design. When the tire loading was simplified as uniform contact stress in rectangular area, the maximum tensile strains in the asphalt layer varied with the width/length ratio of contact area. As the width/length ratio was taken as 1.3, the uniform loading pattern produced closer results with the ones calculated using the realistic nonuniform contact stresses. On the other hand, the sinusoid-shape contact stresses resulted in similar pavement responses as the realistic nonuniform contact stresses as the asphalt layer thickness was relatively thick. The pavement strain responses decreased as the speed increased, and this effect was not affected by asphalt layer thickness.The analysis results indicate that viscoelastic responses of asphalt pavement under moving load and nonuniform contact stress can be calculated using the proposed 2.5-D FEM in a fast manner. The 2.5-D FEM has major advantage of time saving as compared to traditional 3-D FEM. This provides an effective forward model for mechanistic-empirical pavement design and analysis. More importantly, the 2.5-D FEM has the potential for backcalculation of pavement layer modulus from surface measured deflections due to the flexibility of algorithm.



---

*Source: 1029089-2020-01-18.xml*