Advances in Civil Engineering
Volume 2018 (2018), Article ID 9045820, 15 pages
https://doi.org/10.1155/2018/9045820

Research Article

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

Pengfei Liu,¹ Qinyan Xing,² Dawei Wang,1,3 and Markus Oeser¹

¹Institute of Highway Engineering, RWTH Aachen University, 52074 Aachen, Germany
²Department of Civil Engineering, Tsinghua University, Beijing 100084, China
³School of Transportation Science and Engineering, Harbin Institute of Technology, Harbin 150090, China

Correspondence should be addressed to Dawei Wang; wang@isac.rwth-aachen.de

Received 1 August 2017; Accepted 9 October 2017; Published 8 February 2018

Academic Editor: Ghassan Chehab

Copyright © 2018 Pengfei Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 Figure 1, 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 Figure 2 [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:

where E₁, E₂ and η₁, η₂ 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:

The corresponding relaxation function is the following:

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 which z ranges between zero and a, as shown in Figure 4 [18–27, 39].

where l identifies the term of the Fourier series and L is the number of terms considered; and 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 the z-direction is given a similar form as the displacements [19]:

where and represent the external load function.

The pavement is assumed to be held at z = 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 can be written as follows [39]:

where , , and 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]

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 matrix is defined as follows:

By using the principle of minimum potential energy, a typical submatrix of the element stiffness matrix is [39]

A typical term for the force vector becomes

From (9) and (10), the stiffness matrix of one element includes [39]

The integrals exhibit orthogonal properties which ensure that

Only when l and m are both odd or even numbers, the first integral I₁ is zero. Due to the special structure of the Bl matrix, all terms that include I₁ vanish (become zero). This means that the matrix 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]:

The last equation shows that the large system of equations splits up into L separate problems.

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

where G and K are the shear modulus and bulk modulus, respectively. is the deviatoric strain, and is the spherical strain. is the deviatoric stress, and is the spherical stress. 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:

where G₁ and G₂ are the shear modulus of the elastic element in Maxwell and Kelvin elements, respectively. K₁ and K₂ are the bulk modulus of the elastic element in Maxwell and Kelvin elements, respectively. and are the viscosity coefficient of the viscous element for deviatoric strain response in Maxwell and Kelvin elements, respectively. and 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 ratio v 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:

where and are Poisson’s ratio of the viscous element for deviatoric strain response in Maxwell and Kelvin elements, respectively. and 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:

For implementation in the SAFEM algorithm, the incremental deviatoric and spherical strains are written in terms of history integrals as follows:

The variables 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 Figure 5.

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 Table 1. 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 Tables 2 and 3
0.35

Binder course
60
See Tables 2 and 3
0.35

Asphalt base course
70
See Tables 2 and 3
0.35

Road base course
250
1500
0.49

Subbase
250
1200
0.49

Subgrade
5000
50
0.49

Table 2: Material parameters of the Burgers model used in SAFEM.

Layer
E₁ (MPa)
η₁ (MPa·s)
E₂ (MPa)
η₂ (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.1

Table 3: Prony series parameters used in ABAQUS.

Layer
g₁
g₂
τ₁ (s)
τ₂ (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.16

The meshes generated from the SAFEM and ABAQUS are shown in Figure 6. A full-size model was created in the SAFEM, while a half symmetrical model was in the ABAQUS.

(a)

(b)

Figure 6: Mesh automatically generated from (a) SAFEM and (b) ABAQUS in static analysis with linear viscoelastic material property.

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⁻⁸ 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 Figure 7. 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.

(a)

(b)

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.

The computational vertical stresses and strains at the center of the pavement surface are shown in Figure 8. Both computational results are in good agreement with one another. The stresses increase from 0 to 0.7 MPa in 1 × 10⁻⁸ 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.

(a)

(b)

Figure 8: Comparison of (a) vertical stresses and (b) vertical strains at the center of the pavement surface between ABAQUS and SAFEM.

In Table 4, 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⁻⁸
−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.797

The 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 Table 5. 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 Table 6, 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 Table 2
0.35

Binder course
50
See Table 2
0.35

Asphalt base course
110
See Table 2
0.35

Road base course
150
157.8
0.49

Subbase
570
125.7
0.49

Subgrade
2000
98.9
0.49

A standard single axle with dual wheels was utilized in this study, and its configuration can be seen in Figure 9.

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 Table 7.

Table 7: Loading parameters of the asphalt pavement [48].

Axle load, P (kg)
Contact pressure,p (MPa)
Contact area, A (mm²)
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
1600

When 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).

(a)

(b)

Figure 10: The sketch of two loading modes. (a) Pulse loading mode and (b) equivalent loading mode.

The loading time for each pulse loading is

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.

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 Figure 11.

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 × 10⁶ 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 Figure 13(a). For each axle load, the vertical displacement increases with the increase of loading cycles linearly.

(a)

(b)

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 direction Z of the pavement.

The variation of the vertical displacement after loading cycles of 1 × 10⁶ 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 × 10⁶ 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 × 10⁶.

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.84

The 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]:

where Nin situ is the number of load cycles under the axle load of 10,000 kg until the macrocrack initiates; is the shift factor, for this case  = 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 Table 9 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 × 10⁶.

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

Figure 14: The fatigue life of the asphalt pavement subjected to different axle loads after loading cycles of 1 × 10⁶.

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 Figure 15(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 × 10⁶ along the transverse direction Z of the pavement is shown in Figure 15(b).

(a)

(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.

The maximum surface deflections after loading cycles of 1 × 10⁶ 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 × 10⁶.

Driving speed (km/h)
60
80
100
120
140

Maximum deflection (mm)
−3.38
−2.84
−2.52
−2.30
−2.15

The tensile strain and fatigue life of the asphalt pavement are shown in Table 11 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 × 10⁶.

Driving speed (km/h)
60
80
100
120
140

Tensile strain (µm/m)
6200
5200
4500
4100
3800

Fatigue life

Figure 16: The fatigue life of the asphalt pavement subjected to different driving speeds after loading cycles of 1 × 10⁶.

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.

Disclosure

The authors are solely responsible for the content.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This paper is based on parts of the research project carried out at the request of the Federal Ministry of Transport and Digital Infrastructure, requested by the Federal Highway Research Institute, under research project no. 04.0259/2012/NGB, as well as parts of the research project carried out at the request of the German Research Foundation, under research project no. FOR2089.

References

A. E. Zeevaert, “Finite element formulation for the analysis of interfaces, nonlinear and large displacement problems in geotechnical engineering,” Georgia Institute of Technology, Atlanta, GA, USA, 1980, Ph.D. Dissertation.

R. D. Barksdale, Q. Robnett, J. Lai, and A. E. Zeevaert, “Experimental and theoretical behavior of geotextile reinforced aggregate soil systems,” in Proceedings of 2nd International Conference of Geotextiles, vol. 2, pp. 375–380, Las Vegas, NV, USA, August 1982.

L. Raad and J. L. Figueroa, “Load response of transportation support systems,” Transportation Engineering Journal of ASCE, vol. 16, no. 1, pp. 111–128, 1980.

R. S. Harichandran, M. S. Yeh, and G. Y. Baladi, “MICH-PAVE: a nonlinear finite element program for analysis of pavements,” TRB, National Research Council, Washington, DC, USA, 1989, Transportation Research Record 1286.

J. M. Brunton and J. R. De Almeida, “Modeling material nonlinearity in a pavement backcalculation procedure,” TRB, National Research Council, Washington, DC, USA, 1992, Transportation Research Record 1377.

S. F. Brown, “The characterisation of cohesive soils for flexible pavement design,” in Proceedings of 7th European Conference on Soil Mechanics and Foundation Engineering, vol. 2, pp. 15–22, Brighton, UK, September 1979.

S. F. Brown, S. C. Loach, and M. P. O’Reilly, “Repeated loading of fine grained soils,” U.K. Transport and Road Research Laboratory, Crowthome, Berkshire, England, 1987, Contractors Report CR 72.

R. G. Hicks and C. L. Monismith, “Factors influencing the resilient response of granular materials,” TRB, National Research Council, Washington, DC, USA, 1971, Highway Research Record 345.

A. Gonzalez, D. Bodin, D. Jameson, M. Oeser, and B. Vuong, Development of a Nonlinear Finite Element Pavement Response to Load Model, Austroads Publication No. AP–T199–12, Sydney, Australia, 2012.

A. Scarpas, Ed., “CAPA-3D: a mechanics based computational platform for pavement engineering,” Delft University of Technology, Delft, Netherlands, 2005, Ph.D. dissertation.

R. L. Lytton, F. L. Tsai, S. I. Lee, R. Luo, S. Hu, and F. Zhou, “Models for predicting reflection cracking of hot-mix asphalt overlays,” Transportation Research Board, Washington, DC, USA, 2010, NCHRP Report 669.

CESAR-LCPC, “Cesar-LCPC finite element package,” June 2016, http://www.cesar-lcpc.com/.

E. L. Wilson, “Structural analysis of axisymmetric solids,” AIAA Journal, vol. 3, no. 12, pp. 2269–2274, 1965.

H. E. Meissner, “Laterally loaded pipe pile in cohesionless soil,” in Proceedings of 2nd International Conference on Numeric Methods in Geomechanics, vol. 1, Virginia Polytechnical Institute and State University, Blacksburg, VA, USA, June 1976.

L. A. Winnicki and O. C. Zienkiewicz, “Plastic (or visco-plastic) behaviour of axisymmetric bodies subjected to non-symmetric loading-semi-analytical finite element solution,” International Journal for Numerical Methods in Engineering, vol. 14, no. 9, pp. 1399–1412, 1979.

J. P. Carter and J. R. Booker, “Consolidation of axi-symmetric bodies subjected to non axi-symmetric loading,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 7, no. 2, pp. 273–281, 1983.

J. Y. Lai and J. R. Booker, “Application of discrete Fourier series to the finite element stress analysis of axi-symmetric solids,” International Journal for Numerical Methods in Engineering, vol. 31, no. 4, pp. 619–647, 1991.

J. J. Fritz, “Flexible pavement response evaluation using the semi-analytical finite element method,” International Journal of Materials and Pavement Design, vol. 3, no. 2, pp. 211–225, 2002.

S. Hu, X. Hu, and F. Zhou, “Using Semi-Analytical Finite Element Method to Evaluate Stress Intensity Factors in Pavement Structure,” Pavement Cracking, CRC Press, Boca Raton, FL, USA, 2008.

P. Liu, D. Wang, and M. Oeser, Eds., “Leistungsfähige semi-analytische methoden zur berechnung von asphaltbefestigungen,” 2013, Tangungsband, 3. Dresdner Asphalttagein German.

P. Liu, D. Wang, M. Oeser, and X. Chen, “Einsatz der semi-analytischen finite-elemente-methode zur beanspruchungszustände von asphaltbefestigungen,” Bauingenieur, vol. 89, no. 7-8, pp. 333–339, 2014, in German.

P. Liu, D. Wang, F. Otto, J. Hu, and M. Oeser, “Application of semi-analytical finite element method to evaluate asphalt pavement bearing capacity,” International Journal of Pavement Engineering, pp. 1–10, 2016.

P. Liu, D. Wang, J. Hu, and M. Oeser, “SAFEM–software with graphical user interface for fast and accurate finite element analysis of asphalt pavements,” Journal of Testing and Evaluation, vol. 45, no. 4, pp. 1301–1315, 2017.

P. Liu, F. Otto, D. Wang, H. Balck, and M. Oeser, “Measurement and evaluation on deterioration of asphalt pavements by geophones,” Measurement, vol. 109, pp. 223–232, 2017.

P. Liu, Q. Xing, Y. Dong, D. Wang, M. Oeser, and S. Yuan, “Application of finite layer method in pavement structural analysis,” Applied Sciences, vol. 7, no. 6, p. 611, 2017.

P. Liu, Q. Xing, D. Wang, and M. Oeser, “Application of dynamic analysis in semi-analytical finite element method,” Materials, vol. 10, no. 9, p. 1010, 2017.

P. Liu, D. Wang, and M. Oeser, “Application of semi-analytical finite element method coupled with infinite element for analysis of asphalt pavement structural response,” Journal of Traffic and Transportation Engineering, vol. 2, no. 1, pp. 48–58, 2015.

E. V. Dave, G. H. Paulino, and W. G. Buttlar, “Viscoelastic functionally graded finite element method with recursive time integration and applications to flexible pavements,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 36, no. 9, pp. 1194–1219, 2012.

H. H. Hilton and S. Yi, “Anisotropic viscoelastic finite element analysis of mechanically and hygrothermally loaded composites,” Composites Engineering, vol. 3, no. 2, pp. 123–135, 1993.

S. Mukherjee and G. H. Paulino, “The elastic-viscoelastic correspondence principle for functionally graded materials, revisited,” Journal of Applied Mechanics, vol. 70, no. 3, pp. 359–363, 2003.

K. R. Rajagopal and A. S. Wineman, “A quasi-correspondence principle for quasi-linear viscoelastic solids,” Mechanics of Time-Dependent Materials, vol. 12, no. 1, pp. 1–14, 2008.

O. C. Zienkiewicz, M. Watson, and I. P. King, “A numerical method of visco-elastic stress analysis,” International Journal of Mechanical Sciences, vol. 10, no. 10, pp. 807–827, 1968.

R. L. Taylor, K. S. Pister, and G. L. Goudreau, “Thermomechanical analysis of viscoelastic solids,” International Journal for Numerical Methods in Engineering, vol. 2, no. 1, pp. 45–59, 1970.

A. Muliana and K. A. Khan, “A time-integration algorithm for thermo-rheologically complex polymers,” Computational Materials Science, vol. 41, no. 4, pp. 576–588, 2008.

M. A. Zocher, S. E. Groves, and D. H. Allen, “A three-dimensional finite element formulation for thermoviscoelastic orthotropic media,” International Journal for Numerical Methods in Engineering, vol. 40, no. 12, pp. 2267–2288, 1997.

S. Yi and H. H. Hilton, “Dynamic finite element analysis of viscoelastic composite plates in the time domain,” International Journal for Numerical Methods in Engineering, vol. 37, no. 23, pp. 4081–4096, 1994.

M. Arabani and N. Kamboozia, “The linear visco-elastic behaviour of glasphalt mixture under dynamic loading conditions,” Construction and Building Materials, vol. 41, pp. 594–601, 2013.

M. Yazarov, One-Dimensional Viscoelastic Simulation of Ice Behaviour in Relation to Dynamic Ice Action, TU Delft, Delft University of Technology, Delft, Netherlands, 2012.

O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics, Elsevier Butterworth-Heinemann, Oxford, UK, 6th edition, 2005.

ASTM, Ed., Standard Test Method for Indirect Tension Test for Resilient Modulus of Bituminous Mixtures, Test Procedure D 4123, American Society for Testing Materials, West Conshohocken, PA, USA, 1995.

H. Di Benedetto, N. Mondher, C. Sauzeat, and F. Olard, “Three-dimensional thermo-viscoplastic behavior of bituminous materials,” Road Material and Pavement Design, vol. 8, no. 2, pp. 285–316, 2007.

M. Oeser, “Non-linear numerical simulation models for pavement structures under consideration of mechanical, thermal and hydraulic influences/nichtlineare numerische simulationsmodelle für verkehrswegebefestigungen unter berücksichtigung von mechanischen, thermischen und hydraulischen einwirkungen habilitation,” TU Dresden, Veröffentlichungen des Lehrstuhls für Statik, Dresden, Germany, 2010, Professorial-Thesis, Habilitation.

F. Wang, “A 3D FEM analysis of creep deformation law of asphalt pavement based on Burgers viscoelastic model,” Industrial Construction, vol. 44, no. 3, pp. 104–109, 2014, in Chinese.

J. Kong and J. Yuan, “Application of linear viscoelastic differential constitutive equation in ABAQUS,” in Proceedings of the International Conference on Computer Design and Applications (ICCDA 2010), vol. 5, pp. 152–156, Qinhuangdao, China, June 2010.

L. A. Al-Khateeb, A. Saoud, and M. F. Al-Msouti, “Rutting prediction of flexible pavements using finite element modeling,” Jordan Journal of Civil Engineering, vol. 5, no. 2, pp. 173–190, 2011.

A. Szydlo and P. Mackiewicz, “Model for forecasting ruts in rutting tester,” in Proceedings of the 6th International RILEM Symposium Performance Testing and Evaluation of Bituminous Materials (PTEBM’03), pp. 14–16, Zurich, Switzerland, April 2003.

FGSV, Guidelines for the Standardization of the Upper Structure of Traffic Areas, Research Society for Road and Transportation, FGSV Publisher, Cologne, Germany, 2001 edition, 2002, in German.

L. Yu, J. Yang, L. Bao, and F. Yan, “Analysis on the mechanical responds of asphalt pavement with subbase used fly-ash-flushed by seawater under heavy traffic,” Journal of Shenyang Jianzhu University, vol. 29, no. 1, pp. 110–115, 2013, in Chinese.

R. K. A. Al-Rub, M. K. Darabi, C.-W. Huang, E. A. Masad, and D. N. Little, “Comparing finite element and constitutive modelling techniques for predicting rutting of asphalt pavements,” International Journal of Pavement Engineering, vol. 13, no. 4, pp. 322–338, 2012.

R. Leutner, H. Lorenzl, K. Schmoeckel et al., Stoffmodelle zur Voraussage des Verformungswiderstandes und Ermüdungsverhaltens von Asphaltbefestigungen, Berichte der Bundesanstalt Fuer Strassenwesen. Unterreihe Strassenbau, Bergisch Gladbach, Germany, 2006, in German.

FGSV 2009, Arbeitsanleitung zur Bestimmung des Steifigkeits- und Ermüdungsverhaltens von Asphalten mit dem Spaltzug-Schwellversuch als Eingangsgröße in die Dimensionierung, AL Sp-Asphalt 09. FGSV Verlag, Köln, Germany, 2009, in German.