MatFEA: Efficient Finite Element Framework for Analyzing Pavement Structures With Nonlinear Unbound Materials

Abstract

This paper reports the details of a computationally efficient finite element analysis code (called MatFEA), designed specifically for flexible pavements, where stress-dependent characteristics of unbound layers are considered. The stress-dependent nonlinearity of the unbound materials has been modeled in MatFEA through an iterative approach. MatFEA was mainly developed to improve the Mechanistic-Empirical Asphalt Pavement Analysis (MEAPA) web application with capability of analysis/design of nonlinear pavement structures. A specific-purpose mesh generation algorithm was developed for the MatFEA in which the mesh density was dynamically adjustable based on the expected stress magnitude. To verify the results of the MatFEA, 12,000 nonlinear pavement structures with different structural and loading properties were generated and analyzed with MatFEA. The same pavement structures were modeled in a general-purpose FE-based ABAQUS™ program. The results of this study showed that the critical pavement structural responses from the MatFEA were almost identical with those modeled by ABAQUS™ program. The average runtime of MatFEA for analyzing each pavement structure was about 1.13 s. Finally, it was concluded that application of the dynamically adjustable mesh generation algorithm in the MatFEA helped achieve a reasonably high accuracy as well as efficient runtime, which brings the advantage of implementing the MatFEA as a pavement structural response model in the mechanistic-empirical pavement analysis/design procedures.

Keywords

pavement analysis nonlinear unbound material finite element analysis stress-dependency mechanistic-empirical

Flexible pavements are multi-layer structures mainly consisting of asphalt concrete (AC) in the top layers, unbound granular aggregates in intermediate layers and unbound fine-grain soils in lower layers or subgrade. The response of pavement structure to different traffic tire loads is a key factor in the mechanistic-empirical (ME) pavement design procedures. Two common pavement structural analysis methods are: layered elastic analysis (LEA) and numerical finite element (FE) analysis.

The current procedure of ME pavement design uses LEA-based software for calculation of pavement responses ( 1 , 2 ). Thus, although the AC layer at low strain levels is well known to act as a linear viscoelastic material, for pavement design purposes, it is still considered as a linear elastic layer. However, the viscoelastic behavior of AC is indirectly taken into account by computing an equivalent elastic modulus from the complex (dynamic) modulus (E^*) master curve of the AC layers. In this procedure, AC layers are divided into several sub-layers, where the elastic modulus of each asphaltic sub-layer is estimated from its E^* master curves based on duration of applied traffic load and temperature gradient, calculated by enhanced integrated climatic model (EICM).

On the other hand, almost all LEA methods assume a linear elastic behavior for unbound granular aggregates and fine-grained soils in non-asphaltic layers. For example, Michigan Department of Transportation (MDOT) practice is to assume an elastic modulus of 33,000 pounds per square inch (psi) for unbound granular base and 20,000 psi for sand subbase layers during pavement design using the AASHTOWare Pavement ME software ( 3 ). However, it is well known that the behavior of unbound granular materials and fine-grained soils are nonlinear (stress-dependent) under traffic tire loads ( 4 – 6 ). This nonlinearity significantly affects the calculation of critical pavement responses ( 7 , 8 ). The main reason for not including nonlinear behavior in pavement design and analysis software is that, without a finite element analysis (FEA) core engine, nonlinearity of base, subbase and subgrade layers cannot be accurately modeled. Since the existing FEA algorithms are computationally inefficient (compared with LEA algorithms), the pavement research community has historically accepted the simplification of unbound material nonlinear behavior into an elastic material. There is a need for computationally efficient FEA engines so that the true nonlinear behavior of unbound layers can accurately be considered in pavement design and analysis applications.

The objective of this study was to develop a reasonably accurate and computationally efficient FE-based pavement analysis framework, called MatFEA, in which the nonlinearity of unbound materials is included. The ultimate goal for this effort was to enhance the Mechanistic-Empirical Asphalt Pavement Analysis (MEAPA) web application (https://paveapps.com/meapaapp3/) so that MatFEA is used as the main pavement structural response model when nonlinearity of unbound materials needs to be considered ( 1 , 9 ).

Modeling Methods for Unbound Material Nonlinearity

Cyclic loading tests are typically used to characterize the behavior of unbound materials ( 10 ). For an unbound material, it is well known that plastic deformation occurs during initial loading cycles, which tends to stabilize at a certain level as the number of cycles grow. This phenomenon is called shakedown, which can take place in the construction and initial traffic passing stages. After shakedown occurs, the unbound material exhibits steady state and elastic strain response by which the resilient modulus (M_R) at the applied load level is defined:

M_{R} = \frac{σ_{d}}{ε_{r}}

(1)

where M_R is the resilient modulus of unbound material (psi), σ_d is the applied deviator stress (psi) in a triaxial test and ε_r is the recoverable strain after shakedown. The M_R value of unbound materials depends on the applied stress. Unbound granular materials mainly exhibit stress-hardening behavior ( 4 , 6 , 11 ), where the M_R values increase as the load level increases, while fine-grained soils often exhibit stress-softening behavior ( 5 , 11 , 12 ) where M_R decreases with increasing stress state. Hicks and Monismith ( 13 ) developed the so-called K-θ model for unbound granular materials by which the M_R is a nonlinear function of bulk stress, as shown in Equation 2.

M_{R} = K_{1} \cdot θ^{K_{2}}

(2)

where K₁ (psi) and K₂ are model coefficients and θ is the bulk stress (psi), which is calculated as:

θ = (σ_{r} + K_{0} \cdot γ \cdot z) + (σ_{z} + γ . z) + (σ_{θ} + K_{0} \cdot γ \cdot z)

(3)

where σ_r, σ_z, σ_θ are radial, vertical and tangential stresses (psi) at the point of interest; K₀ is the coefficient of earth pressure at rest; γ is the unit weight of the sample (pounds per cubic inch) and z is the depth from the top of sample (inches). It should be noted that tensile stress components (after adding the effect of earth pressure) are assumed to be zero as they do not cause any confinement. As the K-θ model is a single stress component model based on the bulk stress, it cannot incorporate the effect of shear stress on the predicted M_R value. This could cause limitation for pavement structures with relatively large shear stresses.

Witczak and Uzan ( 14 ) introduced a universal model by incorporation of the effect of deviatoric and shear stresses via octahedral shear stress (τ_oct), as shown in Equation 4. Since τ_oct covers the material stress states in all directions, this model has gained more attention for use in the three-dimensional (3D) analysis of nonlinear pavement structures.

M_{R} = K_{1} \cdot p_{a} {(\frac{θ}{p_{a}})}^{K_{2}} {(\frac{τ_{oct}}{p_{a}})}^{K_{3}}

(4)

where p_a is the atmospheric air pressure (≈14.706 psi), K₁, K₂ and K₃ are the model coefficients and τ_oct is the octahedral shear stress (psi), which is calculated as:

τ_{oct} = \sqrt{\frac{2}{3} J_{2}}

(5)

where J₂ is the second invariant of stress (psi²) and is calculated for an axisymmetric geometry using Equation 6.

J_{2} = \frac{1}{6} ({[(σ_{r} + K_{0} \cdot γ \cdot z) - (σ_{z} + γ \cdot z)]}^{2} + {[(σ_{r} + K_{0} \cdot γ \cdot z) - (σ_{θ} + K_{0} \cdot γ \cdot z)]}^{2} + {[(σ_{θ} + K_{0} \cdot γ \cdot z) - (σ_{z} + γ \cdot z)]}^{2} + 6 σ_{rz}^{2})

(6)

where σ_rz is the shear stress (psi). Similar to the calculation of bulk stress, the tensile stress components (after adding the effect of earth pressure) are assumed to be zero, since unbound materials do not take tension. Later on, the universal model was modified during the NCHRP 1-37A project, as shown in Equation 7 ( 2 ). The modification of adding “+1” to the shear softening term of universal model was added to avoid the effect of zero deviatoric shear stresses, as the K₃ exponent is negative.

M_{R} = K_{1} \cdot p_{a} {(\frac{θ}{p_{a}})}^{K_{2}} {(\frac{τ_{oct}}{p_{a}} + 1)}^{K_{3}}

(7)

where all the parameters are previously defined. The NCHRP 1-37A model for nonlinear unbound granular materials was planned to be used in MEPDG, but never implemented in the AASHTOWare Pavement ME software, because of the computational inefficiency of the FE engines available at the time ( 2 ).

The M_R of fine-grained soils is often believed to have stress-softening behavior as a function of deviator stress (σ_d). Several mathematical models have been proposed for representing the M_R of subgrade soils ( 15 ), among which the bilinear model ( 16 ), which is another single stress component model based on deviatoric stress, has been mostly utilized in the pavement community, as shown in Equation 8.

{\begin{matrix} M_{R} = K_{1} + K_{3} (K_{2} - σ_{d}) < M_{R, \max} when σ_{d} \leq K_{2} \\ M_{R} = K_{1} - K_{3} (σ_{d} - K_{2}) > M_{R, \min} when σ_{d} \geq K_{2} \end{matrix}

(8)

where M_R,max and M_R,min are the maximum and minimum M_R values (psi), and K₁, K₂, K₃ and K₄ are the model regression coefficients. Several studies have investigated the existence of a correlation between the M_R model coefficients and soil properties such as gradation, density, moisture content, porosity and so forth ( 4 , 11 , 12 , 17 –20).

Although the M_R and Young’s modulus (E) are different in nature, the use of M_R instead of E for linear elastic layers is a common practice during pavement structural analysis. The reason is that an unbound material experiences its shakedown stages during construction and initial traffic, which means the behavior of unbound materials under regular traffic loading can reasonably be assumed to be purely elastic.

Pavement Structure Analysis Programs

Burmister ( 21 ) developed a theory of layered elastic solution for analysis of two- and three-layered pavement structures in which pavement layers are assumed to be homogeneous, isotropic and linear elastic. This theory was later extended for an arbitrary number of pavement layers, where each pavement layer is represented by its elastic modulus (E), Poisson’s ratio (ν) and thickness. This theory can only model a single circular tire loading with a uniform contact pressure.

Several LEA-based pavement analysis softwares have been developed including: ELSYM5 ( 22 ), WESLEA ( 23 ), VESYS ( 24 ), KENLAYER ( 25 ) and others. The repeated analysis of pavement structure in the ME pavement design procedure has led to several improvements to increase the computational efficiency of LEA-based software: MnLayer ( 26 ), ERAPAVE ( 27 ) and MatLEACANGG ( 28 ) can be listed in this regard.

Although the layered elastic theory was developed based on linear elastic behavior, the nonlinearity of materials may also be approximated through an iterative process. Huang ( 25 ) incorporated the nonlinearity concept into the KENLAYER software by iteratively adjusting the elastic modulus of the nonlinear layer based on the stress states of a specific stress point within that layer. Huang suggested using a point at the mid-depth of each nonlinear layer, called the stress point. To increase the accuracy of this approach, the nonlinear layer was divided into several sub-layers of about 2 in. thickness. The same approach is used by ERAPAVE and MatLEACANGG. As pavement layers are assumed to be homogeneous in both horizontal and vertical directions, the calculated nonlinear pavement responses are expected to diverge from actual responses at locations radially further away from the stress point. Therefore, Huang ( 25 ) recommended locating the stress point under the tire load when the maximum stresses and displacements are required or using his proposed formulation for radial distances for calculation of average responses. Varma and Kutay ( 29 ) also suggested the use of 3.5 times the loading contact radius for the radial distance of stress point from the axis of symmetry. The authors of the current paper conducted extensive runs of LEA-based approximate nonlinear methods (e.g., Huang [ 25 ], Varma and Kutay [ 29 ]) and compared the results with ABAQUS™ nonlinear runs. It was observed that although the LEA-based approximate nonlinear methods are reasonable in some cases, they diverged from the results of nonlinear modeling in ABAQUS™ in many other cases.

Unlike the LEA approach, the FE approach is a mechanistic solution for analyzing the pavement structures by which the nonlinear and/or anisotropic material behavior can be considered. In the FE approach, the pavement structure is divided into a finite number of elements where any specific constitutive model behavior can be defined for each of them. Considering the pavement structure and circular load geometry, two-dimensional (2D) axisymmetric FE solutions have been a common practice for analysis of pavement structures. Such 2D axisymmetric FE solutions are simpler to model and more computationally efficient compared with 3D solutions. ILLI-Pave ( 30 ), MICH-Pave ( 31 ), GT-PAVE ( 32 ) and C-FLEX ( 33 ) are some of the well-known FE-based pavement analysis softwares which are based on 2D axisymmetric solution.

Two-dimensional axisymmetric FE solutions can only model a single tire loading. Although the superposition concept has been extensively used to consider the effect of multiple tire loadings with LEA solutions, it is not used with 2D axisymmetric FE models. This is while full 3D FE models have been extensively used to model the nonlinearity of pavement structure under multiple tire loading more realistically. Therefore, several specific-purpose 3D FE pavement analysis programs have been developed, including SAFEM ( 34 ), CAPA-3D ( 35 ) and others. On the other hand, several studies have used the general-purpose FE softwares, such as ABAQUS, ANSYS and ADINA, for 3D modeling of pavement structures with different material behaviors and various loading conditions ( 7 , 8 , 29 , 36 , 37 ). It has been reported that the computational runtime for analyzing nonlinear pavement structure significantly increases with the use of 3D FE solutions compared with the 2D axisymmetric ones ( 7 , 38 ), which makes the implementation of 3D FE solutions for the pavement design procedures impractical.

Development of MatFEA

MatFEA has been developed as a modularized MATLAB code to efficiently analyze pavement structures with nonlinear unbound layers through the FE approach. It can model the nonlinearity of unbound materials based on the models discussed above. Figure 1 shows the analysis logic of MatFEA in a flowchart. Analysis formulations, boundary condition, mesh generation algorithm and nonlinearity iterations are explained in this section.

Figure 1.

Flowchart of pavement analysis in MatFEA.

Strain-Displacement Relationship

For an axisymmetric FE solution, the generic nodal displacement vector of each node within the r-z plane includes radial (u) and vertical (w) components. The differential strain-displacement relations for the axisymmetric solution are shown Equation 9.

{\begin{matrix} \begin{matrix} ε_{r} \\ ε_{z} \end{matrix} \\ \begin{matrix} ε_{θ} \\ ε_{rz} \end{matrix} \end{matrix}} = [\begin{matrix} \begin{matrix} \frac{\partial}{\partial r} & 0 \\ 0 & \frac{\partial}{\partial z} \end{matrix} \\ \begin{matrix} \frac{1}{r} & 0 \\ \frac{\partial}{\partial z} & \frac{\partial}{\partial r} \end{matrix} \end{matrix}] \times {\begin{matrix} u \\ w \end{matrix}}

(9)

where ε_r is the radial strain, ε_z is the vertical strain, ε_θ is the tangential strain, ε_rz is the shear strain, r is the radial coordinate and z is the vertical coordinate. The stress vector is calculated using the strain vector using Equation 10.

{\begin{matrix} \begin{matrix} σ_{r} \\ σ_{z} \end{matrix} \\ \begin{matrix} σ_{θ} \\ σ_{rz} \end{matrix} \end{matrix}} = D_{4 \times 4} \times {\begin{matrix} \begin{matrix} ε_{r} \\ ε_{z} \end{matrix} \\ \begin{matrix} ε_{θ} \\ ε_{rz} \end{matrix} \end{matrix}} = D_{4 \times 4} \times ε_{4 \times 1}

(10)

where σ_r is the radial stress (psi), σ_z is the vertical stress (psi), σ_θ is the tangential stress (psi), σ_rz is the shear stress (psi), ε is the strain vector and D is the rigidity matrix which includes the constitutive relationship between the strain and stresses. It is noted that the symbol “×” in Equation 10 represents the matrix multiplication, not the cross product. The rigidity matrix can take different forms depending on the material behavior; for example, according to the generalized Hook’s law for an isotropic axisymmetric element, the rigidity matrix is shown in Equation 11.

D_{4 \times 4} = \frac{E}{(1 + ν) (1 + 2 ν)} [\begin{matrix} 1 - ν & ν & ν & 0 \\ ν & 1 - ν & ν & 0 \\ ν & ν & 1 - ν & 0 \\ 0 & 0 & 0 & \frac{1 - 2 ν}{2} \end{matrix}]

(11)

where E is the modulus of elasticity (psi) and ν is the Poisson’s ratio.

MatFEA uses four-node bilinear elements in which each node has two degrees of freedom for radial and vertical displacements. Therefore, the transfer matrix ( B ) for each element in the local coordinate system (ζ-η) is shown in Equation 12.

B_{4 \times 8} = [\begin{matrix} \frac{\partial N_{1}}{\partial ξ} & 0 & \frac{\partial N_{2}}{\partial ξ} & 0 & \frac{\partial N_{3}}{\partial ξ} & 0 & \frac{\partial N_{4}}{\partial ξ} & 0 \\ 0 & \frac{\partial N_{1}}{\partial η} & 0 & \frac{\partial N_{2}}{\partial η} & 0 & \frac{\partial N_{3}}{\partial η} & 0 & \frac{\partial N_{4}}{\partial η} \\ \frac{N_{1}}{R} & 0 & \frac{N_{2}}{R} & 0 & \frac{N_{3}}{R} & 0 & \frac{N_{4}}{R} & 0 \\ \frac{\partial N_{1}}{\partial η} & \frac{\partial N_{1}}{\partial ξ} & \frac{\partial N_{2}}{\partial η} & \frac{\partial N_{2}}{\partial ξ} & \frac{\partial N_{3}}{\partial η} & \frac{\partial N_{3}}{\partial ξ} & \frac{\partial N_{4}}{\partial η} & \frac{\partial N_{4}}{\partial ξ} \end{matrix}]

(12)

where N_i is the shape function of each node in the local coordinate system, ζ and η are the radial and vertical components of the local coordinate system and R is radial distance from the axis of symmetry, as shown in Equation 13.

R = \sum_{i = 1}^{4} N_{i} . r_{i}

(13)

where r_i is the global radial coordinate of the ith node in each element (inches). The elemental stiffness matrix ( k ) was calculated using Equation 14.

k_{4 \times 4} = 2 π \int \int_{A} B^{T} \times D \times B \times r . dr . dz

(14)

An isoparametric approach was then used to convert the global coordinate (r-z) to the local coordinate (ζ-η) system. Thus, the elemental stiffness matrix is represented as shown in Equation 15.

k_{4 \times 4} = 2 π \int_{- 1}^{1} \int_{- 1}^{1} B^{T} \times D \times B \times r \times | J | . d ξ . d η

(15)

where J is the Jacobian matrix and is defined as shown in Equation 16.

J_{2 \times 2} = [\begin{matrix} \frac{\partial r}{\partial ξ} & \frac{\partial z}{\partial ξ} \\ \frac{\partial r}{\partial η} & \frac{\partial z}{\partial η} \end{matrix}] = \sum_{i = 1}^{4} [\begin{matrix} r_{i} \frac{\partial N_{i}}{\partial ξ} & z_{i} \frac{\partial N_{i}}{\partial ξ} \\ r_{i} \frac{\partial N_{i}}{\partial η} & z_{i} \frac{\partial N_{i}}{\partial η} \end{matrix}]

(16)

The Gaussian integration method with three Gaussian points was used to numerically calculate the elemental stiffness matrix, which was then combined together to form the global stiffness matrix ( K ).

As the tire loading is the only external load to the pavement structure, the magnitudes of nodal forces at the surface nodes within the loading area were calculated by multiplying the contact tire pressure to the corresponding loading strip area, as shown in Equation 17.

F_{i} = - CP . π . (r_{out}^{2} - r_{in}^{2})

(17)

where i is counter of nodes on the pavement surface and within the tire loading area, F_i is the vertical nodal force applied to the i^th node (pounds), CP is the contact tire pressure (psi), and r_out and r_in are the outer and inner radius of loading strip (inches), such that r_out = (r_i+1+r_i) / 2 and r_in = (r_i−1+r_i) / 2. Finally, the nodal displacements were calculated by applying boundary conditions and solving Equation 18.

F = K \times d

(18)

where F is the nodal external forces matrix and d is the nodal displacement matrix.

Boundary Condition and Domain Size

Layered elastic theory uses the semi-infinite half-space assumptions in which the pavement structure extends to infinity in both vertical and horizontal directions. This infinite assumption is not implementable in FE solutions. Therefore, the boundary conditions and domain sizes of the FE model should be selected carefully to maintain the accuracy of the solution and avoid boundary edge errors.

According to the 2D axisymmetric geometry of the FE model, the structural roller support should be used on the inner vertical edge of the model. This will restrict the horizontal displacements but vertical displacement is allowed ( 39 ). Structural roller support was also used at the outer vertical edge of the model. The horizontal bottom edge of the model was supported by the structural fixed support, which restricts both vertical and horizontal displacements. The same boundary conditions have been used in the literature ( 7 , 8 , 29 , 30 , 33 , 36 ). Figure 2 shows a conceptual scheme of the boundary conditions in MatFEA.

Figure 2.

Boundary condition in MatFEA.

An appropriate domain size of the pavement structure in the FE model can minimize the effect of boundary conditions on pavement responses. Dehlen ( 40 ) recommended vertical and radial domain sizes of 50 and 12 times the loading radius. The same vertical domain size is also used in ILLI-PAVE ( 30 ). The vertical and radial domain sizes of 140 and 20 times the loading radius have been recommended based on comparison with the linear analytical solutions ( 7 ). Saad et al. ( 36 ) performed a domain size analysis by monitoring the maximum vertical strain at the top of the subgrade and maximum tensile strain at the bottom of the AC layer, while keeping the same mesh design. The same procedure was performed in the current study and the domain size of 150 and 40 times the loading radius was selected for the depth and radius of the pavement structure model in MatFEA.

Mesh Generation Algorithm

Mesh design is one of the most important factors affecting the efficiency of a FE solution. Application of coarse mesh reduces the number of elements which brings reduced computational runtime but less accuracy in pavement responses. Refining the mesh could also improve the solution accuracy while increasing the required computational effort.

The application of a 2D axisymmetric model for analyzing the pavement structure provides a simple rectangular layered geometry. Therefore, a recursive specific-purpose unstructured mesh generation algorithm was developed to generate four-node quadrilateral elements. In this algorithm, the mesh density is dynamically adjusted based on the expected stress state, that is, it generates fine mesh at the areas where high stress is expected close to the load centerline and gradually increases the mesh size as it approaches the bottom and outer boundaries. It is noted that this algorithm does not change the mesh from one iteration to another (see Figure 1) and the mesh density is pre-defined based on the distance from the pavement surface and load centerline. The algorithm also prevented the formation of irregular or rotated (diamond like) quadrilateral elements, as well as elements with aspect ratio of greater than 10. It has been reported that the calculated stresses from the quadrilateral elements with the aspect ratios of less than five are reasonably accurate ( 31 ). The main purpose for developing this mesh generation algorithm was to maintain the solution accuracy and computational efficiency simultaneously. An example of the mesh generated for a typical pavement structure is shown in Figure 3.

Figure 3.

Mesh design in MatFEA.

Nonlinear Iterative Procedure

MatFEA uses an iterative approach for capturing the nonlinear behavior of materials. For this purpose, elemental stress responses were calculated by assuming an initial M_R value for the elements in the nonlinear layers with the same materials. At each iteration, the M_R of each element in the nonlinear layer was updated after running FE analysis based on the calculated stress responses and the corresponding nonlinear M_R model. The relative difference in M_R value (ΔM_R) was then calculated for each element using Equation 19.

Δ M_{R} = \frac{| M_{R, new} - M_{R, old} |}{M_{R, new}} \times 100

(19)

where M_R,new is the updated elemental resilient modulus (psi) and M_R,old is the resilient modulus of the previous iteration (psi). The elemental ΔM_R values of all elements in a given nonlinear layer were statistically evaluated against the convergence criterion. This iterative procedure (see Figure 1) continued until the convergence criterion was satisfied.

The convergence criterion was defined as a mean ΔM_R of all elements in the nonlinear layers at a specific reliability level, as shown in Equation 20.

ψ = μ - Z_{R} . Std

(20)

where ψ is the mean value of ΔM_R values at the specific reliability level (%), Z_R is the normal deviate for the reliability level, μ and Std are the mean and standard deviation of the elemental ΔM_R values for all elements in a nonlinearity layer (%), respectively. MatFEA uses the 95% reliability level (Z_R = −1.645) and 2% convergence criterion under a maximum of 10 iterations.

Verification of MatFEA

ABAQUS™, a general-purpose FE analysis tool, was used to verify the accuracy of the MatFEA implementation. For this purpose, a database of 12,000 nonlinear pavement structures was randomly generated, which includes 3,000 pavement structures with nonlinear bases (K-θ model in Equation 2), 3,000 pavement structures with nonlinear bases (NCHRP 1-37A model in Equation 7), 3,000 pavement structures with nonlinear subgrades (bilinear model in Equation 8), and 3,000 pavement structures with both nonlinear bases (NCHRP 1-37A model in Equation 7) and nonlinear subgrades (bilinear model in Equation 8). Each pavement structure/loading property in this database was generated by selecting a random number within a pre-determined interval for the property of interest, as shown in Table 1. These pavement sections were then modeled through axisymmetric FE method using ABAQUS™ and MatFEA programs.

Table 1.

Variation Interval of Three-Layer Pavement Structure and Loading Properties

Property	Variation interval
Property	Minimum value	Maximum value
Load magnitude (lb)	1,000	16,000
Load contact pressure (psi)	120.0	120.0
AC thickness (in.)	4.0	10.0
AC modulus (psi)	100,000	4,000,000
AC Poisson’s ratio	0.25	0.35
AC unit weight (pcf)	140	150
Base thickness (in.)	6.0	16.0
Base M_R (psi)^a	10,000	30,000
Base Poisson’s ratio	0.30	0.45
Base unit weight (pcf)	115	130
Base K₁ in K-θ model (psi)^b	2,000	10,000
Base K₂ in K-θ model^b	0.4	0.7
Base K₁ in NCHRP 1-37A model^c	639.1	1,015.0
Base K₂ in NCHRP 1-37A model^c	0.2570	0.6335
Base K₃ in NCHRP 1-37A model^c	−1.6649	−0.0948
Subgrade M_R (psi)^d	5,000	15,000
Subgrade Poisson’s ratio	0.35	0.45
Subgrade unit weight (pcf)	100	130
Subgrade K₁ in bilinear model (psi)^e	1,000	12,500
Subgrade K₂ in bilinear model (psi)^e	5.8	6.6
Subgrade K₃ in bilinear model^e	1,020	1,200
Subgrade K₄ in bilinear model^e	160	200
Coefficient of earth’s pressure at rest	0.60	0.60

Note: AC = asphalt concrete; pcf = pounds per cubic foot; psi = pounds per square inch.

Base resilient modulus interval for the cases with linear elastic base.

Based on Huang ( 25 ).

Based on Yau and Von Quintus ( 11 ).

Subgrade resilient modulus interval for the cases with linear elastic subgrade.

Based on Huang ( 25 ).

To be consistent with the algorithm of the MatFEA, a Python script was developed to model each pavement section in the ABAQUS™ program and perform the nonlinearity iterations (instead of utilizing the user material [UMAT] Fortran subroutine). The ABAQUS™ models use a free mesh based on four-node first-order elements with reduced integration (CAX4R). The developed Python script consists of the following steps:

Model the pavement structure and applied pressure loading in an axisymmetric scheme;

Mesh the model such that the mesh is denser at the high stress areas;

Run the analysis for the initial model;

Perform iteration by adjusting the modulus of the elements in the nonlinear layers and re-running the model until the convergence criterion is met.

Four different critical pavement responses were compared against ABAQUS™ for verification of the MatFEA: maximum surface deflection, radial strain at the bottom of the AC layer, vertical strain at the middle of the base layer and vertical strain at the top of the subgrade. Commonly used statistical metrics were utilized for comparison purposes, including coefficient of determination (R²) (Equation 21), mean absolute error (MAE) (Equation 22), root mean square error (RMSE) (Equation 23) and bias (Equation 24).

R^{2} = 1 - \frac{\sum {(Resp_MatFE A_{i} - Resp_ABAQU S_{i})}^{2}}{\sum {(Resp_MatFE A_{i} - Resp_MatFE A_{avg})}^{2}}

(21)

MAE = \frac{\sum | Resp_MatFE A_{i} - Resp_ABAQU S_{i} |}{n}

(22)

RMSE = \sqrt{\frac{\sum {(Resp_MatFE A_{i} - Resp_ABAQU S_{i})}^{2}}{n}}

(23)

Bias = \frac{\sum (Resp_MatFE A_{i} - Resp_ABAQU S_{i})}{n}

(24)

where Resp_MatFEA_i and Resp_ABAQUS_i are the structural responses for each pavement section using the MatFEA and ABAQUS™ programs and n is the number of pavement sections used for verification.

As the ABAQUS™ models are used for verification of the results of MatFEA models, a relatively denser mesh was used in ABAQUS™ models to assure higher accuracy of the pavement structural responses. In this regard, Figure 4 shows the distribution of the number of elements used for modeling each of the pavement structures in the database using the MatFEA and ABAQUS™ programs. As this figure shows, the mesh of ABAQUS™ models was about six times finer than in the MatFEA models.

Figure 4.

Distribution of number of elements used in modeling of pavement structures.

The equality graphs and calculated statistical metrics for comparing the results of MatFEA versus ABAQUS™ models for pavement structures in the database are shown in Figure 5. As this figure shows, almost all of the data points for all critical structural responses were clustered around the line of equality. This observation was supported with the statistical metrics, such as high R², and considerably low MAE, RMSE and bias.

Figure 5.

Comparison of pavement critical responses from MatFEA and ABAQUS™ models at: (a) maximum surface deflection, (b) radial strain at bottom of asphalt concrete layer, (c) vertical strain at middle of base depth, and (d) vertical strain at top of subgrade.

Although the pavement structures were modeled in ABAQUS™ in accordance with the main algorithm of the MatFEA, Figure 5 shows that the structural responses from MatFEA did not exactly match those from ABAQUS™. This can be attributed to the significantly finer mesh of the model in the ABAQUS™ program with which higher accuracy can be achieved. Moreover, a sensitivity analysis revealed that the best match between the ABAQUS™ and MatFEA models was achieved for the pavement structures with nonlinear bases (either the K-θ or NCHRP 1-37A models), while some divergences were observed in vertical strains at the top of the subgrade for pavement structures with nonlinear subgrade (bilinear model).

It is noted that the average runtime of the MatFEA for 12,000 pavement structures in the database was about 1.1 s, while the average runtime of each ABAQUS™ model was about 227.8 s. The significantly higher runtimes of the ABAQUS™ models were the result of their finer mesh, as well as application of the iterative approach using the Python scripts. To better compare the runtime of MatFEA against ABAQUS™, one of the pavement structures with nonlinear base layer were modeled using MatFEA, ABAQUS™ with iterative approach (Python script), and ABAQUS™ with incremental approach using the user material subroutine (UMAT). The MatFEA model had 683 elements. Therefore, the mesh setting of ABAQUS™ models were modified accordingly so that they include 640 elements. The results showed that the MatFEA model ran in 0.80 s, while the ABAQUS™ models with iterative approach (Python script) and incremental approach (using UMAT) ran in 65.34 s and 11.07 s, respectively. These results indicated that the ABAQUS™ model was about six times faster by using the incremental approach instead of iterative one, as the incremental approach uses UMAT files in Fortran language (compared with the Python language in the iterative approach) and it is more efficiently imbedded into the ABAQUS™ software. Also, it is noted that the MatFEA model was still significantly faster than the ABAQUS™ models. This does not necessarily show faster computation in MatFEA, but rather the long runtime of ABAQUS™ can be attributed to the overhead time and several checks before running the model. Finally, the average runtime of 1.1 s makes MatFEA a practical pavement structural response model to implement into ME pavement design procedures.

Summary and Conclusion

In this study, a 2D axisymmetric FE-based MATLAB code, called MatFEA, was developed to analyze multi-layered pavement structures which is capable of modeling the stress-dependent nonlinear behavior of unbound materials. The main goal of developing the MatFEA was to be reasonably accurate and computationally efficient such that it can be implemented into the ME pavement design/analysis procedures, for example, the MEAPA web application.

The radius and depth of the pavement structure model in MatFEA were selected to be 40 and 150 times the loading tire radius, respectively. A customized mesh generation algorithm was developed for MatFEA in which the mesh density is dynamically adjustable. This algorithm generates four-node quadrilateral elements which are finer at the expected highly stressed areas close to the load centerline and pavement surface, while the element size increases gradually as it approaches the bottom and outer boundaries.

To verify the results of the MatFEA, 12,000 pavement structures with different structural and loading properties were generated and analyzed with MatFEA, as well as the general-purpose FEA ABAQUS™ program. The pavement structure models in MatFEA consisted of about six times fewer elements on average compared with those modeled in ABAQUS™, which resulted in a significant computational advantage for MatFEA, where the average runtime of each pavement structure in MatFEA was 1.13 s. This is while the results showed that the critical pavement structural responses were almost identical between both MatFEA and ABAQUS™ models. In conclusion, the results showed that implementation of the dynamically adjustable mesh generation algorithm in MatFEA enables it to be used in ME-based pavement design/analysis procedures through its reasonable accuracy and significantly fast runtime.

Future Work

In this study, a framework for nonlinear (stress-dependent) pavement structural analysis, called MatFEA, was developed. The results of MatFEA were numerically verified against ABAQUS™ software, as a well-accepted general-purpose FEA software. In the next step, the authors will continue working on “validation” of the results of MatFEA against the data from the full-scale pavement structure simulator available at Michigan State University. Moreover, the LEA engine of the MEAPA web application will be replaced with MatFEA to study the effect of nonlinearity (stress-dependency) of unbound materials on the predicted pavement performance.

Footnotes

Acknowledgements

The authors greatly acknowledge the Institute for Cyber-Enabled Research (ICER) at Michigan State University for the computational resources provided.

Author Contributions

The authors confirm contribution to the paper as follows: study conception and design: S. F. Abdollahi, M. E. Kutay; data collection: S. F. Abdollahi; analysis and interpretation of results: S. F. Abdollahi, M. Ghazavi., M. E. Kutay; draft manuscript preparation: S. F. Abdollahi, M. E. Kutay. All authors reviewed the results and approved the final version of the manuscript.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Seyed Farhad Abdollahi

M. Emin Kutay

Mahdi Ghazavi

References

Kutay

M. E.

Lanotte

Formulations of the Pavement Performance Prediction Models in the Mechanistic-Empirical Asphalt Pavement Analysis (MEAPA) Web Application. Michigan State University, 2020.

NCHRP 1-37A. Guide for Mechanistic–Empirical Design of New and Rehabilitated Pavement Structures. Publication NCHRP Project 1-37A. Transportation Research Board, National Research Council, Washington, D.C., 2004.

MDOT. Michigan DOT User Guide for Mechanistic-Empirical Pavement Design. Michigan Department of Transportation, Lansing, MI, 2021.

Rada

Witczak

M. W.

Comprehensive Evaluation of Laboratory Resilient Moduli Results for Granular Material. Transportation Research Record: Journal of Transportation Research Board, 1981. 810: 23–33.

Thompson

M. R.

Elliot

R. P.

ILLI-PAVE-Based Response Algorithms for Design of Conventional Flexible Pavements. Transportation Research Record: Journal of Transportation Research Board, 1985. 1043: 50–57.

Uzan

Characterization of Granular Materials. Transportation Research Record: Journal of Transportation Research Board, 1985. 1022: 52–59.

Kim

Tutumluer

Kwon

Nonlinear Pavement Foundation Modeling for Three-Dimensional Finite-Element Analysis of Flexible Pavements. International Journal of Geomechanics, Vol. 9, No. 5, 2009, pp. 195–208. https://doi.org/10.1061/(ASCE)1532-3641(2009)9:5(195).

Kuo

C. M.

Huang

C. W.

Three-Dimensional Pavement Analysis With Nonlinear Subgrade Materials. Journal of Materials in Civil Engineering, Vol. 18, No. 4, 2006, pp. 537–544. https://doi.org/10.1061/(ASCE)0899-1561(2006)18:4(537).

Ghazavi

Abdollahi

S. F.

Kutay

M. E.

Implementation of NCHRP 9-44A Fatigue Endurance Limit Prediction Model in Mechanistic-Empirical Asphalt Pavement Analysis Web Application. Transportation Research Record: Journal of the Transportation Research Board, 2022. 2676: 696–706.

10.

AASHTO T 307. Standard Method of Test for Determining the Resilient Modulus of Soils and Aggregate Materials. American Association of State Highway and Transportation Officials (AASHTO), Washington, D.C., 2013.

11.

Yau

Von Quintus

H. L.

Study of LTPP Laboratory Resilient Modulus Test Data and Response Characteristics. Report No. FHWA-RD-02-051. Federal Highway Administration, McLean, VA, 2002.

12.

Santha

B. L.

Resilient Modulus of Subgrade Soils: Comparison of Two Constitutive Equations. Transportation Research Record: Journal of Transportation Research Board, 1994. 1462: 79–90.

13.

Hicks

R. G.

Monismith

C. L.

Factors Influencing the Resilient Response of Granular Materials. Highway Research Record, Vol. 345, 1971, pp. 15–31.

14.

Witczak

M. W.

Uzan

The Universal Airport Pavement Design System Report I of IV: Granular Material Characterization. University of Maryland, College Park, 1988.

15.

Selig

E. T.

Resilient Modulus for Fine-Grained Subgrade Soils. Journal of Geotechnical Engineering, Vol. 120, No. 6, 1994, pp. 939–957. https://doi.org/10.1061/(ASCE)0733-9410(1994)120:6(939).

16.

Thompson

M. R.

Robnett

Q. L.

Resilient Properties of Subgrade Soils. Transportation Engineering Journal of ASCE, Vol. 105, No. 1, 1979, https://doi.org/10.1061/TPEJAN.0000772.

17.

Dai

Zollars

Resilient Modulus of Minnesota Road Research Project Subgrade Soil. Transportation Research Record: Journal of Transportation Research Board, 2002. 1786: 20–28.

18.

Archilla

A. R.

Ooi

P. S.

Sandefur

K. G.

Estimation of a Resilient Modulus Model for Cohesive Soils Using Joint Estimation and Mixed Effects. Journal of Geotechnical and Geoenvironmental Engineering, Vol. 133, No. 8, 2007, pp. 984–994. https://doi.org/10.1061/(ASCE)1090-0241(2007)133:8(984).

19.

Nazzal

M. D.

Mohammad

L. N.

Estimation of Resilient Modulus of Subgrade Soils for Design of Pavement Structures. Journal of Materials in Civil Engineering, Vol. 22, No.7, 2010, pp. 726–734. https://doi.org/10.1061/(ASCE)MT.1943-5533.0000073.

20.

Saha

Luo

Lytton

R. L.

Use of an Artificial Neural Network Approach for the Prediction of Resilient Modulus for Unbound Granular Material. Transportation Research Record: Journal of Transportation Research Board, 2018. 2672: 23–33.

21.

Burmister

D. M.

The General Theory of Stresses and Displacements in Layered Soil Systems. III. Journal of Applied Physics, Vol. 16, No. 5, 1945, pp. 296–302. https://doi.org/10.1063/1.1707590.

22.

Kopperman

Tiller

Tseng

ELSYM 5: Interactive Microcomputer Version: User’s Manual: IBM-PC and Compatible Microcomputer. Report No. FHWA-TS-87-206. Federal Highway Administration, McLean, VA, 1986.

23.

Van Cauwelaert

Alexander

White

Barker

Multilayer Elastic Program for Backcalculating Layer Moduli in Pavement Evaluation. ASTM STP 1026. ASTM International, Philadelphia, PA, 1989, pp. 171–188.

24.

Jordahl

P. R.

Rauhut

J. B.

Flexible Pavement Model VESYS IV-B. Report No. FHWA-RD-84-021. Federal Highway Administration, McLean, VA, 1984.

25.

Huang

Y. H.

Pavement Analysis and Design. Pearson, Upper Saddle River, NJ, 2004.

26.

Khazanovich

Wang

Q. C.

MnLayer: A High-Performance Layered Elastic Analysis Program. Transportation Research Record: Journal of Transportation Research Board, 2007. 2037: 63–75.

27.

Erlingsson

Ahmed

Fast Layered Elastic Response Program for the Analysis of Flexible Pavement Structures. Road Materials and Pavement Design, Vol. 14, No. 1, 2013, pp. 196–210. https://doi.org/10.1080/14680629.2012.757558.

28.

Kutay

M. E.

Hasnat

Levenberg

Layered Nonlinear Cross Anisotropic Model for Pavements With Geogrids. Proc., 2nd International Conference on Advances in Materials and Pavement Performance Prediction (AM3P), San Antonio, TX, 2020.

29.

Varma

Kutay

M. E.

Viscoelastic Nonlinear Multilayered Model for Asphalt Pavements. Journal of Engineering Mechanics, Vol. 142, No. 7, 2016, p. 04016044. https://doi.org/10.1061/(ASCE)EM.1943-7889.0001095.

30.

Raad

Figueroa

J. L.

Load Response of Transportation Support Systems. Transportation Engineering Journal of ASCE, Vol. 106, No. 1, 1980, pp. 111–128.

31.

Harichandran

R. S.

Yeh

M.-S.

Baladi

G. Y.

MICH-PAVE: A Nonlinear Finite Element Program for Analysis of Flexible Pavements. Transportation Research Record: Journal of Transportation Research Board, 1990. 1286: 123–131.

32.

Tutumluer

Predicting Behavior of Flexible Pavements with Granular Bases. PhD thesis. Georgia Institute of Technology, Atlanta, GA, 1995.

33.

Lou

Huang

Qamhia

I. I.

Tutumluer

Tingle

J. S.

C-FLEX Advanced Finite Element Analysis Program for Flexible Pavement Analysis and Design. Transportation Research Record: Journal of Transportation Research Board, 2021. 2675: 1238–1249.

34.

Liu

Wang

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, 2017, p. 20150456. https://doi.org/10.1520/JTE20150456.

35.

Scarpas

CAPA-3D: A Mechanistic Based Computational Platform for Pavement Engineering. PhD thesis. Delft University of Technology, The Netherlands, 2005.

36.

Saad

Mitri

Poorooshasb

Three-Dimensional Dynamic Analysis of Flexible Conventional Pavement Foundation. Journal of Transportation Engineering, Vol. 131, No. 6, 2005, pp. 460–469. https://doi.org/10.1061/(ASCE)0733-947X(2005)131:6(460).

37.

Wang

Al-Qadi

I. L.

Importance of Nonlinear Anisotropic Modeling of Granular Base for Predicting Maximum Viscoelastic Pavement Responses Under Moving Vehicular Loading. Journal of Engineering Mechanics, Vol. 139, No. 1, 2013, pp. 29–38. https://doi.org/10.1061/(ASCE)EM.1943-7889.0000465.

38.

Schwartz

C. W.

Effect of Stress-Dependent Base Layer on the Superposition of Flexible Pavement Solutions. International Journal of Geomechanics, Vol. 2, No. 3, 2002, pp. 331–352. https://doi.org/10.1061/(ASCE)1532-3641(2002)2:3(331).

39.

Logan

D. L.

A First Course in the Finite Element Method. CL Engineering, Stamford, CT, 2011.

40.

Dehlen

G. L.

The Effect of Non-Linear Material Response on the Behaviour of Pavements Subjected to Traffic Loads. PhD thesis. University of California, Berkeley, 1969.