# Application of Linear Viscoelastic Properties in Semianalytical Finite Element Method with Recursive Time Integration to Analyze Asphalt Pavement Structure

Authors: Pengfei Liu; Qinyan Xing; Dawei Wang; Markus Oeser
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Traditionally, asphalt pavements are considered as linear elastic materials in finite element (FE) method to save computational time for engineering design. However, asphalt mixture exhibits linear viscoelasticity at small strain and low temperature. Therefore, the results derived from the elastic analysis will inevitably lead to discrepancies from reality. Currently, several FE programs have already adopted viscoelasticity, but the high hardware demands and long execution times render them suitable primarily for research purposes. Semianalytical finite element method (SAFEM) was proposed to solve the abovementioned problem. The SAFEM is a three-dimensional FE algorithm that only requires a two-dimensional mesh by incorporating the Fourier series in the third dimension, which can significantly reduce the computational time. This paper describes the development of SAFEM to capture the viscoelastic property of asphalt pavements by using a recursive formulation. The formulation is verified by comparison with the commercial FE software ABAQUS. An application example is presented for simulations of creep deformation of the asphalt pavement. The investigation shows that the SAFEM is an efficient tool for pavement engineers to fast and reliably predict asphalt pavement responses; furthermore, the SAFEM provides a flexible, robust platform for the future development in the numerical simulation of asphalt pavements.


---

## Body

## 1. Introduction

The analysis of stress states is of considerable importance for the design, construction, maintenance, and rehabilitation of asphalt pavements in practice. In the past decades, many computer software applications have been developed and were increasingly used in many industrial fields as well as in the routine pavement design and assessment process. Finite element (FE) modeling is a commonly used numerical approach to solve for a layered asphalt pavement system. A variety of specialized FE codes were developed for analyzing asphalt pavements.Zeevaert [1] and Barksdale et al. [2] developed one of the most comprehensive finite element programs for the analysis of asphalt pavements, which is named “GAPPS7.” It is a geotechnical axisymmetric finite element problem solver having several essential capabilities to analyze the mechanical response of asphalt pavements. ILLI-PAVE is a commonly used finite element program developed at the University of Illinois [3], and the MICH-PAVE program was developed at the Michigan State University [4] for the analysis of asphalt pavements. Both programs modeled the pavement as an axisymmetric solid of revolution and used the K-θ model for granular materials, the bilinear approximation for fine-grained subgrade soils. The program FENLAP (Finite Element Nonlinear Analysis of Pavements) was developed by Brunton and De Almeida [5], which performs a finite element calculation of an axisymmetric solid and is designed for the structural analysis of pavements. Several material models have been adopted, including linear elastic, Brown’s [6] and Brown et al.’s [7] nonlinear models for fine-grained soils, and the K-θ model for granular materials [8]. APADS 2D was developed from 2008 under the Austroads project named “Developments of Pavement Design Models.” It adopts a two-dimensional (2D) axisymmetric concept and considers the effect of multiple loads through superposition [9].The specialized FE codes mentioned above all adopt axisymmetric concept to reduce the computational time. Due to inherent limitations of axisymmetric models, it is difficult to simulate a pavement model with a determined scale and complex loading condition. And most of them lack consideration of viscosity of asphalt mixtures. Three-dimensional (3D) finite element analysis has been increasingly considered as the best approach to capture realistic behavior of the pavement structure because it eliminates many shortcomings of the existing 2D models. CAPA-3D is a linear/nonlinear, static/dynamic finite element system for the solution of very large-scale 3D solid models such as those typically encountered in pavement and soil engineering. It includes several constitutive material models, for example, linear elasticity, hyperelasticity, elastoplasticity, and viscoelasticity. As such, it can simulate a very broad range of engineering materials under various loading conditions [10]. Unfortunately, due to its 3D characteristics, the high hardware demands and long execution times render it suitable primarily for research purposes [11]. The CESAR-LCPC was developed by Laboratoire Central des Ponts et Chaussées (LCPC) and is able to perform static and dynamic analyses with a 2D or 3D model. However, it does not include a viscoelastic constitutive model which is extremely important for the simulation of creep in asphalt layers [12]. The other general purpose FE software, such as ABAQUS, ANSYS, and ADINA, are being widely used in structural analysis of asphalt pavements. They may provide more powerful capabilities to simulate the response of asphalt pavements to a certain extent. However, the expensive costs to get valid licenses and the time-consuming training process often render it impractical to be used by a pavement engineer.FE code specifically used for analyses of asphalt pavements is still necessary to be developed. The main objective is the development of a fast and accurate FE code with the capability of producing a structural model of asphalt pavement, like the one in Figure1, using specific constitutive models and implementing realistic boundary conditions without requiring overlong computation time. One method was proposed to fulfill these requirements, which assumes the displacements in the geometrical z-direction and can be represented using Fourier series. Exploiting its orthogonal properties, the problem of such a class can thus be simplified into a series of 2D solutions. This method is the so-called semianalytical FE method and was first developed in linear analysis by Wilson [13]. Meissner [14] extended Wilson’s work to an axisymmetric elastoplastic body. Winnicki and Zienkiewicz [15] used a viscoplastic formulation to tackle material nonlinearity. Carter and Booker [16] provided an efficient analysis of the consolidation of axisymmetric elastic bodies subjected to nonsymmetric loading by using the continuous Fourier series. Lai and Booker [17] successfully applied a discrete Fourier technique to analyze the nonlinear behavior of axisymmetric solids under 3D loading conditions. Further development was taken by Fritz [18] and Hu et al. [19], who programmed specific FE codes to apply this method in analyzing the full 3D problems of asphalt pavements. However, their FE codes are relatively simple; for example, only static analysis with linear elastic material properties and a total bond between pavement layers can be performed. Recently, an FE code named “SAFEM” was developed by Liu et al. [20–27]. The partial bond condition between pavement layers can be considered, and infinite elements were coupled to reduce the influence of the boundary on the computational results.Figure 1
Pavement structure geometry and load mode.As mentioned above, asphalt pavements are inherently viscoelastic structures [28], and it is of great importance for pavement engineers to assess the pavement condition when rutting is considered. Extensive research has been conducted in the area of numerical analysis of viscoelastic problems. The literature can be divided into two categories: (1) transformation-based or correspondence principle-based methods and (2) time-integration methods. Examples of the first approach include works by Hilton and Yi [29], which describe a transformation-based approach for anisotropic viscoelastic composites. This approach relies on the use of integral transformations, and therefore, it is imposed with the limitations of the correspondence principle as described by Mukherjee and Paulino [30] and Rajagopal and Wineman [31]. The second approach, time-integration schemes, has been explored for viscoelastic FEs such as incremental schemes [32, 33] and recursive schemes [34–36]. Incremental schemes are based on the determination of response increment in each time step, whereas recursive schemes rely on keeping track of the material history effect and updating it at each time step. The time-integration approaches are not imposed with limitations on material model description or boundary conditions as required by the correspondence principle.In this study, the specific recursive time-integration scheme in conjunction with the viscoelastic constitutive model for SAFEM is proposed. The outstanding advantage of this method is the significant improvement of the computational speed compared with conventional 3D FE simulation with viscoelasticity. In the following sections, the linear viscoelastic constitutive model and mathematical basis of the SAFEM are introduced, which is followed by a recursive time-integration scheme. Whereafter, verification through comparison with numerical solution from ABAQUS is carried out. The latter portion of this paper describes one application example: the impact of axle load and driving speed on creep deformation of asphalt pavements is discussed. Finally, a brief summary and conclusions are provided.


## 2. Application of Linear Viscoelastic Model in Semianalytical Finite Element Method

### 2.1. Linear Viscoelastic Material Constitutive Model

Linear viscoelasticity is usually applicable for small deformations of materials. Analogous to most viscoelastoplastic materials, the creep curve of asphalt mastic can be generally divided into three stages: decelerated creep, equivelocity creep, and accelerated creep. These stages of viscoelastoplastic behavior of asphaltic mixture are shown in Figure2 [37].Figure 2
A typical creep curve of asphaltic mixture at a temperature of 50°C and stress of 0.25 MPa.The Burgers model, which exhibits linear viscoelastic behavior, is a Maxwell model in series with a Kelvin model. When undergoing a creep test, the strain response of a material, described by the Burgers model, is a combination of elastic, delayed elastic, and viscous strain responses [38], as shown in Figure 3. It is suitable to describe the first two stages of creep, namely, decelerated and equivelocity creep stages. Viscoelastic materials, such as the asphalt mixtures, can be modeled by the Burgers model in order to determine their stress or strain interactions as well as their temporal dependencies [37].Figure 3
The Burgers model and its creep response.This model has the following response equation:(1)σ+p1σ˙+p2σ¨=q1ε˙+q2ε¨p1=η1E1+η1E2+η2E1E1E2,p2=η1η2E1E2,q1=η1,q2=η1η2E2,where E1, E2 and η1, η2 are the modulus of elasticity and the viscosity of the material in Maxwell and Kelvin models, respectively; σ and ε are stress and strain, respectively; σ˙, ε˙ and σ¨, ε¨ are the first and second derivative of stress and strain with respect to time t, respectively.When the material exerts a constant stress, the creep function can be expressed as follows:(2)Jt=1E1+tη1+1E21−exp−E2tη2.The corresponding relaxation function is the following:(3)Et=q2p2α−βα−q1q2exp−αt+q1q2−βexp−βtα,β=p1∓p12−4p22p2.


### 2.2. Formulation of Semianalytical Finite Element Method

The first step in the FE formulation is to express the element coordinates and element displacements in the form of interpolations using the natural coordinate system of the element. By using SAFEM, the general form of the shape functions defining the variation of displacements can be written as a Fourier series in whichz ranges between zero and a, as shown in Figure 4 [18–27, 39].(4)u=∑k=16Nkx,y,zuk=∑l=1L∑k=16Nk¯x,ycoslπza+Nk¯¯x,ysinlπzaukl,where l identifies the term of the Fourier series and L is the number of terms considered; Nk¯ and Nk¯¯ are the shape functions of the node in the XY plane which are the same as those for displacement approximation used in 2D problems.Figure 4
Schematic representation of an SAFEM situation.The loading function defining the variation of load along thez-direction is given a similar form as the displacements [19]:(5)f=∑l=1Lp¯x,ycoslπza+p¯¯x,ysinlπza,where p¯x,y and p¯¯x,y represent the external load function.The pavement is assumed to be held atz = 0 and z = a in a manner preventing all displacements in the XY plane but permitting unrestricted motion in the z-direction. The Fourier series expansion should meet this requirement of the boundary condition, and the displacement functions with the three components u,v,andw can be written as follows [39]:(6)U=uvw=∑l=1L∑k=16Nksinlπza000sinlπza000coslπzauklvklwkl=∑l=1LNl·Ul,where ukl, vkl, and wkl are the displacements of the node at the term of the Fourier series along x-, y-, and z-directions, respectively.Similarly, the loading function for the pavement analysis can be simplified as [19](7)f=∑l=1Lsinlπza000sinlπza000coslπza0px,y0=∑l=1Lpl,(8)px,y=∑t=1n2PtlπcoslπaZt1−coslπaZt2,where Pt is the tire load pressure; Zt1 is the z-coordinate where the tire load starts; Zt2 is the z-coordinate where the tire load ends. Equation (8) was used to simulate the rectangular/square contact shape in previous works [18, 19]. Actually, if Zt1 and Zt2 vary with the x-coordinates of different contact elements, different contact shapes can be captured such as circular, ellipsoid, and even irregular shapes [20–27].After determining the element displacement, geometrical and physical equations can be used to obtain the strain and stress of one element. The strain-displacement matrixBkl is defined as follows:(9)Bkl=∂Nk∂xsinlπza000∂Nk∂ysinlπza000−lπaNksinlπza∂Nk∂ysinlπza∂Nk∂xsinlπza00lπaNkcoslπza∂Nk∂ycoslπzalπaNkcoslπza0∂Nk∂xcoslπza.By using the principle of minimum potential energy, a typical submatrix of the element stiffness matrixKlme is [39](10)Klme=∭volBlTDBmdxdydz.A typical term for the force vector becomes(11)Fle=∭volNlTpldxdydz.From (9) and (10), the stiffness matrix of one element includes [39](12)I1=∫0asinlπza·cosmπza·dzI2=∫0asinlπza·sinmπza·dzI3=∫0acoslπza·cosmπza·dz.The integrals exhibit orthogonal properties which ensure that(13)I2=I3=12a,forl=m0,forl≠m.Only whenl and m are both odd or even numbers, the first integral I1 is zero. Due to the special structure of the Bl matrix, all terms that include I1 vanish (become zero). This means that the matrix Klme becomes a diagonal one. In other words, the nonzero values are only located in the diagonal area where l = m. Thus, the stiffness matrix can be reduced, and the final assembled equations have the following form [39]:(14)K11K22⋱KLLU1U2⋮UL+F1F2⋮FL=0.The last equation shows that the large system of equations splits up intoL separate problems.(15)KllUl+Fl=0.According to (15), the Fourier expansion of the loading factors involves only one term for a particular harmonic, so only one set of simultaneous equations needs to be solved. This solution is just like a 2D plane problem except for three degrees of freedom in this case but only two in the case of a 2D plane. The subdisplacement vector calculated from each term of the Fourier series only needs to be assembled into a global vector.


### 2.3. Recursive Time-Integration Scheme

A recursive time-integration scheme is generated for the linear viscoelastic response based on the Burgers model and implemented in the SAFEM code by a recursive-iterative algorithm.The strain response for isotropic materials can be decoupled into deviatoric and spherical parts. It can be presented as(16)εij=eij+13εkkδij=12GSij+σkk9KδijSij=σij−13σkkδijσkk=σ1+σ2+σ3δij=0fori≠j1fori=j,where G and K are the shear modulus and bulk modulus, respectively. eij is the deviatoric strain, and εkk is the spherical strain. Sij is the deviatoric stress, and σkk is the spherical stress. δij is the Kronecker delta.When asphalt mixtures are applied by a constant load, the deviatoric and spherical strains can be expressed as follows by using the constitutive relationship of the Burgers model:(17)eijt=1G1+t2η1G+1G21−exp−G2tη2GSijεkkt=1K1+t2η1K+1K21−exp−K2tη2Kσkk,where G1 and G2 are the shear modulus of the elastic element in Maxwell and Kelvin elements, respectively. K1 and K2 are the bulk modulus of the elastic element in Maxwell and Kelvin elements, respectively. η1G and η2G are the viscosity coefficient of the viscous element for deviatoric strain response in Maxwell and Kelvin elements, respectively. η1K and η2K are the viscosity coefficient of the viscous element for spherical strain response in Maxwell and Kelvin elements, respectively.The experimental measurements showed that the Poisson’s ratiov varies only slightly for wide ranges of temperatures and loading rates [40, 41], and ν is therefore assumed to be time independent leading to the following expressions of the abovementioned parameters:(18)G1=E121+ν1,G2=E221+ν2,K1=E131−2ν1,K2=E231−2ν2,η1G=η121+ν1G,η2G=η221+ν2G,η1K=η131−2ν1K,η2K=η231−2ν2K,where ν1G and ν2G are Poisson’s ratio of the viscous element for deviatoric strain response in Maxwell and Kelvin elements, respectively. ν1K and ν2K are the Poisson’s ratio of the viscous element for spherical strain response in Maxwell and Kelvin elements, respectively.If the load applied on the asphalt mixtures is not constant, the deviatoric and spherical strains should be written by integral convolution:(19)eijt=∫0t1G1+t−τ2η1G+1G21−exp−G2t−τη2GS˙ijdτεkkt=∫0t1K1+t−τ2η1K+1K21−exp−K2t−τη2Kσ˙kkdτ.For implementation in the SAFEM algorithm, the incremental deviatoric and spherical strains are written in terms of history integrals as follows:(20)Δeijti=1G1+Δti2η1G+1G21−η2GG2Δti1−exp−G2Δtiη2GΔSiji+1G21−exp−G2Δtiη2Gɛ2Gi−1+Δtiη1Gɛ1Gi−1Δɛkkti=1K1+Δti2η1K+1K21−η2KK2Δti1−exp−K2Δtiη2KΔσkki+1K21−exp−K2Δtiη2Kɛ2Ki−1+Δtiη1Kɛ1Ki−1ɛ2Gi−1=η2GG2Δti−11−exp−G2Δti−1η2GΔSiji−1+exp−G2Δti−1η2Gɛ2Gi−2,ɛ1Gi−1=ΔSiji−1+ɛ1Gi−2,ɛ2Ki−1=η2KK2Δti−11−exp−K2Δti−1η2KΔσkki−1+exp−K2Δti−1η2Kɛ2Ki−2,ɛ1Ki−1=Δσkki−1+ɛ1Ki−2,ɛ1G0=ɛ2G0=ɛ1K0=ɛ2K0=0.The variablesɛ1Gi,check spaceɛ2Gi,ɛ1Ki,andɛ2Ki are the deviatoric and spherical history variables, respectively, for every increment term i. The history variables are updated at the end of every converged time increment, which will be used for the next time increment.The deviatoric and spherical stress increments can be determined from (20), provided that the strains are given. This algorithm is suitably implemented in the displacement-based SAFEM framework, in which strains are the given variables. The current deviatoric and spherical stresses and the current material stiffness cannot be determined directly because the material stiffness is dependent on the current stress and vice versa. Hence, the modified Newton-Raphson iterative algorithm is adopted to solve for the current stress state, in which material stiffness is assumed at the beginning of each time increment. The force residual criterion is used as the convergence criterion for this iterative algorithm. The force residual vector is defined as the difference between internal forces computed from the FE solution and material constitutive relation. Convergence is assumed to be achieved when the sum of all the elements in the force residual vector is less than 3% of the sum of them at the beginning of each time increment [42].The flowchart of this algorithm is shown in Figure5.Figure 5
The flowchart of the recursive-iterative algorithm.


## 2.1. Linear Viscoelastic Material Constitutive Model

Linear viscoelasticity is usually applicable for small deformations of materials. Analogous to most viscoelastoplastic materials, the creep curve of asphalt mastic can be generally divided into three stages: decelerated creep, equivelocity creep, and accelerated creep. These stages of viscoelastoplastic behavior of asphaltic mixture are shown in Figure2 [37].Figure 2
A typical creep curve of asphaltic mixture at a temperature of 50°C and stress of 0.25 MPa.The Burgers model, which exhibits linear viscoelastic behavior, is a Maxwell model in series with a Kelvin model. When undergoing a creep test, the strain response of a material, described by the Burgers model, is a combination of elastic, delayed elastic, and viscous strain responses [38], as shown in Figure 3. It is suitable to describe the first two stages of creep, namely, decelerated and equivelocity creep stages. Viscoelastic materials, such as the asphalt mixtures, can be modeled by the Burgers model in order to determine their stress or strain interactions as well as their temporal dependencies [37].Figure 3
The Burgers model and its creep response.This model has the following response equation:(1)σ+p1σ˙+p2σ¨=q1ε˙+q2ε¨p1=η1E1+η1E2+η2E1E1E2,p2=η1η2E1E2,q1=η1,q2=η1η2E2,where E1, E2 and η1, η2 are the modulus of elasticity and the viscosity of the material in Maxwell and Kelvin models, respectively; σ and ε are stress and strain, respectively; σ˙, ε˙ and σ¨, ε¨ are the first and second derivative of stress and strain with respect to time t, respectively.When the material exerts a constant stress, the creep function can be expressed as follows:(2)Jt=1E1+tη1+1E21−exp−E2tη2.The corresponding relaxation function is the following:(3)Et=q2p2α−βα−q1q2exp−αt+q1q2−βexp−βtα,β=p1∓p12−4p22p2.


## 2.2. Formulation of Semianalytical Finite Element Method

The first step in the FE formulation is to express the element coordinates and element displacements in the form of interpolations using the natural coordinate system of the element. By using SAFEM, the general form of the shape functions defining the variation of displacements can be written as a Fourier series in whichz ranges between zero and a, as shown in Figure 4 [18–27, 39].(4)u=∑k=16Nkx,y,zuk=∑l=1L∑k=16Nk¯x,ycoslπza+Nk¯¯x,ysinlπzaukl,where l identifies the term of the Fourier series and L is the number of terms considered; Nk¯ and Nk¯¯ are the shape functions of the node in the XY plane which are the same as those for displacement approximation used in 2D problems.Figure 4
Schematic representation of an SAFEM situation.The loading function defining the variation of load along thez-direction is given a similar form as the displacements [19]:(5)f=∑l=1Lp¯x,ycoslπza+p¯¯x,ysinlπza,where p¯x,y and p¯¯x,y represent the external load function.The pavement is assumed to be held atz = 0 and z = a in a manner preventing all displacements in the XY plane but permitting unrestricted motion in the z-direction. The Fourier series expansion should meet this requirement of the boundary condition, and the displacement functions with the three components u,v,andw can be written as follows [39]:(6)U=uvw=∑l=1L∑k=16Nksinlπza000sinlπza000coslπzauklvklwkl=∑l=1LNl·Ul,where ukl, vkl, and wkl are the displacements of the node at the term of the Fourier series along x-, y-, and z-directions, respectively.Similarly, the loading function for the pavement analysis can be simplified as [19](7)f=∑l=1Lsinlπza000sinlπza000coslπza0px,y0=∑l=1Lpl,(8)px,y=∑t=1n2PtlπcoslπaZt1−coslπaZt2,where Pt is the tire load pressure; Zt1 is the z-coordinate where the tire load starts; Zt2 is the z-coordinate where the tire load ends. Equation (8) was used to simulate the rectangular/square contact shape in previous works [18, 19]. Actually, if Zt1 and Zt2 vary with the x-coordinates of different contact elements, different contact shapes can be captured such as circular, ellipsoid, and even irregular shapes [20–27].After determining the element displacement, geometrical and physical equations can be used to obtain the strain and stress of one element. The strain-displacement matrixBkl is defined as follows:(9)Bkl=∂Nk∂xsinlπza000∂Nk∂ysinlπza000−lπaNksinlπza∂Nk∂ysinlπza∂Nk∂xsinlπza00lπaNkcoslπza∂Nk∂ycoslπzalπaNkcoslπza0∂Nk∂xcoslπza.By using the principle of minimum potential energy, a typical submatrix of the element stiffness matrixKlme is [39](10)Klme=∭volBlTDBmdxdydz.A typical term for the force vector becomes(11)Fle=∭volNlTpldxdydz.From (9) and (10), the stiffness matrix of one element includes [39](12)I1=∫0asinlπza·cosmπza·dzI2=∫0asinlπza·sinmπza·dzI3=∫0acoslπza·cosmπza·dz.The integrals exhibit orthogonal properties which ensure that(13)I2=I3=12a,forl=m0,forl≠m.Only whenl and m are both odd or even numbers, the first integral I1 is zero. Due to the special structure of the Bl matrix, all terms that include I1 vanish (become zero). This means that the matrix Klme becomes a diagonal one. In other words, the nonzero values are only located in the diagonal area where l = m. Thus, the stiffness matrix can be reduced, and the final assembled equations have the following form [39]:(14)K11K22⋱KLLU1U2⋮UL+F1F2⋮FL=0.The last equation shows that the large system of equations splits up intoL separate problems.(15)KllUl+Fl=0.According to (15), the Fourier expansion of the loading factors involves only one term for a particular harmonic, so only one set of simultaneous equations needs to be solved. This solution is just like a 2D plane problem except for three degrees of freedom in this case but only two in the case of a 2D plane. The subdisplacement vector calculated from each term of the Fourier series only needs to be assembled into a global vector.


## 2.3. Recursive Time-Integration Scheme

A recursive time-integration scheme is generated for the linear viscoelastic response based on the Burgers model and implemented in the SAFEM code by a recursive-iterative algorithm.The strain response for isotropic materials can be decoupled into deviatoric and spherical parts. It can be presented as(16)εij=eij+13εkkδij=12GSij+σkk9KδijSij=σij−13σkkδijσkk=σ1+σ2+σ3δij=0fori≠j1fori=j,where G and K are the shear modulus and bulk modulus, respectively. eij is the deviatoric strain, and εkk is the spherical strain. Sij is the deviatoric stress, and σkk is the spherical stress. δij is the Kronecker delta.When asphalt mixtures are applied by a constant load, the deviatoric and spherical strains can be expressed as follows by using the constitutive relationship of the Burgers model:(17)eijt=1G1+t2η1G+1G21−exp−G2tη2GSijεkkt=1K1+t2η1K+1K21−exp−K2tη2Kσkk,where G1 and G2 are the shear modulus of the elastic element in Maxwell and Kelvin elements, respectively. K1 and K2 are the bulk modulus of the elastic element in Maxwell and Kelvin elements, respectively. η1G and η2G are the viscosity coefficient of the viscous element for deviatoric strain response in Maxwell and Kelvin elements, respectively. η1K and η2K are the viscosity coefficient of the viscous element for spherical strain response in Maxwell and Kelvin elements, respectively.The experimental measurements showed that the Poisson’s ratiov varies only slightly for wide ranges of temperatures and loading rates [40, 41], and ν is therefore assumed to be time independent leading to the following expressions of the abovementioned parameters:(18)G1=E121+ν1,G2=E221+ν2,K1=E131−2ν1,K2=E231−2ν2,η1G=η121+ν1G,η2G=η221+ν2G,η1K=η131−2ν1K,η2K=η231−2ν2K,where ν1G and ν2G are Poisson’s ratio of the viscous element for deviatoric strain response in Maxwell and Kelvin elements, respectively. ν1K and ν2K are the Poisson’s ratio of the viscous element for spherical strain response in Maxwell and Kelvin elements, respectively.If the load applied on the asphalt mixtures is not constant, the deviatoric and spherical strains should be written by integral convolution:(19)eijt=∫0t1G1+t−τ2η1G+1G21−exp−G2t−τη2GS˙ijdτεkkt=∫0t1K1+t−τ2η1K+1K21−exp−K2t−τη2Kσ˙kkdτ.For implementation in the SAFEM algorithm, the incremental deviatoric and spherical strains are written in terms of history integrals as follows:(20)Δeijti=1G1+Δti2η1G+1G21−η2GG2Δti1−exp−G2Δtiη2GΔSiji+1G21−exp−G2Δtiη2Gɛ2Gi−1+Δtiη1Gɛ1Gi−1Δɛkkti=1K1+Δti2η1K+1K21−η2KK2Δti1−exp−K2Δtiη2KΔσkki+1K21−exp−K2Δtiη2Kɛ2Ki−1+Δtiη1Kɛ1Ki−1ɛ2Gi−1=η2GG2Δti−11−exp−G2Δti−1η2GΔSiji−1+exp−G2Δti−1η2Gɛ2Gi−2,ɛ1Gi−1=ΔSiji−1+ɛ1Gi−2,ɛ2Ki−1=η2KK2Δti−11−exp−K2Δti−1η2KΔσkki−1+exp−K2Δti−1η2Kɛ2Ki−2,ɛ1Ki−1=Δσkki−1+ɛ1Ki−2,ɛ1G0=ɛ2G0=ɛ1K0=ɛ2K0=0.The variablesɛ1Gi,check spaceɛ2Gi,ɛ1Ki,andɛ2Ki are the deviatoric and spherical history variables, respectively, for every increment term i. The history variables are updated at the end of every converged time increment, which will be used for the next time increment.The deviatoric and spherical stress increments can be determined from (20), provided that the strains are given. This algorithm is suitably implemented in the displacement-based SAFEM framework, in which strains are the given variables. The current deviatoric and spherical stresses and the current material stiffness cannot be determined directly because the material stiffness is dependent on the current stress and vice versa. Hence, the modified Newton-Raphson iterative algorithm is adopted to solve for the current stress state, in which material stiffness is assumed at the beginning of each time increment. The force residual criterion is used as the convergence criterion for this iterative algorithm. The force residual vector is defined as the difference between internal forces computed from the FE solution and material constitutive relation. Convergence is assumed to be achieved when the sum of all the elements in the force residual vector is less than 3% of the sum of them at the beginning of each time increment [42].The flowchart of this algorithm is shown in Figure5.Figure 5
The flowchart of the recursive-iterative algorithm.


## 3. Verification of the SAFEM

Analytical verification was carried out to check the accuracy of SAFEM with linear viscoelastic material properties. The results were compared with ABAQUS. The geometry and identical material properties of the pavement used in both programs are listed in Table1. According to preinvestigations by authors, the length and width of the pavement were finally defined as 6000 mm to minimize the influence of the pavement edge effect on the results and also to limit the time required for the mesh generation and the following computational calculation. In both analyses, the deformation of viscoelastic bulk was not considered. The Burgers model was applied in SAFEM to simulate the viscoelastic material property of asphalt layers. Creep tests at a defined temperature of 20°C were carried out, and the material parameters used for the Burgers model were derived by using nonlinear least squares regression based on the creep curve [43], as listed in Table 2. The integral constitutive equations of relaxation shear moduli were adopted in the viscoelastic model of ABAQUS and expanded by a Prony series. The parameters in the Burgers model were converted into relaxation shear moduli of the Prony series form in ABAQUS using the Laplace transform and other equations [44], as listed in Table 3. Once the temperature was determined, the material parameters were thus determined, and the temperature did not explicitly appear in the formula of the materials stiffness matrix in our algorithm.Table 1
Geometry and identical material properties of the pavement in both programs.
Layer
Thickness (mm)
E (MPa)
µ
Surface course
40
See Tables2 and 3
0.35
Binder course
60
See Tables2 and 3
0.35
Asphalt base course
70
See Tables2 and 3
0.35
Road base course
250
1500
0.49
Subbase
250
1200
0.49
Subgrade
5000
50
0.49Table 2
Material parameters of the Burgers model used in SAFEM.
Layer
E
1 (MPa)
η
1 (MPa·s)
E
2 (MPa)
η
2 (MPa·s)
Surface course
229.8
49,196.4
563.8
30,108.4
Binder course
235.8
492,835
582.3
30,187.0
Asphalt base course
233.6
492,955
585.3
31,002.1Table 3
Prony series parameters used in ABAQUS.
Layer
g
1
g
2
τ
1 (s)
τ
2 (s)
Surface course
0.42522
0.57477
33.3445
285.719
Binder course
0.33593
0.66406
34.6577
2605.26
Asphalt base course
0.33285
0.66714
35.5770
2618.16The meshes generated from the SAFEM and ABAQUS are shown in Figure6. A full-size model was created in the SAFEM, while a half symmetrical model was in the ABAQUS.Figure 6
Mesh automatically generated from (a) SAFEM and (b) ABAQUS in static analysis with linear viscoelastic material property.
(a)
(b)The load in SAFEM and ABAQUS was assumed as a square load with the side length of 300 mm. The uniformly distributed contact pressure of 0.7 MPa was fully applied at the first increment with time of 1 × 10−8 s and then kept constant for 1000 s which is divided into 20 equal increments; afterwards, the pressure was removed at the 22nd increment with time of 20 s which was followed by 49 increments with the same duration for each. The bottom nodes of the mesh representing the subgrade in both models were fixed in all directions. According to the theory of SAFEM, the displacements on both edges (z = 0 and z = a) are restricted to zero in the x- and y-directions [19]. The equivalent boundary conditions were also used in the ABAQUS model. The three asphalt layers were totally bound; the two contact layers among the asphalt base course, road base course, subbase, and subgrade were defined as being partially bound, which means the interface between layers keeps the same displacement in the vertical direction and allows different displacements in horizontal directions.The computational vertical displacement at the end of the 21st increment from both programs is compared in Figure7. The cross section is also through the centroid of the full-size pavement model and along the traffic direction. The distribution of the displacement magnitude and the deformation shapes from both FE programs are consistent with one another.Figure 7
The computational vertical displacement from (a) SAFEM and (b) ABAQUS with a scale factor of 50 in static analysis with linear viscoelastic material property.
(a)
(b)The computational vertical stresses and strains at the center of the pavement surface are shown in Figure8. Both computational results are in good agreement with one another. The stresses increase from 0 to 0.7 MPa in 1 × 10−8 s, remain nearly constant for 1000 s, decrease to zero in 20 s, and remain as zero at the end of simulation, which is consistent with the evolution of the applied pressure. The strains increase linearly at the first increment and then increase along a degressively steep path, while the stresses are constant, which is the so-called creep phenomenon. When the loads are removed, the strains decrease gradually and finally are infinitely close to a constant level.Figure 8
Comparison of (a) vertical stresses and (b) vertical strains at the center of the pavement surface between ABAQUS and SAFEM.
(a)
(b)In Table4, the vertical strains at the center of the tire pavement contact area derived from 6 specific increments are considered. The differences are extremely small, which further proves that the results from both programs have a high correlation.Table 4
Comparison between ABAQUS and SAFEM regarding the determined vertical strains (−) at the center of the tire pavement contact area.
Time (s)
SAFEM
ABAQUS
Difference
1 × 10−8
−0.00220
−0.00219
0.457%
500
−0.00450
−0.00444
1.351%
1000
−0.00580
−0.00580
0
1020
−0.00350
−0.00352
−0.797
1500
−0.00250
−0.00254
−1.823
2000
−0.00220
−0.00228
−3.797The SAFEM needs significantly fewer elements and nodes compared with ABAQUS, particularly due to the use of a Fourier series in the transverse direction; this results in a far shorter computational time for the SAFEM as compared to ABAQUS. Both FE analyses were run on a computer with an Intel Core Duo 3.4 GHz, 32 GB RAM. The computational time required by the half symmetrical model in ABAQUS is 3.1 hours, whereas the full-size model in SAFEM requires only 1.8 hours, as shown in Table5. The computational time of the SAFEM is much shorter than that of the ABAQUS, especially when the model is more complicated and full-size model is required. With code optimization, the computation time of the SAFEM will be reduced further.Table 5
Comparison of the elements, nodes, and computational time between SAFEM and ABAQUS in static analysis with linear viscoelastic material property.
Full-size model in SAFEM
Half symmetrical model in ABAQUS
Elements
1321
101459
Nodes
2738
184391
Computational time (h)
1.8
3.1


## 4. Analysis of Creep Deformation of Asphalt Pavements Using SAFEM

With the increase of the vehicle number and overload phenomenon, the creep deformation of asphalt pavements has become more and more serious and caused several asphalt pavement distresses, such as rutting, which has an adverse impact on ride comfort and road safety. It is of great importance for pavement engineers to research the creep deformation of asphalt pavements in order to improve the pavement service life [45, 46]. In this section, 3D finite element models of asphalt pavements using SAFEM were applied to research the impact of different loading conditions on the creep deformation.


### 4.1. Definition of the Models

The thicknesses and the material properties of the pavement layers are listed in Table6, which is a typical pavement structure adopted in Germany [47]. The material properties are derived from the defined temperature of 20°C. The length and width of all layers were also set to 6000 mm. The upper three asphalt layers were totally bound. The two contact layers among the asphalt base course, gravel base layer, frost protection layer, and subgrade were defined as being partially bound.Table 6
Geometry and material properties of the pavement.
Layer
Thickness (mm)
E (MPa)
µ
Surface course
40
See Table2
0.35
Binder course
50
See Table2
0.35
Asphalt base course
110
See Table2
0.35
Road base course
150
157.8
0.49
Subbase
570
125.7
0.49
Subgrade
2000
98.9
0.49A standard single axle with dual wheels was utilized in this study, and its configuration can be seen in Figure9.Figure 9
The configuration of the single axle with dual wheels.Considering different working conditions, five different levels of axle loads were applied. The corresponding loading parameters are shown in Table7.Table 7
Loading parameters of the asphalt pavement [48].
Axle load,P (kg)
Contact pressure,p (MPa)
Contact area,A (mm2)
Wheel width,W (mm)
Wheel length,L (mm)
Distance between two wheel centers
d (mm)
D (mm)
10,000
0.70
35,256
156
226
320
1600
12,500
0.76
40,172
166
242
320
1600
15,000
0.82
45,056
176
256
320
1600
17,500
0.85
50,220
186
270
320
1600
20,000
0.89
55,185
195
283
320
1600When moving loads pass for a large number of times, the pavement can be assumed to undertake a series of pulse loads [49], as shown in Figure 10(a). However, when the viscoelasticity is considered in the material property, the pulse loads are extremely time-consuming due to large numbers of increments and/or iterations used in the computation. Former investigation [49] showed the model in which an equivalent wheel load can reasonably predict the creep deformation of the pulse loads while its computational time is significantly shorter. As a result, the equivalent loading mode was applied in this study, as shown in Figure 10(b).Figure 10
The sketch of two loading modes. (a) Pulse loading mode and (b) equivalent loading mode.
(a)
(b)The loading timeTl for each pulse loading is(21)Tl=Wv,where W is the wheel width and v is the driving speed.Equivalent loading is obtained by substituting pulse loading by a one-step equivalent loading whose duration is equal to the summation of loading times in pulse loading.(22)T=N·Tl,where N is the number of loading cycles.The contact pressure of 0.70 MPa with five different driving speeds, such as 60, 80, 100, 120, and 140 km/h, was applied in this study to investigate the impact of vehicle speed on the creep deformation of asphalt pavements.The automatically generated mesh from SAFEM is shown in Figure11.Figure 11
Mesh automatically generated from SAFEM.


### 4.2. Impact of the Axle Loads

When the driving speed is 60 km/h and the contact pressure is 0.70 MPa, the computational vertical displacement after loading cycles of 1 × 106 is shown in Figure 12. This cross section is across the centers of loading areas. The creep deformation can be visualized from this figure, most of which occurs close to the loading areas and decreases rapidly with the increase of distance from the loading areas. As a result, it is meaningful to investigate the near-surface pavement response with more comprehensive model in the future study, which highlights the importance of high efficiency of SAFEM.Figure 12
The computational vertical displacement (mm) with a scale factor of 10.The impact of different axle loads was investigated. The vertical displacement was derived from the center of the loading area 1 on the pavement surface. The variation of the vertical displacement with the increase of loading cycles is shown in Figure13(a). For each axle load, the vertical displacement increases with the increase of loading cycles linearly.Figure 13
Comparison of the vertical displacement under different axle loads. (a) Variation with the increase of loading cycles and (b) variation along the transverse directionZ of the pavement.
(a)
(b)The variation of the vertical displacement after loading cycles of 1 × 106 along the transverse direction Z of the pavement is shown in Figure 13(b). The most deformation occurs below the axle loading area and increases with the increase of axle loads.The maximum surface deflections after loading cycles of 1 × 106 are listed in Table 8. The deflection increases to 143% when the axle load of 10,000 kg increases to 20,000 kg; that is, the surface deflection is sensitive to the change of the axle load. If the asphalt pavement is subjected to the heavy traffic load more, the thicknesses and stiffness of the pavement structural layers should be increased adequately in order to support the surface deflection; otherwise, the destruction of the asphalt pavement will be brought about easier.Table 8
The maximum surface deflections due to different axle loads after loading cycles of 1 × 106.
Axle load (kg)
10,000
12,500
15,000
17,500
20,000
Maximum deflection (mm)
−3.38
−3.73
−4.12
−4.45
−4.84The fatigue life of the asphalt pavement subjected to the different axle loads can be calculated according to (23) [50], which was derived from the indirect tensile test on drill cores extracted from the asphalt base layers based on the German rules and regulations [51]:(23)Nin  situ=Nf·2.8283·ε−4.194,where Nin situ is the number of load cycles under the axle load of 10,000 kg until the macrocrack initiates; Nf is the shift factor, for this case Nf = 1540; ε is the computational tensile strain at the bottom of the asphalt base course.The tensile strain and fatigue life of the asphalt pavement are shown in Table9 and Figure 14. The fatigue life decreases to 48% when the axle load increases from 10,000 kg to 12,500 kg and to 9.4% when the axle load increases from 10,000 kg to 20,000 kg; that is, the increase of the axle load will lead to the destruction of the asphalt pavement extremely easily, especially at the beginning of the axle load increases.Table 9
The tensile strain and fatigue life of the asphalt pavement subjected to different axle loads after loading cycles of 1 × 106.
Axle load (kg)
10,000
12,500
15,000
17,500
20,000
Tensile strain (µm/m)
6200
7400
8700
9700
10,900
Fatigue life
7.90
×
10
12
3.76
×
10
12
1.91
×
10
12
1.21
×
10
12
7.41
×
10
11Figure 14
The fatigue life of the asphalt pavement subjected to different axle loads after loading cycles of 1 × 106.


### 4.3. Impact of the Driving Speed

The impact of different driving speeds was investigated in this section. The vertical displacement was also derived from the center of the loading area 1 on the pavement surface. The variation of the vertical displacement with the increase of driving speed is shown in Figure15(a). For each driving speed, the vertical displacement increases with the increase of loading cycles linearly; with the increase of driving speed, the vertical displacement decreases, and the downtrend becomes less. The variation of the vertical displacement after loading cycles of 1 × 106 along the transverse direction Z of the pavement is shown in Figure 15(b).Figure 15
Comparison of the vertical displacement under different driving speeds. (a) Variation under different loading cycles and (b) variation along the transverse direction Z of the pavement.
(a)
(b)The maximum surface deflections after loading cycles of 1 × 106 are listed in Table 10. The deflection decreases to 84%, 75%, 68%, and 64% when the driving speed of 60 km/h increases to 80 km/h, 100 km/h, 120 km/h, and 140 km/h; that is, the creep deformation is sensitive to the change of relatively lower driving speed, and when the driving speed reaches a certain high level, the creep deformation rate is approaching constant. Common phenomena in daily life are that the rutting is seldom seen on highway, but near the bus stop, the rutting is extremely serious.Table 10
The maximum surface deflections due to different driving speeds after loading cycles of 1 × 106.
Driving speed (km/h)
60
80
100
120
140
Maximum deflection (mm)
−3.38
−2.84
−2.52
−2.30
−2.15The tensile strain and fatigue life of the asphalt pavement are shown in Table11 and Figure 16. The fatigue life increases linearly with the driving speed increases from 60 km/h to 140 km/h.Table 11
The tensile strain and fatigue life of the asphalt pavement subjected to different driving speeds after loading cycles of 1 × 106.
Driving speed (km/h)
60
80
100
120
140
Tensile strain (µm/m)
6200
5200
4500
4100
3800
Fatigue life
7.90
×
10
12
1.65
×
10
13
3.03
×
10
13
4.48
×
10
13
6.16
×
10
13Figure 16
The fatigue life of the asphalt pavement subjected to different driving speeds after loading cycles of 1 × 106.


## 4.1. Definition of the Models

The thicknesses and the material properties of the pavement layers are listed in Table6, which is a typical pavement structure adopted in Germany [47]. The material properties are derived from the defined temperature of 20°C. The length and width of all layers were also set to 6000 mm. The upper three asphalt layers were totally bound. The two contact layers among the asphalt base course, gravel base layer, frost protection layer, and subgrade were defined as being partially bound.Table 6
Geometry and material properties of the pavement.
Layer
Thickness (mm)
E (MPa)
µ
Surface course
40
See Table2
0.35
Binder course
50
See Table2
0.35
Asphalt base course
110
See Table2
0.35
Road base course
150
157.8
0.49
Subbase
570
125.7
0.49
Subgrade
2000
98.9
0.49A standard single axle with dual wheels was utilized in this study, and its configuration can be seen in Figure9.Figure 9
The configuration of the single axle with dual wheels.Considering different working conditions, five different levels of axle loads were applied. The corresponding loading parameters are shown in Table7.Table 7
Loading parameters of the asphalt pavement [48].
Axle load,P (kg)
Contact pressure,p (MPa)
Contact area,A (mm2)
Wheel width,W (mm)
Wheel length,L (mm)
Distance between two wheel centers
d (mm)
D (mm)
10,000
0.70
35,256
156
226
320
1600
12,500
0.76
40,172
166
242
320
1600
15,000
0.82
45,056
176
256
320
1600
17,500
0.85
50,220
186
270
320
1600
20,000
0.89
55,185
195
283
320
1600When moving loads pass for a large number of times, the pavement can be assumed to undertake a series of pulse loads [49], as shown in Figure 10(a). However, when the viscoelasticity is considered in the material property, the pulse loads are extremely time-consuming due to large numbers of increments and/or iterations used in the computation. Former investigation [49] showed the model in which an equivalent wheel load can reasonably predict the creep deformation of the pulse loads while its computational time is significantly shorter. As a result, the equivalent loading mode was applied in this study, as shown in Figure 10(b).Figure 10
The sketch of two loading modes. (a) Pulse loading mode and (b) equivalent loading mode.
(a)
(b)The loading timeTl for each pulse loading is(21)Tl=Wv,where W is the wheel width and v is the driving speed.Equivalent loading is obtained by substituting pulse loading by a one-step equivalent loading whose duration is equal to the summation of loading times in pulse loading.(22)T=N·Tl,where N is the number of loading cycles.The contact pressure of 0.70 MPa with five different driving speeds, such as 60, 80, 100, 120, and 140 km/h, was applied in this study to investigate the impact of vehicle speed on the creep deformation of asphalt pavements.The automatically generated mesh from SAFEM is shown in Figure11.Figure 11
Mesh automatically generated from SAFEM.


## 4.2. Impact of the Axle Loads

When the driving speed is 60 km/h and the contact pressure is 0.70 MPa, the computational vertical displacement after loading cycles of 1 × 106 is shown in Figure 12. This cross section is across the centers of loading areas. The creep deformation can be visualized from this figure, most of which occurs close to the loading areas and decreases rapidly with the increase of distance from the loading areas. As a result, it is meaningful to investigate the near-surface pavement response with more comprehensive model in the future study, which highlights the importance of high efficiency of SAFEM.Figure 12
The computational vertical displacement (mm) with a scale factor of 10.The impact of different axle loads was investigated. The vertical displacement was derived from the center of the loading area 1 on the pavement surface. The variation of the vertical displacement with the increase of loading cycles is shown in Figure13(a). For each axle load, the vertical displacement increases with the increase of loading cycles linearly.Figure 13
Comparison of the vertical displacement under different axle loads. (a) Variation with the increase of loading cycles and (b) variation along the transverse directionZ of the pavement.
(a)
(b)The variation of the vertical displacement after loading cycles of 1 × 106 along the transverse direction Z of the pavement is shown in Figure 13(b). The most deformation occurs below the axle loading area and increases with the increase of axle loads.The maximum surface deflections after loading cycles of 1 × 106 are listed in Table 8. The deflection increases to 143% when the axle load of 10,000 kg increases to 20,000 kg; that is, the surface deflection is sensitive to the change of the axle load. If the asphalt pavement is subjected to the heavy traffic load more, the thicknesses and stiffness of the pavement structural layers should be increased adequately in order to support the surface deflection; otherwise, the destruction of the asphalt pavement will be brought about easier.Table 8
The maximum surface deflections due to different axle loads after loading cycles of 1 × 106.
Axle load (kg)
10,000
12,500
15,000
17,500
20,000
Maximum deflection (mm)
−3.38
−3.73
−4.12
−4.45
−4.84The fatigue life of the asphalt pavement subjected to the different axle loads can be calculated according to (23) [50], which was derived from the indirect tensile test on drill cores extracted from the asphalt base layers based on the German rules and regulations [51]:(23)Nin  situ=Nf·2.8283·ε−4.194,where Nin situ is the number of load cycles under the axle load of 10,000 kg until the macrocrack initiates; Nf is the shift factor, for this case Nf = 1540; ε is the computational tensile strain at the bottom of the asphalt base course.The tensile strain and fatigue life of the asphalt pavement are shown in Table9 and Figure 14. The fatigue life decreases to 48% when the axle load increases from 10,000 kg to 12,500 kg and to 9.4% when the axle load increases from 10,000 kg to 20,000 kg; that is, the increase of the axle load will lead to the destruction of the asphalt pavement extremely easily, especially at the beginning of the axle load increases.Table 9
The tensile strain and fatigue life of the asphalt pavement subjected to different axle loads after loading cycles of 1 × 106.
Axle load (kg)
10,000
12,500
15,000
17,500
20,000
Tensile strain (µm/m)
6200
7400
8700
9700
10,900
Fatigue life
7.90
×
10
12
3.76
×
10
12
1.91
×
10
12
1.21
×
10
12
7.41
×
10
11Figure 14
The fatigue life of the asphalt pavement subjected to different axle loads after loading cycles of 1 × 106.


## 4.3. Impact of the Driving Speed

The impact of different driving speeds was investigated in this section. The vertical displacement was also derived from the center of the loading area 1 on the pavement surface. The variation of the vertical displacement with the increase of driving speed is shown in Figure15(a). For each driving speed, the vertical displacement increases with the increase of loading cycles linearly; with the increase of driving speed, the vertical displacement decreases, and the downtrend becomes less. The variation of the vertical displacement after loading cycles of 1 × 106 along the transverse direction Z of the pavement is shown in Figure 15(b).Figure 15
Comparison of the vertical displacement under different driving speeds. (a) Variation under different loading cycles and (b) variation along the transverse direction Z of the pavement.
(a)
(b)The maximum surface deflections after loading cycles of 1 × 106 are listed in Table 10. The deflection decreases to 84%, 75%, 68%, and 64% when the driving speed of 60 km/h increases to 80 km/h, 100 km/h, 120 km/h, and 140 km/h; that is, the creep deformation is sensitive to the change of relatively lower driving speed, and when the driving speed reaches a certain high level, the creep deformation rate is approaching constant. Common phenomena in daily life are that the rutting is seldom seen on highway, but near the bus stop, the rutting is extremely serious.Table 10
The maximum surface deflections due to different driving speeds after loading cycles of 1 × 106.
Driving speed (km/h)
60
80
100
120
140
Maximum deflection (mm)
−3.38
−2.84
−2.52
−2.30
−2.15The tensile strain and fatigue life of the asphalt pavement are shown in Table11 and Figure 16. The fatigue life increases linearly with the driving speed increases from 60 km/h to 140 km/h.Table 11
The tensile strain and fatigue life of the asphalt pavement subjected to different driving speeds after loading cycles of 1 × 106.
Driving speed (km/h)
60
80
100
120
140
Tensile strain (µm/m)
6200
5200
4500
4100
3800
Fatigue life
7.90
×
10
12
1.65
×
10
13
3.03
×
10
13
4.48
×
10
13
6.16
×
10
13Figure 16
The fatigue life of the asphalt pavement subjected to different driving speeds after loading cycles of 1 × 106.


## 5. Summary and Conclusion

The linear viscoelasticity is applied in the SAFEM by using a recursive formulation to predict the mechanical responses of asphalt pavements. The accuracy of the program is verified by a comparison with ABAQUS. The analytical verification showed that the pavement responses derived from SAFEM and ABAQUS are in accordance with one another. It should be emphasized that the computational time of the SAFEM is much shorter than that of the ABAQUS. The analysis of creep deformation of asphalt pavements is carried out, and some laws of the impact of different axle loads and driving speeds are revealed.The current code of SAFEM provides a flexible, robust base for further development. For the next step of the development, more various material properties, such as nonlinear elasticity for the subbase of the pavement, would be applied. In order to investigate the crack initiation and propagation, the cohesive elements would be coupled with the finite element model. With these improvements, the SAFEM should be more suited to predict the mechanical performan ces of the asphalt pavement.



---

*Source: 9045820-2018-02-08.xml*