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: FHWAHRT08073 Date: September 2009 
Chapter 5. MVEPCD Characterization and Verification in CompressionThe permanent deformation in HMA is affected by several mixture factors, such as the resistance of the binder to flow, aggregate angularity and gradation, the amount of asphalt, the air void content, etc. A significant amount of research has been conducted to develop laboratory test methods, analysis techniques, and models to study the permanent deformation growth of HMA. The nature of the permanent deformation models available in the literature ranges from empirical to mechanisticempirical to completely mechanistic. Further attempts have been made in recent years to develop a mechanistic permanent deformation model that involves fundamental material characterization.^{(20,52,53)} The objective of this research was to develop a constitutive model of HMA in compression that could predict HMA behavior under loading conditions and temperatures encountered in the field. The following sections present theories used in the modeling, the experimental program, the characterization process, and predictions for various loading histories. For this effort only, the ALF Control mixture has been used, and it should be understood that all figures in this chapter are related to this mixture. 5.1. Engineering Behavior of Asphalt Concrete in Compression5.1.1. Constant Rate CompressionFigure 85 and figure 86 show the ratedependent stressstrain curves of HMA at the confining pressures of 0 and 500 kPa, respectively. As observed from figure 85 and figure 86, the gross trends in stressstrain curves with 500 kPa confinement were quite similar to those with 0 kPa confining pressure; the overall strength of the material decreased with increasing temperature and with decreasing strain rates. Directly comparing the curves at different temperatures (figure 87 through figure 90), it was found that at 5 and 25 °C, the response was not sensitive to confining pressure. The strengths began to deviate during the 40 °C tests and were significantly different at 55 °C, with the confined tests showing higher strengths. This trend was similar to that observed in the confined dynamic modulus tests, which was described for tensioncompression loading in subsection 4.2.2 and for compression loading in section 5.2. This finding suggests that when the material was very stiff, the effects of confinement were negligible.
Figure 85. Graph. Stressstrain curves for unconfined constant strainrate tests.
Figure 86. Graph. Stressstrain curves for 500 kPa confinement constant strainrate tests. Figure 87. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 5 °C. Figure 88. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 25 °C. Figure 89. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 40 °C. Figure 90. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 55 °C. 5.1.2. Repetitive Creep and Recovery Test5.1.2.1. Creep and Recovery Test with VLFigure 91 through figure 93 present the viscoplastic strains for the given VL stress histories. The viscoplastic strains caused by the first loading in each group, except for the first loading group, were close to zero or were at least very small because the deviatoric stress that caused the viscoplastic strains was quite small compared to that of previous loadings. Figure 91. Graph. Viscoplastic strain versus cumulative loading time (unconfined VL). Figure 92. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VL). Figure 93. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement VL). 5.1.2.2. Creep and Recovery Tests with VT and RVTFigure 94 and figure 95 present the viscoplastic strains measured at the end of each rest period in unconfined VT and 500 kPa confinement VT, respectively. In both unconfined and confined VT testing, the overall reproducibility was sufficiently adequate to identify the effects of load level and loading sequence. In RVT testing, less viscoplastic strain was observed than in VT testing despite the fact that the test conditions were the same for both tests, with the exception of the sequence of loading, as described previously. This difference in viscoplastic strain indicated that the sequence of loading plays was an important role in viscoplastic strain development. In addition, a change of slope between the two groups of viscoplastic strain, which could not be explained by the concepts inherent of the conventional viscoplastic model, was observed. It seems the characteristic behavior of HMA was affected by viscoelastic relaxation. Although quantifying the variations in the viscoplastic strain rate under a given loading condition is necessary for rigorous modeling work, no test protocol was available that could capture only the viscoplastic strain rate because HMA showed timedependent viscoelastic strain, too. However, trends for viscoplastic strain rates developed in repetitive creep and recovery testing could be estimated by analyzing the VT test results. Figure 96 presents incremental viscoplastic strain rates (i.e., incremental viscoplastic strain divided by pulse time). As shown in figure 96, at 1percent viscoplastic strain, viscoplastic strain rates from 0.05s loadings (2.0 x 10^{3}) were much greater than those from 6.4s loadings (8.0 x 10^{5}). These results indicated that most of the viscoplastic strain developed at the beginning of the loading period and that the viscoplastic strain that developed during the remainder of the loading period (i.e., at 6.4 to 0.05 s) was relatively small. Because the calculated viscoplastic strain rate was the average of the viscoplastic strain rates during loading, the actual viscoplastic strain rates at the end of the loading were much smaller than 8.0 x 10^{5}. This is another important behavior of HMA, along with the softening concept presented in the following section. Mathematically, this behavior can be formulated in a viscoplastic constitutive model with either increasing viscosity or increasing yield stress due to aggregate interlocking. Figure 94. Graph. Viscoplastic strain versus cumulative loading time (unconfined VT testing). Figure 95. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement for VT and RVT testing). Figure 96. Graph. Incremental viscoplastic strain rate versus viscoplastic strain (500 kPa confinement, 1,600 kPa deviatoric). 5.1.2.3. Creep and Recovery Tests with Constant Load Level and CLTFigure 97 presents the viscoplastic strain history for each loading condition with a confining pressure of 500 kPa. A smaller viscoplastic strain was observed as the pulse time increased. Even considering the ramp time of 0.005 s, which was not taken into account in the cumulative loading time, the difference in the viscoplastic strains was quite significant. As shown in figure 97, more viscoplastic strain was observed in the CLT tests that consisted of shorter loading times at a given cumulative loading time. The differences were significant. For example, at 150 s of cumulative loading time, the viscoplastic strain in 1,800 kPa CLT testing with a 0.4s loading time was over 3 percent, whereas it was around 1.5 percent in 1,800 kPa CLT testing with a 6.4s loading time and was 2.2 percent with a 1.6s loading time. One reason for this behavior could be related to the dynamic effects associated with the ramp to the target load because one difference in the tests, as they are presented in figure 97, is the number of load applications. Alternatively, because the tests are exposed to different total rest times at given cumulative loading times, the differences could be related to material softening. Furthermore, it is not beyond reason to suppose that this softening behavior could be rate dependent. For further study of this issue, two additional VT tests with 0.1 and 0.05s rest periods were performed at a confining pressure of 140 kPa. These test results were compared with the results from VT testing with 200 s of rest. In these tests, the conditions were identical (i.e., the number of loadings and rest periods were the same for each), except for the length of the rest period. Figure 98 presents the viscoplastic strains measured at the end of the rest periods. For VT tests with 0.05 and 0.1 s of rest, 200 s of rest was allowed at the end of the testing to measure the pure viscoplastic strain because it was not possible to measure pure viscoplastic strain immediately after 0.05 or 0.1 s of rest. The deviatoric stresses were 827 and 552 kPa, and the confining pressure was 140 kPa. Figure 97. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement CLT).
Figure 98. Graph. Viscoplastic strain versus cumulative loading time As seen in figure 98, even though the loading histories were identical except for the length of the rest period, a smaller viscoplastic strain developed as the rest period became shorter. Furthermore, the effect of the rest period on viscoplastic development was significant considering the amount of viscoplastic strain when using 552 kPa deviatoric stress and 200 s of rest. This experimental observation did not demonstrate the effects of the dynamic loading ramp, but it did demonstrate the significant effect of ratedependent softening during unloading. Therefore, the modeling effort focused on developing a viscoplastic model that accounted for ratedependent hardening and/or softening. 5.1.3. tTS with Growing Damage in CompressionThe principle of tTS is one of the fundamental and most important concepts for HMA in tension modeling because it provides a strong mechanical background and significantly reduces the experimental efforts. In order to verify the principle for the compression stress state, stress characteristic curves were constructed from constant strainrate test results by utilizing shift factors determined from the dynamic modulus tests. Further verification was also performed using repetitive creep and recovery tests. For the verification of the principle for monotonic loading, a wide range of eight referencestrain values were chosen according to the results of both the uniaxial and triaxial compressive constant strainrate tests (figure 85 and figure 86). According to the procedure shown schematically in figure 55 through figure 57 and discussed in detail elsewhere, the stress and time values were determined for all of these tests at fixed strain levels.^{(13,20)} These plots of stress versus time are shown in figure 99 through figure 114. Then, shift factors obtained from smallstrain LVE testing were applied to determine the reduced time that corresponded to each physical time. If the tTS principle was valid with growing damage, the resulting plots of stress and reduced time would appear continuous at all strain levels. This behavior was indeed observed for the compression tests under both the confined and unconfined conditions (figure 115 and figure 116). These results verified that the tTS concept held true for mixtures subjected to compressive loading as well as to tensile loading, even if there was severe damage and viscoplastic strain. However, to verify that the principle held for the physical mechanisms behind the behavior of repetitive creep and recovery tests, more rigorous verification was neededthis verification compared VT test results at 40 and 55 °C with the same reduced time histories. For this verification, VT testing was first performed at 55 °C (the 200 s rest period results are used here). Next, the time history was used with the tTS shift factors to compute the equivalent reduced time history at 40 °C. However, because the testing time was estimated to take several days (the equivalent time to 200 s at 55 °C is approximately 3,265 s at 40 °C), the following analysis was performed to finish the VT testing within a reasonable time. In this analysis, the measured strain history during the unloading portion of several load pulses was used to compute the strain rate, which was plotted against the rest period time in figure 117. To avoid issues related to the initial loading of a test, the 0.05 s data were taken from the second load block, whereas the other pulse times were taken from the first loading block of the VT test. As shown in figure 117, most of the strain rates became quite small after around 40 s, except for the rest periods following the 1.6 and 6.4 s load pulses. For this reason, 40 s was used to compute the reduced time for pulse times less than 1.6 s (653 s at 40 °C); 50 s was used for a pulse time of 1.6 s (816 s at 40 °C); and 60 s was used for a pulse time of 6.4 s (980 s at 40 °C). Note that strain rates reached an asymptotic value of zero more quickly as the strain level increased, and thus, it was conservative to consider the times used in the first loading block as the reference times. The results of these two tests are shown figure 118. As the figure shows, viscoplastic strains measured at the end of rest periods were well matched to each other. This agreement confirmed that the tTS principle was applicable regardless of loading sequence and the amount of damage and viscoplastic strain in asphalt concrete. Figure 99. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.0001 strain level under uniaxial conditions.
Figure 100. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.0005 strain level under uniaxial conditions.
Figure 101. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.001 strain level under uniaxial conditions.
Figure 102. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.003 strain level under uniaxial conditions.
Figure 103. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.005 strain level under uniaxial conditions.
Figure 104. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.01 strain level under uniaxial conditions.
Figure 105. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.015 strain level under uniaxial conditions.
Figure 106. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.02 strain level under uniaxial conditions.
Figure 107. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.0001 strain level under 500 kPa conditions.
Figure 108. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.0005 strain level under 500 kPa conditions. Figure 109. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.001 strain level under 500 kPa conditions.
Figure 110. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.003 strain level under 500 kPa conditions.
Figure 111. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.005 strain level under 500 kPa conditions.
Figure 112. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.01 strain level under 500 kPa conditions.
Figure 113. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.015 strain level under 500 kPa conditions.
Figure 114. Graph. Stresstime curves for the Control mixture before the application of timetemperature shift factors at a 0.02 strain level under 500 kPa conditions. Figure 115. Graph. Stress mastercurves for the Control mixture under uniaxial conditions. Figure 116. Graph. Stress mastercurves for the Control mixture under triaxial conditions (500 kPa).
Figure 117. Graph. Variation of strain rate during unloading.
Figure 118. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VT at 40 and 55 °C). 5.2. MVECD Characterization in Compression5.2.1. Linear Viscoelastic CharacterizationAs presented in chapter 2, the first stage in MVECD characterization was the determination of the LVE properties of the material. Following the test protocols presented in chapter 3, frequencytemperature sweep tests were conducted in both unconfined and confined states. The results of this characterization under zeromaximum deviatoric stress conditions are shown for the unconfined stress state in table 16, for a confining pressure of 140 kPa in table 17 and for a confining pressure of 500 kPa in table 18. The shift factor function coefficients that resulted as part of the mastercurve construction process for each of the tests in these tables are shown in table 19. Table 16 through table 18 show that a trend similar to that observed for the zeromean deviatoric conditions was also observed for the zeromaximum deviatoric tests. Specifically, under confinement, the material modulus increased as a function of the confinement level. Additionally, under conditions where a given confinement state diverged from the unconfined state, the mastercurve moved toward a lower temperature and higher frequency at the higher confined stress. Additionally, a higher confining stress generally resulted in a lower phase angle (i.e., higher elasticity). Although it may be argued that increased aggregate interlocking was the cause of these behaviors, it was assumed that such an interpretation was flawed due to the strain magnitudes typically induced during the testing. The net effect of the increased dynamic modulus and phase angle was an increase in the relaxation modulus, which is shown for the different confining pressures in figure 119. These curves were obtained by the method outlined in subsection 2.2.1. Figure 119 shows that the relaxation modulus began to diverge at around 0.01 to 1 s in reduced time. These times generally corresponded to a physical time of between Table 16. Linear viscoelastic characterization and variation for the Control mixture in unconfined compression state at selected frequencies and temperatures.
Table 17. Linear viscoelastic characterization and variation for the Control mixture in 140 kPa confined compression state at selected frequencies and temperatures.
Table 18. Linear viscoelastic characterization and variation for the Control mixture in 500 kPa confined compression state at selected frequencies and temperatures.
Table 19. Effect of confining pressure on shift factor function coefficients for the Control mixture in compression state.
Figure 119. Graph. Confining stress effect on the relaxation modulus. 5.2.2. Comparison of ZeroMean and ZeroMaximum Deviatoric Stress ResultsFigure 120 through figure 123 present the LVE characteristics determined from the zeromean and zeromaximum deviatoric stress states. Only the 0 and 500 kPa results are shown in the figures because they were consistent for the two test methods. Overall, the comparisons were very favorable, and little difference was seen between the zeromean and zeromaximum deviatoric stress states. However, a noticeable difference was seen in the phase angle results, but these results were subject to higher variability. Further, the phase angle calculation was less precise due to limitations of the measurement instrumentation. These results indicated that at small strains, asphalt concrete was not bimodal, and it was only after some threshold strain had been exceeded that such behavior occurred. Such findings are consistent with the work presented elsewhere when care is taken in performing the experiments at sufficiently small strains (50 to 70 με).^{(54,55)} When such care is not taken, conflicting data are found in the literature.^{(28,56)} For modeling then, it is important that the LVE characterization is limited to the 50 to 70με range; in this range, either the zeromean or zero‑maximum method may be used. As a final comparison of the two test conditions, the model developed in subsection 4.2.2 was applied to the zeromaximum deviatoric stress state tests. The results of this analysis are shown in figure 124, and good agreement was observed between the model and measured data. Figure 120. Graph. Comparison of zeromean and zeromaximum deviatoric stress dynamic modulus mastercurves in semilogarithmic scale. Figure 121. Graph. Comparison of zeromean and zeromaximum deviatoric stress dynamic modulus mastercurves in logarithmic scale. Figure 122. Graph. Effect of test method on shift factor functions. Figure 123. Graph. Comparison of zeromean and zeromaximum deviatoric stress phase angle mastercurves. Figure 124. Graph. Application of stress statedependent model to zeromaximum deviatoric stress tests. 5.2.3. MVECD Damage Function CharacterizationCharacterization of the VECD model in compression was identical in analysis and general experimental requirements to that in tension. Constant rate compression tests were performed under both unconfined and confined conditions, the results of which were then used to calculate the damage functions C_{11}, C_{12}, and C_{22}. The process was discussed previously in subsection 4.2.3 and is shown schematically in the flow diagram shown in figure 125. 5.2.3.1. Characterization of C_{11}(S)Uniaxial constant rate tests were used to characterize the C_{11} damage function. Test results were used with equation 89 and equation 90 to compute the material function and damage parameter, respectively. Note that stresshardening (equation 85) was used for this analysis. As with the tension damage function, this relationship was refined following the NCHRP 919 methodology.^{(22)} This relationship is shown for the Control mixture in figure 126. This figure is constructed at a reference temperature of 5 °C and is only for the material under compressive loading. For this mixture, a_{1} = 7.679 x 10^{6}, a_{2} =1.6971 x 10^{6}, a_{3} = 7.6793 x 10^{6}, a_{4} = 4.6192 x 10^{7}, a_{5} = 1.5128 x 10^{7}, and a_{6} = 4.1017 x 10^{8}. Figure 125. Diagram. MVECD model characterization. Figure 126. Graph. C_{11} versus S for compression for Control mixture (5 °C reference). 5.2.3.2. Characterization of C_{12}(S)In addition to using uniaxial constant rate tests to characterize the C_{11} damage function, these test results were also used to characterize the C_{12} damage function. As with the tension tests, complications arose due to the timedependent nature of Poisson's ratio in asphalt concrete. Based on the analysis presented in subsection 4.2.3.2, an initial value of approximately 0.179 was assumed for C_{12}. This assumption was not unreasonable based on the results of the characterization tests. A slight difference between the initial values of the compression and tension data shown previously was observed due to the use of slightly different mixtures for the two sets of tests (Control2006 for tension and Control for compression). Tension data for the Control mixture will be shown in subsequent sections of this report. The C_{12} material parameter was defined from equation 91, and the damage parameter was calculated and refined as before. The functional form taken for the compressive C_{12} damage function is shown along with the data in figure 127. It was observed that at early damage stages, Poisson's ratio changed very little, but, in general, as damage increased, Poisson's ratio reduced, as indicated by an increase in C_{12}. For this mixture, it was found that b_{1} = 3.5598 x 10^{8}, 1.6982 x 10^{8}, and 3.2756 x10^{10}. Figure 127. Graph. C_{12} versus S for compression for Control mixture (5 °C reference). 5.2.3.3. Characterization of C_{22}(S)Confined constant rate tests, along with the other two damage functions, were used to characterize the final C_{22} function. Following the study conducted for the tension characterization methodology, a visual basic macro was created in Microsoft Excel^{®} to determine the damage function by optimization. The characterized C_{22} function is shown in figure 128 along with the model function. The parameters of the model were H_{1} = 0.504959 and H_{4} = 1.16247 x 10^{13}. Figure 128. Graph. C_{22} versus S for compression for Control mixture (5 °C reference). 5.2.4. Comparison of MVECD Behavior in Compression and TensionAsphalt concrete, like other composite materials, is known to exhibit bimodal behavior. In particular, the strength of the material it is in and its ductility are much higher in the compression direction than they are in the tensile direction. Comparisons of the MVECD damage functions revealed some interesting behaviors. First, when comparing the primary axial modulus damage function, C_{11}, between tension and compression, the material was substantially more resistant to damage in compression than it was tension. This behavior is shown in the same scale in figure 129. This behavior was explained by the notion that the damage parameter, S, represented some kind of crack density or volumeaveraged cracked area. In tension, these cracks were oriented perpendicular to the principal loading direction, whereas, in compression, the cracks ran parallel to the principal loading direction. In constant rate loading, cracks only grew when they were under a tensile stress. Because tension was only induced (at the microscale level) in compression tests and not directly applied, the material was more resistant to a given cracking volume. This physical interpretation for damage was also supported by the behavior of the second damage function, C_{12}. These damage functions are shown for both tension and compression in figure 130. This figure has both the C_{12} function and Poisson's ratio between the symmetric axis and perpendicular axis (_{ν3132}) for the convenience of the discussion. From the comparison shown in figure 130, it was observed that the primary Poisson's ratio showed an overall decreasing trend in tension, whereas the opposite trend occurred in compression. As with the C_{11} damage function, the tensile behavior of C_{12} was observed to be more sensitive to damage. A possible explanation of this behavior depended on understanding that Poisson's ratio indicated the degree to which the radius changed relative to a unit change in length. Applying the physical interpretation of damage in tension and compression, the physical implications of the patterns in both damage curves were identical. The tension curve was increasing, thus indicating that with higher levels of damage, the direction parallel to the primary cracking direction (radial strain) was less sensitive to changes in strain perpendicular to the primary cracking direction (vertical strain). For compression, the primary cracking direction was parallel to the loading direction and perpendicular to the radial direction, so with the same interpretation (of damage in tension and compression), an increase in Poisson's ratio was observed. In the most general terms, the increased opening of a microcrack within a body was not necessarily accompanied by a relative increase in the parallel to the crack direction dimension of the body. In fact, the results indicated that at higher levels of cracking or damage the material became less likely to change dimension as much when loaded perpendicular to the damage orientation. Figure 129. Graph. Comparison of tension and compression of C_{11} damage function. The physical interpretation of the C_{22} damage function was not as straightforward as the other two damage functions. As seen in equation 45, the C_{11} and C_{12} damage functions can be directly translated to elements in the stiffness matrix. However, the C_{22} damage function entered in the denominator of multiple elements of the material stiffness matrix. Nevertheless, this damage function can be considered as a sort of volumetric compliance, which has little meaning in a transverse isotropic problem. From figure 131, though, it is seen that the material was more volumetrically compliant in compression than it was in tension. Figure 130. Graph. Comparison of tension and compression of C_{12} damage function. Figure 131. Graph. Comparison of tension and compression of C_{22} damage function. 5.3. Viscoplastic Modeling of Asphalt Concrete in Compression5.3.1. A Phenomenological Model Considering Pulse Time EffectsAs a first step toward developing a mechanistic material model for the behavior of HMA in compression, a series of analyses was performed on VT and VL test data, and a phenomenological model was developed. The modeling approach adopted in the phenomenological model was based on the existing strainhardening model presented for describing the tensile behavior in section 4.3 and shown in the general form as follows: (158) Where: = Viscoplastic strain rate. σ = Stress. ε_{vp }= Viscoplastic strain. As shown in equation 158, the viscoplastic strain rate was represented by the combination of two functions, f(_{σ}) and g(_{εvp}), which allowed the stress rate dependency and strain hardening to be taken into consideration in the model. Equation 158 can be generalized as equation 159, which accounts for the effect of the pulse time. (159) Where: t_{p} = The loading time. The exact form of the function, F, is presented along with experimental data in the following sections. 5.3.1.1. Tests Performed in This StudyThree types of repetitive creep and recovery tests were performed for the phenomenological model development, including the creep and recovery test with VL, the creep and recovery tests with VT, and the creep and recovery test with a constant load level and CLT. All the tests were conducted at 55 °C under the confining pressure of 500 kPa. Experimental details for these tests are given in the following sections because this particular study uses different test conditions than those outlined in chapter 3. 5.3.1.1.1. Creep and Recovery Test with VL:The creep and recovery test with a VL test was performed with 200 kPa as the starting load level. An incremental factor of 1.0245 was used for the subsequent load levels to increase the load level until the complete failure of the specimen. The loading time and rest period were 0.1 and 10 s, respectively. 5.3.1.1.2. Creep and Recovery Tests with VT:The creep and recovery test with a VT test was performed with a loading block consisting of 30 different loading times. The loading time varied from 0.005 to 2.0 s with an incremental factor of 1.1356. The rest period for each load cycle was 30 times that of each loading time. The VT tests were performed at three different load levels, 1,800, 2,000, and 2,200 kPa. 5.3.1.1.3. Creep and Recovery Test with CLT:In this test, a constant load level and constant loading time were used for each test. Load levels and loading times were changed between tests. Three load levels of 1,800, 2,000, and 2,200 kPa were used with the loading time of 1.6 s. For the 2,000 kPa load level, the creep and recovery tests were conducted with three additional loading times of 0.1, 0.4, and 6.4 s. The VL and VT tests were used to understand the effects of load level and loading time on the permanent deformation behavior of HMA and to calibrate the phenomenological model. The CLT tests were used to verify the calibrated model. 5.3.1.2. Model CharacterizationBy observing the viscoplastic strain rate versus the viscoplastic strain VT in figure 132, a constitutive relationship between these viscoplastic media was defined, as shown in equation 160. (160) Where: a(t_{p}) = Material function of loading time. D(t_{p},_{ σ}) = Intercept of the curve on viscoplastic strain rate axis. Figure 132. Graph. Incremental viscoplastic strain rate versus viscoplastic strain In equation 160, a was a function of the loading time, and D was a function of the load level and loading time. Equation 160 can be represented in equation 161, which was a generalized form of equation 158. (161) Where: k(t_{p}) = Function of loading time. Equation 160 required the determination of a(t_{p}) and D(t_{p},_{ σ}) to calculate the viscoplastic strain rate for a given viscoplastic strain. Values of a(t_{p}) and D(t_{p},_{ σ}) for given loading times could be found by fitting log functions against each viscoplastic strain rate versus viscoplastic strain curve corresponding to given loading time. At this time, the values of a(t_{p}) can be represented by the second logarithm function, as shown in equation 162. (162) Where: a_{1} and a_{2} = Materialdependent constants. In order to determine the form of D(t_{p},_{ σ}), it was assumed that D(t_{p},_{ σ}) could be represented by the summation of the loading time term and the load level term, as shown in equation 163. (163) Where: The function c(σ) was given by c_{1}log(σ) when the stress was less than 1,000 kPa and by c_{2σ} when the stress was greater than 1,000 kPa. The coefficient b was determined by fitting equation 163 against D(t_{p},_{ σ}) from the VT test. At this time, it was assumed that c(σ) and d constitute one constant that accounts for the effect of load level. The function c(_{σ}) was then fitted by the logarithmic function and the linear function. Finally, coefficient d was determined by fitting equation 163 for the viscoplastic strain rate versus the viscoplastic strain curve of the VL and VT tests. Figure 133 and figure 134 show the fitting results and the coefficients that were determined. Figure 135 and figure 136 presents predictions for the VT and VL tests, which were then used for the characterization process. Figure 133. Graph. Determined fitting results and coefficients of function a(t_{p}). Figure 134. Graph. Determined fitting results and coefficients of function D(t_{p}_{ σ}).
Figure 135. Graph. VT predictions.
Figure 136. Graph. VL predictions. 5.3.1.3. Verification of the ModelAs shown in figure 137 through figure 142, predictions were made for the CLT tests. Although the model was able to account for differences in the viscoplastic development for various loading times, overall predictions were not as accurate. Several causes for the discrepancy between viscoplastic strain predictions and measurements could be suggested. However, the inability of the model to consider strain history was a highly probable cause for this discrepancy, given that the fitting results for the VT and VL tests were acceptable.
Figure 137. Graph. CLT predictions (2.0 MPa deviatoric stress0.1s pulse time).
Figure 138. Graph. CLT predictions (2.0 MPa deviatoric stress0.4s pulse time).
Figure 139. Graph. CLT predictions (2.0 MPa deviatoric stress1.6s pulse time).
Figure 140. Graph. CLT predictions (2.0 MPa deviatoric stress6.4s pulse time).
Figure 141. Graph. CLT predictions (1.8 MPa deviatoric stress1.6s pulse time).
Figure 142. Graph. CLT predictions (2.2 MPa deviatoric stress1.6s pulse time). 5.3.2. HISSPerzyna ModelThe HISSPerzyna model, suggested by the Delft University of Technology and the University of Maryland, was investigated with respect to the data set obtained from experimental tests.^{(26,28,36)} A prediction using the Delft University of Technology model was not be made because of numerical problems. However, the characterization process using the tTS principle and coefficients of the model are described in subsection 5.3.2.1. 5.3.2.1. Delft University of Technology ModelThe model suggested by researchers at the Delft University of Technology required the development of several relationships between the material parameters and the strain rates obtained from constant strainrate tests. In this research, the tTS principle was utilized to simplify the modeling effort and to reduce the number of relationships required. With the assumption that the yield stress in deviatoric stress space was presented as circular, (β=0), equation 72 was be reduced to equation 164. (164) Where: Φ = Softening parameter. α = Hardening parameter. R = Tensile strength of material when deviatoric stress is 0. n = Parameter that determines shape of yield stress. Figure 143 shows peak stresses for a series of compressive and tensile constant strainrate tests; the strain rates are listed in table 4. These peak stresses were used as fundamental quantities to develop relationships between the material parameters and the reduced strain rates. R and _{γ0} could be determined as functions of the reduced strain rate by plotting the compressive and tensile peak stresses obtained from the constant strainrate tests under a certain strain rate and then taking the slope and xintercept of the line. In the model, R and γ0 represented the tensile strength for hydrostatic stress and the softening of the material in the postpeak region, respectively. The parameter n governs the overall shape and size of the yield function and was related to the dilation of the material. The beginning of dilation was defined as the stress at the minimum plastic volumetric strain because the plastic deviatoric strain and elastic strain (or viscoelastic strain) was assumed to not be associated with the volumetric change of a material. In addition, because HMA specimens dilated after a little compression as the compressive stress increased, the dilation stress could be determined. Once the dilation stresses were determined for several strain rates, the value of n could be determined using equation 165. (165) Where: I_{1,diation }= I_{1} at beginning of dilation. J_{2,dilation }= J_{2} at beginning of dilation. Once n is determined, α_{0 }can be readily determined using equation 166. (166) The sigmoidal function was used to represent relationships between the reduced strain rate and the material parameters. The form of the function and the coefficients determined for each parameter are listed in equation 167 and table 20. Figure 144 through figure 147 show a comparison of measured values versus predicted values at various reduced strain rates. (167)
Figure 143. Graph. Compressive and tensile peak stress in SQRT(J_{2})  I_{1} space. Table 20. Delft material model coefficients functions.
Figure 148 shows the strain ratedependent yield surface that was constructed using the characterized parameters when the viscoplastic strain was equal to zero (i.e., the initial yield surface). It was observed that the initial yield surface increased as the temperature decreased; the reduced strain rate increased. This behavior coincided with observations from constant strainrate tests in which more viscoplastic strains were developed under a small, reduced strain rate (or higher temperature). As shown in equation 164, the second term in the square root always had to be smaller than the first term in order to construct a valid yield surface. However, because of the numerical errors involved in the characterization process of α and n, the prediction program was often required to calculate the square root of a negative number during analysis. This situation was encountered without tTS, as mentioned by others.^{(36)} Figure 144. Graph. Determined γ_{0} parameter function. Figure 145. Graph. Determined R parameter function. Figure 146. Graph. Determined n parameter function. Figure 147. Graph. Determined α_{0} parameter function. Figure 148. Graph. Ratedependent initial yield surface. 5.3.2.2. University of Maryland ModelAs shown in equation 73, a simplified HISSPerzyna model was suggested by researchers at the University of Maryland.^{(28)} Equation 168 represents a general form of the hardening function used for the suggested model. (168) Where: α_{0} and Κ = Material constants. However, the observation made in subsection 5.1.2.3 indicates that a single hardening function, equation 168, was not sufficient to represent the characteristic behavior of the material, such as softening during unloading. Therefore, one more variable, the viscoplastic strain increment during loading, was introduced into equation 168. Equation 169 represents the modified functionα to incorporate the variation of the viscoplastic strain rate during the pulse time in the existing model. α_{1} and α_{2} in equation 169 describe general variations of α in terms of viscoplastic strain and a local variation of α in terms of incremental viscoplastic strain in a pulse, respectively. (169) Where: (170) (171) Figure 149 presents the variation of α determined by using a modified alphaviscoplastic relationship for five 6.4s pulses with 1,800 kPa of load level. As shown, α was no longer a simple decreasing function of the viscoplastic strain, but had multiple decreasing functions of which independent variables were incremental viscoplastic strain during each load pulse and overall viscoplastic strain. The incremental viscoplastic strain was reset to zero each time the material was unloaded. Figure 150 through figure 152 present measured and predicted viscoplastic strains by using a modified hardening function. The model was able to describe viscoplastic strain development for various loading conditions, such as VT, RVT, and CLT; this capability was not possible in the existing HISSPerzyna model. However, even though incremental viscoplastic strain in a pulse described multiple hardening rates at certain viscoplastic strains, it was more reasonable to assume that the multiple hardening rates were caused by the viscoelastic property of the material, given the rate dependency of the softening. Therefore, a viscoplastic model with a ratedependent hardeningsoftening function was developed and is presented in subsection 5.3.3. Figure 149. Graph. Variation of α for 1,800 kPa CLT loading (500 kPa confinement). Figure 150. Graph. Viscoplastic strain predictions for VT tests (500 kPa confinement). Figure 151. Graph. Viscoplastic strain predictions for CLT tests (500 kPa confinement). Figure 152. Graph. Viscoplastic strain predictions for RVT tests (500 kPa confinement). 5.3.3. Development of a Viscoplastic Model Using RateDependent Yield StressThe model developed in this research was capable of capturing both additional hardening that was due to aggregate interlocking and ratedependent softening due to viscoelastic relaxation. Viscosity in Perzyna's evolution law is separated into a constant term and a viscoplastic straindependent term that together represent the change of viscosity in viscoplastic flow. The yield stress function that takes into account ratedependent hardening and softening is also described in subsection 5.3.3.1. 5.3.3.1. Flow Rule and Yield Function for Developed Viscoplastic ModelAs an expansion of equation 55, a general flow rule for materials exhibiting kinematic and isotropic hardening is represented in equation 172. m amplified or reduced the stress rate dependency of the model, and D determined the viscosity in the viscoplastic flow. When D was a constant, it was assumed that the effect of the change in viscosity on the response of the material was taken into account by the yield stress function. However, when variations of yield stress were affected by the viscoelastic property of the material, it seemed reasonable to consider the viscosity in the viscoplastic flow as a function that was not subjected to yield stress. (172) Where: α = Kinematic hardening function. γ = Isotropic hardening function. D = Viscosity parameter. m = Ratedependency parameter. Therefore, a flow rule that takes into consideration the variations of viscosity in the viscoplastic flow is suggested in equation 173 by incorporating Perzyna's flow rule and Von Mises' yield criterion. In equation 173, D was the viscosity and represented the scalar hardening and softening as described above. The anisotropic behavior of the material was also integrated by using D_{ij}. Meanwhile, G_{ij} represented the orientationdependent isotropic hardening function that reflected the viscoplastic and viscoelastic property of the material. Because the material was subjected only to compressive stress, the kinematic hardening rule was not introduced in this model. The viscosity (D) was related to aggregate interlocking and was represented as a function of the viscoplastic strain, as shown in equation 174. Because the function could represent both increasing and decreasing viscosity according to the viscoplastic strain, it had the potential to represent the behavior of HMA mixtures in the tertiary region as well as in the primary and secondary regions. (173) Where: D = Viscosity related to aggregate interlocking, equation 174. η_{0}, m, α, β, γ = Material constants. D_{0 }= Initial viscosity. s_{ij }= Deviatoric stress tensor. g_{ij }= Deviatoric back stress tensor. σ_{ij} = Stress tensor. G_{ij} = Yield stress tensor. J(_{σij} G_{ij}) = The second invariant of (σ_{ij}  G_{ij}). (174) However, as discussed in subsection 5.1.2.3 , ratedependent softening, which implied the possibility of a multiple state of the material at certain viscoplastic strains, was observed when HMA was subjected to repetitive loading. In order to introduce the characteristic behavior of HMA into the viscoplastic constitutive model, equation 175, which was one of the simplest forms, was suggested as the hardeningsoftening function. (175) In equation 175, the yield stress increased as the viscoplastic strain and viscoplastic strain rate increased during loading, whereas it decreased during unloading (when the viscoplastic strain rate was zero). Figure 153 presents a schematic concept of the variation of yield stress subjected to a creep and recovery loading condition. The remaining yield stress, which was yield stress at the asymptote, was governed by the viscoplastic strain at the end of the loading condition, _{ε0}^{vp}, and material constants E_{1} and E_{2}. The decreasing yield stress during unloading allowed a multiple viscoplastic strain rate at a certain viscoplastic strain. Figure 153. Illustration. Variation of yield stress (Standard Linear Solid model). Because of their simplicity, equation 176 was suggested as a hardeningsoftening function to confirm the characteristics of the model for an arbitrary stress history, and equation 177 was suggested to predict the actual behavior of the HMA mixture. As shown in equation 177, the hardeningsoftening function was represented as the convolution integral, including the relaxation modulus. Material constants A and B were introduced into the relaxation modulus to develop a relationship between the relaxation modulus, which was the LVE property, and the viscoplastic yield stress. Additionally, by utilizing the material constants, the number of parameters needed to calibrate the model could be reduced. In calculating the yield stress, the state variable approach was used to reduce computational time, as shown in equation 178. (176) Where: E_{1}, E_{2}, _{η}_{1}, _{η}_{2 }= Material constants. (177) Where: A, B = Material constants. E_{0}, E_{i }= Prony coefficients determined for the relaxation modulus. ρi = Relaxation time. Equation 177 was solved using the state variable approach to predict strains and calculate pseudo strains. This approach is shown mathematically in equation 178. (178) Where: (179) (180) 5.3.3.2. Characteristics of the Developed Model for Arbitrary Stress HistoryIn order to confirm the characteristics of the developed viscoplastic model, the following predictions were made for the arbitrary stress histories by using equation 173 and equation 176. In this study, D was considered a constant to simplify the calibration and prediction processes. Two sets of stress histories were generated, as shown in figure 154 and figure 157. Table 21 shows the material constants used in this analysis. Table 21. Material coefficients used for the developed model analysis.
5.3.3.2.1. Effect of Rest Period:Figure 154 shows two different stress histories used to check the sensitivity of the viscoplastic model to the effects of rest periods. For both stress histories, stress levels and the cumulative loading time were fixed to 2,000 units less stress and 160 s, respectively. However, for the first stress history, 8.0 s of rest between the loading pulses were allowed, whereas only 1 s of rest was allowed for the second loading history. Figure 155 presents the variation of yield stress for each stress history, and as expected, the model showed different yield stress developments depending upon the rest period. Figure 156 presents the viscoplastic strain developed by each stress history. It shows more viscoplastic strain for a longer rest period. This result corresponded to the experimental observations made from figure 98. Figure 154. Graph. Stress histories for rest period analysis. Figure 155. Graph. Yield stress versus cumulative loading time (rest period analysis). Figure 156. Graph. Viscoplastic strain versus cumulative loading time (rest period analysis). 5.3.3.2.2. Effect of Loading Time:Figure 157 presents another set of stress histories used to check the effects of loading time. For these stress histories, the load level, rest periods, and cumulative loading time were fixed to 2,000 units, 4 and 66 s, respectively. However, the first loading history consisted of 6 pulses at 11 s long, and the second loading history consists of 22 pulses at 3 s long. The analysis results for the given stress histories are shown in figure 158. The loading history with shorter individual loading times showed more viscoplastic strain, which was identical to the CLT test results. As shown, the viscoplastic model that incorporated the softening rule appeared to account for the pulse time effect. As shown in figure 156 and figure 158, the viscoplastic model with the ratedependent hardeningsoftening capability could account for the effects of loading time. Figure 157. Graph. Stress history for loading time analysis. Figure 158. Graph. Viscoplastic strain versus cumulative loading time (loading time analysis). 5.4. Characterization and Verification of the Viscoplastic Model5.4.1. CalibrationPrior to the calibration process, data points acquired during the unloading period were filtered to reduce the computational time. Strains measured at the end of the rest periods were defined as the objective function. The nonlinear optimization function (lsqnonlin) in Matlab^{TM} was utilized to minimize errors between the measured and predicted viscoplastic strains. Based on the model calibrated for the VT and VL tests, the viscoplastic strains of the other loading conditions, such as CLT and VLT, could be predicted. Table 22 shows the coefficients determined from the calibration process for 140 and 500 kPa confining pressures. Table 22. Compression viscoplastic material model coefficients.
Figure 159 and figure 160 present the calibration results for VT and VL testing at 140 kPa confining pressure, and figure 161 and figure 162 present the calibration results for VT and VL testing at 500 kPa confining pressure. In general, the predicted and measured viscoplastic strains match very well, although there was a slight discrepancy in the VT and VL 500 kPa confining pressure results. Figure 159. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VT). Figure 160. Graph. Viscoplastic strain versus cumulative loading time Figure 161. Graph. Viscoplastic strain versus cumulative loading time Figure 162. Graph. Viscoplastic strain versus cumulative loading time 5.4.2. VerificationViscoplastic strain predictions for the 140kPa confining pressure tests are presented in figure 163 through figure 166. Figure 163 shows the ability of the developed model to consider the effects of rest periods on the viscoplastic strain development even though the viscoplastic strains were slightly underpredicted. Figure 164 and figure 165 show the predictions for the VLT tests and a low load level VL test, respectively; these predictions were quite good. Figure 166 presents predictions for complex loading histories, which were a combination of VT test results and flow number test results. Up to 0.5percent strain, the prediction of the viscoplastic strain matched well with the measured viscoplastic strain; however, the viscoplastic strain was underpredicted for the last strain level. This discrepancy could indicate a need to refine the softening function. Figure 167 to figure 172 present the viscoplastic strain predictions made for 500 kPa confining pressure tests. The overall prediction was good considering the complexity of the loading history. Figure 163. Graph. Viscoplastic strain versus cumulative loading time Figure 164. Graph. Viscoplastic strain versus cumulative loading time Figure 165. Graph. Viscoplastic strain versus cumulative loading time Figure 166. Graph. Viscoplastic strain versus cumulative loading time Figure 167. Graph. Viscoplastic strain versus cumulative loading time Figure 168. Graph. Viscoplastic strain versus cumulative loading time Figure 169. Graph. Viscoplastic strain versus cumulative loading time Figure 170. Graph. Viscoplastic strain versus cumulative loading time Figure 171. Graph. Viscoplastic strain versus cumulative loading time Figure 172. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement VLT). 
