Top-Down Cracking of Hot-Mix Asphalt Layers: Models for Initiation and Propagation (2010)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

12 3. FINDINGS: THE VECD-BASED MODEL The primary feature of the VECD-based crack initiation model was to account for effects of damage zones on response prior to cracking. For this reason, several important material property sub-models, including aging, healing, failure criteria, viscoplasticity, and thermal stress models, were developed, modified, and/or investigated, and then incorporated into the existing VECD model. These material sub-models were then converted into and/or combined with the structural sub-models. The integrated model was implemented into the VECD-FEP++, and an extrapolation method was developed for predicting top-down cracking initiation in HMA pavements. 3.1 Framework The overall framework guiding the VECD-FEP++ analysis is shown in Figure 3-1. The analysis is divided into five sub-modules: the input module, the material properties sub-models, the analytical sub-models, the performance prediction module, and the output module. These modules provide an analytical/computational method for identifying the location and time of crack initiation in the pavement structure. Each sub-model requires specific parameters to be obtained and/or determined. Table 3-1 summarizes the inputs for the VECD-based model along with the quantity of interest and test method for each input (a detailed description of the laboratory specimen fabrication and test method is given in Appendix A).

13 Figure 3-1. VECD-based model framework 3.1.1 Inputs Module Preprocessors have been developed to facilitate easy and rapid analysis of pavement systems using the FEP++. Specifically, the preprocessor helps in the rapid development of input models for analysis and also helps in making consistent changes for repeated analysis. This tool is ANSI-compliant and developed with portable libraries, thus making it easy to transfer to other platforms. Inputs Module CHAPTER 1Preprocessor Material Property Sub-models CHAPTER 2Linear Viscoelastic (LVE) CHAPTER 3VECD CHAPTER 4Aging CHAPTER 5Failure Criteria CHAPTER 6Healing CHAPTER 7Viscoplasticity CHAPTER 8Thermal (TSRST) Analytical Sub-models CHAPTER 9EICM CHAPTER 10Structural Aging CHAPTER 11Damage Factor Performance Prediction Module (VECD-FEP++) CHAPTER 12Traffic CHAPTER 13Extrapolation Outputs Module CHAPTER 14Damage Contour CHAPTER 15Crack Initiation: time,

14 Table 3-1. Inputs Required for the VECD-based Model Input Quantity of Interest Test Method LVE Material coefficients ( E∞ , iE , iρ ) Temperature/frequency sweep test VECD Material coefficients (a, b) Monotonic direct tension test at 5°C, Constant crosshead cyclic test at 19°C Aging LVE and VECD with time Same as above and long-term oven aging Failure Criteria Critical pseudo stiffness Constant crosshead cyclic test Healing Change in pseudo stiffness with rest duration Random loading and healing period cyclic tests Viscoplasticity Material coefficients (p, q, Y) Monotonic direct tension test at 40°C Thermal Thermal expansion coefficient Thermal expansion test, or from literature EICM Temperature distribution within a pavement structure -- Structural Aging Effective age of pavement materials with depth Binder viscosity or |G*| (A and VTS) Damage Factor Viscoplastic sensitivity -- Traffic Number of Equivalent Single Axle Loads (ESAL count) -- Extrapolation Number of analysis segments -- The preprocessor provides • Graphical forms for user input; • Pavement loading with support for multiple moving loads— a. Material data with support for elastic, viscoelastic, and viscoelastic continuum damage (VECD) models, and b. Analysis parameters and finite element (FE) mesh customization; • Pavement temperature profiles generated from the EICM; • Data validation and consistency checks; • Visualization of the finite element model (FEM), consisting of a. A graphical view of the model and mesh; b. Standard graphical operations, such as zoom, pan, rotate, etc.; and c. Mesh information, such as node and cell numbers; • Data persistence using a database for storing and retrieving input data; and • One-click installer and user documentation.

15 3.1.2 Material Property Sub-models A brief review of the concepts of the LVE model and the VECD model and their most important formulations are given in Appendixes A.1.2.1 and A.1.2.2, respectively; a review of both the viscoplastic and thermal stress models is given in Appendix A.1.2.5 and A.1.2.6. However, because the aging, healing, and failure criteria sub-models were developed in this project, the details of each are given below. 3.1.2.1 Aging model The effect of aging in asphalt binders on the performance of asphalt aggregate mixtures is well recognized. Traditionally, the effects of binder aging on mixture performance have been investigated using two approaches. The first is to subject different asphalt binders to various aging conditions and measure the resultant changes in physical properties to assess the aging potential of the different binders. However, this approach does not account for the effect of aggregate particles and does not yield accurate and realistic information about the performance of the mixtures. The other approach is to subject the asphalt aggregate mixtures to various aging conditions and then measure the physical properties of the aged mixtures. Although this method is more realistic because the aging process is conducted directly on the mixtures, little work has been done to determine the best method for incorporating these changes in mixture properties due to aging into the framework of constitutive modeling or fatigue performance modeling (15). The original VECD uniaxial constitutive model is based on the assumption that the material is a non-aging system. This study focuses on establishing an analytical/experimental methodology for incorporating the effects of aging in the current constitutive model. 3.1.2.1.1 Aging method for laboratory mixtures NCHRP Project 9-23 has found that the Strategic Highway Research Program (SHRP) protocol is not sufficient to simulate field-aging behavior in the laboratory because of its inability

16 to account for variables such as field-aging conditions and mix properties (16). Nevertheless, the SHRP method 1. Is simple to implement; 2. Provides a general relationship between laboratory and field-aging behavior; and 3. Has been used successfully for previous studies conducted at NCSU (15). In addition, because this part of the study focuses on developing a place-holder relationship for the aging effects in the current VECD model rather than on matching the laboratory aging to that of the field, it was decided to use the SHRP aging method in this study. Four levels of asphalt mixture aging were simulated, as follows (17): (1) Short-term aging (STA): The loose, uncompacted mixture is conditioned at 135°C for 4 hours and then compacted. Specimens are cored and cut for testing. (2) Long-term aging, Level 1 (LTA1): The aging procedure is the same as for STA, except the specimens are conditioned at 85°C for 2 days before testing after coring and cutting. (3) Long-term aging, Level 2 (LTA2): The aging procedure is the same as for STA, except the specimens are conditioned at 85°C for 4 days before testing after coring and cutting. (4) Long-term aging, Level 3 (LTA3): The aging procedure is the same as for STA, except the specimens are conditioned at 85°C for 8 days before testing after coring and cutting. These aging processes for asphalt mixtures follow the AASHTO R30 specification (18), with the exception that three different long-term aging times are used. To minimize slump in the specimens during the oven aging procedure, the method suggested by the NCHRP 9-23 project, whereby the specimens are wrapped in wire mesh that is held in place by three steel clamps, was adopted. Diameter, height, weight, and air void percentage of all specimens were measured before and after oven aging. The differences in the dimensions were negligible; the dry weights were reduced by 0 to 0.5 g, and the air voids (%) increased by 0 to 0.5%, depending on the level of aging. The maximum specific gravity (Gmm) value may also change with aging; however, a

17 single Gmm value for the STA mixture was used for all the aged specimens. Therefore, it was concluded that no apparent damage occurred due to aging. 3.1.2.1.2 Materials and specimen preparation for laboratory testing Two different mixtures were used for the aging study: A mix and AL mix. The A mix, one of the most common mixture types used in North Carolina, uses granite aggregate that was obtained from the Martin-Marietta quarry in Garner, NC. The aggregate structure is a fine 9.5 mm nominal maximum aggregate size (NMAS) mixture composed of 36% #78M stone, 25% dry screenings, 38% washed screenings, and 1% baghouse fines. The blended gradation is shown in Figure 3-2. An unmodified asphalt binder (PG 64-22) from Citgo in Wilmington, NC, was also used for this study. The asphalt content for the A mix is 5.7% by weight of the total mix. After a series of laboratory tests, it was determined that this mix was insensitive to aging and thus did not meet the study goals. A weak trend was found between the VECD characterization parameters and the aging times. As a result, the AL mix was used because it was known to be prone to aging. The AL mix consisted of AAD-1 asphalt binder and local limestone aggregate. The binder (PG 58-28) was obtained from the Materials Reference Library (MRL), and the local limestone aggregate was obtained from the Martin-Marietta Castle Hayne quarry. The aggregate structure is a fine 9.5 mm Superpave mixture composed of 50% # 78M stone, 33% 2S-sand, 14.5% dry screenings, and 2.5% baghouse fines. The blended gradation is shown in Figure 3-2. A Superpave mix design procedure was performed, and an optimal asphalt content of 6.2 % by weight of total mix was chosen. The mixing and compaction temperatures used for the study were 147°C and 135°C, respectively.

18 0 20 40 60 80 100 Sieve Size (mm) % P as si ng A mix AL mix Control Points 0. 07 5 12 .59. 5 4. 75 2. 36 0. 60 0 0. 15 0 1. 18 0. 30 0 Figure 3-2. Mixture gradation 3.1.2.1.3 VECD characterization of aged mixtures The A mix The material characteristics of primary importance in VECD modeling are the LVE and damage characteristics. Figure 3-3 (a) and (b) present the replicate averaged dynamic modulus (|E*|) master curves for four different aging levels of the A mix in both semi-log and log-log scales. In general, the stiffness increases with aging time, but the increases are small. This finding suggests that the short-time aging kinetics are more active than the long-time aging kinetics. Figure 3-3 also shows that the STA mixture exhibits similar stiffness characteristics at higher reduced frequencies (physically representing cooler temperatures and/or faster loading frequencies) than the LTA mixtures, but deviates from the LTA mixtures at approximately 1 Hz. Below this reduced frequency, the STA mixture shows lower modulus values than the LTA mixtures. Such a trend is expected because as temperature increases or frequency reduces the modulus gradient between the asphalt binder and aggregate particles becomes larger. As a result,

19 the effect of the asphalt binder properties on the total mixture behavior may become more noticeable. Similar trends can also be seen in the phase angle master curve graph, shown in Figure 3-3 (c); the three LTA mixtures exhibit similar phase angle values for all frequencies, but the STA mixture shows a different phase angle in the lower frequency range. With regard to the damage characteristics, Figure 3-3 (d) indicates that the LTA3 mixture exhibits the most favorable damage characteristic curve, followed by the LTA2 and LTA1 mixtures and, finally, by the STA mixture. The differences among the aged mixtures are not significant and cannot be used directly to assess fatigue performance because the curves shown Figure 3-3 (d) do not account for stiffness effects, boundary conditions, or failure characteristics. The importance of these characteristics, as well as a method to account for them, is explained in the next section. The AL mix Figure 3-4 presents the averaged dynamic modulus mastercurves for the four different levels of aging of the AL mix in both semi-log and log-log scales. Compared to the A mix, the effect of aging on the dynamic modulus mastercurves is more significant over all of the frequency ranges. Moreover, a clearer trend with regard to aging time is observed with the AL mix when compared to the A mix. A similar trend can be seen in the phase angle mastercurves graph, shown in Figure 3-5. To evaluate the statistical differences among the mastercurves for the various aging levels, all replicates results have been analyzed by the step-down bootstrap method (details of the analysis and results are presented in Appendix A.1.2.4). The results reveal that test data for the four aging levels are overall statistically different, except under the most extreme reduced frequencies.

20 0 5000 10000 15000 20000 25000 30000 1E-08 1E-05 1E-02 1E+01 1E+04 Reduced Frequency (Hz) |E *| (M Pa ) STA LTA1 LTA2 LTA3 (a) 100 1000 10000 100000 1E-08 1E-05 1E-02 1E+01 1E+04 Reduced Frequency (Hz) |E *| (M Pa ) STA LTA1 LTA2 LTA3 (b) 0 5 10 15 20 25 30 35 40 45 1E-08 1E-05 1E-02 1E+01 1E+04 Reduced Frequency (Hz) |E *| (M Pa ) STA LTA1 LTA2 LTA3 0.0 0.2 0.4 0.6 0.8 1.0 0.E+00 1.E+05 2.E+05 3.E+05 S C STA LTA1 LTA2 LTA3 (a) (b) Ph as e An gl e (d eg ) (d)(c) Figure 3-3. VECD characterization summary of the A mix: (a) dynamic modulus mastercurves in semi-log space, (b) dynamic modulus mastercurves in log-log space, (c) phase angle mastercurves, and (d) damage characteristic curves With regard to the damage characteristics, the trend is the same as for the A mix where the damage characteristic curve for the LTA3 mixture is positioned the highest of all, followed by the LTA2 and LTA1 mixtures and, finally, by the STA mixture, as shown in Figure 3-6. However, the differences among the aged mixtures are more significant than those of the A mix. Care should be taken in concluding from these observations that aging affects the mixture properties relative to performance. The aging effects must be quantified by considering both resistance to deformation (stiffness) and resistance to damage. Hence, the aging effects on pavement performance will be discussed in the next section along with the results of the FEP++ simulations.

21 0 5000 10000 15000 20000 25000 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 Reduced Frequency (Hz) |E *| (M Pa ) AL-STA AL-LTA1 AL-LTA2 AL-LTA3 (a) 100 1000 10000 100000 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 Reduced Frequency (Hz) |E *| (M Pa ) AL-STA AL-LTA1 AL-LTA2 AL-LTA3 (b) Figure 3-4. Dynamic modulus mastercurves for mixtures of the AL mix in: (a) semi-log space, and (b) log-log space

22 0 5 10 15 20 25 30 35 40 45 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 Reduced Frequency (Hz) Ph as e A ng le (d eg ) AL-STA AL-LTA1 AL-LTA2 AL-LTA3 Figure 3-5. Phase angle mastercurves for mixtures of the AL mix 0.0 0.2 0.4 0.6 0.8 1.0 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 S C AL-STA AL-LTA1 AL-LTA2 AL-LTA3 Figure 3-6. Damage characteristic curves for mixtures of the AL mix 3.1.2.1.4 Incorporation of aging effects into the VECD model To incorporate the effects of aging into the VECD model, two major parameters must be evaluated: the viscoelastic properties and damage characteristics. It is analytically possible to

23 include the effects of aging in the formulation of the current constitutive model using another time variable that accounts for the aging time, as shown in Equations 3-1 and 3-2. The major advantage of this approach is that the interaction between loading and aging can be modeled realistically, thus allowing a more accurate evaluation of the effects of aged binders on mixture properties and the performance of the mixtures. Accounting for the effects of aging on the damage growth of asphalt aggregate mixtures using fundamental damage mechanics principles is probably the best approach to realistically simulate the interaction between damage and aging that occurs simultaneously in actual pavements. The pseudo strain (εR) of an aged material under uniaxial conditions is represented by the hereditary integral (Equation 3-1). ( ) 0 1 ,R R dE t d E d ξ εε ξ ξ ξ = ∫ (3-1) where εR is a particular reference modulus included for dimensional compatibility and typically taken as one. Two time variables are used in this integral: the aging time, t , and the loading time, ξ (i.e., the time that has elapsed since the specimen was fabricated and ε is uniaxial strain). From this equation, the relaxation modulus (E) accounts for the aging effects; it is a function of aging time, loading time, and temperature, and is given as ( , , ) ( , )E E t t T E tτ ξ= − = (3-1) where ( ) / Tt aξ τ= − , and Ta is the time-temperature (t-T) shift factor. Theoretically, Equation 3-2 results in the mastercurves of different aging times at a specific reference temperature. However, the relaxation mastercurves of the different aged mixtures could not simply be shifted horizontally or vertically to construct a single mastercurve that includes the

24 effects of aging. Moreover, solving Equation 3-1 with the two-dimensional relaxation modulus (Equation 3-2) is analytically difficult. Thus, an alternative, quasi-static approach was followed that takes advantage of the relatively long time scale of aging as compared to the time of loading. As a result of an analytical investigation, it was found that the coefficients for the dynamic modulus mastercurve and damage characteristic curve, expressed by Equations 3-3 through 3-6, vary as either a power law or exponential model with aging time. To develop a final set of relationships between these coefficients and in-service aging time, the lab aging-to-field aging times suggested by SHRP (17) were selected. This relationship is as follows: 1 year for the STA mix, 4 years for the LTA1 mix, 7.5 years for the LTA2 mix, and 18 years for the LTA3 mix. The relationships between aging time (te) and the LVE and VECD model coefficients are expressed by Equations 3-7 through 3-10. These equations show that the models have been formulated by using the ratio of aged values to original (un-aged) values so that the final function could be applied universally to other mixtures to simulate the aging effects. The model considers ways that the LVE and damage characteristics of a material vary with age. As expected, the model suggests that AC becomes stiffer and less time-dependent with age. ( )log log * 1 1 Rd g f bE a e + = + + (3-2) where |E*| is dynamic modulus, and fR is reduced frequency. ( ) 21 2 3log Ta T Tα α α= + + (3-3) where aT is time-temperature shift factor, and T is temperature. 1 1 u α = + (3-5) where u is the material constant related to material time dependence.

25 nmSC e= (3-4) where C is normalized pseudo secant modulus, and S is damage parameter. 0 0 0 0 1 2 3 4, , ,e e e et t t ta a a ae e e e t t t t a b d g t t t t a b d g = = = = (3-5) where a1, a2, a3, and a4 are regression constants. 0 0 0 12 22 32 1 2 3 , , 1 2 3 e e et t t e e e t t t t t tα α α α α α α α α = = = (3-6) where α12, α22, and α32 are regression constants. 0 max 2 31(1 ( ) ) et k kt e k t α α α = + (3-7) where αmax, k1, k2, and k3 are regression constants. 0 0 ( 3 4* ) ( 3 4* ) 1 1 , 2 2 1 1 1 1 e e e e t t t t m m t n n t m n m nm nm n e e+ + = = + + + + (3-8) where m1, m2, m3, m4, n1, n2, n3, and n4 are regression constants. The relationships between aging time and the LVE and VECD model coefficients are shown in Figure 3-7 through Figure 3-9. Generally, the sigmoidal coefficients (a, b, d, g), shift factor coefficients (α1, α2, α3), and damage evolution rate (α) vary as a function of aging time in power relationships. The damage coefficients (m, n) change with aging time by following the sigmoidal type of function shown in Figure 3-9. The best fit relationship between the experimental results and the sigmoidal function was selected to represent the VECD model and, as a result, the relationship appears to be somewhat arbitrary. It cannot be stated definitively that rapid transition between years 7 and 12, as suggested by the model, will occur in reality.

26 However, a more thorough investigation into the correlation between the oven aging and field aging times for damage characteristics is beyond the scope of the current research. 0.996 1.000 1.004 1.008 1.012 1.016 1.020 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S ig m oi da l F un ct io n (a ) (a) 0.997 0.998 0.999 1.000 1.001 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S ig m oi da l F un ct io n (b ) (b) 1.00 1.05 1.10 1.15 1.20 1.25 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S ig m oi da l F un ct io n (d ) (c) 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S ig m oi da l F un ct io n (g ) (d) Figure 3-7. Relationship of sigmoidal functions to aging time

27 0.996 0.997 0.998 0.999 1.000 1.001 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S hi ft Fa ct or (a 1) (a) 0.99 1.00 1.01 1.02 1.03 1.04 1.05 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S hi ft Fa ct or (a 2) (b) 0.00 0.20 0.40 0.60 0.80 1.00 1.20 0 5 10 15 20 Aging Time (yrs) N or m al iz ed S hi ft Fa ct or (a 3) (c) 0.98 1.00 1.02 1.04 1.06 1.08 1.10 0 5 10 15 20 Aging Time (yrs) N or m al iz ed a lp ha (d) Figure 3-8. Relationship of shift factor coefficients and alpha to aging time 0.00 0.20 0.40 0.60 0.80 1.00 1.20 0 5 10 15 20 Aging Time (yrs) N or m al iz ed D am ag e C oe ffi ci en ts (m ) (a) 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12 1.14 0 5 10 15 20 Aging Time (yrs) N or m al iz ed D am ag e C oe ffi ci en ts (n ) (b) Figure 3-9. Relationship of damage coefficients to aging time 3.1.2.2 Failure criteria as functions of aging and temperature Aging effects on the performance of asphalt mixtures must be evaluated in terms of the combined effects of stiffness, damage, and failure criteria. A series of cyclic fatigue tests has been conducted to investigate the effects of aging on these characteristics directly and closely. One important finding from this experimental work is that the failure of the specimen occurs at different degrees of damage, depending on the aging level of the mixture and the testing

28 temperature. This subsection compares the fatigue performance for mixtures aged at different levels, then presents a new method to analyze cyclic fatigue tests, and finally, presents the failure criteria used to interpret the FEP++ results. Detailed information on this new method for cyclic test analysis, the so-called Simplified VECD Model (S-VECD), is provided in Appendix A.1.2.2.4 along with material model fatigue life predictions that include a more general failure criterion to determine the capability and accuracy of the S-VECD (see Appendix A.1.2.3). 3.1.2.2.1 Cyclic fatigue tests Cyclic fatigue tests using the AL mix were conducted for all four aged mixtures under controlled stress (CS) and controlled crosshead (CX) cyclic conditions. Two amplitudes, high and low, were applied at two temperatures, 5°C and 19°C. The test conditions and some important results, especially the number of cycles to failure (Nf), are summarized in Table 3-2. To define the number of cycles to failure, the phase angle criterion suggested by Reese (19) is utilized because it seems to work well under both CX and CS test conditions. The phase angle quickly drops in the CX tests and suddenly increases in the CS tests once the specimen fails. One example from the test results for the AL-STA mixture is shown in Figure 3-10. 0 1000 2000 3000 4000 0.E+00 1.E+03 2.E+03 3.E+03 4.E+03 5.E+03 Number of Cycles |E *| 20 30 40 50 Ph as e A gn le |E*| Phase Angle (a) 0 2000 4000 6000 8000 10000 12000 14000 0.E+00 1.E+03 2.E+03 3.E+03 4.E+03 5.E+03 Number of Cycles |E *| 10 20 30 Ph as e A gn le |E*| Phase Angle (b) Figure 3-10. Dynamic modulus and phase angle versus load cycle for AL-STA mix: (a) the CX test (AL229), and (b) the CS test (AL248) Overall, the dynamic modulus decreases and the phase angle increases according to the number of cycles until the specimen fails, which is when the localization of on-specimen strain starts.

29 Using the cyclic test results, comparisons can be made with regard to four conditions: (1) magnitude of the input, (2) testing mode, (3) temperature, and (4) aging level. However, care must be taken in comparing these results, because the conditions to be compared are not always the same (e.g., there may be a difference in on-specimen strain for each aged mixture). Nonetheless, the comparisons led to the following conclusions: 1. Regardless of testing mode, temperature, and aging level, fatigue resistance decreases as the magnitude of the input increases; 2. By comparing the CX test results at a similar initial strain magnitude, as the temperature decreases or the aging time increases, the resulting initial stress magnitude increases and the number of cycles to failure (Nf) decreases; and 3. By comparing the CX test results at a similar initial stress magnitude, as the temperature decreases or the aging time increases, the resulting initial strain magnitude decreases and the number of cycles to failure (Nf) increases. Overall, the different testing modes result in the opposite fatigue performance (e.g., the STA mix shows better performance in the CX tests, but the LTA3 mix shows better performance in the CS tests). It is known, based on energy principles, that stiff materials tend to perform better in CS test protocols, whereas soft materials yield better performance in CX tests, all other factors being equal. In this case, however, all the factors are not equal because aging has occurred, and has somewhat embrittled the LTA mixes. The experimental data suggest that this embrittlement may not be as significant for the CS tests. As noted above, soft materials should perform better in the CX test protocols, and because the softest mixture is also the least embrittled, i.e., the STA, the CX test results cannot be used reliably to determine the exact effects of this embrittlement process.

30 Table 3-2. Cyclic Fatigue Test Summary (Frequency of 10 Hz) Mix ID Specimen ID Test Designation Initial Stress Amplitude in Tension (kPa) Initial Peak-to-Peak Strain (Microstrain) Nf AL-STA AL229 19a-CXb-Hc 1,700 582 3,091 AL232 19-CX-H 1,400 415 8,908 AL231 19-CX-L 610 162 > 150,000 AL246 19-CS-H 750 152 1,610 AL238 19-CS-L 250 50 46,200 AL233 5-CX-H 2,870 274 3,498 AL239 5-CX-L 980 87 > 280,000 AL248 5-CS-H 1,500 120 4,000 AL247 5-CS-L 900 76 44,700 AL-LTA1 AL228 19-CX-H 1,790 466 3,613 AL242 5-CX-H 2,900 259 3,100 AL244 5-CX-L 920 74 > 280,000 AL243 5-CS-H 1,500 107 4,590 AL245 5-CS-L 900 66 77,000 AL-LTA2 AL225 19-CX-H 2,020 483 2,753 AL250 5-CX-H 3,210 264 713 AL241 5-CX-L 1,180 84 > 280,000 AL249 5-CS-H 1,500 99 4,200 AL240 5-CS-L 900 64 97,200 AL-LTA3 AL220 19-CX-H 2,030 419 1,383 AL234 5-CX-H 3,220 253 800 AL235 5-CX-L 1,810 134 > 280,000 AL237 5-CS-H 1,500 95 8,900 AL236 5-CS-L 900 52 167,000 a Test temperature in degrees Celsius; b Test mode; c Relative magnitude: H=high, L=low Damage characteristics from cyclic tests Work by Underwood et al. (20) concluded that the damage characteristic curves obtained from CS, CX, and monotonic testing collapse into a single curve, indicating that the property represented by these curves is fundamental and independent of temperature and test type. Initially, the cyclic data for the aged mixtures were analyzed following the method by Underwood et al. (20) (details are provided in Appendix A.1.2.2.4), but the damage characteristic curves obtained from the different modes of testing did not collapse well, although the collapse in each test mode was acceptable. This analysis showed that the samples conditioned for the LTA mixes produced poorer results than those conditioned only for the STA mix. Thus,

31 some trials were conducted, and the resultant damage characteristic curves collapsed reasonably using two options taken from the originally suggested analysis method: (1) the  function is defined as (1/u + 1) for the monotonic test and (1/u) for the CX and CS tests, and (2) the correction factor, I , is set to one for all the cyclic tests. These differences in the current and earlier work for analyzing cyclic test results still need to be investigated. Because the definition of  in the original model is based on theoretical arguments and has been verified using STA materials only, its universality with regard to materials tested at other aging conditions was not verified in the original work (20). Figure 3-11 shows the damage characteristic curves of the AL- STA mixture obtained from the monotonic, CX, and CS tests. All the cyclic test data are included in this plot and collapse well with each other. The damage characteristic curves for all four aged mixtures are shown in Figure 3-12 for only the tests that have been performed under the same testing condition (i.e., 19-CX-H). The collapses are good except for the discrepancy between the monotonic and cyclic results of the LTA1 mixture. 3.1.2.2.2 Failure criterion As shown in Appendix A.1.2.3.1, the general failure criterion developed in previous work is based on reduced frequency and does not include the aging effect on the failure. Because the failure criterion does not include aging time and because the reduced frequency basis approach is not amenable for use with the VECD-FEP++, a temperature-only based failure criterion has been developed that includes the aging effect in the failure definition.

32 0.0 0.2 0.4 0.6 0.8 1.0 0E+00 1E+05 2E+05 3E+05 4E+05 5E+05 S (Method D) C * 19-CX-H 19-CX-H 19-CX-L 19-CS-H 19-CS-L 5-CX-H 5-CX-L 5-CS-H 5-CS-L Monotonic (α=1/u+1) Monotonic (α=1/u) STA Figure 3-11. Damage characteristic curves of AL-STA mix (all cyclic test data included) 0.0 0.2 0.4 0.6 0.8 1.0 0E+00 1E+05 2E+05 3E+05 4E+05 5E+05 S (Method D) C * 19-CX-H 5-CX-H 5-CX-L 5-CS-H 5-CS-L Monotonic STA (a) 0.0 0.2 0.4 0.6 0.8 1.0 0E+00 1E+05 2E+05 3E+05 4E+05 5E+05 S (Method D) C * 19-CX-H 5-CX-H 5-CX-L 5-CS-H 5-CS-L Monotonic LTA1 (b) 0.0 0.2 0.4 0.6 0.8 1.0 0E+00 1E+05 2E+05 3E+05 4E+05 5E+05 S (Method D) C * 19-CX-H 5-CX-H 5-CX-L 5-CS-H 5-CS-L Monotonic LTA2 (c) 0.0 0.2 0.4 0.6 0.8 1.0 0E+00 1E+05 2E+05 3E+05 4E+05 5E+05 S (Method D) C * 19-CX-H 5-CX-H 5-CX-L 5-CS-H 5-CS-L Monotonic LTA3 (d) Figure 3-12. Damage characteristic curves of: (a) STA, (b) LTA1, (c) LTA2, and (d) LTA3 aged mixtures of AL mix

33 In Figure 3-13, the pseudo stiffness (C) values at failure are plotted against the test temperatures for the CX tests for all aged mixtures. The C values at failure decrease as the temperature increases for all aging levels. The aging effect on failure is most noticeable at 19°C. Such a trend can be expected because of the aforementioned tendency of AC to embrittle with age. On the other hand, all the aging mixes failed at similar C values at 5°C. This result is due mainly to the fact that AC is brittle enough at low temperatures that the aging effect is not significant. Because of limited test data, it is assumed that the failure criterion varies linearly between 5°C and 19°C and is constant beyond this range. All data at 5°C are averaged, and that the LTA1 and LTA2 data are averaged because they show similar values, as shown in Figure 3-14. Therefore, the piecewise function is used for the failure criterion, as shown in Equation 3-11. Like the aging model, the failure criterion is formulated by the ratio of the aged value to the original (un-aged) value so that it can be universally applied to other mixtures according to the aging effect. 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 25 Temperature (°C) C * a t f ai lu re STA LTA1 LTA2 LTA3

34 Figure 3-13. Variation of C at failure from cyclic fatigue tests 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 Temperature (°C) C * a t f ai lu re STA (C* = -0.0266xTemp + 0.7276) LTA1&2 (C* = -0.0205xTemp + 0.6972) LTA3 (C* = -0.0084xTemp + 0.6367) (a) Figure 3-14. (a) Failure criterion λ = -0.0396xTime + 1.0224 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 5 10 15 20 Aging Time (yrs) C o ef fic ie n t ( λ) (b) κ = -0.0072xTime + 1.0041 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 0 5 10 15 20 Aging Time (yrs) C o ef fic ie n t ( κ) (c) Figure 3-14. (b) and (c) variation of coefficients λ and κ as functions of temperature and aging, respectively 0 0 0 0 0 0 ( )*5 ( ) 5 ( )* ( ) 5 19 ( )*19 ( ) 19 e e e e e e t t t t t t f t t t t t t T C C T C T C T C λ κ λ κ λ κ λ κ λ κ λ κ  + ≤ °   = + ° < < °    + ≥ °  (3-9)

35 where 0 et t λ λ = * e1 t 2λ λ+ , 0 et t κ κ = * e1 t 2κ κ+ , and , , ,1 2 1 2λ λ κ κ = constant. 3.1.2.3 Healing model As a placeholder for a fully mechanistic, potentially micromechanically motivated healing model, a simplified version of the phenomenological mechanical model, which was derived, characterized, and verified by Lee and Kim (21), has been applied. A simplified version of this model is used here because (1) the previous formulation, of which the healing model is a part, has certain shortcomings in its rigor, and (2) modifying the material model for the FEP++ would require significant time. The primary disadvantage of this model is the lack of characterization and verification data under various conditions needed for the FEP++ simulations. For this reason, an empirical adjustment has been made to improve the engineering reasonableness of the simplified model. 3.1.2.3.1 Previous formulation The healing model formulation is encapsulated in the general framework of the VECD model, as it is coded in the FEP++. However, the healing model includes two additional damage parameters (S2 and S3) and pseudo stiffnesses (C2 and C3) that the current model does not include. The VECD model, shown in Equation 3-12, includes only a single damage parameter, S1, and material integrity term, C1. The additional terms are included mathematically in the healing function shown in Equation 3-13 where they are shown to physically represent the

36 increase in pseudo stiffness during the rest periods, C2(S2), and the reduction in pseudo stiffness as the healed material is redamaged, C3(S3). These material functions are independent of the amount of damage in the virgin material, i.e., C1(S1). ( )1 1 RC Sσ ε= (3-10) [ ]1 RC Hσ ε= + (3-11) where ( ) ( ) ( ) 1 , 2 2, 3 3, 1 1 , , 1 i R R R B i i i B j C j j H S C S C S C S S S − =  = + − − −  ∑ , if both the healed and virgin materials are damaged (see Region 1 in Figure 3-15); or 1 , , 1 i R R B j C j j H S S − = = −∑ , if only the virgin material is damaged (see Region II in Figure 3-15). Figure 3-15 schematically illustrates the model behavior in terms of pseudo stiffness for a single rest and reloading period. The evolution laws for the S2 and S3 damage parameters are given in Equations 3-14 and 3-15, respectively. The relationships between these damage parameters and their respective material integrity parameters are given by Equations 3-16 and 3- 17. The coefficients for these equations, as characterized by Lee and Kim (21), are summarized in Table 3-3.

37 Ps eu do S tif fn es s (k Pa ) No. of Cycles B B’ C D’ D Region II Region I ∆Nf Ps eu do S tif fn es s (k Pa ) ASRA Ps eu do S tif fn es s (k Pa ) Ps eu do S tif fn es s (k Pa ) Figure 3-15. Effect of rest period on pseudo stiffness ( ) 2 2 2 2 20 1 2 rpt R CS dt S α ε  ∂ =  ∂  ∫ (3-12) ( ) 3 2 3 3 32 R RA CSS S α ε  ∂ = − ∂   (3-13) ( ) ( ) 222 2 20 21 2 CC S C C S= + (3-14) ( ) ( ) 323 3 30 31 3 CC S C C S= − (3-15)

38 Table 3-3. Healing Model Formulation Coefficients (21) Coefficient AAM Value AAD Value Function C1 C10 0.9900 1.1000 C11 0.0065 0.0764 C12 0.4400 0.2000 α1 2.6950 2.7860 Function C2 C20 -0.3750 0.0015 C21 0.0941 0.0004 C22 0.2000 0.6200 α2 1.7000 1.5000 Function C3 C30 1.0100 1.0100 C31 0.0019 0.0114 C32 0.5500 0.3000 α3 2.1000 2.0000 t-T Shift Factor α1 0.0009 0.0012 α2 -0.1602 -0.1843 α3 3.3709 3.9471 TR 25 25 3.1.2.3.2 Simplified formulation The basic trends found in the above formulation and numerous experimental studies (15, 22, 23) are shown in Figure 3-16. In short, the more time that elapses between load applications, the more the pseudo stiffness will recover, and hence, the longer the fatigue life will be extended. The amount of recovery is most sensitive to the rest period duration, but the model also suggests that recovery is a function of the previous load intensity and the total cumulative amount of damage that exists in the material. The previous formulation includes these factors by including multiple parameters and damage functions. For the simplified model, an attempt is made to use only a single material function, C1(S1), to accomplish the same objective. A closer examination of the previous formulation shows that it is not possible to directly quantify the increase in

39 pseudo stiffness healing with a single function because, as shown in Figure 3-15, the healed material is more sensitive to loading than the virgin material (i.e., the damage rate of the healed material is greater than that of the virgin material in Region I, while it is the same for the healed and virgin materials in Region II). C C ∆ restξ S1, ∆Si S2, ∆Si S2, ∆Sj Figure 3-16. Effect of rest period on material healing The basic concept of the simplified model is diagrammed schematically in Figure 3-17 where damage characteristic curves are shown for the first three cycles of loading. The points evenly labeled (0, 2, and 4) represent the pseudo stiffness values at the beginning of load pulses 1, 2 and 3. The points labeled with the odd numbers (1, 3, and 5) represent the pseudo stiffness values at the end of load pulses 1, 2 and 3. The single damage parameter model will predict that during the first cycle, damage will grow from point 0 to point 1. The simplified model would then suggest that during some rest period the material would heal back (i.e., recover) to point 2. Then, when loading begins for the second cycle, the damage would grow to point 3. Upon resting after the end of the second cycle, healing would increase the pseudo stiffness to point 4. Finally, after the third loading cycle, the damage would grow to point 5. In total, then, there are three different segments where damage is assumed to occur in the virgin material: from point 0

40 to point 1; from point 1 to point 3 (second cycle); and finally from point 3 to point 5 (third cycle). Similarly, there are only two segments where damage is assumed to occur in the healed material: from point 2 to point 1 (second cycle) and from point 4 to point 3 (third cycle). Between each of these points, though, a single damage function, C1(S1), is used to compute the damage. Ps eu do S tif fn es s (k Pa ) Damage, S 0 4 2 1 3 5 Lo ad in g Time 0 1 3 52 4 Ps eu do S tif fn es s (k Pa ) Ps eu do S tif fn es s (k Pa ) Lo ad in g Ps eu do S tif fn es s (k Pa ) Lo ad in g Ps eu do S tif fn es s (k Pa ) Ps eu do S tif fn es s (k Pa ) Lo ad in g Figure 3-17. Conceptual schematic diagram for the simplified model In order to use the previous formulation, a systematic factorial data set was created for the different rest periods (x 7), damage levels (x 9), and energy inputs, ∆S (x 7), which have values that could be expected in any given FEP++ simulation. For each of these combinations, the previous formulation was used to predict the pseudo stiffness at the end of the cycle after the rest period, i.e., point 3. Then, SH in Equation 3-18 was determined through optimization to give the same pseudo stiffness value after a load pulse; e.g., damage was found at point 2, so when calculating damage growth, the resulting pseudo stiffness was equal to that at point 3. The pseudo stiffness at this healed damage level was then calculated using Equation 3-19, and finally used to compute the value of the change in pseudo stiffness (∆C), which is the difference

41 between the pseudo stiffness values at point 1 and point 2. In order to apply the outcomes universally, the damage level, S, and energy input, ∆S, were normalized to the damage level at failure, Sf, and the change in pseudo stiffness was normalized to the pseudo stiffness level. ( ) 1 2 1 1 2 R i H CS S S α ε+ ∂ = + − ∂  (3-16) ( ) ( ) 1210 11 C H H HC S C C S= − (3-17) After calculating the change in pseudo stiffness for the multiple conditions, a regression model was characterized for easy use. The basic model is shown in Equation 3-20 where the coefficients are found to depend on the damage level and energy input, as seen in Equations 3-21 through 3-29. Note that the model is presented in terms of reduced rest time, ξrest, which is defined in Equation 3-30 by combining the physical rest period, trest, and the t-T shift factor function, which is a second order polynomial, at some reference temperature, TR (25°C). 1 Heal rest C C δγ γ κ β ξ ∆ =     +      (3-18) 2 1 f S S κ κ κ  ∆ =      (3-19) * 1 f Sb Saeκ      = (3-20) 2 d f Sc Sκ  = −     (3-21) 2 1 f S S β β β  ∆ =      (3-22)

42 * * 1 h f f S Sg iS Sfeβ    +       = (3-23) * 2 f Sk Sjeβ      = − (3-24) 1γ γ= (3-25) * 1 f Sz Syeγ      = (3-28) 2 log 1 f Sp q fS o Sn S e δ δ  +       ∆ = +      +  (3-26) ( ) ( )2 21 210 R R rest rest rest T T T T T t t a α α ξ − + − = = (3-27) The coefficients for this model are summarized in for AAD and AAM mixtures (21). These two mixtures represent, respectively, a light healing and a heavy healing material. Because the model is normalized, i.e., it predicts the ratio of change in pseudo stiffness to current pseudo stiffness, ∆C/C, it can be applied universally to other mixtures to simulate heavy healing and light healing. The strength of the characterized model is shown for the AAD mixture in Figure 3-18.

43 Table 3-4. Simplified Healing Model Coefficients Coefficient AAM Value AAD Value Function Κ a 0.1679 0.1657 b 2.2695 2.8865 c 0.0548 0.0261 d 0.2572 0.2379 Function β f 27.3913 1.5144 g -1.5887 -4.7659 h 30.1058 5.0763 i -5.8228 -3.3953 j 0.5070 0.5290 k -8.1700 -3.0063 Function γ y 1.4707 1.2448 z -1.6387 -1.6425 Function δ n -0.3593 1.0526 o 30.1626 7.2289 p 2.8853 1.7666 q -0.2850 -2.0453 δ2 -0.0843 -0.1597

44 0.E+00 1.E+00 2.E+00 3.E+00 4.E+00 0.E+00 1.E+00 2.E+00 3.E+00 4.E+00 Observed dC/C Pr ed ic te d dC /C (a) 1.E-09 1.E-06 1.E-03 1.E+00 1.E-09 1.E-06 1.E-03 1.E+00 Observed dC/C Pr ed ic te d dC /C (b) R2 = 0.962 Se/Sy = 0.195 R2 = 0.970 Se/Sy = 0.202 Figure 3-18. Strength of simplified model for AAD mixture in: (a) arithmetic and (b) logarithmic scales 3.1.2.3.3 Healing potential factors One shortcoming of the simplified model, which has transferred from the previous formulation, is its handling of rest periods introduced before the healed material has been entirely redamaged (e.g., the rest periods in Region I). This issue is represented schematically in

45 Figure 3-19, but was not encountered during the development of the more robust healing model because the earlier work (21) did not include experiments in which the subsequent healing and redamaging occurred in Region I of the previous healing cycle. This situation is important in the context of FEP++ simulations because each load pulse is followed by a rest period and only a single load pulse is applied after each rest period. It is believed that the behavior suggested by the previous formulation, shown in Figure 3-19, is unrealistic, because the repeated redamaging of the healed material should actually limit the healing potential; thus, the response depicted in Figure 3-20 is deemed more appropriate. P se ud o S tif fn es s (k P a) No. of Cycles R BS R undamagedS A A A P se ud o S tif fn es s (k P a) P se ud o S tif fn es s (k P a) Figure 3-19. Effect of damage and healing in Region I using previous formulation

46 P se ud o S tif fn es s (k P a) No. of Cycles R BS R undamagedS A <A P se ud o S tif fn es s (k P a) Figure 3-20. Representation of expected material behavior Because not all damage occurs in the healed material, a correction factor was formulated to be a function of the total summed damage that occurs in the material. This correction factor is referred to as healing potential factor-C. Damage related to the redamage of healed material is considered to be negative damage, and the damage related to virgin material is considered to be positive damage. Referring to Figure 3-17, the damage that occurs from point 0 to point 1, point 1 to point 3 (second cycle), and finally point 3 to point 5 (third cycle) is calculated as a positive value. Then, the damage that occurs from point 2 to point 1 (second cycle) and from point 4 to point 3 (third cycle) is calculated as a negative value. If there is a net positive value, then the full healing potential suggested by Equation (3-18) is achieved. If, however, there is a net negative value, then the full healing potential is not reached, and the healing potential factor-C is less than one. The exact mathematical function used is shown in Equation 3-32. Another correction factor, healing potential factor-O, is applied with the consideration that the healing potential generally reduces as the material becomes more damaged. In other words, even when the healing in a cycle occurs at newly damaged sites, i.e., when the total net damage

47 is positive, the healing potential should nonetheless decrease as the damage grows. This situation is believed to be related to the fact that high damage levels represent both an increase in crack density and an increase in the average size of the cracks. As the size of the cracks increases, the overall potential for healing should be adjusted. The mathematical formulation for this healing potential factor is shown in Equation 3-33. Because the value of healing potential factor-O is related to the maximum damage level ( S ) it always reflects an overall drop in healing potential in damaged materials, regardless of the type of damage, new or rehealed. Figure 3-21 shows the healing potential factors, and Figure 3-22 represents each healing potential factor that will be used in the FEP++ simulation. The healing potential factor-C and healing potential factor-O are more important, respectively, in a single cycle and in overall periods within the simulation; therefore, both factors should be considered together, as shown in Equation 3-34. Healing Potential Factor 0 Total Net Damage (S_total) Damage Growth 1 Healing Potential Factor-C Healing Potential Factor-O (Positive)(Negative) Figure 3-21. Effect of net damage growth on healing potential correction factor

48 1E-02 1E-01 1E+00 1.E-01 1.E+01 1.E+03 1.E+05 Total "S" Healed * -1 H ea l P ot en tia l F ac to r - C (a) 1E-16 1E-12 1E-08 1E-04 1E+00 0.0 0.1 1.0 (Smax/Sf) H ea l P ot en tia l F ac to r - O (b) Figure 3-22. Healing potential factors *Adjusted HealC C Healing Potential Factor∆ = ∆ (3-28) ( )y1+z1*log -1* 1 0 1 0 1+e t t tS S Healing Potential Factor C S ≥ − =  < (3-29) y2+z2*( / ) 1 1+e max fS S Healing Potential Factor O− = (3-30) ( ) y2+z2*( / ) y2+z2*( / )y1+z1*log -1* 1 0 1+e 1 1 * 0 1+e1+e max f max ft tS S tS SS S Healing Potential Factor S  ≥=   <  (3-31) where maxS = maximum S in history, fS = S at failure, tS = total net damage, and 1, 2, 1, 2y y z z = constant. 3.1.3 Analytical Sub-models To predict pavement performance, the material sub-models must be converted to, or implemented into, structural models to consider the different structures, boundary conditions, climate conditions, etc. Consequently, three analytical sub-models have been developed: (1) a

49 structural aging model, (2) a damage correction factor (DCF) model, and (3) a temperature variation model. 3.1.3.1 Structural aging model Details of the development of a material level aging model are provided in Section 3.2.2.1. The time scale used for the material level aging model corresponds physically only to that used for the top layer of a real pavement cross-section. To apply this model to different depths, the age of each sub-layer relative to that of the surface must be found. This goal is achieved by coupling the principles of the Global Aging System (GAS), first proposed by Mirza and Witczak (24), with an effective time concept. A 10-year-old pavement can serve as an example of the effective time concept whereby after 10 years of service, the surface layer has aged 10 years, but the material at a depth of 3 inches may behave as the surface layer behaved when the pavement was only 4 years old. In this example then, the effective time of the sub-layer 3 inches from the surface 10 years after construction is 4 years. To compute the material properties of this sub- layer at year ten, material aging models described in Section 3.2.2.1.4 (Figure 3-7 through Figure 3-9) can be used to find the value of the coefficients at 4 years. The GAS model was used to determine the effective time of a given pavement structure. The GAS model predicts the viscosity of the asphalt binder as a function of depth, mean annual air temperature (MAAT) representing the effect of geographical location, and rolling thin film oven (RTFO) binder viscosity. The effective time is determined by finding, for some physical time and depth, the time that gives the same viscosity for the binder at the pavement surface. A flow chart of the structural aging model, including the equiviscosity concept as well as the plot of effective time versus depth for a typical simulation, is shown in Figure 3-23.

50 A & VTS, MAAT Viscosity (η) 1 yr D e p t h 20 yrs η(t1,d1) η(t2,d2) η(t3,d3) || || Equi-Viscosity GAS 0, (4 ) (1 4 ) 4(1 ) aged t aged z A A z Az          η(t, d) = η(te, d0) |E*| (ti, di) = |E*| (te, d0) = ratio function (|E*|STA) alpha (ti, di) = alpha (te, d0) = ratio function (alphaSTA) C(S) (ti, di) = C(S) (te, d0) = ratio function (C(S)STA) te = f (t, d) VECD change along aging 0 2 4 6 8 10 12 14 16 0 5 10 15 20 Effective Time (te) D ep th (c m ) Actual Time Figure 3-23. Structural aging model

51 As the GAS model is a function of the MAAT, the effective time varies for different climatic regions. To represent the wide range of conditions encountered in the United States, five regions are included in the study. These regions were selected based on the AASHTO climatic classification (25) and consist of (1) wet, freeze-thaw cycling (Washington, D.C.); (2) wet, hard freeze, spring thaw (Chicago, Illinois); (3) wet, no freeze (Tallahassee, Florida); (4) dry, freeze- thaw cycling (Dallas, Texas); and (5) dry, hard freeze, spring thaw (Laramie, Wyoming). Figure 3-24 is the plot of effective time at 20 years versus pavement depth for the five regions. The temperature in the legend is the MAAT for each region. The higher the MAAT, the higher the effective time at the same depth, which indicates that the aging effect is more significant for a higher MAAT region. 0 2 4 6 8 10 12 14 0 5 10 15 20 Effective Time (te) D ep th (c m ) DC (12.2°C) FL (19.4°C) IL (9.4°C) TX (18.3°C) WY (7.2°C) Figure 3-24. Effective time contours for various regions 3.1.3.2 Damage correction factor The current framework uses the VECD material model, which does not account for viscoplasticity and results in an overprediction of stress and increase in damage.

52 The viscoplasticity of a given material and structure is considered by applying the damage correction factor (DCF). The DCF is calculated by combining the simple strain-hardening viscoplastic material model (discussed in Appendix A.1.2.5) with the LVE stress and strain responses of the pavement under consideration. These responses are predicted by using the FEP++ software and turning off the damage functionality. The linear strain responses are adjusted to account for viscoplasticity, and then these adjusted strain responses are used to find an effective viscoelastoplastic modulus. A second run of the FEP++ was performed using the effective modulus values. Comparisons made between the responses of the first and second FEP++ runs revealed the sensitivity of the given pavement structure and materials to variations in viscoplastic strain. This sensitivity is reflected by a damage factor, which varies from 0 to 1. A value of 0 reflects a condition where the material is either very sensitive to structural viscoplastic effects and/or materially sensitive; conversely, a value of 1 reflects a condition where both structural and material viscoplastic effects are not very important (e.g., during winter). The idea behind this scheme is that a pavement simulation is neither in controlled stress nor controlled strain mode but somewhere in between. In a stress-controlled situation, the strain computed in the simulation corresponds to the viscoelastic strain and fully contributes to damage. However, in a strain-controlled situation, the viscoelastic strain that contributes to damage is significantly weaker than the strain computed in the simulation. Thus, depending on which end of the spectrum the analysis lies, the amount of strain that contributes to damage varies significantly. To account for this unpredictability, the DCF has been introduced and is defined as ( )1 dmg ve fem DCF ε β β ε   = − +     , (3-32) where

53 femε = the strain computed by an FEP++ simulation, dmg veε = the strain contributing to damage, and β = a factor that is 0 for control strain and 1 for control stress. Thus, the problem of significant damage due to total strain can be handled by using the DCF to scale the strain used in the damage calculation. The method used to calculate the DCF is as follows. First, the stress and strain in the pavement are computed using viscoelastic analysis. Because the main interest is vertical cracking, the radial stress and strain (denoted as veσ and veε , respectively) are used in this analysis. Now, the viscoplastic strain, vpε , can be approximated using the simple strain-hardening material model, Equation 3-36, and the viscoelastic stress computed in the previous step. Then, to approximate the structural effects of this viscoplastic strain, an effective viscoelastoplastic relaxation modulus, effiE , is computed by scaling the Prony coefficients in Equation 3-37 using Equation 3-38, as follows: 11 11 0 1 t pp q vp p dt Y ε σ ++  + =        ∫ . (3-36) where vpε = the viscoplastic strain, , ,p q Y = the material coefficients. ( ) 1 i tm i i E t E E e ρ − ∞ = = + ∑ (3-37) where ( )E t = the relaxation modulus, E∞ = the elastic modulus, iE = the modulus of the i th Maxwell element (fitting coefficient), and

54 iρ = the relaxation time (fitting coefficient). eff ve i i ve vp E Eε ε ε   =   +  (3-38) A second analysis was performed using the effective viscoelastoplastic Prony coefficients instead of the viscoelastic Prony coefficients. The stress and strain obtained from this analysis are referred to as the effective stress and strain and denoted as effveσ and eff veε , respectively. The mode-of-loading factor β is then obtained by comparing the effective stress and strain ( effveσ and eff veε ) to the original stress and strain ( veσ and veε ): (1 ) (1 ) eff eff ve ve ve eff ve eff E E σ β σ βσ ε β ε βε = − + = − + (3-39) where β is a factor that is 0 for controlled strain and 1 for controlled stress. Equation 3-39 gives two values of β , and hence their average is taken as the degree of controlled stress: 1.0 1.0 0.5 1.0 . 1.0 1.0 eff ve ve eff ve ve avg i i eff eff i i E E E E σ ε σ ε β   − −     = − +   − −     (3-40) Finally, the DCF is given by 1.0 1.0 0.5 1.0 1.0 1.0 eff ve ve eff ve ve i i eff eff i i DCF E E E E σ ε σ ε    − −     = + −   − −      . (3-41) Because the strain-hardening viscoplastic model is valid only for tension, the DCF is taken to be 1 for the region in which the radial stress is compressive.

55 3.1.3.3 Temperature variation The variation of temperature in a pavement has two effects: a change in stiffness of the AC and a change in the thermal stress due to thermal expansion of the material. Thermal stress is generated in the pavement depending on the boundary conditions. These two effects of temperature have been implemented in the FEP++. The actual temperature variation used for the pavement performance prediction is generated from the EICM. Temperature profiles generated from the EICM are input directly into the FEP++ preprocessor. 3.1.3.3.1 Temperature profile Five regions have been selected based on the six climatic zones identified by the AASHTO Guide for Design of Pavement Structures (25) to represent the wide range of temperature variation in the United States. The regions selected for this analysis are the same as those used in evaluating aging. They are (1) wet, freeze-thaw cycling (Washington, D.C.); (2) wet, hard freeze, spring thaw (Chicago, Illinois); (3) wet, no freeze (Tallahassee, Florida); (4) dry, freeze-thaw cycling (Dallas, Texas); and (5) dry, hard freeze, spring thaw (Laramie, Wyoming). Figure 3-25 and Figure 3-26 show the variation of temperature and temperature rates in 2001 for each of the five regions; the numbers in the legend signify the month. Each month’s data are the averaged data for that month at 0.7 cm below the surface of the pavement. Although a detailed method for pavement performance prediction is explained later, it has been determined that analysis on a monthly basis with three segments of a day might be appropriate for long-term simulations, such as 20 years. Based on temperature and rate profiles, three segments of a day are indicated by vertical divisions in the plots, as follows: 5:00 AM – 2:00 PM (5:00 – 14:00) for the heating cycle; 2:00 PM – 9:00 PM (14:00 – 21:00) for the cooling cycle; and 9:00 PM – 5:00 AM (21:00 – 5:00) for the constant cycle. Temperature profiles according to depth can also be obtained for

56 the same time segments and will be input for the simulation. Further, these segments generally follow daily traffic distributions. -10 0 10 20 30 40 50 0 5 10 15 20 25 Time (hr) Te m pe ra tu re (° C ) 1 2 3 4 5 6 7 8 9 10 11 12DC -10 0 10 20 30 40 50 0 5 10 15 20 25 Time (hr) Te m pe ra tu re (° C ) 1 2 3 4 5 6 7 8 9 10 11 12FL -10 0 10 20 30 40 50 0 5 10 15 20 25 Time (hr) Te m pe ra tu re (° C ) 1 2 3 4 5 6 7 8 9 10 11 12IL -10 0 10 20 30 40 50 0 5 10 15 20 25 Time (hr) Te m pe ra tu re (° C ) 1 2 3 4 5 6 7 8 9 10 11 12TX -10 0 10 20 30 40 50 0 5 10 15 20 25 Time (hr) Te m pe ra tu re (° C ) 1 2 3 4 5 6 7 8 9 10 11 12WY Figure 3-25. Temperature variation (at 0.7 cm below pavement surface)

57 -4 -3 -2 -1 0 1 2 3 4 5 0 5 10 15 20 25 Time (hr) R at e (°C /h r) 1 2 3 4 5 6 7 8 9 10 11 12DC -4 -3 -2 -1 0 1 2 3 4 5 0 5 10 15 20 25 Time (hr) R at e (°C /h r) 1 2 3 4 5 6 7 8 9 10 11 12FL -4 -3 -2 -1 0 1 2 3 4 5 0 5 10 15 20 25 Time (hr) R at e (°C /h r) 1 2 3 4 5 6 7 8 9 10 11 12IL -4 -3 -2 -1 0 1 2 3 4 5 0 5 10 15 20 25 Time (hr) R at e (°C /h r) 1 2 3 4 5 6 7 8 9 10 11 12TX -4 -3 -2 -1 0 1 2 3 4 5 0 5 10 15 20 25 Time (hr) R at e (°C /h r) 1 2 3 4 5 6 7 8 9 10 11 12WY Figure 3-26. Temperature rate variation (at 0.7 cm below pavement surface) 3.1.3.3.2 Stiffness change as a function of temperature For material level modeling, the change in stiffness in AC due to temperature is taken into account using the concept of reduced time (see Equation 3-42). The t-T shift factor is obtained from the characterization of the relaxation modulus of the material using dynamic modulus tests at different frequencies and temperatures. The reduced time is calculated with respect to a reference temperature and captures the history of temperature variations in the material.

58 ( ) ( )0 1t T t d a ξ τ τ = ∫ , (3-42) where ( )tξ = the reduced time, Ta = the time-temperature shift factor. In the computational implementation, all analyses were performed in terms of physical time, but the relaxation modulus was adjusted to the temperature of interest. This adjustment approach is similar in concept to the reduced time approach and is possible due to the use of a state variable formulation for stress and strain (cf. reference 26) and can be performed efficiently by adjusting the relaxation times in Equation 3-37 by Equation 3-43. ( ) ( ) ˆ ii T t a t ρ ρ = . (3-43) 3.1.3.3.3 Thermal stress Thermal stress in the material is incorporated by defining the mechanical strain, ( )m tε , as follows: ( ){ } ( ){ } { }( )0m Et t T Tε ε α= − − , (3-44) where ( ){ }tε = the strain in the material, { }Eα = the coefficient of isotropic thermal expansion for the material in three dimensions, T = the current temperature, and 0T = the reference temperature at which there is no initial stress. The constitutive law for a viscoelastic material undergoing damage is given as { } [ ]{ },RmDσ ε= (3-45)

59 where { }σ = the stress vector, [ ]D = the stiffness matrix that is a function of damage, and { }Rmε = the ve ctor of mechanical ps eudo s train, c omputed us ing t he s tate variable solution to pseudo strain (26). 3.1.4 Performance Prediction Module The performance of a pavement can be characterized by predicting the damage accumulation. Using a suitable failure criterion that relates the damage parameter to the onset of cracking in a uniaxial test, the crack initiation points in the pavement under fatigue cracking can be predicted. However, when the period of interest is in the order of years, it becomes computationally impossible to simulate pavement performance. Consequently, an extrapolation method is introduced so that the simulations run in a more time-efficient manner. 3.1.4.1 VECD FEP++ The VECD FEP++ is a finite-element implementation of the VECD model (27). The model assumes that a material is isotropic when undamaged and that damage growth under loading leads to local transverse isotropy (i.e., the material has a local axis of symmetry oriented along the maximum principal stress direction, which is perpendicular to the crack direction). This model is currently used in modeling tensile damage and offers capabilities for modeling the damage progression in a viscoelastic material and the three-dimensional state of stress and strain. It also accounts for the effects of anisotropy in the material as a result of cracking in one direction. The current framework of the VECD FEP++ is formulated for the axisymmetric case, but can be extended easily to three dimensions. The finite-element model (FEM) implementation enables the redistribution of stress in pavements that are being damaged, which is an important

60 phenomenon to be captured. The FEM implementation is also useful in studying the effects of layers with different stiffness values. By running the FEM simulation of a pavement for numerous cycles of fatigue loading, the progression of damage in the pavement over time may be observed, and critical points in the pavement that show early degradation can be noted. Using a suitable failure criterion that relates the damage parameter to the onset of cracking in a uniaxial test, the crack initiation points in the pavement under fatigue cracking can be predicted, and the life of the pavement can be better understood. Additionally, the program can be used in ranking the performance of pavements constructed from different asphalt mixtures. The performance of a pavement can be characterized by predicting the damage accumulation and using a failure criterion to judge the state of the pavement. This evaluation can be achieved by using an FEM simulation with the VECD model over the period of interest. However, when the period of interest is in the order of years, it becomes computationally impossible to simulate the pavement performance. In these cases, approximate schemes are introduced so that the simulations run in a more time-efficient manner. The FEP++ uses an extrapolation scheme that can significantly reduce the running time while still capturing the essential characteristics. The current scheme computes the damage caused by load and thermal variations at representative times in a day. These data are then extrapolated using a nonlinear scheme to obtain the total damage accumulation in a given month. This damage is then applied to the pavement as the initial condition for the next month’s simulation. This process is continued for the entire simulation period. 3.1.4.2 Traffic Traffic loading is a complicated characteristic to quantify and include in pavement simulations. Complexities arise due to the statistical distribution of the numerous different axle types available, each with its own load level distribution. Predicting the exact distribution of the

61 various vehicle loadings and their accompanying load level distributions is an imprecise science at best, and predicting the order of the loadings is impossible. Furthermore, predicting the sequential effects of different loading configurations becomes computationally time consuming. Although, at any instant in time, the exact traffic history will determine the pavement performance, over time an averaging effect occurs such that the exact sequence of loading is not that critical. In light of this statistical effect and in order to simulate pavement performance in a reasonable amount of time, the concept of equivalent single-axle loads (ESALs) is used. For all of the subsequent pavement simulations, loading consists of 10 million ESALs equally distributed over the three daily simulation periods (as defined in Section 3.2.3.3.1) and for a 20-year period. This distribution translates to approximately 41,700 loadings each month or 450 for each of the three daily simulation periods. Traffic growth factors and loading distributions other than uniform distribution could have been used, but these factors were not investigated in this study. Each load is assumed to last for 0.1 second and, as a consequence of the uniform loading distribution, the rest period between each loading is constant and equal to approximately 62.1 seconds. These loading and rest times ensure compatibility between the total simulation time and the physical time. The representative load consists of a single axle with an 80 kN (18 kip) load and tires inflated to 689 kPa (100 psi). It was assumed that the axle would be wide enough that the damage induced by each tire would be independent of that of the other. This assumption was verified through examinations of the damage contours that revealed little damage growth at distances of 0.5 m (1.6 ft) from the wheel center. Due to this damage independence and due to the axisymmetric nature of the current work, the simulations can take advantage of symmetry

62 and analyze a single side of the tire load. Currently, the effect of wheel wander is not considered in the simulations. 3.1.4.3 Extrapolation The basic idea behind extrapolation is to use the damage growth law to obtain an approximate expression for the damage growth rate. The growth law is given by 2( ) 2 RRS W C t S S α α ε ∂ ∂ ∂ = − = −   ∂ ∂ ∂   , (3-46) where RW = pseudo strain energy density function, ( )12 R RW Cε= , Rε = pseudo strain, t = total duration of a load cycle, C = pseudo stiffness, S = damage, and α = damage evolution rate. Using a change of variables and Equation 3-46, the growth rate of the pseudo stiffness can be written as 1 2( ) 2 RC S C C t t S S α ε+∂ ∂ ∂ ∂ = =  ∂ ∂ ∂ ∂  , (3-47) or equivalently as 1 2( ) 2 t R o C C dt N S α ε + ∂ ∂ =  ∂ ∂ ∫ , (3-48) where N is a load cycle. The pseudo stiffness and damage parameters, C and S , are related by ln( ) nC mS= , (3-49) which, upon differentiation, yields

63 ( ) 1 11 lnn n C m nC C S −∂ = ∂ . (3-50) Finally, assuming the pseudo strain, Rε , remains constant during a cycle, ( ) ( )( )( ) 11 11 ln n C AC C N αα + −+∂ = ∂ , (3-51) where A is the constant over the duration of the cycle. The above relationship can be solved numerically using the results from the FEM simulation as the initial condition to obtain an approximate value of the pseudo stiffness and, hence, the damage, at the end of the extrapolation period. The equations for the numerical scheme to solve Equation 3-51 are given by ( ) ( )( ) ( ) ( )( ) 1( 1)(1 )( 1) 1 1 1 1 1 1( 1)(1 )( 1) ln , , , ln nend begin begin begin i i i i begin end heal i i i end begin i i nbegin begin i i C C A C C C C C C CA C C αα αα + −+ + + + + + + −+ = + − = + ∆ − = × − (3-52) where begin iC = the pseudo stiffness at the beginning of cycle i , end iC = the pseudo stiffness at the end of cycle i , and heal iC∆ = the recovered stiffness due to healing in cycle i . 3.1.4.4 Damage calculation summary A flow chart of the total damage calculation process, including extrapolation, is shown in Figure 3-27. This flow chart outlines the combination of elements discussed previously: traffic loading, thermal considerations, damage factor, daily divisions for damage calculation, healing, etc. To verify the accuracy of some of the simplifications, particularly the extrapolation scheme, a simulation was performed using the FEP++ software without any extrapolation. In a real pavement simulation, the number of cycles simulated in an analysis step would be closer to

64 13,000, but predicting this many cycles without extrapolation would require a significant amount of time. Instead, the simulation was carried out for only about 1,300 cycles, and comparisons were made at this number of cycles. Initial C Average temperature of 30 days for 5:00 ~ 14:00 (10°C~30°C) Average temperature of 30 days for 14:00 ~ 21:00 (30°C~10°C) Average temperature of 30 days for 21:00 ~ 5:00 (10°C~10°C) ΔC1 ΔC2 ΔC3 ΔC= ΔC1+ ΔC2+ ΔC3 Extrapolated C Analysis period = 20 yrs Design ESALs = 10 million Duration of analysis segment ≈ 8 hrs Number of analysis segments per day = 3 Number of ESALs per segment ≈ 13889 Number of days per block = 30 Cycle jump = Max. 13889 * 30 Actual EICM Thermal Stress +One Loading +Thermal Damage at Rest Period Extrapolation for 13889 cycles Damage Factor Thermal Stress +One Loading +Thermal Damage at Rest Period Thermal Stress +One Loading +Thermal Damage at Rest Period Healing Figure 3-27. Extrapolation flow chart

65 The results of two different cases, (a) a healing-dominant case and (b) a damage-dominant case, are shown in Figure 3-28, along with the results based on the extrapolation. To compute the extrapolated value, the damage growth in the selected cell for the first cycle was extracted from the VECD-FEP++ analysis and input into a spreadsheet-implemented version of the extrapolation algorithm. The spreadsheet-based algorithm is used here for convenience only. In a full simulation, the extrapolation process is completed automatically in the VECD-FEP++ software. Figure 3-28 shows that the extrapolation algorithm effectively matches the output from the complete analysis at 1,300 cycles. The figure also shows that the agreement between the reference simulation and the extrapolated simulation decreases as the number of cycles increases, together with a tendency for the pseudo stiffness to approach an asymptotic value. Thus, it is expected that the differences observed in the extrapolation to 1,300 cycles would be similar to the differences observed at 13,000 cycles. Direct proof of this hypothesis was not performed due to the time necessary to create a reference simulation for 13,000 cycles. 0.59 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0 500 1000 1500 Number of Cycles C Extrapolation VECD-FEP++ (a) 0.66 0.67 0.68 0.69 0.70 0.71 0.72 0.73 0 500 1000 1500 Number of Cycles C Extrapolation VECD-FEP++ (b) Figure 3-28. Effect of extrapolation calculation on predicted damage growth: (a) healing- dominant case and (b) damage-dominant case 3.1.4.5 Thermal stress calculation In the current simulation extrapolation scheme, the temperature for a given month is averaged for each of the time segments identified in Section 3.2.3.3.1. This method is acceptable for computing the stress and damage due to mechanical loading, but could yield erroneous

66 thermal damage results because it effectively smoothes some of the more extreme events. The thermal stress-related damage calculations are much more sensitive to changes in temperature because these damage calculations are directly proportional to the temperature raised to the VECD damage power, α . This proportionality can be observed by substituting the thermal strain, Equation 3-44, into the pseudo strain calculation, Equation 3-53, and then combining the result with the derived evolution law, Equation 3-54. Conversely, the temperature dependence for mechanical damage enters indirectly into the constitutive model by changes in reduced time in the calculation of the pseudo strain. This indirect entry does not mean that the mechanical damage is not sensitive to temperature, only that it is less sensitive to temperature than is thermal damage. ( ) 0 1 tR R dE t d E d εε τ τ τ = −∫ , (3-53) ( ) ( ) 1 12 1 1 2 RdS C α α αε ξ + + = − ∆ × ∆    (3-54) To compute the thermal damage more accurately, the theoretically observed proportionality between thermal damage and temperature raised to the VECD α was used. Equation 3-55 shows the averaging function that was applied to generate the weighted average temperature. Figure 3-29 provides a comparison of the average and weighted average temperature values computed for a daily temperature cycle. In this figure, the lines represent the daily temperature cycles for the 31 days of December for Laramie, Wyoming. The circles represent the mathematical average temperature, and the squares represent the weighted average temperature. Figure 3-29 shows that the overall average temperatures are not that different from each other.

67 1 1 N i WA T T N α α =      =       ∑ (3-55) -15 -10 -5 0 5 10 0 5 10 15 20 Time (hr) Te m pe ra tu re (° C ) Average Weighted Average -6 -4 -2 0 2 4 6 8 0 5 10 15 20 Time (hr) Te m pe ra tu re R at e (° C ) Figure 3-29. Original versus weighted average temperature and temperature rate for WY-Dec To assess the consequences of these differences, the thermally induced damage was computed for each case and used in the extrapolation scheme. Coincident with this prediction, a rigorous time-wise analysis was performed using the daily temperature variations. Comparisons between the rigorous calculation and the average temperature calculations are shown in Figure 3-30. In this figure, the weighted average and arithmetic average predictions are shown only at the end of the month because they have been included in the extrapolation process. From Figure 3-30 it is clear that, even though the two averaging techniques yield only a slight difference according to Figure 3-29, the weighted average technique yields the more accurate solution for thermal-related damage. For computational efficiency, the FEP++ uses only a single thermal gradient to compute the damage due to both thermal and mechanical loading. Thus, it is important that the method that yields acceptably accurate calculations for the thermal case does not yield unrealistic answers for the mechanical loading. The effect of using the weighted average technique in lieu of the arithmetic average is assessed by performing a limited simulation using both methods. The

68 results of this simulation are shown in Figure 3-31. In this figure, the x-axis shows the damage calculated via the arithmetic average temperature, and the y-axis shows the damage calculated via the weighted average temperature. These calculations were performed for each of the three daily time increments. Because both approaches yield basically the same amount of mechanical- related damage, and because the arithmetic average technique results in accuracy issues for thermal damage, the weighted average temperature was chosen for further analysis. 0.793 0.794 0.795 0.796 0.797 0.798 0.799 0.800 0 5 10 15 20 25 30 Time (day) C Actual Average Weighted Average Figure 3-30. Damage accumulation at loading edge of pavement surface for WY-Dec

69 1E-06 1E-05 1E-04 1E-03 1E-02 1E-06 1E-05 1E-04 1E-03 1E-02 Reduction of C - Arithmetic Average R ed uc tio n of C - W ei gh te d A ve ra ge Figure 3-31. Arithmetic versus weighted average temperature and temperature rate for WY-Dec 3.1.5 Output Module The output module consists of the tools and techniques necessary to view and interpret the VECD-FEP++ performance predictions. It creates a single file that can be opened, processed, and manipulated to view visual interpretations (contours) of the predicted damage, stress distribution, or other quantities of interest. This file can also be processed to extract different indices to quantify the visual observations. 3.1.5.1 Output in the form of a contour plot Example contour plots and typical values are shown in Figure 3-32 a, b, c, and d for radial stress, radial strain, pseudo stiffness, and the condition index, respectively. Through qualitative observations of these contour plots, the damage distribution and evolution in the overall pavement structure can be captured easily. More importantly, because the VECD-FEP++ does not assume the location where the damage is concentrated, and usually the severe damage occurs in small areas (i.e., a few structural elements as compared to several thousands of elements within the whole pavement structure), the damage concentration needs to be tracked for the

70 entire structure. The contour plots in Figure 3-32 show the grids that represent the elements that build the AC layers in the pavement structure. For the sensitivity analysis, the element size of 1.4 cm by 1.4 cm was determined for the AC layers. Moreover, significant damage in an element of this size could be considered to be the initiation of cracking, based on field observations. 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 Srr 600 400 200 0 -200 -400 -600 -800 -1000 -1200 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 Err 2.0E-04 1.5E-04 1.0E-04 5.0E-05 0.0E+00 -5.0E-05 -1.0E-04 -1.5E-04 -2.0E-04 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 C 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 CI 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (a) (b) (c) (d) Figure 3-32. Examples of contour plots for: (a) radial stress, (b) radial strain, (c) pseudo stiffness, and (d) the condition index 3.1.5.2 Quantification of cracking prediction A scaled parameter, referred to as the condition index (CI) and defined by Equation 3-56, was computed for each element in the pavement structure to interpret the loss in material integrity. The CI is a variable that ranges from 0 to 1, with 0 being failure according to the material level criteria discussed in Subsection 3.2.2.2 (see Figure 3-13), and 1 representing a completely intact body. , , i f i intact f i C C CI C C − = − (3-56) where CI = the condition index, intactC = the intact pseudo stiffness value, iC = pseudo stiffness at instant i , and f,iC = failure pseudo stiffness at instant i .

71 3.1.5.3 Output in the form of a stacked bar graph Examining the damage contours has some advantages, but it is not the most efficient analytical method for evaluating the many cases that are simulated in this work. To better quantify the simulation results, a quantity termed the condition index area is employed. This quantity is defined as the percentage of the total investigation area with CI values of 0-0.2, 0.2- 0.4, 0.4-0.6, 0.6-0.8, and 0.8-1. To appropriately represent the pavement condition, the area from the load center to 0.5 meters in the horizontal direction and from the pavement surface to the bottom of the AC layer in the vertical direction has been chosen for the analysis (see Figure 3-32). The respective effects of top-down and bottom-up cracking can be determined using this quantity by separately calculating the CI for the upper and lower halves of the pavement structure. Once computed for each block (0-0.2, 0.2-0.4, etc.) the CI area (%) can be shown graphically using a stacked bar graph. Figure 3-33 shows an example stacked bar graph for the CI area that was calculated from the CI contour plot in Figure 3-32 (d). CI<0.2 0.2<CI<0.4 0.4<CI<0.6 0.6<CI<0.8 0.8<CI<1 0 20 40 60 80 100 Example Example C I A re a (% TOP BOTTOM C I A re a (% ) Figure 3-33. Example of condition index area

72 3.1.6 Model Integration Reasonableness and sensitivity studies were carried out by a series of FEP++ simulations to show whether or not the framework development and integration into the VECD-FEP++ were performed correctly. First, a reasonableness check was performed to verify the proper implementation of the VECD model, and hence, the analytical sub-models (healing, aging, thermal damage, etc.) were not included. Second, a sensitivity check was conducted to investigate the effect of each analytical sub-model on pavement performance. From the sensitivity studies it was revealed that the sub-models affect pavement performance in different ways and to varying degrees (details of these studies are illustrated in Appendix A). Two pavement structures, as shown in Figure 3-34, were investigated: thin (127 mm or 5 in.) and thick (304.8 mm or 12 in.). Two AC mixes, the Control and SBS mixes used in the FHWA ALF research, were selected for the simulation because all the necessary model data for both mixes were available, and it has been reported that both mixes show significantly different performance behavior (28). The aggregate gradation, asphalt content, and target air void levels are the same for both mixtures, as shown in Table 3-5. AC Layer (127mm) Base Layer (203mm) Subgrade (infinite) AC Layer (305mm) Subgrade (infinite) CL CL Thin Pavement Thick Pavement Figure 3-34. Pavement structures for FEP++ simulation

73 Results from the material characterization show large differences in the material responses (see Figure 3-35). The base layer for the thin pavement structure was 203 mm (8 in.) thick, and the subgrade for both structures was considered to be semi-infinite. The AC layers were modeled using the VECD model and additional analytical models to include such factors as aging, healing, damage factors, and thermal stress. The unbound material layers were assumed to be linear elastic at three different levels of modulus values (weak, medium, strong), as presented in Table 3-6. The moduli of both the base layer and subgrade materials were varied to yield either a weak, moderate, or strong support condition. Depending on the purpose of the simulation, the appropriate pavement temperature was selected from the five regional temperature profiles generated from the EICM. The temperatures varied in accordance with time and depth within a given AC layer. To simulate the variations in temperature with depth, each row of elements in the FEM mesh was assigned a different temperature consistent with the EICM output at the average depth of the row. A moving load was simulated by applying a 0.1-second haversine loading pulse with a magnitude of 40 kN (9 kip) and contact pressure of 689 kPa (100 psi) on the pavement surface, followed by 62.1 seconds of rest.

74 Table 3-5. Volumetric Properties of Laboratory Mixtures for ALF Pavements Mix Designation Control SBS Binder Type Unmodified Styrene Butadiene Styrene Binder Grade PG 70-22 PG 70-28 Binder Content 5.30% NMSA* 12.5 mm Target Air Voids 4% Sieve Size Gradation, % Passing 37.5 mm 100.0 25.0 mm 100.0 19.0 mm 100.0 12.5 mm 93.8 9.5 mm 85.2 4.75 mm 56.0 2.36 mm 35.6 1.18 mm 25.1 0.600 mm 18.4 0.300 mm 13.1 0.150 mm 9.3 0.075 mm 6.7 * Nominal Maximum Size Aggregate 0 5000 10000 15000 20000 25000 30000 35000 40000 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 Reduced Frequency (Hz) Control SBS (a) 100 1000 10000 100000 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 Reduced Frequency (Hz) (b) 0 5 10 15 20 25 30 35 40 45 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 Reduced Frequency (Hz) (c) 0.0 0.2 0.4 0.6 0.8 1.0 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 S C (d) Ph as e An gl e (d eg ) |E *| (M Pa ) |E *| (M Pa ) Figure 3-35. Characteristics of Control and SBS mixtures for FEP++ simulation: (a) |E*| mastercurve in semi-log scale; (b) |E*| mastercurve in log-log scale; (c) phase angle mastercurve; and (d) damage characteristic curve

75 Table 3-6. Elastic Modulus of Unbound Layers Layer Elastic Modulus, MPa (ksi) Weak Moderate Strong Base Layer 138 (20) 276 (40) 552 (80) Subgrade 41 (6) 83 (12) 166 (24) To save computational time for the example simulation and parametric study, the simulations were performed until one element reached a pseudo stiffness value of 0.25 or until 10 years had been simulated. In the case of the thin pavement, the length of time was somewhat longer than a full year, but the thick pavement was simulated for the full 10 years. The pseudo stiffness value of 0.25 has been found to describe the failure point of monotonic direct tension tests and is a conservative value for the failure point of cyclic fatigue tests. Based on experience and some limited trials, it was found that when one element reached a pseudo stiffness value of 0.25, the CI was zero or close to zero for many elements (including those neighboring the element with the pseudo stiffness value of 0.25). The failure criterion for any one element, e.g., that used in defining the CI, was always based on the criterion described in Section 3.2.2.2. 3.2 Example Simulation 3.2.1 Introduction An example simulation of the FEP++ using the VECD model and all accompanying analytical sub-models was carried out to demonstrate the capabilities of the modeling approach. For this example problem, a single region – wet, freeze-thaw cycling for Washington, D.C. – was chosen. Only a single support condition, moderate, was used for this example, and all other factors were fixed, except the pavement thickness, which was evaluated at the 127 mm (5 in.) and 305 mm (12 in.) levels. All of the analytical sub-models were active for these simulations. The simulation details are summarized in Table 3-7.

76 Table 3-7. Simulation Details for Example Simulation Item Number of Cases Details Region 1 DC Structure 2 Thin (127 mm or 5 in.), Thick (304.8 mm or 12 in.) AC Material 1 ALF Control Support Condition 1 Moderate Aging 1 Yes Healing 1 High Thermal Stress 1 High Thermal Coefficient Damage Factor 1 5-yr Average EICM Climate 1 5-yr EICM Repeated Load Level 1 40 kN (9 kip) Contact Pressure 1 689 kPa (100 psi) Total Number of Cases 2 3.2.2 Simulation Results and Discussion Results from the simulation are presented by stacked bar graphs for the CI areas shown for the thin and thick pavements in Figure 3-36 and Figure 3-37, respectively (contour plots for the same simulation results can be found in Appendix A). In Figure 3-36 and Figure 3-37, the simulation month number (all simulations began in July) and the abbreviated month that corresponds to the simulation are shown on the x-axis. Figure 3-36 shows the results for the months of October through July of the first simulation year for the thin pavement. Figure 3-37 presents the final nine months of the thick pavement simulation. A few observations can be made from Figure 3-36 and Figure 3-37. First, it is evident that the damage in the thin pavement section is much larger than the damage in the thick pavement at the equivalent time. This result is to be expected because both simulations were conducted for the equivalent of a 10 million ESAL. It is also observed from examining the damage evolution between the months of April and July that pavement healing may constitute a major component of a pavement’s total fatigue performance. However, it is seen, through the mechanisms associated with repeated loading, healing and aging, that the thick pavement shows an area of concentrated damage at the pavement surface near the wheel load edge. Therefore, it appears that

77 healing is more effective near the pavement surface of a thin pavement than that of a thick pavement (i.e., the thick pavement appears to heal better at the bottom of the pavement than does the thin pavement). CI<0.2 0.2<CI<0.4 0.4<CI<0.6 0.6<CI<0.8 0.8<CI<1 0 20 40 60 80 100 4 (O ct ) 5 (N ov ) 6 (D ec ) 7 (J an ) 8 (F eb ) 9 (M ar ) 10 (A pr ) 11 (M ay ) 12 (J un ) 13 (J ul ) Cumulative Month C I A re a (% ) TOP 0 20 40 60 80 100 4 (O ct ) 5 (N ov ) 6 (D ec ) 7 (J an ) 8 (F eb ) 9 (M ar ) 10 (A pr ) 11 (M ay ) 12 (J un ) 13 (J ul ) Cumulative Month C I A re a (% ) BOTTOM Figure 3-36. Condition index area (%) for thin pavement simulation CI<0.2 0.2<CI<0.4 0.4<CI<0.6 0.6<CI<0.8 0.8<CI<1 0 20 40 60 80 100 11 2 (O ct ) 11 3 (N ov ) 11 4 (D ec ) 11 5 (J an ) 11 6 (F eb ) 11 7 (M ar ) 11 8 (A pr ) 11 9 (M ay ) 12 0 (J un ) Cumulative Month C I A re a (% ) TOP 0 20 40 60 80 100 11 2 (O ct ) 11 3 (N ov ) 11 4 (D ec ) 11 5 (J an ) 11 6 (F eb ) 11 7 (M ar ) 11 8 (A pr ) 11 9 (M ay ) 12 0 (J un ) Cumulative Month C I A re a (% ) BOTTOM Figure 3-37. Condition index area (%) for thick pavement simulation 3.3 Model Evaluation: Parametric Study 3.3.1 Introduction A systematic evaluation using the VECD-FEP++ code was carried out to gain insight into the parameters most responsible for the development and/or dominance of top-down cracking. Combinations of three regions (DC, FL, and WY), two structures (thin and thick pavements), two materials (Control ALF and SBS ALF) and two support conditions (weak and strong) were investigated so that simple and interactive variables could be studied. The factorial levels are

78 summarized in Table 3-8. The most important factors identified, for which all the necessary inputs were available, are material type (modified versus unmodified), pavement thickness (thick versus thin), climate (wet, no freeze versus wet, freeze-thaw cycling versus dry, hard freeze, spring thaw), and support condition (weak versus strong). The CI area graphs are shown for the thin and thick pavements in Figure 3-38 and 3-39, respectively. The results for the thin pavements are all at the end of the ninth month (March), whereas the results for the thick pavements are also for March, but at year 10 (contour plots for the same simulation results are provided in Appendix A to show overall trends). The naming convention used for the different simulations is by letter and in the following order: pavement type, thick (T) or thin (t); material type, ALF control (C) or ALF SBS (S); and support condition, weak (W) or strong (S). Table 3-8. Simulation Details for the Parametric Study Item Number of Cases Details Region 3 DC, FL, WY Structure 2 Thin (127 mm or 5 in.), Thick (304.8 mm or 12 in.) AC Material 2 ALF Control vs. ALF SBS Support Condition 2 Weak vs. Strong Aging 1 Yes Healing 1 High Thermal Stress 1 High Thermal Coefficient Damage Factor 1 5-yr Average EICM Climate 1 5-yr EICM Repeated Load Level 1 40 kN (9 kip) Contact Pressure 1 689 kPa (100 psi) Total Number of Cases 24 Combination of all cases The stacked bar graphs for the CI areas were calculated for March to best demonstrate the observed trends. The magnitude of these trends varies somewhat from month to month due to the relative effects of the different analytical sub-models. Healing, for example, may substantially reduce the amount of accumulated damage during the summer months by slightly different amounts in Wyoming and in Florida. However, it was observed that even with these differences,

79 the major trends do not change substantially, and that these trends are observed most clearly in March. CI<0.2 0.2<CI<0.4 0.4<CI<0.6 0.6<CI<0.8 0.8<CI<1 0 20 40 60 80 100 t-C-W- FL t-C-W- DC t-C-W- WY t-C-W- FL t-C-W- DC t-C-W- WY C I A re a (% ) TOP BOTTOM 0 20 40 60 80 100 t-C-S- FL t-C-S- DC t-C-S- WY t-C-S- FL t-C-S- DC t-C-S- WY C I A re a (% ) TOP BOTTOM 0 20 40 60 80 100 t-S-W- FL t-S-W- DC t-S-W- WY t-S-W- FL t-S-W- DC t-S-W- WY C I A re a (% ) TOP BOTTOM 0 20 40 60 80 100 t-S-S- FL t-S-S- DC t-S-S- WY t-S-S- FL t-S-S- DC t-S-S- WY C I A re a (% ) TOP BOTTOM Figure 3-38. Condition index area (%) for thin pavement

80 CI<0.2 0.2<CI<0.4 0.4<CI<0.6 0.6<CI<0.8 0.8<CI<1 0 20 40 60 80 100 T-C-W- FL T-C-W- DC T-C-W- WY T-C-W- FL T-C-W- DC T-C-W- WY C I A re a (% ) TOP BOTTOM 0 20 40 60 80 100 T-C-S- FL T-C-S- DC T-C-S- WY T-C-S- FL T-C-S- DC T-C-S- WY C I A re a (% ) TOP BOTTOM 0 20 40 60 80 100 T-S-W- FL T-S-W- DC T-S-W- WY T-S-W- FL T-S-W- DC T-S-W- WY C I A re a (% ) TOP BOTTOM 0 20 40 60 80 100 T-S-S- FL T-S-S- DC T-S-S- WY T-S-S- FL T-S-S- DC T-S-S- WY C I A re a (% ) TOP BOTTOM Figure 3-39. Condition index area (%) for thick pavement 3.3.2 Effect of Region A comparison of the results at a consistent thickness and time for pavements constructed with similar materials indicates that the more extreme and cooler climates, as compared to less extreme and warmer climates, tend to show a higher concentration of highly damaged areas at both the top and bottom of the pavement (e.g., damage accumulates much more quickly in Wyoming than in Florida) because of the relatively stiff PG 70-22 binder used in the mixtures. Thus, mixtures with these binders would likely not be used in a climate such as that of Wyoming. The key factors that may lead to this increased damage accumulation are (1) increased thermal damage potential and (2) decreased healing potential. Both factors worsen due to extreme temperature variations and generally cool temperatures. Longitudinal thermal stress was used to calculate the associated thermal damage, which tends to overestimate the thermal effect. In the future, transverse thermal stress should be used, which is limited by the maximum

81 frictional resistance that can develop between the HMA surface and base layers The SBS material tends to soften the regional differences slightly as does the use of a strong support condition. The support condition effect is more prevalent for the thin pavement. Care should be taken when considering these support conditions, though, because they do not account for differences due to support layer moisture or freeze/thaw action. 3.3.3 Effect of Structure Structural effects are found to be similar to those shown previously in Figure 3-36 and Figure 3-37. Specifically, for the same traffic conditions, a thin pavement shows more significant damage earlier than a thick pavement. Also, thin pavements tend to show high damage zones at both the top and the bottom of the pavement structure that start to grow from the early stages of loading. Thick pavements, on the other hand, tend to develop damage at a much slower rate than thin pavements and also tend to show higher concentrations of damage near the pavement surface. These differences are exacerbated by cooler and more extreme temperature conditions. It should be noted that the designation of a pavement as thin or thick is somewhat arbitrary; it is defined as thin or thick relative to traffic loadings. For example, a 127 mm (5 in.) pavement may be considered thin for the case of 10 million ESALs, but would be considered thick for a significantly less loading (e.g., 100,000 ESALs). 3.3.4 Effect of Unbound Layer Property The support layer stiffness has an important effect on the distribution of damage in the pavement structure. In general, weak support results in more extensive and more severe damage growth than strong support. This effect is more pronounced in thin pavements, but tends to have a more significant effect on damage at the top of the pavement in thick pavements. The effect is also more pronounced in cooler and more extreme climatic conditions.

82 3.3.5 Effect of Asphalt Mixture Properties From the thin pavement simulation, it appears that the SBS material results in improved performance with regard to damage at the bottom of the pavement but slightly poorer performance in terms of damage at the top of the pavement for the one month considered. The damage progression from October through June is shown in Figure 3-40. From examining the behavior at the end of June, it is clear that the SBS material shows an overall benefit in the minimum CI levels for both top-down and bottom-up related damage. This difference occurs because the damage in SBS pavements tends to slow considerably during the summer months due to the overall lower modulus of the SBS material relative to the Control pavement and the tendency of the SBS pavement to exhibit more viscoplastic behavior at high temperatures and thus to have a more favorable DCF than the Control pavement. However, the healing potential function applied in this simulation is the same for both the Control and SBS materials. The same basic mechanisms that are active in the thin pavement are also active in the thick pavement, although the effects are somewhat less noticeable due to the less extreme temperature variations at the bottom of the pavement structure in the thick pavement. The effect of the SBS material is also seen more clearly in the thick pavement, irrespective of potential differences in healing and the DCF. In the early stages (year 1), the SBS pavements tend to show a lower minimum CI value than do the Control materials. However, a transition begins to occur around year 5 when the SBS material tends to show a larger damaged area, but this damaged area has a higher minimum CI value than the corresponding Control section. This effect appears to be slightly more pronounced for the Wyoming climatic condition.

83 0 20 40 60 80 100 4 (O ct ) 6 (D ec ) 8 (F eb ) 10 (A pr ) 12 (J un ) 4 (O ct ) 6 (D ec ) 8 (F eb ) 10 (A pr ) 12 (J un ) Cumulative Month C I A re a (% ) TOP BOTTOMt-C-W-DC 0 20 40 60 80 100 4 (O ct ) 6 (D ec ) 8 (F eb ) 10 (A pr ) 12 (J un ) 4 (O ct ) 6 (D ec ) 8 (F eb ) 10 (A pr ) 12 (J un ) Cumulative Month C I A re a (% ) TOP BOTTOMt-S-W-DC CI<0.2 0.2<CI<0.4 0.4<CI<0.6 0.6<CI<0.8 0.8<CI<1 Figure 3-40. Damage progression and healing in SBS and Control pavements 3.4 Summary of Findings An enhanced VECD-based model for predicting the initiation of top-down cracking in HMA layers has been established in this portion of the study. This effort was accomplished by developing, modifying, and/or investigating several important material property models, such as an aging model, healing model, failure criteria, viscoplasticity, and thermal stress, and then incorporating these sub-models into the existing VECD model. The material models were converted to and/or combined with the structural models. These sub-models were implemented into the VECD-FEP++, and an extrapolation method for predicting top-down cracking initiation was developed. Reasonableness and sensitivity studies were undertaken to verify that the framework development and implementation into the VECD-FEP++ were performed correctly. A sensitivity study to investigate the effect of each sub-model on pavement performance revealed that each sub-model affects pavement performance in different ways and to varying degrees. To demonstrate the full capabilities of the VECD-FEP++, two example simulations were carried out, and the results indicated that the interactions among the sub-models and overall trends in terms of pavement behavior were captured successfully.

84 A systematic evaluation using the VECD-FEP++ was carried out to gain insight into the parameters most responsible for the development and/or dominance of top-down cracking. In general, more damage at both the top and bottom of the pavement was observed in the case of a cold climate, thin structure, weak support layer, and the use of the Control mixture. However, because the longitudinal thermal stress was used to calculate the associated thermal damage, the thermal effect in cold climates may have been overestimated. It was also found that thin pavements tend to show more damage at the bottom of the pavement, whereas thick pavements tend to show a higher concentration of damage at the top of the pavement.