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. [2027]. 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.

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 [3436]. 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].

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

This model has the following response equation: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:

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 [1827, 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.

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 [2027].

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

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.

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.

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

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

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.

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.

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.

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

Considering different working conditions, five different levels of axle loads were applied. The corresponding loading parameters are shown in Table 7.

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

The loading time for each pulse loading iswhere 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.

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.

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.

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.

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.

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

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.

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.