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

Publication Number: FHWAHRT08073 Date: September 2009 
PDF Version (3.85 MB)
PDF files can be viewed with the Acrobat® Reader®
The constitutive modeling of asphalt concrete behavior is a topic that has gained national importance in the past few years. Such modeling efforts have the explicit goal of providing for better design and analysis of asphalt pavement structures to resist failure and/or better predict when failure will occur. These efforts should thus provide the tools necessary to better utilize available resources and/or to gain maximum results from limited resources. One such modeling effort that encompasses the two main forms of pavement distress, cracking and permanent deformation, is the multiaxial viscoelastoplastic continuum damage (MVEPCD) model and finite element package, finite element program (FEP++). The MVEPCD model combines elements of viscoelasticity, continuum damage mechanics, and viscoplasticity to model the material behavior, and FEP++ is used to model the interaction of material and structure.
The MVEPCD model has been characterized and verified using asphalt concrete mixtures tested at the Federal Highway Administration's Accelerated Load Facility in McLean, VA. A novel approach to modeling this process is suggested and verified in this work. In light of practical concerns related to constant rate tests using the Asphalt Mixture Performance Tester and due to the complexities of performing true timedependent analysis of cyclic fatigue tests, a refined and simplified viscoelastoplastic continuum damage model is presented. A robust FEP++ has been developed to account for the effects of loading and boundary conditions. Analysis can be performed in either twodimensional or threedimensional configurations. The resulting predictions are deemed reasonable and, thus, a reliable simulation of pavement response.
Cheryl Allen Richter
Acting Director, Office of Infrastructure
Research and Development
This document is disseminated under the sponsorship of the U.S. Department of Transportation in the interest of information exchange. The U.S. Government assumes no liability for the use of the information contained in this document.
The U.S. Government does not endorse products or manufacturers. Trademarks or manufacturers' names appear in this report only because they are considered essential to the objective of the document.
The Federal Highway Administration (FHWA) provides highquality information to serve Government, industry, and the public in a manner that promotes public understanding. Standards and policies are used to ensure and maximize the quality, objectivity, utility, and integrity of its information. FHWA periodically reviews quality issues and adjusts its programs and processes to ensure continuous quality improvement.
1. Report No. FHWAHRT08073 
2. Government Accession No. 
3. Recipient's Catalog No. 

4. Title and Subtitle Development of a Multiaxial Viscoelastoplastic Continuum Damage Model for Asphalt Mixtures 
5. Report Date September 2009 

6. Performing Organization Code 

7. Author(s) Y. Richard Kim, Ph.D. P.E., M.N. Guddati, Ph.D., B.S. Underwood, T.Y. Yun, V. Subramanian, S. Savadatti 
8. Performing Organization Report No. 

9. Performing Organization Name and Address North Carolina State University 
10. Work Unit No. (TRAIS) 

11. Contract or Grant No. DTFH6105H00019 

12. Sponsoring Agency Name and Address Office of Research and Technology Services 
13. Type of Report and Period Covered 

14. Sponsoring Agency Code FHWA 

15. Supplementary Notes Project performed under FHWA DTFH6105H00019; FHWA AOTR: Katherine Petros 

16. Abstract This report highlights findings from the FHWA DTFH6105H00019 project, which focused on the development of the multiaxial viscoelastoplastic continuum damage model for asphalt concrete in both compression and tension. Asphalt concrete pavement, one of the largest infrastructure components in the United States, is a complex system that involves multiple layers of different materials, various combinations of irregular traffic loading, and various environmental conditions. The performance of this structure is closely related to the performance of asphalt concrete. To predict the performance of asphalt concrete with reasonable accuracy, a better understanding of its deformation behavior under realistic conditions is urgently needed. Over the past decade, the authors have been successful in developing uniaxial material models that can accurately capture various critical phenomena such as microcrackinduced damage that is critical in fatigue modeling, strainrate temperature interdependence, and viscoplastic flow, which is critical for high temperature modeling. The resulting model is termed the viscoelastoplastic continuum damage model. However, to consider the complicated nature of inservice stress states, a multidimensional model is needed. To predict the performance of the real pavement structures, it is also important to incorporate the material model in a pavement model that considers the vehicle and climatic loads as well as the boundary conditions; the inhouse finite element package has been developed for this purpose. 

17. Key Words Constitutive modeling, Asphalt concrete, Viscoelastic, Viscoplastic, Multiaxial, VECD, Damage mechanics 
18. Distribution Statement No restrictions. This document is available through the National Technical Information Service, Springfield, VA 22161. 

19. Security Classif. (of this report) Unclassified 
20. Security Classif. (of this page) Unclassified 
21. No. of Pages 264 
22. Price 
Form DOT F 1700.7 (872) Reproduction of completed page authorized
SI (Modern Metric) Conversion Factors
Figure 1. Graph. Schematic representation of dynamic modulus shifting process with unshifted data
Figure 2. Graph. Schematic representation of dynamic modulus shifting process with shifted data
Figure 4. Graph. Constant crosshead test results in stressstrain space
Figure 5. Graph. Constant crosshead test results in stresspseudo strain space
Figure 6. Graph. Comparison of refined and approximate damage calculation techniques
Figure 7. Graph. Comparison of refined and approximate damage characteristic relationship
Figure 8. Illustration. Mechanical analog for the viscoplastic model
Figure 9. Illustration. Isotropic hardening diagram
Figure 10. Illustration. Kinematic hardening diagram
Figure 11. Illustration. Strain decomposition from creep and recovery testing
Figure 12. Illustration. Typical yield surface of HISS model
Figure 13. Graph. Mixture gradation chart
Figure 14. Graph. Comparison of test lane and laboratory gradations
Figure 15. Graph. Stress history of VL testing (unconfined and 140 kPa confinement VL)
Figure 16. Graph. Stress history of VL testing (500 kPa confinement VL)
Figure 17. Graph. Stress history of VLT testing (140 kPa confinement)
Figure 18. Graph. Stress history of VLT testing (500 kPa confinement)
Figure 19. Graph. Effect of 500 kPa confining pressure on the dynamic modulus in semilog space
Figure 20. Graph. Effect of 500 kPa confining pressure on the dynamic modulus in loglog space
Figure 21. Graph. Effect of 500 kPa confining pressure on observed elasticity in the Control mixture
Figure 22. Graph. Effect of 500 kPa confining pressure on the log shift factor function
Figure 23. Graph. Effect of confining pressure on the relaxation spectrum
Figure 24. Graph. Use of the uniaxial relaxation spectrum for multiaxial test results
Figure 26. Graph. Multiaxial equilibrium characterization results
Figure 27. Graph. Multiaxial dynamic modulus model strength
Figure 28. Graph. C_{11} versus S for tension for Control2006 mixture (5 °C reference)
Figure 29. Graph. Effect of timedependent Poisson's ratio on C_{12} calculation
Figure 30. Graph. C_{12} characteristic curve for Control2006 mixture
Figure 31. Graph. C_{12} characteristic curve for Control2006 mixture in semilogarithmic space
Figure 33. Graph. C_{12 }versus S for tensile loading for Control2006 mixture (5 °C reference)
Figure 46. Graph. C_{22} as a function of S calculated by different methodologies for 51T
Figure 47. Graph. C_{22} as a function of S calculated by different methodologies for 53T
Figure 48. Graph. C_{22} as a function of S calculated by different methodologies for 54T
Figure 49. Graph. C_{22} as a function of S calculated by different methodologies for 55T
Figure 51. Graph. Effect of confining pressure on coefficient Y
Figure 59. Graph. tTS with growing damage under confinement verification at a 0.0001_{ε} level
Figure 60. Graph. tTS with growing damage under confinement verification at a 0.0005_{ε} level
Figure 61. Graph. tTS with growing damage under confinement verification at a 0.001_{ε} level
Figure 62. Graph. tTS with growing damage under confinement verification at a 0.0022_{ε} level
Figure 63. Graph. tTS with growing damage under confinement verification at a 0.004_{ε} level
Figure 64. Graph. tTS with growing damage under confinement verification at a 0.007_{ε}
Figure 65. Graph. Effect of 500 kPa confining pressure on strength mastercurves
Figure 66. Graph. Effect of 500 kPa confining pressure on ductility in constant crosshead rate tests
Figure 67. Illustration. A schematic representation of the concept of average dC/dt
Figure 69. Illustration. Mathematical equivalence of the formulation used by Lee, Daniel, and Kim
Figure 85. Graph. Stressstrain curves for unconfined constant strainrate tests
Figure 86. Graph. Stressstrain curves for 500 kPa confinement constant strainrate tests
Figure 87. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 5 °C
Figure 88. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 25 °C
Figure 89. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 40 °C
Figure 90. Graph. Comparison of 500 kPa confinement and unconfined constant rate tests for 55 °C
Figure 91. Graph. Viscoplastic strain versus cumulative loading time (unconfined VL)
Figure 92. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VL)
Figure 93. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement VL)
Figure 94. Graph. Viscoplastic strain versus cumulative loading time (unconfined VT testing)
Figure 97. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement CLT)
Figure 98. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VT)
Figure 115. Graph. Stress mastercurves for the Control mixture under uniaxial conditions
Figure 116. Graph. Stress mastercurves for the Control mixture under triaxial conditions (500 kPa)
Figure 117. Graph. Variation of strain rate during unloading
Figure 119. Graph. Confining stress effect on the relaxation modulus
Figure 122. Graph. Effect of test method on shift factor functions
Figure 125. Diagram. MVECD model characterization
Figure 126. Graph. C_{11} versus S for compression for Control mixture (5 °C reference)
Figure 127. Graph. C_{12} versus S for compression for Control mixture (5 °C reference)
Figure 128. Graph. C_{22} versus S for compression for Control mixture (5 °C reference)
Figure 129. Graph. Comparison of tension and compression of C_{11} damage function
Figure 130. Graph. Comparison of tension and compression of C_{12} damage function
Figure 131. Graph. Comparison of tension and compression of C_{22} damage function
Figure 133. Graph. Determined fitting results and coefficients of function a(t_{p})
Figure 134. Graph. Determined fitting results and coefficients of function D(t_{p} _{σ})
Figure 135. Graph. VT predictions
Figure 136. Graph. VL predictions
Figure 137. Graph. CLT predictions (2.0 MPa deviatoric stress0.1s pulse time)
Figure 138. Graph. CLT predictions (2.0 MPa deviatoric stress0.4s pulse time)
Figure 139. Graph. CLT predictions (2.0 MPa deviatoric stress1.6s pulse time)
Figure 140. Graph. CLT predictions (2.0 MPa deviatoric stress6.4s pulse time)
Figure 141. Graph. CLT predictions (1.8 MPa deviatoric stress1.6s pulse time)
Figure 142. Graph. CLT predictions (2.2 MPa deviatoric stress1.6s pulse time)
Figure 143. Graph. Compressive and tensile peak stress in SQRT(J_{2})  I_{1} space
Figure 144. Graph. Determined γ_{0} parameter function
Figure 145. Graph. Determined R parameter function
Figure 146. Graph. Determined n parameter function
Figure 147. Graph. Determined α_{0} parameter function
Figure 148. Graph. Ratedependent initial yield surface
Figure 149. Graph. Variation of α for 1,800 kPa CLT loading (500 kPa confinement)
Figure 150. Graph. Viscoplastic strain predictions for VT tests (500 kPa confinement)
Figure 151. Graph. Viscoplastic strain predictions for CLT tests (500 kPa confinement)
Figure 152. Graph. Viscoplastic strain predictions for RVT tests (500 kPa confinement)
Figure 153. Illustration. Variation of yield stress (Standard Linear Solid model)
Figure 154. Graph. Stress histories for rest period analysis
Figure 155. Graph. Yield stress versus cumulative loading time (rest period analysis)
Figure 156. Graph. Viscoplastic strain versus cumulative loading time (rest period analysis)
Figure 157. Graph. Stress history for loading time analysis
Figure 158. Graph. Viscoplastic strain versus cumulative loading time (loading time analysis)
Figure 159. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VT)
Figure 160. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VL)
Figure 161. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement VT)
Figure 162. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement VL)
Figure 163. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VT)
Figure 164. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VLT)
Figure 165. Graph. Viscoplastic strain versus cumulative loading time (140 kPa confinement VT)
Figure 172. Graph. Viscoplastic strain versus cumulative loading time (500 kPa confinement VLT)
Figure 174. Graph. Layout of the numerical experiment specimen
Figure 177. Graph. Evolution of damage parameter, S, for the monotonic test
Figure 178. Graph. Plot of function C(S) with time
Figure 179. Diagram. Domain module
Figure 180. Illustration. Infinite elastic layer on a rigid base
Figure 183. Diagram. Analysis module
Figure 184. Screen capture. Main window of the preprocessor
Figure 185. Screen capture. Control panel
Figure 186. Screen capture. Mesh discretization
Figure 187. Screen capture. Zoom operation on a mesh
Figure 188. Diagram. Stages involved in an FEP++ analysis
Figure 189. Screen capture. Sample data entry window
Figure 190. Screen capture. Error dialog for a semantic error
Figure 191. Screen capture. Analysis runtime window
Figure 192. Graph. Temperature variations used for simulations
Figure 193. Illustration. Vertical strains in winter
Figure 194. Illustration. Vertical strains in summer
Figure 195. Illustration. Longitudinal strains in winter
Figure 196. Illustration. Longitudinal strains in summer
Figure 197. Illustration. Transverse strains in winter
Figure 198. Illustration. Transverse strains in summer
Figure 199. Illustration. Vertical strains for Control mixture
Figure 200. Illustration. Vertical strains for SBS mixture
Figure 201. Illustration. Longitudinal strains for Control mixture
Figure 202. Illustration. Longitudinal strains for SBS mixture
Figure 203. Illustration. Transverse strains for Control mixture
Figure 204. Illustration. Transverse strains for SBS mixture
Figure 205. Illustration. Vertical strains for a wheel speed of 13.41 m/s
Figure 206. Illustration. Vertical strains for a wheel speed of 26.82 m/s
Figure 207. Illustration. Longitudinal strains for a wheel speed of 13.41 m/s
Figure 208. Illustration. Longitudinal strains for a wheel speed of 26.82 m/s
Figure 209. Illustration. Transverse strains for a wheel speed of 13.41 m/s
Figure 210. Illustration. Transverse strains for a wheel speed of 26.82 m/s
Figure 211. Screen capture. General Information dialog
Figure 212. Screen capture. Material Properties dialog
Figure 213. Screen capture. Elastic material properties dialog
Figure 214. Screen capture. Viscoelastic material properties dialog
Figure 215. Screen capture. Prony Coefficients dialog
Figure 216. Screen capture. Layer properties dialog
Figure 217. Screen capture. Mesh properties dialog
Figure 218. Screen capture. Load properties dialog
Figure 219. Screen capture. Analysis Parameters dialog
Figure 220. Screen capture. Summary dialog
Figure 221. Screen capture. FEP++ analysis in progress
Table 1. Relevant asphalt binder information
Table 2. Summary of constructed lanes' air void and asphalt content
Table 3. Controlled crosshead testing matrix for Control2006 in tension
Table 4. Controlled crosshead testing matrix for Control in compression
Table 5. Creep and recovery testing matrix for Control mixture in compression
Table 6. Test conditions for the VT and RVT tests
Table 7. Loading times for VLT test
Table 12. Cyclic tests performed
Table 13. Summary of cyclic correction factors
Table 14. Summary of αvalues for refined model
Table 15. Cyclic tests performed for 9.5mm Fine mixture
Table 20. Delft material model coefficients functions
Table 21. Material coefficients used for the developed model analysis
Table 22. Compression viscoplastic material model coefficients
Table 23. Properties of pavement
Table 24. Properties of moving wheel load
2D 
Twodimensional 
3D 
Threedimensional 
ALF 
Accelerated Load Facility 
AMPT 
Asphalt Mixture Performance Tester 
CLT 
Constant loading level and time test 
CRTB 
Crumb Rubber Terminal Blend 
CS 
Controlled stress 
CX 
Controlled crosshead 
EICM 
Enhanced Integrated Climatic Model 
FEM 
Finite element method 
FEP++ 
Finite element package (proper name) 
FHWA 
Federal Highway Administration 
GUI 
Graphical user interface 
HISS 
Hierarchical single surface 
HMA 
Hot mix asphalt 
IDT 
Indirect tension test 
LVDT 
Linear variable displacement transducers 
LVE 
Linear viscoelastic 
MVECD 
Multiaxial viscoelastic continuum damage 
MVEPCD 
Multiaxial viscoelastoplastic continuum damage 
NCHRP 
National Cooperative Highway Research Program 
NCSU 
North Carolina State University 
NMSA 
Nominal maximum size aggregate 
RVT 
Reversed variable loading time test 
SBS 
Styrenebutadienestyrene 
SHRP 
Strategic Highway Research Program 
SPT 
Simple Performance Tester 
tTS 
Timetemperature superposition 
TFHRC 
TurnerFairbank Highway Research Center 
TRS 
Thermorheologically simple 
UVECD 
Uniaxial viscoelastic continuum damage 
UVEPCD 
Uniaxial viscoelastoplastic continuum damage 
VECD 
Viscoelastic continuum damage 
VECDFEP++ 
Viscoelastic continuum damage model included in finite element package 
VEPCD 
Viscoelastoplastic continuum damage 
VP 
Viscoplastic 
VL 
Variable loading level test 
VLT 
Variable loading level and time test 
VT 
Variable loading time test 
VTK 
Visualization Toolkit 
ε^{e}_{if}  elastic strain 
ε^{p}_{if}  plastic strain 
ε^{c}_{if}  creep strain 
ε^{v}_{if}  viscoplastic strain 
ε^{vp}_{if}  elastic plus linear viscoelastic strain due to damage 
ε_{ve}  elastic plus linear viscoelastic strain due to damage 
∫_{R}  reduced physical frequency (Hz) 
ω_{R}  reduced angular frequency (=2π∫RΔt) 
α_{t}  timetemperature shift factor at specific temperature, T 
π  3.141593 
σ  stress 
ε  strain 
ε_{x }  strain in direction x 
με  microstrain 
Ε(t)  relaxation modulus 
D(t)  creep compliance 
τ  dummy integration variable 
Ρ subscript i  relaxation time (fitting coefficient) 
τ_{j}  retardation time (fitting coefficient) 
Ε_{∞}  relaxation modulus at infinite time 
Ε_{i}  Prony coefficient for relaxation modulus 
D_{ε}  glassy compliance 
Dj  Prony coefficient for creep compliance 
Ε'  storage modulus 
 Ε*   dynamic modulus 
Φ  phase angle 
t  time 
ξ  reduced time 
∂  partial derivative 
Δ_{X}  delta, finite difference in X 
∂_{X}  delta, finite difference in X 
ε^{R}  pseudo strain 
Ε_{R}  reference modulus 
C  normalized pseudo secant modulus 
I  normalization parameter 
S  damage 
W  strain energy density 
W _{d}^{R}  pseudo energy density 
α  viscoelastic damage growth rate 
C_{ikl}  general stiffness matrix 
ν  Poisson's ratio 
Ε  elastic elongational modulus 
G  elastic shear modulus 
Z_{ij}  stiffness matrix for transversely isotropic material 
Ε_{3}  modulus along axis of symmetry 
Ε  modulus perpendicular to axis of symmetry 
v_{3132}  Poisson's ratio between axis of symmetry and perpendicular plane 
v_{1323}  Poisson's ratio between perpendicular plane and axis of symmetry 
v_{12}  Poisson's ratio on perpendicular plane 
Υ  Poisson's ratio term 
S_{i}  compliance matrix for transversely isotropic material 
A_{j}  alterative stiffness matrix terms for transversely isotropic material 
e_{v}  volumetric strain (dilation) 
e3  major deviatoric strain 
e_{2}  strain difference 
Ρ  pressure 
C_{11}  first material integrity term 
C_{12}  second material integrity term 
C_{22}  third material integrity term 
υ  dilation 
υ^{R}  pseudo dilation 
Q_{i}  generalized loads 
q_{i}  generalized displacements 
υ_{s}  complementary strain energy 
ε^{vp}  viscoplastic strain rate 
λ  positive scalar factor 
g  plastic potential function (general viscoplasticity) 
Φ  overstress function 
η  viscosity 
Γ  fluidity 
I_{1}  first stress invariant 
J_{2}  second deviatoric stress invariant 
J_{3}  third deviatoric stress invariant 
s_{ij}  deviatoric stress 
δ_{IJ}  Kroneker delta 
Κ  isotropic hardening parameter 
α_{IJ}  viscoplastic kinematic hardening parameter 
g  stress function (strain hardening model specific) 
P_{α}  atmospheric pressure 
R  tensile strength parameter 
n  yield stress shape parameter 
β  parameter determining shape of yield stress in deviatoric stress space 
ε^{in}  inelastic strain rate 
κ  isotropic hardening function 
α  viscoplastic back stress parameter (kinematic hardening) 
D  drag stress 
R  isotropic hardening function 
H  kinematic hardening function 
G  back stress (viscoplastic model) 
ζ  viscoplastic strain trajectory 
θ  angular direction in axisymmetric coordinate system 
α _{1}  timetemperature shift factor function coefficient 1 
α_{2}  timetemperature shift factor function coefficient 2 
α_{3}  timetemperature shift factor function coefficient 3 
θ  bulk stress 
ε^{R}_{pressure}  pseudo strain due to pressure 
υ^{R}_{pressure}  pseudo dilation due to pressure 
ε^{R}e  effective pseudo strain 
ε^{R}s  permanent pseudo strain 
ε^{R}_{m }  total maximum pseudo strain 
ε^{R}_{me}  total effective pseudo strain at peak of loading 
M  time change correction factor 
t_{p}  pulse time 
ζ_{p}  reduced pulse time 
N  number of points in a calculation 
F  pseudo strain slope function 
G  pseudo strain hysterisis function 
H  healing function 
analytical expression for pseudo strain as a function of time  
Q  pseudo strain shape factor 
Z  combined pseudo strain shape and pseudo stiffness time factor 
ε^{R}_{0,ta}  pseudo strain tension amplitude only 
σ_{0,ta}  stress tension amplitude only 
σ_{pp}  peaktopeak stress magnitude 
β  factor quantifying time under tensile loading 
σ_{peak}  maximum value of stress in a cycle 
σ_{valley}  minimum value of stress in a cycle 
ζ_{i }  reduced time within a cycle when tension loading begins 
ζ_{f}  reduced time within a cycle when tension loading ends 
R_{t}  form adjustment factor for characterization 
κ_{l}  form adjustment factor for prediction 
ε_{0,ta }  strain tension amplitude only 
α(^{t}_{p})  phenomenological viscoplastic model slope pulse time function 
D  phenomenological viscoplastic model intercept function 
γ  viscoplastic softening parameter 
I_{l,dilation}  first stress invariant value at the beginning of dilation 
J_{2,dilation}  second deviatoric stress invariant value at the beginning of dilation 
ε_{reduced}  reduced strain rate 
D  viscosity parameter (final viscoplastic model) 
G  orientationdependent isotropic hardening function 
G^{n+1}  value of hardening function at next time step 
G^{n}  value of hardening function at current time step 
ζ_{o}^{n}  value of elastic state variable at current time step 
ζ_{o}^{n+l}  value of elastic state variable at next time step 
ζ_{i}^{n}  value of state variable for element i at current time step 
ζ_{i}^{n+l}  value of state variable for element i at next time step 