Application of Dynamic Analysis in Semi-Analytical Finite Element Method.

Analyses of dynamic responses are significantly important for the design, maintenance and rehabilitation of asphalt pavement. In order to evaluate the dynamic responses of asphalt pavement under moving loads, a specific computational program, SAFEM, was developed based on a semi-analytical finite element method. This method is three-dimensional and only requires a two-dimensional FE discretization by incorporating Fourier series in the third dimension. In this paper, the algorithm to apply the dynamic analysis to SAFEM was introduced in detail. Asphalt pavement models under moving loads were built in the SAFEM and commercial finite element software ABAQUS to verify the accuracy and efficiency of the SAFEM. The verification shows that the computational accuracy of SAFEM is high enough and its computational time is much shorter than ABAQUS. Moreover, experimental verification was carried out and the prediction derived from SAFEM is consistent with the measurement. Therefore, the SAFEM is feasible to reliably predict the dynamic response of asphalt pavement under moving loads, thus proving beneficial to road administration in assessing the pavement's state.

1. Introduction

The analysis of loading in finite element (FE) simulation of asphalt pavement can be modelled in two types: dynamic or static analyses. Static analysis is relatively simpler and thus reduces computational time. This loading type has advantages when elastic materials are investigated due to the exactly same loading and unloading path for elastic materials [1]. The early studies by Duncan et al. [2], Raad and Figueroa [3] and Harichandran et al. [4] applied static analysis on asphalt concrete surface, which were followed by more researchers using static analysis to investigate other aspects of asphalt pavement [5,6,7]. However, the inertial forces and time dependency of materials cannot be taken account in the static analysis, because it ignores the force induced by mass and material damping, which is not the reality of the pavement.

Zaghloul and White [8] conducted one of the first studies of dynamic analysis in the general-purpose finite element method (FEM) program ABAQUS in the year of 1993, in which uniformly distributed pressure varying as a trapezoidal shape in time was used to model a moving load. In this study, it was found that when the speed of loading increased, the surface deflection decreased and there was more deflection derived from static loading than from dynamic loading. A dynamic analysis of flexible pavement was carried out by Desai and Whitenack [9], which was further developed by Desai [10]. The series of studies proposed a new material constitutive model to model the distress induced in pavement during dynamic loading. The stored strain energy was used to consider the effect of material changes caused by cycles of loading. Particularly, a reduced dissipation of energy happened in each cycle, leading to stiffer material behavior in the following cycle. Saad et al. [11] assumed the dynamic loading as a triangular pulse in a 0.1 s period and then performed a dynamic analysis of flexible pavement. The researchers pointed out that the results of surface deflection calculated from dynamic loading may be less than half of the value derived from static loading. The reason was believed to be that the absorption of energy by the damping and mass inertia of the pavement system. A comprehensive review of the previous methods to simulate dynamic loading was made by Beskou and Theodorakopoulos [12]. They categorized the dynamic analysis according to representative models for asphalt pavement. The dynamic analysis of a pavement system can be achieved by analytical method, hybrid method of analytical and numerical methods or purely numerical method. Only numerical methods can solve realistic problems with complex material behavior and geometries; the FEM, the boundary element method (BEM) and the BEM/FEM hybrid scheme are the most popular numerical methods. Considering efficiency, availability, versatility and accuracy, the FEM appears to be the best. Most of the layered pavement systems were analyzed by linear formulation with respect to the material behavior and deformation. The dynamic analysis of these models required the development of an innovative efficient and accurate numerical method and its implementation in computer software.

One method was proposed to meet the requirement of accuracy and efficiency, which is the so-called semi-analytical finite element method (SAFEM) [13,14,15,16,17,18,19,20]. This method was first developed by Wilson [21] for linear analysis and Meissner [22] extended Wilson’s work to an elastoplastic body. Winnicki and Zienkiewicz [23] adopted a viscoplastic formulation to deal with material nonlinearity. An efficient analysis of the consolidation of elastic bodies’ under nonsymmetric loading was provided by Carter and Booker [24] through continuous Fourier series. A discrete Fourier technique was successfully applied by Lai and Booker [25] to analyze the nonlinear behavior of solids under three-dimensional (3D) loading conditions. Fritz [13] and Hu et al. [14] further developed the SAFEM, and programed specific FE codes using this method to analyze the problems of asphalt pavement. However, their FE codes are relatively simple, e.g., only static analysis with linear elastic material properties and a total bond between pavement layers can be performed. As a result, the dynamic analysis should be further applied to the SAFEM. In the following sections, the mathematical basis of the SAFEM and the algorithm of application of the dynamic analysis to SAFEM will be introduced, followed by the verification using the results derived from commercial FE software (ABAQUS) and the field measurement. A brief summary and conclusions are provided at the end.

2. Description of Semi-Analytical Finite Element Method

For the sake of simplification, the FE solution for static analysis is introduced first. The shape functions in a 3D SAFEM model, which are to define the variation of displacements, are written as a two-dimensional (2D) traditional shape function multiplied by Fourier series, as shown in Figure 1: where l identifies the term of the Fourier series and is up to L, which is the number of total terms; and are the shape functions of the node in the XY plane.

Figure 1

Schematic representation of an SAFEM situation.

In order to derive the displacements in three dimensions, three degrees of freedom at each node of the triangular finite element should be considered. The vector of nodal displacements is: where , and are the displacements of the node k in x-, y- and z-directions, respectively.

Displacements at some point inside the element can be determined with the use of nodal displacements and shape functions :

The loading function defining the variation of load along the z-direction is given by [15]: where and represent the pavement load.

The pavement is assumed to be supported at both edges (z = 0 and z = a), particularly, all displacements in the XY plane are prevented, while the motion in the z-direction is unrestricted. In order to meet this requirement, the displacement functions with three components u, v and w are rewritten as follows:

Similarly, the function of loading can be formulated as [13]: where P is the tire load pressure; Z1 is the z coordinate where the tire load area starts; Z2 is the z coordinate where the tire load area ends.

Through displacements at nodal points the strains are determined as:

Matrix is called the strain-displacement matrix. is the strain-displacement matrix of the node k at lth term of the Fourier series. It is obtained by differentiation of displacements in the form of shape functions and nodal displacements:

The in Equation (8) is split into two matrices—each includes only one set of trigonometric terms.

Now, the general form of total potential energy can be expressed through nodal displacements:

Nodal displacements corresponding to the minimum of the functional are determined by the following condition:

The following equilibrium equation is produced by differentiation of in respect to for one element: which can be reformulated as follows:

Here [k] is the stiffness matrix of one element; is the vector of load.

From Equations (8) and (12), the stiffness matrix of one element includes the following integrals due to [B] [14]:

The integrals exhibit orthogonal properties and ensure that:

Only when l and m are both odd or even numbers, the first integral I1 is zero. But due to the special structure of the B matrix, all terms that include I1 vanish (become zero). This means that the matrix becomes a diagonal one, i.e., the non-zero values only exist in the diagonal area where l = m, which reduces the stiffness matrix to [14]: where g, k represent the element nodes, respectively. Area is the area of the element.

A typical term for the force vector becomes:

The global linear system is achieved by assembling the elemental stiffness matrix to the global domain:

The large system of equations can be divided into L separate problems: where K is the global stiffness matrix; U is the global displacement vector; F is the global loading vector.

For analysis of dynamic response in asphalt pavement, the time coordinates are introduced to the FE algorithm. The motion equation can be expressed based on Newton’s second law of motion: where , and are global displacement, velocity and acceleration vectors and = 0, = 0 and = 0. M, C and K represent the global mass, damping and stiffness matrices, respectively. is the force vector.

The M, C and K are assembled by element mass, damping and stiffness matrices , and , which are: where ρ, μ and D are the density, damping factor and elastic matrices, respectively.

As the same with static one, these matrices are diagonal and the global equations are as follows:

Thus, the large system also splits up into L separate problems:

According to Equations (16) and (21), the harmonics of the Fourier series are decoupled in the global equations, which is a benefit to the parallel calculation, and thus the computational time can be significantly reduced compared to a sequential solving procedure [26].

3. Time Discretization for Dynamic Analysis

A modal method or a direct integration method is usually used in FE analysis to undertake dynamic response analysis. No transformation to a special form is required in the direct integration method and the governing equation such as Equation (22) is solved step-by-step in the time domain.

One of the integration methods popularly used in the FE program for the solution of structural dynamic problems for both blast and seismic loading is called the Newmark method, which was proposed by Newmark in 1959 [27]. The Newmark method has been applied to the dynamic analysis of many engineering structures in the past several decades. A brief description of this method specialized for the SAFEM is provided here in two approaches: displacement based and acceleration based approach.

3.1. Displacement Based Approach

Within the time step (t, t + ∆t), the standard form of Newmark’s equations should be modified according to SAFEM as: where and are algorithm parameters which are determined by the requirement of integral accuracy and stability. When = , = and = , = it becomes average acceleration method and linear acceleration method, respectively.

From Equation (24), one can derive the following equation:

It is now substituted into Equation (22) to give the result:

The linear dynamic equilibrium equation used in SAFEM is written in the form of Equation (22). The substitution of Equations (25) and (26) into Equation (22) allows the dynamic equilibrium of the system at time “” to be written in terms of the unknown node displacements for lth harmonic of Fourier series: where is effective dynamic stiffness matrix and is effective load vector for lth harmonic of Fourier series.

The Newmark direct integration algorithm based on displacements is summarized as follows [28]: The calculation for the constants needs to be carried out only once. For linear systems, the effective dynamic stiffness matrix is formed and triangularized only once.

Step 1: Initial Calculation

Form static stiffness matrix K, mass matrix M and damping matrix C.

Specify time step and integration parameters , .

Calculate integration constants .

Form effective stiffness matrix .

Step 2: For Each Time Step

Specify initial conditions , and .

Calculate effective load vector .

Solve for node displacement vector at time t according to Equation (27).

Calculate node velocities and accelerations at time t according to Equations (25) and (26).

Go to Step 2 (b) with .

3.2. Acceleration-Based Approach

The unknown node accelerations for lth harmonic of Fourier series can be derived from substitution of Equations (23) and (24) into Equation (22) directly:

In this approach, the incremental equations of motion are solved for incremental acceleration first, and the incremental velocity and displacement are calculated from this incremental acceleration according to Equations (23) and (24). The numerical integration algorithm based on this scheme approach is similar to that based on displacements.

In the code of SAFEM, the acceleration based approach with the average acceleration method is adopted.

4. Analytical Verification of SAFEM in Dynamic Analyses

The accuracy of dynamic analyses under moving loads using SAFEM was verified by comparing the results with the data derived from ABAQUS (Abaqus 6.14-1, Johnston, RI, USA). The responses from both models were evaluated with the pavement type in Table 1, which is widely used in Germany. The thicknesses of all layers, excluding the subgrade, were derived from the guideline RStO 12 [29]. The thickness of the subgrade was defined as 2000 mm; setting such a large value aims to minimize the influence of the boundary condition on the results. The length and width of all layers were set to 6000 mm for a similar reason. Homogeneous E-moduli in each layer were defined according to the guideline RDO Asphalt 09 [30] for assumed pavement surface temperatures of 27.5 C which is the typical temperature in the summer in Germany.

Table 1

Geometrical data and material properties of the pavement in dynamic analysis under moving loads.

Layer	Thickness (mm)	E (MPa)	µ	Density (t/mm3)
Surface course	40	22,690	0.35	2.377 × 10−9
Binder course	80	27,283	0.35	2.448 × 10−9
Asphalt base course	140	17,853	0.35	2.301 × 10−9
Road base course	150	10,000	0.25	2.400 × 10−9
Sub-base	340	100	0.49	2.400 × 10−9
Subgrade	2000	45	0.49	2.400 × 10−9

The 2D mesh generated from SAFEM was as shown in Figure 2a. A 3D 10-node quadratic tetrahedron element was applied in ABAQUS due to its three-dimensional discretization. In order to save computational time of the analysis, only a half-symmetrical model was created in ABAQUS, whereas a full size model was created in SAFEM. The model with the mesh algorithm in ABAQUS is illustrated in Figure 2b. Due to the different dimensions of both models, the mesh algorithms are different, while the general mesh regulation is the same in both models, i.e., the element size is gradually increased from the loading area to periphery.

Figure 2

Mesh automatically generated from (a) SAFEM and (b) ABAQUS.

The load in SAFEM and ABAQUS was assumed as a square load and its side length is 300 mm and the speed is 52 km/h. The uniformly distributed contact pressure was 0.7 MPa. The load path traverses the center of the pavement, as shown in Figure 2. The element size in ABAQUS is 60 mm in the direction of traffic and 50 mm in the transverse direction, as shown in Figure 2b. The load moved forward one element in each increment by default; therefore there were 100 increments in both analyses. It should be mentioned that currently the loading mode of moving loads in ABAQUS is not offered in its graphical user interface (GUI) [31], i.e., the user has to write a subroutine to realize this function which is difficult for common pavement engineers.

The bottom nodes of the mesh representing the subgrade in both models were fixed in all directions. On both edges (z = 0 and z = a), the displacements are restricted in the x- and y-directions due to the theory of SAFEM. 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, sub-base and subgrade were defined as being partially bound.

The computational vertical displacements obtained from both models when the load is at the center of the pavement (the 50th increment) are illustrated in a Moiré pattern, as shown in Figure 3. The cross-section is through the centroid of the full size pavement model and along the traffic direction. It can be seen that the distribution of the displacement, magnitude and the deformation shapes from both FE programs are in good agreement with one other.

Figure 3

The computational vertical displacement from (a) SAFEM and (b) ABAQUS with a scale factor of 500 in dynamic analysis under moving loads.

The computational results from four response points in the loading history are shown in Figure 4.

Figure 4

Comparison of the results between ABAQUS and SAFEM. (a) Vertical displacement at the top of the surface course; (b) Horizontal stress at the bottom of the asphalt base course; (c) Vertical stress at the top of the sub-base course; (d) Vertical stress at the top of the subgrade.

The corresponding results from Figure 4 at a loading time of 0.2 s are listed in Table 2. As seen in Figure 4 and Table 2, it can be stated that the results obtained from both programs have a high correlation considering that their mesh algorithms are totally different due to the different dimensions of discretization.

Table 2

Comparison between ABAQUS and SAFEM regarding the computational results at critical points when the loading time is 0.2 s.

Result	SAFEM	ABAQUS	Difference
Vertical displacement (mm) at the top of the surface course	−0.335	−0.327	2.44%
Horizontal stress (MPa) at the bottom of the asphalt base course	0.974	0.912	6.79%
Vertical stress (MPa) at the top of the sub-base course	−0.0104	−0.0112	−7.14%
Vertical stress (MPa) at the top of the subgrade	−0.0100	−0.0109	−8.26%

Benefit from the use of a Fourier series in the transverse direction the SAFEM needs significantly fewer elements and nodes compared with ABAQUS; this results in a far shorter computational time for the SAFEM as compared to ABAQUS. A computer with the configuration of an Intel Core Duo 3.4 GHz and 32 GB RAM was used to run the two FE analyses. Generally, the computational time required by the half-symmetrical 3D model in ABAQUS is about 281 min, whereas the full-size pavement model in SAFEM requires 10 min, which is only 3.6% of that used by ABAQUS, as shown in Table 3.

Table 3

Comparison of the efficiency between SAFEM and ABAQUS.

	SAFEM	ABAQUS
Elements	1144	127,095
Nodes	2431	218,333
Computational time	10 min	281 min

5. Experimental Verification of SAFEM in Dynamic Analyses

A test track in the German Federal Highway Research Institute (BASt) was used in this study to experimentally verify the SAFEM, as shown in Figure 5. The test track is comprised of five layers, which are asphalt surface course, asphalt binder course, asphalt base course, frost protection layer and subgrade with the thickness of 40, 50, 130, 680, and 1440 mm, respectively. Strain gauges were placed at the bottom of the asphalt base course during the test track construction, which is to measure strains along the traffic direction for the verification. A passing truck with speed of 30 km/h was used to apply the moving loads in the experiment. The geometry and the distribution of the loads are illustrated in Figure 6. The left wheels of five axles were considered in the verification in order to simplify the model and thus reduce the computational time. The specimens were drilled from the test track to obtain the material properties and then the material parameters were used in SAFEM. More details can be found from [19,32].

Figure 5

The test track in German Federal Highway Research Institute (BASt).

Figure 6

Geometrical and loading data of the truck.

When the second wheel load was exactly above the locations where the sensors were embedded the computational strains from SAFEM were compared with the field measurement, as shown in Table 4. The computational strain at the bottom of the asphalt base course is higher than the measured value with a difference of 5.88%. The measured values can be accepted in a range of 20% due to the uncertainties and fluctuations [19], the computational strains are therefore in a quite good agreement with the measurement.

Table 4

Comparison of the strains derived from measurement and SAFEM.

	Measurement	SAFEM	Difference
Strain along the traffic direction at the bottom of the asphalt base course (10−6)	81.5	86.3	5.88%

6. Summary and Conclusions

This paper proposes to use the SAFEM for predicting the asphalt pavement dynamic responses under moving loads. The accuracy of the program is analytically and experimentally verified by comparison with ABAQUS and field measurements, respectively. The results predicted by SAFEM and ABAQUS are generally consistent with each other. Furthermore, the efficiency of the SAFEM is much higher than that of the ABAQUS. The computational strain derived from SAFEM at the critical point is quite close to the field measurement, which further proves the accuracy of the developed program.

In conclusion, the SAFEM has potential to reliably analyze the dynamic response of asphalt pavement under moving loads. For further investigation, the SAFEM allows the application of various material properties, such as viscoelasticity for asphalt and nonlinear elasticity for the sub-base of the pavement. With these improvements, the SAFEM should be more appropriate to predict the impact of the moving load on the asphalt pavement.