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.


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 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; N k and = N k are the shape functions of the node in the XY plane. 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: , k = 1, 2, 3 . . . 6 (2) where u k , v k and w k are the displacements of the node k in x-, yand z-directions, respectively. at the end.

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:  Displacements at some point inside the element can be determined with the use of nodal displacements {d k } and shape functions N k : The loading function defining the variation of load along the z-direction is given by [15]: where p(x, y) and = p(x, y) 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 t is the tire load pressure; Z t1 is the z coordinate where the tire load area starts; Z t2 is the z coordinate where the tire load area ends. Through displacements at nodal points the strains are determined as: Matrix B is called the strain-displacement matrix. B l k 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 B l k 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 {d} 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 {d} for one element: which can be reformulated as follows: (12) Here [k] is the stiffness matrix of one element; { f } 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 I 1 is zero. But due to the special structure of the B matrix, all terms that include I 1 vanish (become zero). This means that the matrix k lm e 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 U(t), The M, C and K are assembled by element mass, damping and stiffness matrices (M lm ) e , (C lm ) e and K lm e , 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].

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.

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 δ = 1 2 , α = 1 4 and δ = 1 2 , α = 1 6 it becomes average acceleration method and linear acceleration method, respectively.
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 "t + ∆t" to be written in terms of the unknown node displacements U(t + ∆t) l for lth harmonic of Fourier series: where K ll is effective dynamic stiffness matrix and F(t + ∆t) l 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 b i needs to be carried out only once. For linear systems, the effective dynamic stiffness matrix K ll is formed and triangularized only once.

•
Step 1: Initial Calculation (d) Form effective stiffness matrix K ll .
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). (e) Go to Step 2 (b) with t = t + ∆t.

Acceleration-Based Approach
The unknown node accelerations ..
U(t + ∆t) l 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.

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. 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. 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. 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 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 xand 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. 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. The computational results from four response points in the loading history are shown in Figure 4.  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.
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 2. Comparison between ABAQUS and SAFEM regarding the computational results at critical points when the loading time is 0.2 s.

SAFEM ABAQUS Difference
Vertical   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. 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.

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

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.