Investigation of flexible pavement maintenance patching factors using a finite element model

Patching of flexible pavements is one of the most important functions of pavement maintenance. Although finite element modeling has become commonplace in the world of pavement engineering, modeling has not yet been significantly leveraged for maintenance applications which improve safety, ride quality, and pavement service life. The objective of this study was to model viscoelastic properties of pavement and patching materials to determine the effect of various repair factors on pavement performance using the finite element method. Specifically, surface permanent deformation, local shear stress concentration, and horizontal strain distribution were investigated. Two types of models were simulated; the first model applied static loading to a surface layer fixed on a plate and the second model applied cyclic traffic loading to a two-layered flexible pavement system. The results demonstrate the importance of patching using a semi-permanent method. The results also demonstrated the accumulated effect of repeated loading using a time-dependent material response. Results also indicated that a larger patching area resulted in less influence of the shape of the area, while a circular area proved superior to a conventional rectangular patch for sizes near the tire footprint and smaller than it. Different responses were observed depending on the type of patching material modeled, demonstrating the effect of material choice in maintenance applications. Finally, mesh optimization was performed to ensure appropriate mesh sizes are used in future studies to accurately represent the pavement layers and patches.


Introduction
Potholes are one of the most critical and widespread types of damage encountered by flexible highway pavement agencies. The most typical causes of pothole formation are freeze and thaw cycles, although other factors such as construction, material quality, and traffic loading can play roles in all weather. Loading and temperature related stresses cause cracking on the surface of asphalt layer, where moisture is prone to penetrate, freeze, and widen the previous cracks. However, as heating and cooling flux cycles occur, evaporated moisture leave holes on the (2021) 2:30 Page 2 of 16 materials selected. cold mix asphalt (CMA) is a commonly used material in maintenance applications which, like hot mix asphalt (HMA), consists of mineral aggregates and a bituminous binder. Especially for cold climates when hot mix asphalt plants are closed, cold mix often becomes the main alternative. CMA is a material with viscoelastic, or time-dependent properties similar to HMA. However, there are few published studies focusing on the effects of viscoelastic properties of cold mix on pavement repairs. One study from Shanbara et al. modeled fibrereinforced and unreinforced cold mix asphalt to evaluate the effectiveness of reinforced cold mix asphalt in mitigating rutting and expanding pavement service life [3].
Another study by Shabara et al. [4] presented a cold mix design method and used Finite Element Analysis (FEA) to model the designed mixes. Moreover, no matter what materials are employed for patching, stress concentration is a concern at the boundary between existing materials and new patched materials. Ghosh et al. conducted two-dimensional finite element simulations to determine the effect of geometric configuration of pothole on reducing stresses at the boundary between existing pavement and the patching area, considering the pavement linear elastic materials [5] and used a planar shell model [5]. Other studies have also focused on modeling the repair of pavement materials. Chen et al. [6] proposed a 3-D FEM to investigate the mechanical response of patching in partial depth repair of concrete pavements. They found that tensile stress depends on repair depth, damage location along the slab, and material properties. Patching materials were recommended based on the slab location. Pratico et al. [7] focused on the model of macrotexture of flexible pavement. They simulated common surface treatments based on different material property inputs. This was also validated with experimental results. On the side of modeling shear stresses, Ling et al. [8] utilized 3-D FEM to simulate and investigate the top-down crack propagation. And they claimed that pavement top-down cracking depends on various factors. They focused especially on shear stress caused by loading at the edge of or within the tire contact area.
The present study is focused on investigating the effect of vertical tire pressure on patched areas. Static and cyclic loading conditions are applied to a patched pavement section modeled using a finite element engine. Unlike past studies, this study evaluates the surface layer in a layered system from both viscoelastic (time-dependent) and three-dimensional perspectives. In addition to determining the three-dimensional viscoelastic response, a second goal of this study is to understand the effect of the size of the repaired volume. Different materials are evaluated to determine the impact of changes in the interface between existing pavement materials and patching materials, and the effect of using elastic versus viscoelastic patching materials. Finally, the authors focus on developing a mesh optimization framework to ensure that future studies are conducted with the right level of fineness to ensure reliable results while maintaning computational efficiency. This paper represents part of an ongoing effort by the researchers [9] to better leverage finite element modeling for assisting roadway maintenance engineers with decision making. This study builds on the previous study by refining the models to be more realistic, evaluating more combinations including the most critical patched area size, and performing mesh optimization and comparison analysis to mechanistic strain analysis solutions derived via elastic layer theory.

Viscoelastic modeling of asphalt materials
Asphalt binder presents its viscoelastic properties as a combination of an elastic solid and a viscous fluid. The binder's viscoelatic properties are transferred to the mix, resulting in a visocelstic mix which can be represented by mathematical models. The viscous part is represented by a dashpot, whose viscosity is η i . The elastic part is represented by a spring, whose elastic modulus is E i . Among different mathematical models, the generalized Kelvin-Voigt model which can also be called the Voigt model can be used to represent the time-dependent compliance and modulus. The generalized nature of the model means it represents the material with a series of individual Voigt elements which consist of an elastic spring and Newtonian dashpot in parallel, since one single Voigt element is not enough to describe the long period behavior [10]. Another commonly used model is the generalized Maxwell model, which aims at representing the relaxation modulus and represents the material using springs and dashpots in series. Due to its efficiency in modeling linear viscoelastic behavior and its compatibility with Abaqus software, the generalized Maxwell model is applied in this study and is represented by a Prony series [11]. The resulting formulation [12] is Eq. (1).
where E(t) is the relaxation modulus; E ∞ is the equilibrium elastic modulus, E ∞ = lim t→∞ E(t); t is the loading time; ρ i is the relaxation time; Eq. (1) is also written as Eq. (2) in Abaqus.
where G(t) is the shear relaxation modulus; G 0 is instantaneous shear modulus; G 0 = lim t→0 G(0); t is the loading   [11] and fatigue tests [15]. Through Laplace transformation and intercoversion [16][17][18], the creep compliance results can be converted into corresponding relaxation results. Furthermore, if with relaxation tests results, parameters including shear modulus components and relaxation times required to fit a Prony series can be obtained by fitting test data [19,20].

Materials characterization
Since the study focuses on the investigating pavement response with and without patching materials, the surface layer, where maintenance typically occurs, is the main study object. Therefore, the existing surface layer was simulated as viscoelastic HMA materials and a conventional aggregate base layer was simulated as a linear elastic material as shown in Table 1. Material properties for input were obtained from existing literature in which FEA simulations were performed on conventional flexible pavement structures. The effect of temperature was not considered in this study. Therefore, materials properties were adopted at 25°C. For a preliminary evaluation of the effect of patching, CMA was the patching material in most of the following simulations. However, to consider the implications of patching with different repair materials, various patching materials, including CMA, HMA, and epoxy materials (properties shown in Table 2), were compared for a sensitivity analysis, described later in this paper.

Finite element model simulations
The finite element method is one of the most successful techniques to predict and analyze flexible pavement behavior [3]. Finite element models developed in the present study are twofold: 1) A surface layer fixed on a plate and 2) a two-layered system with a surface layer and unbound base layer. The first model was employed to investigate the deformation of patching on the surface layer, with the goal of determining the effect of patching potholes on reducing stresses and further defect formation. This model also assesses the horizontal stress as a comparison with the work of Ghosh et al. [5]. The second model was established based on repaired flexible pavements, simulating a series of loading-and-unloading cycles to determine the effect of accumulated stresses when various factors were changed.

Model geometry and boundary conditions
Based on the nature of pavement structure, continuum 3dimensional 8 node solid elements (C3D8) were selected. Because of the symmetry of loading and model geometry and to reduce computational effort required, the pavement section was simulated as one quarter section applied with symmetrical boundary conditions [23]. For part 1, each quarter section was modeled as a 50 mm thick, 600 mm by 800 mm rectangle. The tire imprint was a 216.4 mm (8.52 in) by 149.0 mm (5.87 in) rectangular area, which is discussed further in the loading condition section. When evaluating the importance of patching, the model compared the structure with a circle patch of 149 mm radius and without any patch, as shown in Fig. 1.
When assessing the horizontal stress as the comparison with previous study, the scale of repairs was the same or similar to the tire imprint at a 149mm radius circle or a 149mm by 216.4 mm rectangle, as shown in Fig. 2.
In the second phase of the study, each quarter section was modeled as a 50 mm thick, 600 mm by 800 mm rectangle with a filled volume of various sizes and shapes as descried below. Two shapes for the repair area were evaluated, rectangular and circular. Small repairs were considered as those smaller than the tire footprint, with dimensions of one quarter size being an 81.6 mm radius circle or an 81.6 mm by 149 mm rectangle. Medium size repairs were the same or similar to the tire imprint at a 149mm radius circle or a 149mm by 216.4 mm rectangle.
Large size repairs were much larger than the tire imprint area. For large repairs, due to meshing problems, the longitudinal dimension instead of transverse dimension was selected as the reference. Therefore, large repairs were a 298mm by 432.8 mm rectangle and a 432.8 mm radius circle. The plan view of this series of model are elaborated in Fig. 3. Corresponding with the description above, green parts in this figure represent the tire imprint, the yellow parts represent the repairs. Therefore, with rectangular shape and circle shape matched with small, medium and large patches, six models were established in this part.
Based on previous studies [2] and manuals [24], semipermanent patching is a common method for repairing potholes. In this patching procedure, workers do not dig deeper material out of the pothole, and patching is placed in a depth consistent with the depth of the holes. Therefore, although it is critical in terms of effect on repair performance, thickness is not evaluated in this study due to the difficulty of controlling it. In this study, to avoid the thickness effect on pacthing, patches are placed at the same depth.
The interface is assumed as a rough interaction between the original surface and patching material. When the model was established, the merge function was used to connect the existing region and patching region, which ensures that the interface is fully bonded. It is important to note that the interface plays an important role in impacting the patch performance (strains, stresses

Loading conditions
Two loading conditions were evaluated. For part 1, the static loading is applied for 60 sec for evaluating the permanent deformation of the surface layer. The loading area is located close to the pothole area, shown in Fig. 1, due to the impossibility of applying loading on a empty pothole in Abaqus. Additional modeling was performed in which the load covers the patched area, as shown in Fig. 2.
In the second part of this study, to enable better future comparison of our results with laboratory tests such as the wheel tracking test (as done by [3]), the applied loading frequency models a wheel passing over the section at a speed of 0.17 m/s. In terms of loading magnitude, based on the work of Huang et al. [10], the standard dual-tandem with 40 kips loads, consisting of 0.7MPa (100 psi) of each tire imprint, is applied on a 216.4 mm (8.52 in) by 149.0 mm (5.87 in) rectangular area. Also, the loading scheme for this series of models is fixed while the patch shapes and sizes change as discussed before. The pressure is applied regularly for 41 cycles, shown in Fig. 4. Previous research has shown that vertical pulse loading is an effective way to represent moving vehicle loads on flexible pavement structures [25]. Within each cycle, which is 1.44 s, loading pressure is applied for 0.17s followed by 1.27s relaxation time.

Simulation results
Static loading simulations Figure 5 shows the maximum deformation of the surface layer exposed to static loading both with and without patches. Figure 6 illustrates the maximum displacement along the bottom edge of the tire imprint area. From the comparison, the pavement layer with the patches lowers the maximum displacement, which happens at the left edge adjacent to the wheel imprint, while the layer with the unfilled pothole deforms two times larger than the one with the patches at the left edge of the pothole. Therefore, the static loading condition model confirms the importance of patching pavement defects. After understanding the necessity of patching, the horizontal stress analysis of the model surface is conducted for patched sections. Figure 7 reflects the effect of different shapes of patches on horizontal stress. Rectangular patches result in higher local stress concentrations around the imprint area than circle patches, which is consistent with the findings of Ghosh et al. [5] based on a 2D tensile elastic model.

Cyclic loading simulations
For small dimensions of holes in Fig. 8, with regard to horizontal strains, the most critical location is at the bottom of the surface layer at the end of the loading period. Among small repairs, a circular shape helps reduce the horizontal tensile strain at the bottom the surface layer and the bottom bonding area, as shown in Table 3.
As cycles increase, the cumulative effect shows that the tensile strain increases after each loading period, while decreasing after each recovery period. For deformation comparisons, the deepest rutting happens at the center of the repaired area due to the higher flexibility of CMA compared with existing aged HMA. Even though the there is less difference in terms of rutting between rectangular and circular repairs, heaving taking placing at the imprint area at the end of all cycles for circular repairs is much smaller than that of rectangular repairs, which is consistent with the findings of Shanbara et al. [3]. Shear stress is also considered an evaluation index since the loading is vertically applied to pavement, which can result in formation of surface defects or top-down cracks. In order to make sure the effect from the tire imprint edge is analogous in both models, shear stress S13 is selected as an evaluation factor such that the horizontal distance of the edge of patching areas to the edge of the tire imprint is the same for both rectangular and circular repairs, as shown in Fig. 8. The results indicate two areas of interest within the structure. One of the most critical areas is located at the middle depth of the surface layer, while the other critical part is around the bonding area. Results in Table 3 show that at the middle part of the surface layer under the right edge of the imprint, geometric shape does not play a big role in influencing shear stress, but at the interface area, the circular shape helps reduce stress at both the surface and bottom. For medium sized patching areas, Fig. 9 and Table 4 show horizontal strain in the surface layer and surface vertical deformations follow the same trends as the smaller patched areas. Meanwhile, the second most critical part is located at the surface interface between patches and the existing pavement, where the circle shape shows more benefits for reducing tensile strains. This is a crucial difference between medium sized repairs and the smaller patches and the reason that surface bonding area is investigated in Table 4. As the number of loading cycles increases, the advantage of the circular shape becomes more prominent, reflecting the importance of using viscoelastic analysis to understand pavement maintenance treatments. From Table 4, at the imprint area, the circular shape reduces the shear stress only after the last cyclic recovery period by 10.47% compared to the rectangular shape, whereas at the bonding area, the circular shape helps to decrease shear stress by 14.23% and 23.39% at the end of the final load cycle and recovery period, respectively. Figure 10 shows that for a large patching area, the horizontal strains at the interface between patching material and existing pavement at the bottom of the surface layer does not show an explicit difference between circular repairs and rectangular repairs. Therefor, horizontal strain comparisons were made only in terms of imprint areas, instead of at bonding areas, in Table 5. In this table, the tension zone refers to the area directly below the loading area. The adjacent area is referred to as the compression zone due to the presence of the compressive strain. In both the compression zone and tension zones, the circular shape reduces the horizontal strain throughout the simulation. Deformation caused by tire loading is independent of the shape of repairs because the loading contact areas are almost identical for both shapes. However, the heaving deformation at the location of bonding caused by circular patches is still smaller than that caused by rectangular patches at the end of the recovery period. As for the shear stress, S23 was selected as the indicator in order to keep the same longitudinal distance between the edge of patching areas to the edge of the wheel imprint for both shapes. Based on the results in Table 5, circular repairs perform better than rectangular repairs for both the surface bonding area and under imprint area. Figure 11 shows the responses of pavements with circular and rectangular repairs at the bonding areas. Regardless of the repair shape, responses at the bonding area, including horizontal strain, surface deformation and vertical shear stress, indicate that the medium size repair results in the worst case scenario, because the tire loading imprint's edge is close to the boundary between existing pavement and repair material, where critical shear can lead to debonding. The largest repair areas generally result in the least shear stress at the interface. However, regardless of the repair shape, large repairs results in the largest level of permanent deformation, followed by medium and then small repairs. Overall, this indicates that recurrent failure is a patch may occur most often for areas of similar size to the tire footprint, but large areas may be the most susceptible to future structural issues.

Mesh sensitivity analysis
Because the mesh size in FEA highly influences the accuracy of the results, to determine the optimum mesh size and distribution, a mesh sensitivity analysis was conducted. There are several methods to analyze mesh sensitivity, including assessing the convergence of simulation results by changing element distributions and comparing the value with other numerical analysis methods. According to the method used by Hernanadez et al. [26], the total strain energy of the model was utilized as a criterion to evaluate the optimized mesh configuration. However, because the present model involved cyclic loading, the total strain energy changed with time. Therefore, in the present study, in order to compare the total strain energy under the same condition, the total strain energy at the time when the energy reaches its maximum value was measured for all mesh configurations, and is presented in Fig. 12. The normalized total strain energy as a function of the number of elements is plotted in Fig. 12. As the number of elements increases from Mesh 0 to Mesh 1 to Mesh 2, the total strain energy tends to be stable at 95% of that of Mesh 0. Moreover, Mesh C was created to represent a coarser mesh distribution scenario than Mesh 0, which only contained 464 elements. The normalized total strain energy claims that a higher and more unstable total strain energy comes from coarser mesh condition as expected. Therefore, the conclusion can be drawn that increasing the elements in materials bonding area appropriately leads to the convergence of the results, but that Mesh 0 is suitable for the purposes of this study.
The other method is to use comparison with the mechanistic computational results of WinJULEA, a well-known and commonly used pavement analysis software. Because the mesh sensitivity analysis was the goal in this section, the former material properties in FEA were all converted to elastic propertiwa in the mesh sensitivity analysis Note: Comparisons indicates the percent improvement of the circle shape compared to the rectangular one. Positive number means circle shaped patch is performs better than rectangular one for the corresponding impact factor. And the larger the positive percentage is, the more noticeable this corresponding improved effect is model to be consistent with the WinJULEA input. Besides materials, the loading area was simulated as a circular area in FEA, which is what is used in WinJULEA. As with the above method of evaluating convergence, three similar types of mesh arrangement were established to observe the trend of the difference of FEA results from that of WinJULEA results. Table 6 presents that Mesh 2, the densest mesh, reduced changes of tensile strain at the bottom of surface layer and shear stress at the edge of the tire decrease to 10% and 6%, separately compared with those of original Mesh 0. This confirms the mesh quality improves as the number of elements are increased. However, when comparing the surface displacement, difference between FEA and WinJULEA increases as the mesh is densified, although the difference is largely negligible. Given that WinJULEA was based on semi-infinite space as boundary condition, while the FEA in this study modeled one quarter of pavement structure to represent patched area combined with existing pavement, some error is expected as a baseline due to the boundary conditions. Likewise, Mesh C was simulated to compare with WinJULEA results to represent a worse scenario. Comparison results indicates significant change in terms of above three indicators. However, it is overall critical to note that densifying the mesh can improve the results based on this analysis. It is the intention of the authors to conduct such analysis using finer meshes to determine if the results of this study are confirmed; however, this analysis is beyond the scope of this study due to limitations on available computational resources.

Asphalt interface model
In the present study, the bonding area was simulated in the way that assumed full bonding, even though in practice, the bonding property relies on the adhesive nature of the materials and good construction practices. However, the effect of usage of the merge function and various tangential interactions between two parts were not yet well-established for pavement FE simulations. The former assumed that repairs and existing pavements were combined in one integrated instance, while the latter assumed that two parts interacted with a specific interaction property. Therefore, to determine the effect of interaction properties on shear stress within pavements, different tangential properties of the interface interaction were investigated, as shown in Table 7. In addition, the effect of materials was investigated here by using HMA with similar properties to the existing pavement as a baseline for comparison.
The results of shear stress at the bonding area revealed almost the same results using the merge function and a rough (very high friction) interaction as the tangential property. Because the Coulomb model was used to describe the friction, if the coefficient friction is 0.2, the shear stress at the bonding area was reduced to nearly half of that of the full bonding condition. According to    Yoo et al. [27], the common friction coefficient should range from 0.2 to 1.0. Additionally, if the interface of repairs and existing pavement was simulated as frictionless bonding, the shear stress at the bonding area was much lower than any other simulations with friction at the interface. The frictionless simulation also led to shear stress to be less than that of pavements patched with the same HMA materials. Therefore, using frictionless bonding as the interaction was unreliable, because the two parts were completely separated apart in the frictionless condition. Patching with the same HMA as existing pavement leads to less shear stress compared to the same rough condition with CMA. This indicates that shear stresses also depend on the patching materials. In addition, as the mesh density increases from Mesh 0 to Mesh 1 to Mesh 2, the shear stress at bonding tends to reach a stable value, which is consistent with the conclusions of the "Mesh sensitivity analysis" section.

Material sensitivity analysis
HMA materials and epoxy materials were applied in similar models to understand the effect of material choice on performance of a patched pavement structure. Epoxy materials were represented as linear elastic to understand the benefits of patching with non-bituminous elastic materials. To minimize the effect of varying sizes and shapes, only circular patches with a smaller size than tire imprint were employed in these simulations. At three critical times, the vertical surface displacements and vertical shear stresses at the bonding area were compared in terms of materials in Fig. 13. Note that horizontal tensile strain at bottom of the surface layer and surface vertical deformation of HMA are less than that of CMA. Horizontal tensile strain at bottom of the surface layer and vertical deformation of epoxy are less than that of HMA. These differences are increased by additional cycles of loading, indicating again the importance of using viscoelastic modeling for CMA and HMA. For Fig. 13c, since HMA is similar to existing pavement materials, HMA results in the least shear stress at the bonding area, while CMA patches lead to the largest shear stress at bonding among three patching materials. Therefore, it can be concluded that the stiffness differential between repair materials and existing pavement materials plays a significant role in influencing bonding quality.

Conclusions and recommendations
In this study, viscoelastic models of pavement sections with and without patches of different types were established using an Abaqus finite element model. Observations of deformations, stresses and strains within the existing pavement and patch result in the following conclusions: 1. Short-term loading from the static loading condition test, which represents the response to moving traffic loads best, indicates large differences in strains between patched and unpatched sections. This indicates the importance of patching pavement defects to protect a pavement structure and surface from further damage. 2. Within cyclic loading simulations, circular patches of the same size as tire imprints or smaller show a reduction in bottom tensile strains, shear stresses at bonding areas, and surface displacement (especially heaving) compared to rectangular repairs. However, circular patches do not lead to an obvious improved influence on shear stresses under the edge of the tire imprint area. Large circular repairs reduce horizontal strain and vertical shear stress, while resulting in less impact on vertical displacements.  3. When a patch size is close to the size of the tire imprint, the patching area is more vulnerable because the effect on the two most critical regions (imprint area and bonding area) overlap, which greatly influences the patch performance. These differences are often compounded by modeling the effects of time dependent material behavior over many loading cycles. This indicates the importance of using a viscoelastic model to capture the material behavior. 4. Stiffer materials may cause higher shear stress at the interface between existing pavement and patches, resulting in increased risk of debonding. However, using a stiffer patching material results in less strain and rutting. Overall, these conclusions have important practical implications. Based on the predicted mechanical response of flexible pavements when exposed to patching, patching in a circular shape should be considered for patches of similar size or smaller than the tire imprint. For a larger patch, the shape of is not as important. In addition, a selection of patching materials with high modulus is recommended when patching. Finally, defects similar to the tire size may be the most critical to patch, since they create high stresses.

Future studies
In future studies, the viscoelastic properties of multiple sources of CMA and HMA obtained from lab tests can be used to construct and show a clearer and more reliable comparison using the developed model. Laboratory tests such as the wheel tacking test should be considered as a comparison to verify the reliability of the simulation in a laboratory setting. Future studies should also focus on optimizing the simulated bonding condition based on a realistic representation of the default bond between different materials and the use of tack coats. Frequency and loading position were also not considered, as they are outside the scope of this study, and evaluating them would require more computational resources than available. These should be a subject of future research in this area. Considering the computational costs and diverse models, the current element sizes and distributions has provided sufficient numbers of elements for reasonable results. However, finer meshes with smaller and more elements are expected to provide more accurate results. Infinite elements can be used in FEA to better simulate the field condition and help promote the comparison with WinJULEA computational results.