Prediction of Rutting in Flexible Pavements using Finite Element Method

In this research study three dimensional (3D) finite element analysis are performed on a flexible pavement section for different material properties, temperature and loading conditions. The main objective of this study is to predict the rut depth under different conditions of temperature, loadings and material properties. Three dimensional finite element model of flexible pavement is developed using ABAQUS to predict rut depth. The pavement system is assumed to be an elastic multi-layers system with each layer being isotropic, homogeneous with specified Resilient Modulus (Mr) and poisson ratio (μ). With the exception of the bottom subgrade layer, each layer is extending to an unlimited horizontal extent and has a finite thickness. The pavement system analyze in this study for a cyclic load of 10000 cycles taken as 0.01sec per cycle. Standard Axle Load (ESAL) of 18 kips (80 kN) loading on an axle with a dual set of tires, the wheel spacing is 13.78 in (350 mm) with a tire contact pressure of 100 psi (0.69 MPa) is used. After performing a series of analysis the results showed that rut depth increases with increase in temperature and loading and decreases by using base stabilizer.


Introduction
Pakistan has national highways with a length over 9555 km. Premature rutting in flexible pavements of highways is one of the major distress being faced in Pakistan. Rutting occurs when persistent deformation occurs in all pavement layers as a result of frequent traffic loadings. Among the various pavement layers' contributions to rut depth, the cumulative permanent deformation in the asphalt pavement's surface course is known to be responsible for a significant amount of the final rut depth measured on the pavement surface. As a result, rutting only occurs in flexible pavements as evidenced by permanent deformation or rut depth along the wheel paths.
Rutting is the most common problem with asphalt surfaces, especially when the weather is hot. Temperature variation is the most significant factor that affects flexible pavement deformation, as it affects the stiffness of the One of the most essential features of pavement design is the rutting potential of asphalt concrete mixes [12]. It's the result of load repetition causing permanent pavement deformation along the wheel path. Rutting can be found deep within the subgrade or only on the surface. Cumulative persistent deformation in the surface is recognized to be responsible for the majority of the rut depth of a pavement surface, among the several layers that contribute to rutting [13]. Rutting has a substantial impact on the performance of asphalt pavements by increasing road roughness and retaining water, resulting in traffic accidents and hydroplaning due to the loss of tire-pavement friction [14]. Excessive rutting also causes structural failures. As a result, rutting not only shortens the life of pavement, but also puts highway safety at risk [15]. One of the goals of road engineers is to build pavement structures with improved characteristics that can provide greater performance and a longer service life [16].
For many years, evaluating asphalt concrete mixes for their rutting potential has been a major research topic. Pavement design requires the capacity to estimate the quantity and growth of permanent deformation or rutting in flexible pavements. Different pavement problems that could not be handled using the simpler multi-layer elastic theory can be successfully simulated utilizing Finite Element (FE) techniques. The structural characteristics of the pavement layers (thickness and material quality), traffic loads, and environmental variables all influence the width and depth of the rut [17].
The Finite Element Method (FEM) is a strong methodology utilized in a variety of engineering disciplines. The goal of this strategy is to come up with a numerical solution to difficult engineering problems. It broke the complicated structure into microscopic little parts, analyzed the small parts, and obtained a general approximation solution to the complex structure by summing the solutions of the small parts. ABAQUS, a three-dimensional, dynamic finite element program [18], has the capability to simulate actual pavement loading conditions. Material properties such as elasticity, viscoelasticity, and plasticity can be represented using a variety of material models. Previous research has detailed certain prediction models that take into account the influence of stress on the resilient modulus [19,20]. Establishing an appropriate resilient modulus prediction model and transplanting it into the finite element (FE) method is an important way for analyzing the mechanical response of flexible pavement structures [21,22]. The first reported use of the FE method in flexible pavement engineering was published in Duncan et al. [23]. Since then, a considerable deal of study has been done using the FE approach to do linear and nonlinear analysis [24,25].
Using 3D finite-element modeling created a better match between predicted and measured results than did 2D modeling. Therefore the present 3D model can predict the rut depth of dense-graded rubberized and conventional mixtures in different loading conditions. To minimize computational cost, a derivation of reduced-order models (model reduction) for engineering applications might be performed using the current simulation approach [26]. Developed models can accurately predict rutting behavior of bituminous mixtures under various loads and at 40 °C temperature [27]. Under combined traffic and temperature conditions, the model has a lesser resilience to rutting failure. By raising the local temperature on the surface of flexible pavement to around 45°C, the maximum number of repetitions required to generate rutting is reduced by roughly three times less than the model of traffic loading alone [28]. Temperature, loading and bituminous material properties have an effect on the depth of rutting of asphalt mixtures [29]. Higher traffic loading result higher stress, strain, and displacement at top of the subgrade with or with no geo-jute inclusion [30].

FEM of Flexible Pavement
The finite element method is a numerical methodology for achieving an approximation solution in the board field of a continuum solid. It is one of the most powerful methods for simulating behavioral reactions to diverse structural engineering problems. Flexible pavements are well-known for being complicated engineering structures with several layers, materials, and boundary conditions. As a result, this study provides a good approach for analyzing flexible pavements using the finite element method, as well as how to obtain the data needed to complete finite element analysis. ABAQUS, a three-dimensional solid structural study of pavement using element C3D8R, is used to analyse and model the pavement system's behavioral response [31]. The C3D8R element is a general-purpose linear brick element having three degrees of freedom at each node translations in the nodal x, y, and z directions, as shown in Figure 1, and reduced integration (1 integration point) situated in the middle of the element. Small elements are required to capture a stress concentration at the structural boundary.

Section Assignment
Flexible pavement is a multi-material structure with a complex structure. Some of these materials are elastomeric, while others are viscoelastic and plastic. As a result, a flawless mathematical model is impossible to achieve. Figure 2 shows the analytical three dimensional model mesh for the pavement system including different layer properties: Asphalt layer (90 mm thick), granular base (110 mm thick), sub base (110 mm thick) and subgrade (150 mm thick), Longitudinal and transverse directions (30,000 mm × 8230 mm) to limit accuracy in the edge effect. In finite element analysis, each pavement layer is believed to be visually continuous. Domains are divided into sub-domains that are visual, minimal and continuous. For each sub-domain divided into finite elements through meshing. These finite elements come together with the problem during the analysis. Common nodes link these elements together. These analyses provide an approximate solution to the problem in different loading types using the boundary conditions [32]. It is expected that the standard axle load of the tire (80 kN) [33] will be applied on the surface of the pavement with a unified contact pressure of one tire (690 kPa) with neglecting the stiffness effect of tire wall, the contact pressure will equal to tire pressure [34]. Measured stresses and strains have been applied to investigate the effect of various adoptive parameters on pavement response. All pavement models in this study simulated the pavement design of Indus Highway (N-55), Pakistan. Selecting the appropriate element type is a critical step in finite element modeling; ABAQUS provides a diverse elements library. An earlier study successfully employed an eight-node linear brick element in a pavement response analysis. To improve the convergence rate, eight-node linear brick elements with reduced integration (i.e., C3D8R) for various thicknesses were selected. A generated three dimensional mesh was designed to provide optimal accuracy for finer elements in regions requiring detailed load modeling and relatively large elements far from the loading area. In this model fine mesh was used in the load region and coarse mesh was used at a distance from the load region. A total of 108630 elements were used, and a mesh convergence study was conducted to determine the optimum number of elements. In this study the purpose is to examine the finite elements of the pavement analysis and study the effect of the isotropic material features on the pavement response to load in which all pavement layers were modeled and analyzed in both isotropic systems.

Materials Assignment
In the isotropic model, the elastic properties of the material, such as the Elasticity modulus and Poisson's ratio, are identical in all directions. This model has been used in numerous pavement analysis models. In the present model, the strain stress formula is used. Material properties assigned to pavement layers are shown in Table 1.

Loading Conditions
Tyre Pavement Contact Area Huang (2004) suggests simplifying the shape, which consists of a rectangular and two half circles. The simplified shape will have a width of 0.6 L and an equivalent area of 0.5227 L2. For the plane strain model, this assumption was made. Because the form of this contact region is not axisymmetric, this assumption cannot be applied in the axisymmetric model. Because some of the models utilized in this study are axisymmetric, a circular contact area with an equivalent size of 0.5227 L2 is used for the three-dimensional finite element studies, whereas a rectangular area with 0.5227 L2 is suggested and used for the two-dimensional finite element analyses. The multilayer elastic theory has previously used this assumption. Although it is not totally accurate, the inaccuracy is thought to be minor. A car wheel is generally represented by a rectangle for our convenience for the length and width shown in Figure 3. This is the tyre of the car whose width has been taken as 2.5 m and being modeled in ABAQUS along the path shown in Figure 4.  Moving load that varies in space and time has been simulated in this research using ABAQUS computer programme, which allows the user to apply several types of load, magnitude, location, and directions through the application of load module to determine the damage of rutting per each pass of wheel load on the flexible pavement. Load defined in AQABUS is converted into Equivalent Single Axle Load (ESAL) first which is 690 kPa and defined in the form of dynamic explicit in step module for which amplitude will be applied. 10000 wheel cycles are applied over the pavement layer to observe plastic analysis in pavement layer for 100 s interval which gives us 0.01 sec for a single pass of equivalent single axle load on the asphalt layer (Load will be applied over asphalt layer but transferred to other layers by applying certain changes in boundary conditions. The path of moving load is the longitudinal distance divided into multiple steps (10 steps complete one wheel cycle) to replicate the moving load on the pavement surface, which can be accomplished using the step module. Loading is applied over the contact area of the wheel and the spacing between both wheels is taken as 2.5 m which is the average width of the vehicle as shown in Figure 5. The edges for each part of the pavement geometry model are allowed to move vertically only and are fixed in the horizontal direction. The boundary conditions for the pavement model, which is fixed at the bottom of the subgrade layer (no horizontal and vertical movement), and the edges for each part of the pavement geometry model are allowed to move vertically only and are fixed in the horizontal direction. For this reason, bottom support is defined as encastre while other layers i.e. Asphalt, Base, Subbase are allowed to move vertically while fixed horizontally. Finally, as shown in Figure 6, the Interaction module of ABAQUS is used to simulate the interaction between distinct pavement layers (asphalt, base, and subgrade) with tie-contact for both (asphalt and base) and (subgrade and base) layers.

Results and Discussion
The finite element approach is used to perform flexible pavement analysis, and a three-dimensional finite element model is created using the ABAQUS computer application. The deformation, stress, strain, and rutting damage of pavements are evaluated using several forms of material characterization for the base and subgrade layers. As in rutting analyses, it was mandatory to get displacement "U2" along the wheel path so it was necessary to specify the output frequency in increments and request output according to a set of time points but instead of getting the only displacement, a request for stresses and strains are also submitted for getting output database.

Effect of Material Properties on Rutting
Three different materials having different properties are analyzed to check the effect of material on rut depth.

Material M1
Section analyzed for strains and displacement along horizontal and vertical directions is addressed. Deformation in the negative y-direction is of sole concern in rutting and is discussed briefly at the end of each analysis for each material property. Material properties assigned to each layer are shown in Table 2. Loads on the pavement's surface cause two stresses that are thought to be significant for design purpose in pavement analysis. These are the horizontal tensile strain at the asphalt layer's bottom and the vertical compressive strain at the subgrade layer's top. The strain induced is shown in Figure 7. The horizontal displacement of asphalt layer is kept ignored that's why it is kept as encastre for the horizontal direction (horizontal movement is restrained) when there will be no horizontal stresses and displacement in asphalt concrete top surface ultimately there will no stress transfer to other layers and no displacement in the x-direction. Figure 8 indicates the magnitude of displacement along y-direction in the pavement layer. It shows approximate rut depth measured in negative y-direction along the wheel path which is found out to be 0.07792 m or 3.07 inches. Figure 9 shows deformations vs no of load repetitions are plotted in which the total number of cycles is to be kept at 10000, In other words each wheel pass is happening at an interval of 0.01 sec which is pretty much high by keeping loading condition of Indus Highway in mind. In other words, 100 passes of standard ESAL are traveling on the 30 m section of the pavement layer in 1 sec. To prevent distortion in the pavement layer for such a greater number of wheel passes boundary conditions play a crucial role as it helped prevent particle dislocation by restraining stresses in the horizontal direction. Rutting is observed more for load 2000 to 4000 cycles and eventually decreases up to 10000 cycles and keeps on increasing afterward. The reason for this rutting increases initially because at first our pavement undergoes elastic deformation up to certain limit states and tries to regain its original position but afterward when No of load repetition increases after 10000 it undergoes permanent plastic deformation. Section analyzed for second material strains and displacement along horizontal and vertical directions is addressed. Deformation in the negative y-direction is of sole concern in rutting. The materials properties of each layer of pavement are shown in Table 3. Loads on the pavement's surface cause two stresses that are thought to be significant for design reasons in pavement analysis. These are the horizontal tensile strain at the asphalt layer's bottom and the vertical compressive strain at the subgrade layer's top. The strain induced is shown in Figure 10. Deformation (mm)

No of Load Repetitions
Rutting in pavement (mm) Figure 11 shows displacement "U" magnitude M2 decreases by approx. 11.65% compared to the displacement incurred in M1. The path selected on successive crests is taken on the same nodes for M2 as taken for M1 and horizontal movement is restrained. Figure 12 indicates the magnitude of displacement along y-direction in the pavement layer. It shows Approximate Rut depth measured in negative y-direction along the wheel path which is found out to be 0.06884 m (2.7 inches). Blue Lined wheel path suggests that there exists the least of stresses along this path because it undergoes plastic deformation and elastic stresses which are marked by red color pavement layers show that there are most stresses on these parts of the pavement.  For M1 max rut depth at 10000 cycles was 80 mm now for M2 rut depth decreases to 70 mm it means there occurred 12.5 % reduction in rut depth by changing base course crushed aggregate.

Material M3
M3 material is taken as the material which is used as a base stabilizer and can have very vast significances. Stabilizer being used can be a base stabilizer as well as a subgrade stabilizer. A base stabilizer is more effective in countering rutting comparatively to subgrade stabilizer. The base stabilizer used can be a cement base stabilizer or silica fume. In M3 silica fume is used as a base stabilizer and would be analyzed to check its effect on strains and ultimately on rutting. The base stabilizer is silica fume, which has certain qualities. The production of silicon metal or ferrosilicon alloys produces silica fume as a byproduct. Elastic Modulus has been calculated from CBR and 15 m 2 /g, and the specific gravity is 2.15 which increases by 75% by adding 10% silica fume in the base which is shown in Table 4.  So it has been observed that by using 10% silica fume the maximum dry density for modified soil has been increased to 11% at the same time it can be seen that CBR for 10 % silica fume increased by 75%. Using silica fumes as soil modifier improves the strength and bearing capacity for all tested samples. On the other hand effect of silica fume on rutting is most obvious as for 10000 cycles where rutting is observed to be at the peak, rutting for using untreated base layer was around 77.92 mm while by using treated based layer it reduces to 7 mm which result in 91 % reduction in rutting by the provision of 10% silica fume as a base stabilizer. The results after analysis are shown in Figure 13.

Effect of Temperature on Rutting
For T1=10 o C Approximate magnitude of displacement along 1, 2, 3 i.e. x, y, and z-direction can be shown in Figure 14.
Magnitude at different panels on the pavement layer can be described by a multi-colored dialogue box. Each color indicates the average magnitude of displacement for its specific color. Color change longitudinally from red to orange and yellow show that stresses are being transferred from the wheel path to the pavement surface. Stresses are transferred due to inter-particle connectivity as well as cohesion. The pavement section defined in Methodology is considered as temperature-dependent. First, the section is analyzed for the 10 ºC field taken as site temperature. So the temperature of the pavement layers is shown in the Table  5.  Figure 15 shows approximate rut depth measured in negative y-direction along the wheel path which is found out to be 0.077 m or 2.98 inches. Blue Lined wheel path suggests that there exists the least of stresses along this path because it undergoes plastic deformation and elastic stresses which are marked by red color pavement layers show that there are most stresses on these parts of the pavement. The results after analysis is shown in graph form as in Figure  16. Rutting is observed more for load 2000 to 4000 cycles and eventually decreases up to 10000 cycles and keeps on increasing afterward. The reason for this rutting increases initially because at first our pavement undergoes elastic deformation up to certain limit states and tries to regain its original position but afterward when number of load repetition increases after 10000 it undergoes permanent plastic deformation. In the temperature-dependent analysis the temperature is increased to 50 ºC for checking the effect of temperature on the rut depth while considering the material properties and loading condition the same as discussed in the above section. In most countries, the temperature varies between 10-50 o C. Approximate magnitude of displacement along 1, 2, 3 i.e. x, y, and z-direction can be shown in Figure 17. Magnitude at different panels on the pavement layer can be described by a multi-colored dialogue box. Each color indicates the average magnitude of displacement for its specific color. Color change longitudinally from red to orange and yellow show that stresses are being transferred from the wheel path to the pavement surface. Stresses are transferred displacement "U" magnitude T2 increases by approx. 3.24 % compared to the displacement incurred in T1. In the temperature-dependent analysis, the temperature is increased to 50 o C for checking the effect of temperature on the rut depth while considering the material properties and loading conditions the same as discussed in the above section. In most countries, the temperature varies between 10-50 ºC. So the temperature of the pavement layers is shown in Table 6.  Figure 18 shows approximate rut depth measured in negative y-direction along the wheel path which is found out to be 0.0798 m (3.2 inches). As from the analysis, it is shown that if the temperature increases from 10oC to 50oC the increase in rutting depth is 3 mm which is 9% increase in rut depth. The increase in rut depth due to temperature is changed but not significantly as the rutting due to applied is already large so that's why temperature effect is decreased. If the applied load on section decreases the effect of temperature increase on the rut depth increase significantly. The results after analysis are shown in Figure 19 and comparison is shown by Figure 20.  The above section defined is first analyzed under a wheel load of 690 kPa for ten thousand cycles. Approximate magnitude of displacement along 1, 2, 3 i.e. x, y, and z-direction can be shown in Figure 21. Magnitude at different panels on the pavement layer can be described by a multi-colored dialogue box. Each color indicates the average magnitude of displacement for its specific color. Color change longitudinally from red to orange and yellow show that stresses are being transferred from the wheel path to the pavement surface. Stresses are transferred due to inter-particle connectivity as well as cohesion. The pavement section defined in methodology is considered as temperature-dependent. First, the section is analyzed for 690 kPa. Figure 22 shows Approximate Rut depth measured in negative y-direction along the wheel path which is found out to be 0.077 m (3 inches). Blue Lined wheel path suggests that there exists the least of stresses along this path because it undergoes plastic deformation and elastic stresses which are marked by red color pavement layers show that there are most stresses on these parts of the pavement. The results after analysis are shown by Figure 23.  Displacement "U" magnitude L2 increases by approx. 5% compared to the displacement incurred in L1. Now the Load is increased to 750 kPa for checking the effect of temperature on the rut depth while considering the material properties and loading condition the same as discussed in the above section. Figure 25 shows Approximate Rut depth measured in negative y-direction along the wheel path which is found out to be 0.081 m or 3.2 inches. The results of analysis are shown in Figure 26.  Approximate rut depth measured in negative y-direction along the wheel path which is found out to be 0.0865 m (3.4 inches). As from the analysis it is shown that if load increases from 750 to 800 kPa the increase in rut depth is 3mm which is 6 % increase in rut depth. The increase in rut depth due load is changed but not significantly.

Conclusions
A 3D dimensional model of flexible pavement with proper material characterizations and bonding conditions is created using the finite element analysis approach. From finite element model of flexible pavement under different material properties, temperature and loading conditions, it is concluded from the results of the analysis of three models:  Rutting is observed more for load 2000 to 4000 cycles and eventually decreases up to 10000 cycles and keeps on increasing afterward. The reason for this rutting increases initially because at first our pavement undergoes elastic deformation up to certain limit state and tries to regain its original position but afterward when number of load repetition increases after 10000 it undergoes permanent plastic deformation;  For M1 approximate rut depth at 10000 cycles is 77.92 mm and for M2 rut depth decreases to 68.84 mm it means there occur 11.65 % reduction in rut depth by changing base course crushed aggregate;  For M1 max rut depth at 10000 cycles was 80 mm now for M2 rut depth decreases to 70 mm it means there occurred 12.5 % reduction in max rut depth by changing base course crushed aggregate;  Effect of silica fume on rutting is most obvious as for 10000 cycles where rutting is observed to be at the peak, rutting for using untreated base layer was around 77.92 mm while by using treated based layer it reduces to 7 mm which result in 91 % reduction in rutting by the provision of 10% silica fume as a base stabilizer;  The rut depth increase with increase in temperature if the loading condition and material properties and other parameters are kept constant. The rut depth is increases as we can see in first analysis it is up to 77 mm and second analysis rutting is 2.5% but in next same temperature interval between higher temperatures the effect is decreased to 1%;  For L1 Approximate Rut depth measured in negative y-direction along the wheel path which is found out to be 0.077 m or 3 inches. For L2 Approximate Rut depth measured in negative y-direction along the wheel path which is found out to be 0.081 m or 3.2 inches;  For L3 Approximate Rut depth measured in negative y-direction along the wheel path which is found out to be 0.0865 m or 3.4 inches. As from the analysis it is shown that if the loading increases from 750 kPa the increase in rut depth is 5.5 mm which is 6 % increase in rut depth. The increase in rut depth due to load is changed but not significantly.

Recommendations
 Instead of using untreated base layer it is most preferably to use treated base layer by mixing cement stabilizer or silica fume because silica fume added 2% , 4 % and onwards results in the reduction of rutting;  Crushed aggregates with more density are generally to be preferred in base layer rather than subgrade layer. Base layer materials affect rutting more than subgrade layer i.e. 65 and 10% respectively;  The temperature has a significant effect on rut depth so the material which effect significantly by temperature like bituminous material care should be taken while using this material to reduce the effect in flexible pavement;  According to NHA (National Highway Authority) Specifications (2005) 75mm rutting is of very high severity so complete overlay may be consider for the already existing pavement;  The load applied is not suitable for existing structure capacity therefore load should be decrease or quality of material used may be increased.

Data Availability Statement
The data used in this research study for different analysis are available in the article.