U.S. Department of Transportation
Federal Highway Administration
1200 New Jersey Avenue, SE
Washington, DC 20590
2023664000
Federal Highway Administration Research and Technology
Coordinating, Developing, and Delivering Highway Transportation Innovations
This report is an archived publication and may contain dated technical, contact, and link information 

Publication Number: FHWAHRT10065 Date: December 2010 
Development of a useful simulation environment to study field compaction of HMA provides a unique challenge because of the complex physical interactions involved. A framework composed of a nonlinear material model coupled with an advanced numerical simulation environment (CAPA 3D) was chosen for this study. The theoretical model employed for the FE model development was adapted from the HMA material model presented in Koneru et al.^{(48)} This model was subsequently utilized for the simulation of the compaction mechanism of the SGC. The results of the study have been documented in Masad et al.^{(49)} Here, the chosen numerical framework is used to extend the work presented in Koneru et al. and Masad et al. to specifically study the material model's behavior when subject to the conditions pertinent to field compaction of HMA.^{(48,49)}
The simulation of the rolling compaction process includes many challenges in terms of representing the complex physical interactions in a numerical model. Some of the chief factors that influence the research related to the field compaction simulation are as follows:
A typical pavement structure used for simulating field compaction is shown in figure 61. It consists of the following three layers: a base layer consisting of aggregate and/or soil approximately 18 inches (457.2 mm) thick, an overlay on top of the base layer consisting of old HMA approximately 2 inches (50.8 mm) thick, and a 3 inch (76.2 mm)thick top layer of freshly laid HMA that will be compacted by rolling.
Figure 61. Illustration. Pavement structure typically employed for studying field compaction.
The general structure of CAPA3D facilitates quick, flexible testing of material models and has been presented in detail by Scarpas.^{(53,54)} For the implementation of the chosen material model, two interface subroutines are input into CAPA3D. The necessary model parameters, boundary conditions, and algorithmic logic necessary for the required stress updates at each material point are provided to CAPA3D through userinterface files for each type of necessary data. The required stress update subroutine was devised from the equations obtained through a numerical discretization of the constitutive relations of the material model chosen. The stressupdate subroutine consists of a system of nonlinear algebraic equations that are iteratively integrated at each materialintegration point and each step in the time domain. The computational sequences devised are then provided as Fortran instructions that can be embedded directly into the CAPA 3D interface. Details on the usage of the computational system are provided by Scarpas.^{(53)}
Some of the important characteristics of the FE mesh used for simulating compaction (shown in figure 62) are as follows:
Figure 62. Illustration. Sectional view of the FE mesh used for setting up field compaction simulations.
Table 4. Typical material properties of pavement structural layers.
Layer 
Model 
Thickness 
Properties 

New compacted asphalt layer 
Compressible nonlinear viscoelastic 
3 inches initially, changes during compaction 
Presented as parameter sets later in this chapter 
Old asphalt layer 
Elastic 
2 inches 
E = 300 ksi, = 0.35 
Base layer 
Elastic 
18 inches 
E = 145 ksi, = 0.20 
Boundary Conditions for Simulating Field Compaction
The boundary conditions most frequently utilized are presented in figure 63. The schematic illustrates that the fixed edge corresponds with the inner edge of the pavement lane while the free edge corresponds with the outer edge of the lane. The displacement on the fixed edge is thereby constrained in the lateral direction (zcoordinate direction for the given mesh orientation), the free edge is allowed to move unconstrained, and a tractionfree boundary condition is employed.
When the middle lane of a highway is to be simulated, simplifying further requires assuming symmetry in boundary conditions and employing a half mesh along the direction of rolling. In this case, the halfsectional mesh represented in figure 62 is sufficient for the simulation.
Figure 63. Illustration. Schematic diagram illustrating the edges of the lane that correspond to fixed and free edges of the mesh in figure 62.
This section discusses the typical response of the material to a moving load in field compaction. The simulations were performed with material parameters that correspond to those in table 5. The same set of parameters in SGC and field compactions do not provide the same level of volumetric reduction because of significantly different boundary conditions such as confinement, temperature changes, uniformity and nonuniformity of applied loads, etc.
Table 5. Parameters used to study the typical material response during field compaction.
Highway Project 
(MPa) 
n_{1} 
λ_{1} 
q_{1} 
(MPa·s) 
n_{2} 
λ_{2} 
q_{2} 

SH21 
620 
4 
0.25 
−15 
1,400 
2.5 
0.25 
−20 
The typical response curve obtained from the simulations employing the previously described mesh subject to the specified boundary conditions is shown in figure 64. The response curve represents the compaction experienced at a point in the roller path as the roller passes over it and returns for a second pass from the opposite direction. The annotation arrows signify the forward motion in the first pass and the reversedirection motion in the second pass. Finally, the load is removed and the material allowed to relax as much as possible. Thereafter, the material finally exhibits its state of permanent deformation, as seen in figure 64.
Figure 64. Chart. Typical displacement curve for a node under the cylindrical load for a cycle with forward and return passes followed by a forward pass with the load removed.
Figure 65 and figure 66 illustrate the relationship between the roller position with a contact area of 3.9 inches (100 mm) and the deflection experienced by a point P near the surface of the pavement. The illustrations are for a mixture that does not have instantaneous elastic response. Due to the viscoelastic nature of the material, the deflection at point P starts when the center of the roller has passed the point. The point continues to deform vertically as the roller moves forward.
When the roller is at point 4, the mixture at P rebounds some of its vertical deflection because it is out of the region of influence of the roller. Also, point 4 marks the time at which the roller starts to move back toward P. Further deflection occurs as the roller again approaches P (point 5). The final deflection is marked by point 6.
Figure 65. Chart. Vertical displacement observed at point P on the path of the roller.
Figure 66. Chart. Roller location from observation point P.
Each pass over the material results in permanent deformation upon unloading, as shown in figure 67. The deformation therefore builds up incrementally with each pass. However, as shown in figure 68, the increments decrease in magnitude due to material densification as the number of passes increases.
Figure 67. Chart. Deflection at the node of interest when a load is applied for a short duration and then removed.
Figure 68. Chart. Permanent deformation prediction of the model in multiplepass loading.
Implementing a ContactArea Algorithm
The objective was to implement an algorithm for load application on flexible pavement. Development of the new algorithm furthers the ability to simulate the compaction loading. The loading algorithm used in this project was based on translating a distributed load across the pavement.
The load was applied in the form of a grid of loading patches independent of the FE mesh. The pressure in each loading patch was uniform. By combining a suitable number of loading patches, any desired distribution of surface loads could be simulated. The loading patches were translated along the surface of the pavement to simulate the moving roller.
A process was developed for determining the pressure of each of the loading patches. The process involves the determination of the contact area between the roller and the pavement (see figure 69). A displacement that mimics the cylindrical profile of the roller is applied on the surface of the pavement. The number of nodes in the contact area is determined such that the summation of their nodal forces is equal to the weight of the roller (see figure 70). Different numbers of elements are maintained in contact to balance the dead weight of the roller.
Figure 69. Illustration. Roller contact geometry for static indentation and during motion.
Figure 70. Chart. Indentation of a cylindrical roller into pavement.
The plots presented in figure 71 represent the nodal reaction forces in the elements in response to the applied dead load of the roller as the contact area (represented by the number of elements in contact) is changed. A1 and A2 represent the edge and middle nodal forces of the same orientation in an element when considering six elements in contact between the roller and the material. Similarly, B1 and B2 represent the respective nodal forces at the edge and middle nodes of an element when four elements are chosen in contact, and C1 and C2 are nodal forces when two elements are chosen in contact. On the basis of the nodal forces in the contact area, the pressure in each loading patch was determined according to FE procedures.^{(56)}
Figure 71. Chart. Change in nodal reaction forces as the load is applied over a smaller area.
Once the compaction load is determined, the load is translated over the pavement, providing a loading effect similar to sliding a roller over pavement. The exact mechanics to establish a rollingtype contact between the roller and the nonlinear pavement material require indepth exploration of the field of nonlinear contact mechanics and has not yet been implemented in the FE model developed for this project.
The loading pattern was used to simulate the compaction of asphalt pavements with a roller travel speed of 3.0 mi/h (4.8 km/h). This speed is in the range of values typically used in the field and reported in table 6.
Table 6. Typical rollerspeed ranges.^{(1)}
Type of Roller 
Breakdown 
Intermediate 
Finish 

Static steel wheel 
2.0–3.5 mi/h 
2.5–4.0 mi/h 
3.0–5.0 mi/h 
Pneumatic 
2.0–3.5 mi/h 
2.5–4.0 mi/h 
4.0‑7.0 mi/h 
Vibratory steel wheel 
2.0–3.0 mi/h 
2.5–3.5 mi/h 
Not used 
The plots in figure 72 though figure 74 are indicative of the differences observed in the displacement and stress response in an element in the direct path of the traversing load using the currently developed contactarea algorithm and the previously available loading algorithm in CAPA3D, where the load is applied uniformly over all the finite elements in contact with the roller. The plots are representative of the response of the material at a node directly under the load (experiencing maximum normal distress) at the beginning of the simulation. The stress response of the pavement while the new cylinder indentation algorithm (nonuniform pressure) is employed is characteristic of the nonuniform stress reactions obtained in a 20noded element, whereas the algorithm previously in use strictly enforced uniform pressure at all nodes under a loaded area at a given instant of time. This, therefore, validates the choice of the use of the cylinder indentation algorithm to better represent the physical contact situation between the roller and the pavement. The simulation results presented in this section utilize the model parameters for the SH21 project simulations for the SGC (see Masad et al.)^{(49)} The material parameters used for these comparisons are presented in table 7.
Figure 72. Chart. Comparison of vertical displacement for the two different patterns of loading.
Figure 73. Chart. Comparison of the normal stress distribution at the pavement's top surface due to different loading patterns.
Figure 74. Chart. Comparison of the shearstress (XY) distribution of the pavement's top surface due to different loading patterns (on the plane of symmetry).
Table 7. Material parameters used to investigate the contact algorithm, interface effects, impedance element effects, and the effect of base stiffness.
Highway Project 
(MPa) 
n_{1} 
λ_{1} 
q_{1} 
(MPa·s) 
n_{2} 
λ_{2} 
q_{2} 

SH21 
2,500 
4 
0.25 
−25 
2,000 
5 
0.22 
−30 
Normal Stresses Predicted by Model
The solution of a cylinder indented into an elastic medium (semiinfinite half space) was compared to the solution from the material model by balancing the selfweight of the roller with the normal reaction from the material. Then, the mean contact pressure used for the elastic solution was calculated using the equation in figure 75.
Figure 75. Equation. Mean contact pressure.
In this equation, P is the load per unit thickness in the zdirection, and a is the radius of contact. Here, the selfweight per unit of pavement width is used for P. The normal stress in the elastic medium is given by the equation in figure 76.
Figure 76. Equation. Normal stress in elastic medium.
This representation of the solution indicates that the normal stress is expressed solely as a function of the position x and the mean contact pressure P_{m}. This relationship is illustrated in figure 77 in comparison to the instantaneous normal stress produced in the HMA material used in the compaction simulations with the same mean normal pressure applied. The comparatively higher stresses in the elastic medium at a given mean contact pressure are due to the lack of any dissipative mechanism in the elastic material. The model employed for representing the HMA develops dissipative responses in a nonlinear manner, and hence, the instantaneous response shows less stress. Also, the low level of stress in the HMA shows the lowviscosity fluidlike precompaction behavior of the model.
Figure 77. Chart. Comparison of the normal stress predicted by the model to the normal stress predicted for an elastic medium (for the same mean contact pressure, P_{m}).
The FE analysis is dependent on the use of interface element layers between the top asphalt layer and the previously compacted asphalt layer as well as between the top newly laid asphalt and the old asphalt sublayer. This is due to the differences in properties of the loose mix and the stiff old asphalt layer. Hence, the specialized interface elements are employed to resolve the contact issues that arise due to the large variation in properties of material on either side of the interface. With the establishment and implementation of a useful contactarea algorithm in CAPA3D, the analysis is performed to obtain insight into the effect of the interface elements. The use of interface elements provides a numerically stable and converging simulation. These interface elements do not contribute to the structural properties of the whole structure, and their applicability has been documented.^{(54)} The interface elements to be used require the specification of two shearstiffness and one normalstiffness moduli. These stiffness values are used in the numerical algorithm to ensure contact is maintained between the layers by artificially attaching the corresponding nodes of the layers through a mathematical expression that governs their relative motion. The shear stiffness is used to govern the relative sliding of the layers, whereas the normal stiffness is used to prevent mesh penetration of the finite elements of a stiffer material into those of a soft material.
The results of FE analysis using the interface layer are presented in figure 78 through figure 81. The two shearstiffness and one normalstiffness properties are varied, and the resulting structural behavior of the HMA material is observed. The shear moduli are individually made stiffer, and a comparison is drawn to a reference set of interfaceelement stiffnesses. The same comparison is made after normal stress is made less stiff. It can be observed in figure 78 through figure 81 that the use of interface elements does not adversely affect the normal and inplane displacement response of the model. The slight deviation observed when zdirection shear is stiffened is attributed to the structure of the mesh. The mesh used is longer in the zdirection, and hence, the effects of free boundaries are more pronounced. The situation is rectified by using a higher stiffness value because, beyond such a value, the interface exhibits similar behavior in the xdirection.
Figure 78. Chart. Comparison of pavement vertical response as the normal stiffness of the interface elements is lowered by an order of magnitude from 1,451 ksi (10,000 MPa) to 145 ksi (1,000 MPa).
Figure 79. Chart. Comparison of pavement response in the rolling direction as the normal stiffness of the interface elements is lowered by an order of magnitude from 1,450 ksi (10,000 MPa) to 145 ksi (1,000 MPa).
Figure 80. Chart. Comparison of pavement response to increasing the shear stiffness in the xdirection.
Figure 81. Chart. Comparison of pavement response to increasing the shear stiffness in the zdirection.
The impedance layers built into the FE mesh at this stage were not exploited for the current analyses since essentially nondynamic loading was studied. However, their impact on the FE model is recorded, for the sake of completeness, through some simple examples. These elements were included as part of the plan to study the effects of truly dynamic loads on the asphalt pavement structure. The use of impedance layers is necessitated by the need to dampen the wave reflections caused at the interface between layers when the pavement is subjected to dynamic loads. The techniques used to overcome the illposed nature of the initial value problem, which causes wave reflections at interfaces in highly viscous and elastic media, involve progressively dampening the displacement waves propagating through the medium through use of numerically cheap impedance layers. These layers are primarily composed of elastic material, with each enclosing layer proportionally less stiff in comparison to the inner one it is in contact with. This proportion, the "acoustic impedance," is maintained to keep an interface between different materials from reflecting waves that can interfere with the solution. Further details of the uses and applications of such special FE layers have been documented by Sluys.^{(55)}
In the mesh employed, three impedance layers are used to encase the pavement structure, each 1 inch (25.4 mm). The maximum value of the innermost layer is finally set at 725 ksi (5,000 MPa), and the densities are calculated such that the ratio of moduli between consecutive layers is kept at 10 (the crosssectional area is fixed by the thickness of the layers). These values eliminate the wave reflections in the asphalt mix and hence are used for further analysis for all impedance layers.
The response plots that have impedance layers around the top fresh mix, old asphalt, and base (as shown in figure 82 and figure 83) are useful in reducing the vertical and horizontal reflections to the point that they are completely eliminated, as indicated in figure 84 and figure 85. Because of this observation, a high enough impedance is sought when use of such layers is warranted, keeping the thickness of layers and the density of the material of the layers fixed.
Figure 82. Chart. Comparison of the dampening effect provided by impedance layers surrounding the structure laterally in the xdirection.
Figure 83. Chart. Comparison of the dampening effect provided by impedance layers surrounding the structure laterally in the ydirection.
Figure 84. Chart. Comparison of the horizontal dampening effect provided by impedance layers surrounding the top layer (pavement and old asphalt) laterally in the xdirection.
Figure 85. Chart. Comparison of the vertical dampening effect provided by impedance layers surrounding the top layer laterally in the xdirection.
FE analyses were performed by assigning properties for the base material in order to study the effect of having a single base layer underneath the pavement. The results indicate that the material response does not seem affected much when the material properties of the single base layer are varied over a wide range with a minimum base stiffness value of 72.5 ksi (500 MPa) and a maximum stiffness value of 290 ksi (2,000 MPa).
The plots of the volumetric response shown in figure 86 and figure 87 indicate that the material response is not very sensitive to changes in stiffness of the base when the response is viewed in an averaged sense. The volumetric response is examined through the evolution of the volumetric component of the viscousdeformation gradient (see Koneru et al. and Masad et al.).^{(48,49)} However, the FE model is responsive to a drastic increase in base modulus (as shown in figure 88 through figure 90) due to increasing rigidity in the base structure leading to interference waves (as a result of structural rigidity and not due to any other dynamic characteristic) reflected internally. This is in contrast to what is observed in the field, where mix designs with stiff base materials produce more compaction in the pavement layer because of the compressive reaction acting on the HMA layer from below. As a result, for further FE analysis, a value between the two extremes (i.e., between the minimum base stiffness value of 72.5 ksi (500 MPa) and the maximum stiffness value of 290 ksi (2,000 MPa), for example 145 ksi (1,000 MPa)), is chosen to conduct the parametric and sensitivity analysis presented later.
Figure 86. Chart. Comparison of the volumetric component of the viscousevolution gradient for soft (72.5 ksi (500 MPa)) versus stiff (290 ksi (2,000 MPa)) bases.
Figure 87. Chart. Comparison of the volumetric component of the viscousevolution gradient at two basestiffness moduli of interest.
Figure 88. Chart. Comparison of the effect on the xdisplacement of a node in the roller path as the base stiffness is varied from 72.5 ksi (500 MPa) (soft base) to 290 ksi (2,000 MPa) (stiff base).
Figure 89. Chart. Comparison of the effect on the ydisplacement of a node in the roller path as the base stiffness is varied from 72.5 ksi (500 MPa) (soft base) to 290 ksi (2,000 MPa) (stiff base).
Figure 90. Chart. Comparison of the deflection for two basestiffness moduli of interest.
Compaction Equipment Characteristics
As a further development to the simulation model, a vibratoryloading algorithm was developed for use in CAPA3D specifically for the compaction project. The vibratory loads take into account the effect of the displacement amplitude and frequency specifications of the machine used for rolling in a field project by calculating the equivalent load to be used to achieve the same level of compaction. Analytical expressions for calculating the contact forces for cylinder indentation into an elastic pavement material are adapted from Johnson.^{(57)} Thereafter, the equivalent loads are calculated for different amplitudes of indentation as necessary. Then, the force necessary for the vibrations is calculated using these analytical expressions at an instant of time, and the pressure is redistributed over the contact area.
To summarize, the algorithm has the following features:
Simulations were set up on a mesh representative of a 12ft (3.66m)wide pavement lane with a roller weight of 27,000 lb (12.258 kg), a diameter of 51 inches (1,295.4 mm), and a width of 7 ft (2.135 m). The results of the simulations (shown in figure 91 through figure 94) show that the material responds according to the mechanics of the loading algorithm implemented. The results are presented at a point in the interior of the material along the width (approximately 5 ft (1.525 m) from the edge). In the results that follow, the term compaction is used as a measure to represent the viscous evolution, det(G)1 (see Masad et al. for clarification of the term).^{(49)} The effect of a cyclic load vibrating at a certain angular frequency on the deformation induced on the material model of choice, with a certain relaxation time λ associated with it, is also presented (see figure 94).
The following observations can be made regarding the response of the material:
Figure 91. Chart. Compaction over a sequence of passes as the amplitude of vibration increases.
Figure 92. Chart. Compaction over a sequence of passes at different frequencies.
Figure 93. Chart. Material response to change in dead load carried by each roller.
Figure 94. Chart. Field compaction response at a constant frequency over multiple passes on a point.