VISCOELASTIC INVERSE ANALYSIS OF FWD DATA USING GENETIC
ALGORITHMS
By
Sudhir Varma

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Civil Engineering−Doctor of Philosophy
2015

ABSTRACT
VISCOELASTIC INVERSE ANALYSIS OF FWD DATA USING GENETIC
ALGORITHMS

By
Sudhir Varma
Dynamic modulus (|E*|) master curve is a fundamental material property for an asphalt
pavement. It is also a key input to Pavement-ME, a pavement design and analysis software that
can predict progression of distresses. Falling Weight Deflectometer (FWD) is a nondestructive
test whose results are typically used for backcalculating layer properties of pavements in situ.
Backcalculation of |E*| master curve of an in-service pavement using FWD data can lead to more
accurate estimation of its remaining service life. Flexible pavements are multilayered structures,
with typically viscoelastic asphalt layer followed by unbound/bound layers. Conventionally,
multilayered elastic analysis is performed to obtain response of flexible pavements for design
and inverse analyses, however, assuming asphalt pavement as a linear elastic material is an
oversimplification of its actual behavior. It is well known that the asphalt pavements’ responses
are both rate and temperature dependent. Appropriate characterization of individual layer
properties is crucial for mechanistic analysis of flexible pavements. Hence, although elastic
analysis is computationally efficient and well accepted in the engineering community, the theory
cannot produce the viscoelastic properties of the asphalt concrete (AC) layer. Backcalculation of
the entire |E*| master curve, including the time-temperature superposition shift factor
coefficients, requires more data than the surface deflection time-histories of FWD drops. In
theory, it should be possible to obtain the |E*| master curve, provided that the data contain time
changing response at different temperature levels.

The specific objectives of the study were to (i) develop a layered viscoelastic pavement
response model in the time domain, (ii) investigate whether the current FWD testing protocol
generates data that is sufficient to backcalculate the |E*| curve using such a model, (iii) if needed,
recommend enhancements to the FWD testing protocol to be able to accurately backcalculate the
|E*| curve as well as the unbound material properties of in-service pavements. Two different
approaches have been proposed to obtain comprehensive behavior of asphalt: (i) using series of
FWD deflection time histories at different temperature levels and (ii) using uneven temperature
profile information existing across the thickness of the asphaltic layer during a single or multiple
FWD drops deflection histories. The models presented can consider the unbound granular
material as both linear elastic as well as nonlinear-stress dependent material. Depending on the
assumed unbound granular material property, and known temperature profile, several
viscoelastic flexible pavement models were developed. The developed forward and
backcalculation models for linear viscoelastic AC and elastic unbound layers have been referred
to as LAVA and BACKLAVA, respectively. The developed forward and backcalculation models
for linear viscoelastic AC and nonlinear elastic unbound layers have been termed as LAVAN and
BACKLAVAN respectively in the study. LAVA and BACKLAVA algorithms assume a
constant temperature along the depth of the AC layer. The algorithms were subsequently
modified for temperature profile in the AC layer and have been referred to as LAVAP and
BACKLAVAP. The major recommendations of the work are the estimated set of temperatures
and number of deflection sensors where FWD tests should be conducted, in order to maximize
the portion of the |E*| curve that can be accurately backcalculated. The results indicate that there
exists a range of temperatures at which the FWD response leads to better inverse solutions. A
genetic algorithm based optimization scheme is offered to search for the pavement properties.

Copyright by
SUDHIR VARMA
2015

TO MY PARENTS

v

ACKNOWLEDGEMENTS
I would like to dedicate my sincere gratitude to my adviser Dr. M. Emin Kutay for giving
me his full support and encouragement throughout my study at Michigan State University. His
enthusiasm and love towards research kept me motivated. Without his advice I would have not
reached this point in my life. The valuable discussions with him gave me an insight into the field
which will remain as a lifelong asset with me.
I would also like to thank Dr. Karim Chatti for his untiring efforts in making valuable
suggestions and ideas. I would also like to thank my committee members, Dr. Neeraj Buch and
Dr. Thomas Pence for their valuable time and comments.
Special thanks to Dr. Gilbert Baladi for his support and guidance. I would also like to
thank Dr. Imen Zaabar for valuable discussions on dynamic backcalculation and parallel
computing.
I would also like to thank my colleagues and friends for their support during my stay at
Michigan State University.
Last but not the least I would like to thank my family for their support and
encouragement.
This study was mainly funded by the Federal Highway Administration (FHWA) grant
number DTFH61-11-C-00026. This support is greatly appreciated. I would like to express my
gratitude towards Dr. Nadarajah Sivanesvaran from the FHWA for his valuable comments.

vi

TABLE OF CONTENTS
LIST OF TABLES .......................................................................................................................... x 
LIST OF FIGURES ...................................................................................................................... xii 
CHAPTER 1 ................................................................................................................................... 1 
INTRODUCTION .......................................................................................................................... 1 
1.1 Objectives ........................................................................................................................... 3 
1.2 Outline................................................................................................................................. 3 
CHAPTER 2 ................................................................................................................................... 5 
LITERATURE REVIEW AND BACKGROUND ........................................................................ 5 
2.1 FWD Data Collection, Analysis and Interpretation ............................................................ 5 
2.2 Forward Analysis Methods to Calculate FWD Response................................................... 7 
2.2.1 Boussinesq’s Solution Method .................................................................................. 7 
2.2.2 Multilayer Elastic Theory .......................................................................................... 9 
2.2.3 Finite Element Analysis ........................................................................................... 10 
2.2.4 Forward Analysis Programs ..................................................................................... 11 
2.3 Backcalculation Approaches ............................................................................................. 12 
2.3.1 Traditional Optimization Techniques ...................................................................... 13 
2.3.2 Non-Traditional Optimization Methods .................................................................. 17 
2.4 Modeling Issues in Backcalculation ................................................................................. 22 
2.5 Characterizing Linear Viscoelastic Properties of HMA in Laboratory ............................ 25 
2.5.1 Interconversion of E(t) to |E*| .................................................................................. 29 
2.5.2 Interconversion D(t) to E(t) .................................................................................... 31 
2.6 Motivation and Scope ....................................................................................................... 32 
CHAPTER 3 ................................................................................................................................. 34 
LAYERED VISCOELASTIC PAVEMENT MODEL ................................................................ 34 
3.1 Introduction ....................................................................................................................... 34 
3.2 Layered Viscoelastic (Forward) Algorithm (LAVA) ....................................................... 36 
3.3 Implementation of Temperature Profile in LAVA ........................................................... 40 
3.3.1 Layered Viscoelastic (Forward) Algorithm for Temperature Profile (LAVAP) ..... 40 
3.3.2 Comparison of LAVAP with FEM Software ABAQUS ......................................... 45 
3.4 Layered Viscoelastic Nonlinear (LAVAN) Pavement Model .......................................... 48 
3.4.1 Nonlinear Multilayer Elastic Solutions.................................................................... 49 
3.4.2 Proposed Constitutive Model for Multilayer Pavement Model under Uniaxial
Loading: Combining Linear Viscoelastic AC and Nonlinear Base .................................. 52 
3.4.3 Proposed Generalized Model for Multilayer Pavement Model: Combining Linear
Viscoelastic AC and Nonlinear Base ................................................................................ 55 
3.4.3.1 Applicability of Existing Theories of Nonlinear Viscoelasticity.......................... 55 
3.5 Proposed Model: LAVAN ................................................................................................ 59 
3.5.1 Validation of the LAVAN ....................................................................................... 61 
3.5.2 Comparison of LAVAN with FEM Software ABAQUS ........................................ 62 
vii

3.6 Chapter Summary ............................................................................................................. 67 
CHAPTER 4 ................................................................................................................................. 68 
BACKCALCULATION USING VISCOELASTICITY BASED ALGORITHM....................... 68 
4.1 Introduction ....................................................................................................................... 68 
4.2 Backcalculation of Relaxation Modulus Master Curve Using Series of FWD Tests Run at
Different Temperatures ........................................................................................................... 74 
4.2.1 Sensitivity of E(t) Backcalculation to the use of Data from Different FWD Sensors
........................................................................................................................................... 76 
4.2.2 Effect of Temperature Range of Different FWD Tests on Backcalculation ............ 80 
4.2.3 Normalization of Error Function (Objective Function) to Evaluate Range of
Temperatures..................................................................................................................... 85 
4.2.4 Backcalculation of Viscoelastic Properties using Various Asphalt Mixtures ......... 91 
4.3 Backcalculation of Relaxation Modulus Master Curve using a Single FWD Test and
Known Pavement Temperature Profile ................................................................................... 94 
4.3.1 Linear Viscoelastic Backcalculation using Single Stage Method............................ 95 
4.3.2 Backcalculation of the viscoelastic Properties of the LTPP Sections using a Single
FWD Test with Known Temperature Profile.................................................................... 99 
4.3.3 Backcalculation of Linear Viscoelastic Pavement Properties using Two-Stage
Method ............................................................................................................................ 112 
4.4 Layered Viscoelastic-Nonlinear Backcalculation (BACKLAVAN) Algorithm ............ 117 
4.4.1 Introduction ............................................................................................................ 117 
4.4.2 Verification of the Nonlinear Backcalculation Procedure ..................................... 122 
4.4.3 Validation of BACKLAVAN using Field FWD Data ........................................... 129 
4.5 Chapter Summary ........................................................................................................... 133 
CHAPTER 5 ............................................................................................................................... 135 
EVOLUTION OF DYNAMIC MODULUS OVER TIME ........................................................ 135 
5.1 Introduction ..................................................................................................................... 135 
5.2 Aging Models.................................................................................................................. 136 
5.2.1 GAS Model ............................................................................................................ 136 
5.2.2 Mixture Conditioning of HMA: AASHTO R30 .................................................... 138 
5.3 Effect of Aging on Viscoelastic Properties of HMA ...................................................... 139 
5.3.1 Effect of Variation of c1 ......................................................................................... 140 
5.3.2 Effect of Variation of c2 ......................................................................................... 140 
5.3.2 Effect of Variation of c3 ......................................................................................... 141 
5.3.2 Effect of Variation of c4 ......................................................................................... 142 
5.3.2 Effect of Variation of a1 ......................................................................................... 142 
5.3.2 Effect of Variation of a2 ......................................................................................... 143 
5.4 Proposed Solution ........................................................................................................... 144 
5.4.1 Results and Discussion of Age Hardening in LTPP Test Section ......................... 147 
5.5 Summary ......................................................................................................................... 158 
CHAPTER 6 ............................................................................................................................... 159 
SENSITIVITY ANALYSIS FOR LAYER THICKNESS ......................................................... 159 
6.1 Introduction ..................................................................................................................... 159 
viii

6.2 Sensitivity Analysis for Layer Thickness: 6 Inches AC Layer ....................................... 163 
6.2.1 Sensitivity Analysis for AC Layer Thickness ........................................................ 164 
6.2.2 Sensitivity Analysis for Base Layer Thickness ..................................................... 167 
6.2.3 Sensitivity Analysis for Total Pavement Thickness .............................................. 171 
6.2.4 Summary and Discussion ....................................................................................... 174 
CHAPTER 7 ............................................................................................................................... 177 
SUMMARY AND CONCLUSION ........................................................................................... 177 
APPENDIX ................................................................................................................................. 182 
BIBLOGRAPHY ........................................................................................................................ 192 

ix

LIST OF TABLES

Table 1: List of commonly used pavement response models ....................................................... 11 
Table 2: List of commonly used backcalculation programs (Smith et. al 2010). ......................... 16 
Table 3: LAVA computation times for different numbers of discrete time steps ........................ 39 
Table 4: Pavement properties used in LAVAP validation with ABAQUS. ................................. 41 
Table 5: Pavement section used in LAVAP validation ................................................................ 46 
Table 6: Peak deflections at temperature profile {19-30}oC and at constant 30oC temperature
using LAVA. ................................................................................................................................. 47 
Table 7: Comparison of ABAQUS and LAVAN: Peak surface deflections for different boundary
conditions. ..................................................................................................................................... 64 
Table 8: Comparison of ABAQUS and LAVAN: Percent error (PEpeak) calculated using the
peaks of the deflections and average percent error (PEavg) calculated using the entire time
history. .......................................................................................................................................... 66 
Table 9: Upper and lower limit values in backcalculation ........................................................... 73 
Table 10: Pavement Properties in viscoelastic backcalculation of optimal number of sensors.... 76 
Table 11: Backcalculation run time for ga-fminsearch seed Runs ............................................... 82 
Table 12: Details of the pavement properties used in single FWD test backcalculation at known
temperature profile ........................................................................................................................ 96 
Table 13: List of LTPP sections used in the analysis ................................................................. 101 
Table 14: Structural properties of the LTPP sections used in the analysis ................................. 102 
Table 15: Depths of stiff layer in each LTPP section estimated using Ullidtz’s method ........... 103 
Table 16: Elastic backcalculation results for LTPP unbound layers .......................................... 103 
Table 17: Viscoelastic backcalculation results for LTPP unbound layers .................................. 104 
Table 18: Variables in two stage linear viscoelastic backcalculation analysis ........................... 113 
Table 19: Pavement properties in two stage linear viscoelastic backcalculation analysis ......... 113 
x

Table 20: Pavement properties and test inputs in staged nonlinear viscoelastic backcalculation.
..................................................................................................................................................... 121 
Table 21: Typical FWD test load levels ..................................................................................... 122 
Table 22: Pavement geometric and material properties in two stage nonlinear viscoelastic
backcalculation. .......................................................................................................................... 123 
Table 23: FWD test data from LTPP section 01-0101 for years 2004-2005. ............................. 130 
Table 24: Nonlinear elastic backcalculation results for section 0101......................................... 130 
Table 25: Comparison of viscoelastic nonlinear viscoelastic backcalculated results obtained at
different stages. ........................................................................................................................... 132 
Table 26: Backcalculation performance in predicting aging of LTPP section. .......................... 154 
Table 27: Pavement properties in sensitivity analysis ................................................................ 164 
Table 28: Error in backcalculated relaxation modulus in structure-A ........................................ 175 
Table 29: Error in backcalculated relaxation modulus in structure-B ........................................ 175 
Table 30: Error in backcalculated base modulus in structure-A ................................................. 176 
Table 31: Error in backcalculated base modulus in structure-B ................................................. 176 

xi

LIST OF FIGURES
Figure 1: Half space notation for Boussinesq’s solution ................................................................ 8 
Figure 2: Typical steps in iterative backcalculation ..................................................................... 14 
Figure 3: Schematic figure of a neuron and ANN network .......................................................... 17 
Figure 4: Typical steps in GA ....................................................................................................... 20 
Figure 5: Figure illustrating GA operations (a) Crossover and (b) Mutation ............................... 20 
Figure 6: Figure illustrating steps in developing |E*| master curve using laboratory measured |E*|
test data :(a) raw |E*| data, (b) master curve development (c) time-temperature shift factors ..... 29 
Figure 7: Generalized Maxwell model ......................................................................................... 30 
Figure 8: Typical flexible pavement geometry for analysis ......................................................... 37 
Figure 9: Discretization of stress history in forward analysis....................................................... 37 
Figure 10: Discretization of the relaxation modulus master curve ............................................... 38 
Figure 11: Deflections calculated under unit stress for points at different distances from the
centerline of the circular load at the surface. ................................................................................ 39 
Figure 12: Schematic diagram of temperature profile .................................................................. 40 
Figure 13: Comparison of response calculated using LAVAP and original LAVA..................... 42 
Figure 14: Comparison of responses calculated using LAVAP at temperature profile {40-3020}oC and original LAVA at constant 40oC temperature. ............................................................ 43 
Figure 15: Comparison of responses calculated using LAVAP at temperature profile {40-3020}oC and original LAVA at constant 30oC temperature. ............................................................ 43 
Figure 16: Comparison of responses calculated using LAVAP at temperature profile {40-3020}oC and original LAVA at constant 20oC temperature. ............................................................ 44 
Figure 17: Region of E(t) master curve (at 19oC reference temperature) used by LAVAP for
calculating response at temperature profile {40-30-20}oC ........................................................... 44 
Figure 18: Relaxation modulus and shift factor master curves at 19oC reference temperature .... 45 
Figure 19: Comparison between LAVAP and ABAQUS at a temperature profile of {19-30}oC
xii

(Terpolymer) ................................................................................................................................. 46 
Figure 20: Comparison between LAVAP and ABAQUS at a temperature profile of {19-30}oC
(SBS 64-40) .................................................................................................................................. 47 
Figure 21: Cross section of multilayer viscoelastic nonlinear system used in uniaxial analysis.. 52 
Figure 22: Flexible pavement cross section. ................................................................................. 58 
Figure 23: Variation of g(σ) with stress and E(t) of AC layer. ..................................................... 58 
Figure 24: Relaxation modulus master curve for LAVAN validation (at 19oC reference
temperature). ................................................................................................................................. 61 
Figure 25: Surface deflection comparison of ABAQUS and LAVAN for the Control mix (*AS1=
ABAQUS sensor 1; *LS1=LAVAN sensor 1). ............................................................................ 66 
Figure 26: Surface deflection comparison of ABAQUS and LAVAN for the CRTB mix (*AS1=
ABAQUS sensor 1; *LS1=LAVAN sensor 1). ............................................................................ 67 
Figure 27: Figure illustrating steps in E(t) master curve development. ........................................ 69 
Figure 28: Backcalculation flow chart in BACKLAVA............................................................... 73 
Figure 29: Schematic of fitness evaluation in BACKLAVA........................................................ 75 
Figure 30: Error in unbound layer modulus in optimal number of sensor analysis ...................... 77 
Figure 31: Backcalculation results using only sensor 1: (a) Backcalculated and actual E(t) master
curve at the reference temperature of 19oC (b) Variation of error,  AC (ti ) . ................................ 79 
Figure 32: Error in unbound layer modulus using FWD data only from farther sensors ............. 80 
Figure 33: Variation of error in backcalculated unbound layer moduli when FWD data run at
different sets of pavement temperatures are used. ........................................................................ 81 
Figure 34: Error in backcalculated E(t) curve in optimal backcalculation temperature set analysis
minimizing percent error............................................................................................................... 82 
Figure 35: Results for backcalculation at {10, 30}oC temperature set: (a) and (b) Only GA is
used, (c) and (d): GA+fminsearch used. ....................................................................................... 83 
Figure 36: Results for backcalculation at {30, 40}oC temperature set ......................................... 84 
Figure 37: Results for backcalculation at {30, 40, 50}oC temperature set ................................... 85 

xiii

Figure 38: Search domain for unconstrained constrain ................................................................ 87 
Figure 39: Search domain for constrained constrain .................................................................... 88 
Figure 40: Backcalculated lE*l master curve using FWD data at temperature set {10, 30}oC
minimizing normalized error ........................................................................................................ 88 
Figure 41: Backcalculated E(t) master curve using FWD data at temperature set {10, 20, 30}oC
minimizing normalized error ........................................................................................................ 89 
avg
Figure 42: Variation of AC
at different FWD temperature sets using modified sigmoidal
variables. ....................................................................................................................................... 90 

Figure 43: Variation of  unbound at different FWD temperature sets using modified sigmoidal
variable.......................................................................................................................................... 90 
Figure 44: Backcalculation results obtained using modified sigmoid variables........................... 91 
Figure 45: Viscoelastic properties of field mix in optimal temperature analysis: (a) Relaxation
modulus at 19oC, (b) Time-temperature shift factor. .................................................................... 92 
Figure 46: Variation of error,  AC (ti ) : (a) ti = 10-5 to ti = 1 sec used in  AC (ti ) , (b) ti = 10-5 to ti
= 10+2 sec used in  AC (ti ) and (c) ti =10-5 to ti = 10+3 sec used in  AC (ti ) computation. ............ 93 
Figure 47: Schematic of fitness evaluation in BACKLAVA........................................................ 95 
Figure 48: Comparison of actual and backcalculated values in backcalculation using temperature
profile. ........................................................................................................................................... 98 
Figure 49: Error in backcalculated E(t) curve for a three-temperature profile ............................. 99 
Figure 50: Time lag in LTPP section 06A805 (shifted based on constant displacement) .......... 101 
Figure 51: LTPP section 01-0101 cross section and temperature profile in AC layer ............... 106 
Figure 52: Comparison of measured and forward calculated deflection histories for section
010101 (a) Viscoelastic backcalculation (b) Forward calculated using E(t) obtained from LTPP
measured D(t) and elastic backcalculated modulus values. ........................................................ 106 
Figure 53: Comparison of measured and backcalculated E(t), aT(T) and |E*| at section 010101.
..................................................................................................................................................... 107 
Figure 54: LTPP section 340802 cross section and temperature profile in AC layer................. 108 
Figure 55: Comparison of measured and forward calculated deflection histories for section
xiv

340802 (a) Viscoelastic backcalculation (b) Forward calculated using E(t) obtained from LTPP
measured D(t) and elastic backcalculated modulus values. ........................................................ 108 
Figure 56: Comparison of measured and backcalculated E(t), aT(T) and |E*| at section 340802.
..................................................................................................................................................... 109 
Figure 57: LTPP section 460804 cross section and temperature profile in AC layer................. 110 
Figure 58: Comparison of measured and forward calculated deflection histories for section
460804 (a) Viscoelastic backcalculation (b) Forward calculated using E(t) obtained from LTPP
measured D(t) and elastic backcalculated modulus values. ........................................................ 111 
Figure 59: Comparison of measured and backcalculated E(t), aT(T) and |E*| at section 460804.
..................................................................................................................................................... 111 
Figure 60: Elastic backcalculation of two-step temperature profile FWD data, assuming AC as a
single layer in two stage backcalculation.................................................................................... 113 
Figure 61: Elastic backcalculation of three-step temperature profile FWD data, assuming AC` as
a single layer in two stage backcalculation ................................................................................. 114 
Figure 62: Elastic backcalculation of two-step temperature profile FWD data, assuming two AC
sublayers in two stage backcalculation ....................................................................................... 115 
Figure 63: Elastic backcalculation of three-step temperature profile FWD data, assuming three
AC sublayers in two stage backcalculation ................................................................................ 115 
Figure 64: Error in backcalculated E(t) curve from two-step temperature profile. .................... 116 
Figure 65: Error in backcalculated E(t) curve from three-step temperature profile. .................. 117 
Figure 66: Nonlinear elastic backcalculated AC modulus (in stage-1) for Control and CRTB
Mixes, using FWD data at different test temperatures. .............................................................. 124 
Figure 67: Nonlinear elastic backcalculated unbound layer properties (in stage-1) for Control and
CRTB Mixes, using FWD data at different test temperature...................................................... 125 
Figure 68: (a) Comparison of backcalculated and actual E(t) master curves (b) average
avg
percentage error,	  AC
for the Control mixture and (c) time-temperature shift factor................ 127 
Figure 69: (a) Comparison of backcalculated and actual E(t) master curves (b) average
avg
percentage error,	  AC
for the CRTB mixture and (c) time-temperature shift factor. ................ 128 
Figure 70: Backcalculated and measured FWD deflection time histories for LTPP section 0101:
(a) test run in 2004 and (b) test run in 2005. ............................................................................... 131 

xv

Figure 71: Nonlinear viscoelastic backcalculation of section 0101 (a) Relaxation modulus (b)
Shift factor at 19oC. .................................................................................................................... 132 
Figure 72: Effect of variation in c1 on |E*| master curve ............................................................ 140 
Figure 73: Effect of variation in c2 on |E*| master curve ............................................................ 141 
Figure 74: Effect of variation in c3 on |E*| master curve ............................................................ 142 
Figure 75: Effect of variation in c4 on |E*| master curve ............................................................ 143 
Figure 76: Effect of variation in a1 on aT(T) master curve ......................................................... 143 
Figure 77: Effect of variation in a2 on aT(T) master curve ......................................................... 144 
Figure 78: Backcalculated |E*| and aT(T) at section 01-0101 for year 1993, 1996 and 2000. ... 148 
Figure 79: Backcalculated |E*| and aT(T) at section 34-0802 for year 1993 and 1998. ............. 148 
Figure 80: Backcalculated |E*| and aT(T) at section 46-0802 for year 1993, 1995, 1999, 2004 and
2011............................................................................................................................................. 149 
Figure 81: Comparison of static and viscoelastic backcalculation (|E*| at 10 Hz) for section 010101............................................................................................................................................. 151 
Figure 82: Comparison of static and viscoelastic backcalculation (|E*| at 10 Hz) for section 340802............................................................................................................................................. 152 
Figure 83: Comparison of static and viscoelastic backcalculation (|E*| at 10 Hz) for section 460802............................................................................................................................................. 153 
Figure 84: Normalized |E*| and aT(T) coefficients with ageing time for section 01-0101 ......... 155 
Figure 85: Normalized |E*| and aT(T) coefficients with aging time for section 34-0802 ........... 156 
Figure 86: Normalized |E*| and aT(T) coefficients with aging time for section 46-0802 ........... 157 
Figure 87: Relaxation modulus master curves for the sensivity analysis ................................... 163 
Figure 88: Error in backcalculated relaxation modulus from variation in 6 inch AC thickness. 165 
Figure 89: Error in backcalculated relaxation modulus from variation in 12 inch AC thickness.
..................................................................................................................................................... 165 
Figure 90: Backcalculated relaxation modulus in 6 inch AC structure: Case-1 sensitivity analysis.
..................................................................................................................................................... 166 
xvi

Figure 91: Backcalculated relaxation modulus in 12 inch AC structure: Case-1 sensitivity
analysis........................................................................................................................................ 167 
Figure 92: Error in backcalculated relaxation modulus from variation in base thickness for 6 inch
AC structure. ............................................................................................................................... 168 
Figure 93: Error in backcalculated relaxation modulus from variation in base thickness for 12
inch AC structure. ....................................................................................................................... 168 
Figure 94: Backcalculated relaxation modulus in 6 inch AC structure, case-2 sensitivity analysis.
..................................................................................................................................................... 169 
Figure 95: Backcalculated relaxation modulus in 12 inch AC structure, case-2 sensitivity
analysis........................................................................................................................................ 170 
Figure 96: Error in backcalculated relaxation modulus in 6 inch AC structure, case-3 sensitivity
analysis........................................................................................................................................ 172 
Figure 97: Error in backcalculated relaxation modulus in 12 inch AC structure, case-3 sensitivity
analysis........................................................................................................................................ 172 
Figure 98: Backcalculated relaxation modulus in 6 inch AC structure, case-3 sensitivity analysis.
..................................................................................................................................................... 173 
Figure 99: Backcalculated relaxation modulus in 12 inch AC structure, case-3 sensitivity
analysis........................................................................................................................................ 174 
Figure 100: Backcalculated relaxation modulus for structure-A: Case-1 sensitivity analysis. .. 183 
Figure 101: Backcalculated relaxation modulus for structure-B: Case-1 sensitivity analysis.... 185 
Figure 102: Backcalculated relaxation modulus for structure-A: case-2 sensitivity analysis. ... 186 
Figure 103: Backcalculated relaxation modulus for structure-B: case-2 sensitivity analysis. ... 188 
Figure 104: Backcalculated relaxation modulus for structure-A, case-3 sensitivity analysis. ... 189 
Figure 105: Backcalculated relaxation modulus for structure-B, case-3 sensitivity analysis..... 191 

xvii

CHAPTER 1

INTRODUCTION
Flexible pavements are multilayered structures with typically viscoelastic asphalt
concrete (AC) as top layer, followed by unbound/bound granular layers. Combined response of
linear viscoelastic and elastic materials that are in perfect bonding is linear viscoelastic.
Assuming there is full bonding between the asphalt layer and the underlying base and subgrade
layers, the overall response of the entire pavement system becomes viscoelastic. Recent
developments in pavement design and rehabilitation drift the traditional empirical analysis of
pavement to a more mechanistic framework. The MEPDG design guide (MEPDG, 2004)
developed under National Cooperative Highway Research Program 1-37A, using mechanisticempirical framework tries to address various shortcomings of previous empirical pavement
design guides (AASHTO, 1993). MEPDG is based on the notion that actual material and loading
condition along with locally calibrated mechanistic models would yield accurate response. Long
Term Pavement Performance (LTPP) is aimed to develop such a database, which can be used to
calibrate MEPDG models. MEPDG uses fundamental properties of pavement materials and
acknowledge that the pavement distresses are the consequence of collective factors such as
environmental, material and loading conditions. Out of the various contributing factors,
temperature distribution in pavement has been given belligerent attention in MEPDG and is
regarded as one of the most critical. One of the significant accomplishments of MEPDG is the
recognition of AC layer as viscoelastic material. The viscoelastic properties vary with time
(frequency) and temperature; and hence require both time and temperature functions to define it.
1

The characteristic mechanistic properties of an isotropic-thermorheologically simple
viscoelastic system are the relaxation modulus E(t) and the creep compliance D(t), which are
function of time, the dynamic modulus |E*|, which is a function of frequency, and the timetemperature shift factors (aT(T)), which is a function of temperature. These characteristic
properties are often expressed at a specific reference temperature, in terms of a ‘master curve’.
It can be shown that if any of the three properties E(t), D(t) or |E*| is known, the other
two can be obtained through an inter-conversion method such as the Prony’s series approach
(Park and Schapery 1999). |E*| is used as the fundamental material property input in MEPDG to
incorporate the viscoelastic properties of Hot Mixed Asphalt (HMA) mixes and for the
temperature function, the time-temperature superposition principal is exploited. Being an
important fundamental material property, knowledge of |E*| master curve of an in-service
pavement using falling weight deflectometer (FWD) data can lead to more accurate estimation of
its remaining life and rehabilitation design. Hence, the viscoelastic properties proposed to be
backcalculated from FWD data include two functions - a time function and a temperature
function. In the present study, the time function refers to the relaxation modulus master curve
,

in which t is physical time and TR is the corresponding (constant) reference temperature.

The temperature function refers to the time-temperature shift factor

, which is a positive

definite dimensionless scalar. Furthermore, AC has been assumed to be thermorheologically
,

simple, which allows applying

for any temperature level

replacing physical time with a reduced time
such that

1 if

; therefore,

.

2

(different than

) by simply

is a function of both

and

,

1.1 Objectives
Appropriate characterization of individual layer properties is crucial for mechanical
analysis of flexible pavements. For meaningful interpretation of FWD data it is important that all
the factors that influence the test are considered in the analysis. The specific objectives of the
study were:
(i)

Develop a layered viscoelastic flexible pavement response model in time domain
which can consider variation in temperature with depth (temperature profile) in AC
layer.

(ii)

Develop a layered viscoelastic flexible pavement model in time domain which can
consider stress sensitive nonlinear unbound layers.

(iii)

Develop robust backcalculation algorithms based on GA which can calculate
viscoelastic properties of AC layer and elastic properties of unbound layers.

(iv)

Investigate whether the current FWD testing protocol generates data that are
sufficient to backcalculate the |E*| master curve using such models.

(v)

If needed, recommend enhancement to the FWD testing protocol (testing temperature,
number of sensors) to be able to accurately backcalculate the |E*| master curve as
well as the unbound material properties of in-service pavements.

1.2 Outline
This dissertation consists of seven chapters including the present one. Brief description of
each chapter is provided below:


The second chapter presents a literature review and an overview of the problem
background.

3



In the third chapter, several algorithms have been developed to implement multilayer
viscoelastic forward solution. The models have been validated using Finite Element
Method (FEM) based solutions.



In chapter four, the forward models developed in chapter three were used to develop
genetic algorithm (GA) based backcalculation schemes for determining E(t) or |E*|
master curve and time-temperature shift factors of AC layers and unbound material
properties of in-service pavements. Further, the effect of FWD test temperatures and
number of surface deflection sensors on backcalculation of |E*| master curve were
studied. Suitable FWD test data requirements have been discussed in the key findings
from the study.



In chapter five, the backcalculation algorithms developed in chapter four were used to
develop a procedure to quantify evolution of |E*| in in-service flexible pavements.



Chapter six presents a sensitivity analysis on viscoelastic backcalculation.



Summary and conclusions are presented in chapter seven.

4

CHAPTER 2

LITERATURE REVIEW AND BACKGROUND
2.1 FWD Data Collection, Analysis and Interpretation
The FWD test involves dropping a load and then measuring surface deflections through
sensors placed at certain pre-designated distances. FWD field data collection usually involves
collecting stress history, deflection history, sensor location and surface temperature. During the
FWD testing, load is released from a given height onto a load plate, where the stress is recorded
through a load cell. The stress and deflection histories are collected at a time step of 0.1-0.2
msec, with the deflection sensors placed at radial distance of 0, 8, 12, 18, 24, 36, 48, 60 and 72
inches. The primary issue with FWD data collection (stress and deflection history) is the error
associated with data acquisition system. This issue is particularly relevant to time-based analysis
of FWD. In time-based analysis, both the load and deflections are assumed to be coincident,
which means any synchronization issue between load and deflection can seriously affect the
results.
Most of the FWD devices are either trailer-mounted or vehicle-mounted. In a trailermounted system (such as KUAB, Dynatest 8000), the FWD device is fixed on a trailer which is
mounted to a towing vehicle, whereas in a vehicle-mounted system, the FWD device is directly
fixed into a van. These vehicles are loaded with infrared temperature gauge to collect pavement
surface and air temperatures. The primary issue with the collected surface and air temperature is
that they may not represent the actual temperature in the asphalt layer. Variation in temperature
along the depth of AC layer may significantly influence flexible pavement response. To negate
5

the effect of temperature FWD deflections are often corrected for temperatures. Temperature
correction in FWD analysis is generally done in two steps. In step one, temperature in AC is
predicted and in step two, either the FWD deflections or backcalculated modulus values are
corrected using the predicted temperature. These temperature prediction equations (Park et al.,
2002; Park et al., 2001 and Lukanen et al., 2000) and modulus correction equations (Chen et al.,
2000; Park et al., 2001 and Lukanen et al., 2000) are developed empirically using data collected
from limited test regions. Further the models are calibrated to predict temperatures at specific
depths (generally the mid depth) and reported to predict temperatures with accuracy of ±2oC to
5oC. In order to predict any profile, temperatures at minimum two depths are needed. Further,
temperature profile in field is typically small and (as discussed later in chapter 3 and chapter 4)
this accuracy may not be sufficient for predicting actual property (i.e. |E*|) of AC layer. During
LTPP FWD field tests, in addition to the surface temperature measurement, the temperature
profile along the depth of AC layer is also collected. The temperature measurements are recorded
at depths, 0mm, 25mm, 50mm, 10mm, 200mm, and 300mm (Schmalzer, 2006). In situ
measurement of temperature profile eliminates the errors associated with temperature prediction
and provides reliable data for mechanistic analysis.
FWD testing is commonly employed in project level analysis for assisting overlay
designs. At network level, FWD testing is done to understand structural capacity of the network
at large. Since the inception of FWD, several methods have been developed for the interpretation
of FWD data. These methods are based on the understanding of pavement response theories, and
techniques to interpret it in a backcalculation scheme. In general, the surface deflection at given
radial offsets are calculated using a suitable response theory. The calculated deflections are then

6

compared with the measured deflections. Various pavement response theories and
backcalculation process available in literature are discussed in the following section.
2.2 Forward Analysis Methods to Calculate FWD Response
In forward analysis, pavement surface deflections are calculated at specific radial
distances. Some of the commonly used mechanistic pavement response theories are


Boussinesq’s solution method



Multilayered elastic theory



Finite element model

2.2.1 Boussinesq’s Solution Method
Bousssinesq developed closed form solutions to compute deflection, strain and stresses in
a homogeneous, isotropic, linear elastic half-space. The solution was developed for an axisymmetric structure under a vertical point load. Vertical deflection at any point A in half-space
can be calculated using equation:

	

2 1

where,

is vertical deflection at any point A, P is point load,

of elasticity of the half space, and R and

and

	

	

	1	

is Poission’s ratio, E is modulus

are geometric distance and angle shown in the Figure

1. Equation 1 implies that surface deflection
be calculated substituting

	

at any radial distance

from the point load can

90
2 1
7

		

	

	

	

	

	2	

where,

is the radial distance between load and point of evaluation. Although the formulation

given in Equation 2 is an easy to implement close-form solution, pavements in reality are not
half spaced and are not subjected to point load. Pavements are multilayered structure with
different material properties in each layer.

Load P
x

R
y

θ 

z
A

z
Figure 1: Half space notation for Boussinesq’s solution
Odemark (Ullidtz, 1988) developed a method to transform multilayer system with
different moduli into a single modulus system that makes Boussinesq’s equations applicable. The
method is based on the assumption that the stiffness of the transformed layers is equivalent to the
stiffness of the original (untransformed) layers. Theoretically, the layer transformation can be
applied to any number of layers using the transformation relationship:

		

∑

where,

	

	

	

is is the equivalent thickness of first n-1 layers, C is correction factor,

	

	3	

and

are the

thickness and modulus of ith layer. The method is commonly known as method of equivalent
thickness (MET). ELMOD3 and BOUSDEF are some of the programs that use MET as forward
8

solution in their backcalculation algorithms. The method is an approximation and requires
appropriate correction factor for predicting pavement response.
2.2.2 Multilayer Elastic Theory
Multilayer elastic theories are much closer to the pavement system and are most
commonly used in pavement analysis. Solution for 2-layer and 3-layer system was first
developed by Burmister (1943, 1945). The theory has been extend to multilayer system by
various researchers (Schiffman, 1962 and Michelow, 1963) and programed in several multilayer
analysis software.
The basic assumptions of the layered elastic solution are:


The pavement system is axisymmetric multilayer structure.



A circular uniformly distributed load is applied normally on the surface of the pavement.



All the layers consist of homogeneous, isotropic and linearly elastic material and follow
Hooke’s law.



All the layers are of homogeneous thickness and extend horizontally to infinity.



The bottom layer is a semi-infinite half-space.

Basic equations in multilayer formulation are:


Equilibrium equations



Compatibility equations

9

A stress function

,

0 can be

that satisfies the governing differential equation

assumed for each layer to derive stresses and displacement (Love, 1923). The vertical
displacement in each layer interms of stress function can be expressed as:

	

2 1

where,

is vertical displacement,

,

is Airy’s stress function,

	

	

	

is Poission’s Ratio and

	4	
is

modulus of elasticity. The main advantage of multilayer solution over the Boussinesq’s solution
is it’s ability to consider distributed surface loading and multiple layers with different properties.
However, the solution cannot be used for nonlinear materials or viscoelastic materials.
EVERCALC, MICHBACK, MODULUS and CHEVDEF are some of the programs that use
linear layered elastic based forward solution in their backcalculation process.
2.2.3 Finite Element Analysis
Finite element analysis (FEA) offers much more flexibility in selecting material
constitutive equations, loading condition and geometric variation. Some of the major advantages
of FEA in pavement analysis are:


Ability to consider nonlinearity in both vertical and horizontal direction



Ability to consider the viscoelastic properties of HMA in analysis



Ability to consider dynamics in analysis.
Researchers have mainly used FEA in pavement analysis using general purpose FEA

software such as ABAQUS and ANSYS. General purpose software are not pavement-specific, as
an example, stress nonlinearity in unbound pavement material are not predefined and need to be
10

defined as a user defined material (UMAT). Researchers have developed programs
(MICHPAVE, CAPA-3D) that have inbuilt pavement-specific models to analyses unbound
nonlinearity, dynamics and viscoelasticity. However, they are not designed for a backcalculation
type analysis and usually slow.
2.2.4 Forward Analysis Programs
With the advancement of computing facilities, several pavement analysis softwares have
been developed. Flexible pavements are multilayered structures, with typically viscoelastic
asphalt layer followed by nonlinear unbound/bound layers. Conventionally, multilayered elastic
analysis is performed to obtain response of flexible pavements for design and inverse analyses,
however, assuming asphalt pavement as a linear elastic material is an oversimplification of its
actual behavior. It is well known that the asphalt pavements’ responses are both rate and
temperature dependent. A list of common pavement analysis computer programs and their
capability to account for nonlinear and viscoelastic behavior of pavement is presented in Table 1.
Table 1: List of commonly used pavement response models
Program Name
BISAR
CAPA-3D
CHEVRON
ELSYMS
Everstress
KENLAYER
MICHPAVE

Developer
Shell Oil Co.
Delf University
of Technology
Chevron
FHWA
University of
Washington
Yang H. Huang
Michigan State
University

Response Model
Layered elastic analysis

Nonlinearity
No

Viscoelasticity
No

3D-FEM

Yes

Yes

Layered elastic analysis
Layered elastic analysis

No
No

No
No

Layered elastic analysis

Yes

No

Layered elastic analysis

Yes

Yes

Axi-symmetric FEM

Yes

No

Most of the softwares for pavement analysis are based on layered elastic theories and do
not consider nonlinearity or viscoelasticity. The FEM based softwares CAPA-3D and
MICHPAVE can consider nonlinearity. CAPA-3D is also capable of accounting for nonlinearity
11

as well as viscoelastic behavior. An approximate nonlinear analysis of pavements can also be
performed using multilayered elastic based solution. Among the layered elastic analysis based
softwares, Everstress and KENLAYER (Huang, 2004) account for nonlinearity through
iteratively adjusting layer moduli. However, since the multilayer elastic theory assumes each
individual layer to be vertically as well as horizontally homogeneous, it can be used to depict
nonlinearity only through approximation. For incorporating variation in modulus with depth,
Huang (Huang, 2004) suggested dividing the nonlinear layer into multiple sublayers.
Furthermore, he suggested choosing a specific location in the nonlinear layers to evaluate
modulus based on the stress state of the point. He showed that when the midpoint under the load
is selected to calculate modulus values, the predicated response near the load is close to the
actual response. However, the difference between actual and predicted response increased at
points away from the loading.
2.3 Backcalculation Approaches
Backcalculation analysis is essentially an inverse analysis, which involves determining
the pavement layer properties using measured pavement surface deflections. The backcalculated
results are sensitive to the backcalculation technique used in the analysis (Chou and Lytton,
1991; and Harichandran et al., 1993). Typically backcalculation programs adopt techniques such
as nonlinear regression, iterative methods, close-form solutions and database search to predict
pavement properties. Regression method and close form solution methods are simple and time
efficient, but are mainly used to predict subgrade modulus (Newcomb, 1986; Horak, 1987; and
AASHTO, 1993). Most of the backcalculation programs employ iterative methods to reach
solution. Backcalcultion using iterative methods are essentially optimization problem. The
problem involves iteratively changing the unknowns till an objective function is minimized
12

under a pre-defined tolerance limit. For instance, in elastic backcalculation this objective
function is generally root mean square error or percent difference between the measured and
predicted deflection basin (peak deflection from each deflection sensor) and the unknowns are
the elastic moduli of the layers. The optimization methods can be classified into two major
categories. The first category is based on the traditional optimization techniques which
repeatedly use a forward analysis method in a direct search or gradient based search technique.
The second category is based on non-traditional optimization techniques which mimics
mechanisms observed in nature.
2.3.1 Traditional Optimization Techniques
In traditional optimization direct search (such as simplex, Hooke-Jeeves) or gradient
search (such as secant method, Newton Raphson, modified Levenberg-Marquardt algorithm, and
modified Powell hybrid algorithm) based methods are used. Figure 2 illustrates the typical steps
involved in a traditional optimization based backcalculation procedure.
An example of direct search method is MODULUS, which uses Hooke-Jeeves algorithm.
The MODULUS program is a multilayer elastic based backcalculation program developed by
Texas Transportation Institute. The program generates a database of deflection basin using linear
elastic program WESLEA (Rohde and Scullion, 1990). The number of deflection basins
generated is based on the number of layers and the expected moduli range. The generated
deflection database is then searched. The program uses Hooke-Jeeves pattern search algorithm to
minimize sum of square difference between the measured and calculated deflections. The
program can also calculate depth to stiff layer and perform temperature correction. The program
(MODULUS 6.0) can handle up to four layers and seven deflection sensors. COMDEF

13

(Anderson, 1989) is another example of backcalculation program which generates a database of
deflection basin. It uses CHEVRON as the forward program to generate deflection database.

Knowns

Measured 
deflection

Pavement 
information 
(Thickness, 
Poission’s 
ratio)

Occasional path
Usual path

Start

Seed 
moduli

Range of 
moduli

Calculated 
deflections

New 
moduli

Error 
<tolerance 
limit

No

Yes
End

Figure 2: Typical steps in iterative backcalculation
MICHBACK, EVERCALC and CHEVDEF are some of the examples of gradient search
method based backcalculation programs. EVERCALC iteratively varies elastic moduli and
compute surface deflection using linear elastic forward program WESLEA. The calculated
deflections are then compared with field measured deflections till one of the following
termination criteria is reached: (a) Deflection tolerance: root mean square error between the
calculated and measured deflection basin. (b) Moduli tolerance: percentage difference in the
predicted moduli values between consecutive iterations. (c) Maximum number of iterations. The
14

software uses a Levenberg-Marquardt algorithm to minimize the error. Although the forward
algorithm in EVERCALC is based on linear layer elastic solution, the software offers nonlinear
characterization of unbound layers using linear regression. The software can also provide
temperature correction and depth to stiff layer.

MICHBACK is another layered elastic

backcalculation program, developed by Michigan State University. It uses linear elastic layered
solution program CHVERONX for forward calculation. The program uses modified NewtonRapson method to minimize root mean square error between the measured and calculated
deflection basin. The program can also predict depth to stiff layer in the analysis. A list of
commonly used backcalculation programs is presented in Table 1. Backcalculation using
traditional iterative methods are popular and most commonly used, but they have the following
issues:


All the conventional backcalculation procedures either require a seed value as an input or
select a seed value within the software. These methods have the tendency to move
towards the closest solution in the vicinity of the initial seed value. If the closest solution
is not the global solution, then a local solution is approached.



Convergence time may vary depending on the initial seed value chosen and the
optimization technique used. Often if the selected seed value is not close to the optimum
solution, the backcalculation would require longer time to reach the solution.
The disadvantages of classical methods listed above do not mean that they cannot be used

in backcalculation procedure. In fact the classical methods can be hybridized along with nontraditional optimization techniques in developing more effective backcalculation procedures.
Non-traditional techniques are much more capable of addressing these issues.
15

Table 2: List of commonly used backcalculation programs (Smith et. al 2010).
Program
Name

Developer

BISDEF

USACE-WES

BOUSDEF

Zhou et al.
(Oregon State)

CHEVDEF

USACE-WES

COMDEF

USACE-WES

DBCONPAS

Tia et al.
(Univ. of
Florida)

ELMOD/EL
CON

Ulidtz
(Dynatest)

ELSDEF

Texas A&M,
USACE-WES

Forward
Calculation
Method
Multi-Layer
Elastic Theory
Method of
Equivalent
Thickness
Multi-Layer
Elastic Theory
Multi-Layer
Elastic Theory
Finite Element
Method of
Equivalent
Thickness
Multi-Layer
Elastic Theory

Forward
Calculation
Sunroutine

Backcalculation
Method

Non-Linear
Analysis

Layer
Interface
Analysis

Maximum Number of
Layers

Seed Moduli

Range of
Acceptable
Modulus

Ability
to Fix
Modulus

Convergence
Scheme

Error
Weighting
Function

BISAR

Iterative

No

Variable

Number of deflections;
best for 3 unknowns

Required

Required

Yes

Sum of squares
of absolute error

Yes

MET

Iterative

Yes

Fixed
(rough)

At least 4

Required

Required

CHEVRON

Iterative

No

Number of deflections;
best for 3 unknowns

Required

Required

BISAR

Database

No

3

No

No

FEACONS
III

Database

Yes?

Yes?

No

No

MET

Iterative

Yes
(subgrade
only)

Fixed
(rough)

4 (exclusive of rigid
layer)

No

No

ELSYM5

Iterative

No

Fixed
(rough)

Number of deflections;
best for 3 unknowns

Required

Fixed
(rough)

3

Required

5

Fixed
(rough)
Fixed
(rough)

Sum of percent
errors
Sum of squares
of absolute error

Yes

Various

No

Yes

Relative error of
5 sensors

No

Required

Yes

Sum of squares
of absolute error

Yes

Required

Yes

Sum of relative
squared error

No

Required

Yes

Sum of absolute
error

Yes

Yes

EMOD

PCS/Law

Multi-Layer
Elastic Theory

CHEVRON

Iterative

Yes
(subgrade
only)

EVERCALC

Mahoney et al.

Multi-Layer
Elastic Theory

WESLEA

Iterative

Yes

Variable

FPEDD1

Uddin

Multi-Layer
Elastic Theory

BASINF?

Iterative

Yes

Fixed
(rough)
Fixed
(rough)

4

Required

Required

Yes

Required (4
and more
layers)
Program
Generated

No

ISSEM4

Ullidtz,
Stubstad

Multi-Layer
Elastic Theory

ELSYM5

Iterative

Yes (finite
cylinder
concept)

MICHBACK

Harichandran
et al.

Multi-Layer
Elastic Theory

CHEVRON
X

Newton method

No

Fixed
(rough)

Number of deflections;
best for 3 unknowns

Required

Required

Yes

MODCOMP
5

Irwin,
Szebenyl

Multi-Layer
Elastic Theory

CHEVRON

Iterative

Yes

Fixed
(rough)

2 to 15 layers; maximum
of 5 unknown layers

Required

Required

Yes

MODULUS

Texas
Transportation
Institute

Multi-Layer
Elastic Theory

WESLEA

Database

Yes?

Fixed?

4 unknown plus stiff
layer

Required

Required

Yes

PADAL

Bown et al.

Multi-Layer
Elastic Theory

Iterative

Yes
(subgrade
only)

Fixed?

Required

RPEDD1

Uddin

Iterative

Yes

Fixed?

Program
Generated

WESDEF

USACE-WES

Multi-Layer
Elastic Theory
Multi-Layer
Elastic Theory

BASINR
WESLEA

Iterative

No

Variable

16

5

Required

Relative
deflection error
Sum of relative
squared error
Relative
deflection error
at sensors
Sum of relative
squared error

No

No
Yes

Sum of relative
squared error
No
Required

Yes

Sum of squares
of absolute error

Yes

2.3.2 Non-Traditional Optimization Methods
The difficulties and limitations of traditional techniques have contributed to the
development of non-traditional methods in backcalculation. Limited literature exists on
application of non-traditional methods in pavement backcalculation. Two methods that have
been used by researchers in pavement backcalculation are genetic algorithm (GA) and artificial
neural network (ANN). In ANN the biological neural networks in brain are mimicked by
network of artificial neurons. Figure 3 shows a schematic figure of an artificial neuron and an
artificial neural network.

Weight

Inputs

Activation
Function

Output

Inputs

Output

(a)
(b)

Figure 3: Schematic figure of a neuron and ANN network
Meier and Rix (1994) developed ANN network using synthetic deflection basin for three
layer flexible pavement. They reported that the ANN performed three fold faster than the
conventional gradient search method. The authors (Meier and Rix, (1995)) also developed an
ANN for dynamic backcalculation of flexible pavement properties. Gopalakrishnan et al. (2006)
trained an ANN for predicting nonlinear layer moduli of airfield pavements. The aim of the study
was to develop networks which can perform backcalculation of airfield pavements subjected to
heavy weight deflectometers. To train the ANN, synthetic deflections were generated using finite
element (FE) based software ILLI-PAVE. Khazanovich and Roesler (1997) developed an ANN
17

based backcalculation program DIPLOBACK to backcalculate elastic properties of composite
pavements. The developed ANN was verified using field FWD data. Similar trend was reported
in the backcalculated results using ANN and the results obtained from WESDEF. One of the
main advantages of ANN in backcalculation is its speed. However, ANN needs to be trained
beforehand using synthetic or field data. Two main issues with ANN based backcalculation
models are:
(1) The quality of the backcalculated results depends on the quality of ANN trained.
(2) Generally, to train ANN, synthetic database of deflections is generated by varying
various parameters such as, number of pavement layer, layer thickness, layer moduli
and depth to stiff layer. The developed ANN is expected to give reasonable results
only for the problems which are within the bounds of the cases considered in the
training of the ANN.
Genetic algorithm is another non-traditional method which mimics biological process in
nature (Bremermann, 1958; Holland, 1975; and Goldberg, 1989). GA comes under the class of
evolutionary algorithms in optimization. There are three major categories of evolutionary
algorithm: genetic algorithm (GA), evolutionary stratergies (ES) (Rechenberg, 1965) and
evolutionary programming (EP) (Fogel, 1962; and Fogel et al., 1966). The differences between
the three methods are very slim and mainly lie in the evolution process.
GA was first introduced by Bermermann in 1985. Holland (Holland, 1975) developed the
much more formal mechanism of GA which was implementable in computers. He developed the
mechanics of GA using schemata theory, which was based on the assumption that every
candidate solution can be encoded as a series of binary numbers, so-called “chromosomes”.
18

Once the candidate solutions were transformed into chromosomes (parents), they were operated
to reproduce new candidate solutions.
In general terms, GA performs the following operations: (a) initialization, (b) selection,
(c) generation of offspring and (d) termination. Figure 4 illustrates the typical steps involved in a
GA based optimization procedure.
In initialization, the GA generates a pool of solutions using a subset of the feasible search
space, so-called ‘population’. Each parameter in a solution is encoded as string which is attached
to form a chromosome. It should be noted that the size of the population plays an important role
in the success of the optimization. A large population size would unnecessarily lead to high
computational time, whereas a small population size may lead to premature convergence at a
local solution.
In selection process, each chromosome (solution) is evaluated using an objective function
and the best fitted solutions are selected to create a mating pool. The parents from the mating
pool are randomly selected to generate next generation population (offspring). This process
mainly involves two operators: crossover and mutation.
In crossover, a new solution is formed by exchanging information between two parent
solutions, which is done by swapping a portion of parent chromosomes (Figure 5a). In mutation,
a new solution is formed by random bit-wise flipping a parent chromosome (Figure 5b). The
newly generated population is evaluated using the objective function. This process is repeated
until a termination criterion is reached. Through guided random search from one ‘generation’ to
another, GA minimizes the desired objective function.

19

Start
Encoding

Initial pool of solution
Fitness evaluation & 
selection
Selected Parents

Variation & mutation

New generation
(offsprings)

No

Error <tolerance limit/
No of generations>max no

Yes
End

Figure 4: Typical steps in GA

Crossover
Parent

offspring

1  0  0  1  1  0  0  1  1
0  0  0  1  0  1  1  0  1

Parent

1 0 0 1 1 0 1 0 1
0 0 0 1 0 1 0 1 1

Mutation
         offspring

1  0  0  1  1  0  0  1  1

(a)

1 0 0 1 1 0 1 0 1
(b)

Figure 5: Figure illustrating GA operations (a) Crossover and (b) Mutation
GAs are much more suited for the following type of problems:


Problems which do not have a close-form solution.



Problems which have complex objective function and gradients are difficult to obtain.
20



Problems with multiple constraints and limits on the variables.



Problems involving multiple optimum solutions.



When the search domain involved is very large.
Limited literature exists on application of GA in flexible pavement backcalculation (Fwa

et al., 1997; Reddy et al., 2004; Tsai et al., 2004; Sangghaleh et al., 2013; and Alkasawneh,
2007). Fwa et al. (1997) developed a GA based backcalculation program NUS-GABACK. The
program was compared with four other backcalculation programs, MICHBACK, MODULUS,
EVERCALC and EVERCALC-Alt. For comparison, synthetic deflection basins were generated
using 3 and 4 layer pavement structures. The population size and generation size used in the
study were 60 and 150 respectively. The program was reported to produce more consistent and
accurate backcalculation results compare to the other programs. Tsai et al. (2004) demonstrated
the applicability of GA in flexible pavement backcalculation. Reddy et al. (2004) studied the
effect of GA parameters on performance of GA in backcalculation. A population size of 60 and
generation size of 60 was recommended for 3 and 4 layer pavement structures. Alkasawneh
(2007) developed a GA-based backcalculation model BackGenetic3d, which can consider any
number of layers in backcalculation. In elastic backcalculation, peak deflection measured at 7-8
sensors are used in backcalculation. This restricts the number of unknowns (preferably ≤ 5) in
most of the traditional optimization based backcalculation programs. However, GA based
optimization do not carry any indeterminacy with number of unknown. Although GA is an
attractive optimization technique, if not used appropriately, it may lead to the following issues:


If the objective function is computationally expensive or require large number of
population and generations, then GA may take long time to reach the solution.
21



The parameters in the algorithm (population size, number of generation, crossover and
mutation probabilities) should be selected such that global solution is reached with
minimum computational effort. Optimal parameters for running GA differ from problem
to problem. Often for any problem a trial and error process is required to choose the best
suited parameter values.



Because of the inherent randomness in the optimization process, GAs produce results
which are near-optimum solution.

2.4 Modeling Issues in Backcalculation
Practical modeling of asphalt pavement systems has been traditionally carried out using
layered elastic theory (Shell 1978, Shook el al. 1982, Monismith 2004); and the FWD data
analyzed using elastic inverse analysis, which assumes pavement to be a multilayered elastic
structure (Harichandran et al., 1993; Fwa et al., 1997; Irvin, 1994; and Bush, 1985). Such an
analysis typically involves identifying layer moduli by matching measured peak displacements
with computed displacements obtained at peak load. However, it is well known that the AC is
linear viscoelastic at low strain levels (Kutay et al., 2011; and Levenberg, 2013). Like any
viscoelastic material it shows properties dependent on time (or frequency) as well as
temperature. In an attempt to provide the engineering community with a better mechanistic
framework, there is a current movement towards treating the asphaltic layers, and hence the
entire pavement system, as layered viscoelastic medium. Such a campaign naturally generates
interest for obtaining viscoelastic pavement properties in situ, nondestructively, by way of
inverse analysis (Kutay et al., 2011). For the asphaltic layers, assuming isotropy, constant
Poisson’s ratio, and thermo-rheological simplicity, the sought after viscoelastic properties

22

include the relaxation modulus E(t), the dynamic modulus |E*|, and the time-temperature shifting
a T (T ) .

In fact, the impulse loading in the FWD test makes it a dynamic test (Uzan, 1994;
Foinquinos et al., 1995; and Roesset et al., 1995). In dynamic backcalculation methods, typically,
damped elastic or elasto-dynamic system is analyzed using finite layer or FE based forward
computations (Uzan, 1994; Al-Khoury et al., 2001; Al-Khoury et al., 2002; and Losa, 2002),
which can be performed either in the frequency domain or in the time domain (Uzan, 1994). In
frequency domain backcalculation, the applied load and deflection response histories are
transformed into the complex domain at various frequencies. This deflection response is then
matched with the deflection obtained from optimized complex moduli of pavement layers. The
aforementioned frequency domain based backcalculation procedure is known to be sensitive to
deflection truncation (Chatti, 2004, Chatti et al., 2006), which is commonly applied about 60
milliseconds after the loading pulse. Recently Lee (2013) developed a time domain multilayer
dynamic forward analysis program, ViscoWave. The solution used in the program is based on
continuous Laplace and Hankel Transforms. Zaabar et al. 2014 implemented the ViscoWave
solution in C++ and used parallel processing to develop a backcalculation program
DYNABACK-VE. Although dynamic backcalculation can consider viscoelastic properties of the
AC layer, it is computationally expensive and difficult to carry out by pavement engineers.
Contribution of dynamics in the FWD test depends mainly on the presence and depth of a
stiff layer (Foinquinos et al., 1995; Roesset et al., 1995 and Uzan, 1994). Foinquins et al. (1995)
compared the static and dynamic response of pavements with various stiff layer depths. They
showed that the difference between the two becomes smaller as the depth of the stiff layer
increases. Lei (2011) showed that the dynamic effects in a multilayered viscoelastic structure are
23

negligible when the stiff layer is located at a depth of more than 5.5 m (18ft) from the pavement
surface. Usually, backcalculation techniques are based on elasto-static analysis; expected to give
solutions comparable to dynamic solutions when there is no stiff layer. Although elastic
backcalculation is simple and computationally efficient, it lacks the inherent ability to reproduce
the observed (time-temperature dependent) system behavior and infer the time and temperature
dependent layer properties.
Kutay et al. (2011) developed a computationally efficient inverse technique based on
layered viscoelastic forward analysis to backcalculate viscoelastic layer properties.
Backcalculation of |E*| from FWD field data needs both load and corresponding deflection time
histories. The peak deflections alone may not be sufficient to provide enough information about
the AC material property (Kutay et al., 2011). The method was used to backcalculate the
relaxation modulus E(t) of an AC layer using FWD time history from a single drop under a
single temperature level (temporally and spatially constant). In the study, a sigmoidal shape was
a-priori assumed for E(t) and the backcalculated relaxation modulus was subsequently converted
to provide the |E*| master curve. A main limitation in the above backcalculation scheme is that it
assumes the entire depth of the AC layer to be at a constant temperature. However, it is likely
that a non-uniform temperature profile will exist in the field (MEPDG, 2004; and Alkasawneh,
2007a); a situation that should be taken into consideration because of the thermo-sensitivity of
the AC layers, especially in comparison with base and subgrade materials (Alkasawneh, 2007a).
The MEPDG accounts for the effects of non-uniform temperature profiles by subdividing the
pavement layers into multiple sublayers and assigning a different modulus value to each on the
basis of the prevailing temperature conditions. Similarly, Alkasawneh et al. (2007b) proposed to
use sublayering to capture modulus variation in the AC layer through elastic layered analysis.
24

However, dividing the AC layer into multiple sublayers increases the number of variables in
backcalculation and hence limits the information that can be obtained regarding the timetemperature properties of the AC layer.
2.5 Characterizing Linear Viscoelastic Properties of HMA in Laboratory
Linear viscoelastic materials possess both elastic as well as viscous characteristics. To
ensure that the properties measured in laboratory are within the linear range, the applied strains
are kept under 100-120 µε. The theory of viscoelasticity states that if one of the following linear
viscoelastic properties is known (i.e. |E*|, D(t) or E(t)), the others can be calculated
mathematically through numerical interconversion procedures (Park and Schapery, 1999; and
Kim 2009). For ease in measurement, typically |E*| is measured in laboratory through cyclic
sinusoidal loading. As shown in equation below the measured strain lags the applied stress by an
angle δ.
	
	
where,

is the stress amplitude,

frequency, is time and

	 	

	

	

	

	

	5	

	

	

	

	

	

	6	

	

2

is the strain amplitude,

is applied angular

is phase angle. The |E*| is calculated as the ratio of stress amplitude to

strain amplitude as shown in the equation (Kim, 2009)
|

∗

	

|

	

	

	

	

	

	7	

A much convenient way to define the viscoelastic behavior is using complex form representation
of stress and stain

25

∗

	 	

∗

The dynamic complex modulus

∗

	

	

	

	

	

	8	

	

	

	

	

	9	

is defined through the stress strain relationship as shown in

the equation.
∗

∗

"			

∗

	

	

	 10 	

where,
cos

and

"

is the ratio of in-phase stress and strain and is called as the storage or elastic modulus. " is
the ratio of out-of-phase stress and strain and is called as the loss or viscoelastic modulus. For a
perfectly elastic material " is expected to be zero, whereas for a perfectly viscoelastic material
is expected to be zero. In complex plain the dynamic modulus is defined as the norm of the

∗

and phase angle is defined as the argument express as:

| E* | E'2 E"2 		 	
  tan1 (E" / E' ) 	 	

	
	

	 11 	

MEPDG considers HMA as a viscoelastic material. It recognizes that |E*| values of HMA
vary with both, the test frequency and test temperature. To incorporate the effect of both the
frequency and the temperature in a single curve, a master curve is generated at a chosen
26

reference temperature. Once a master curve is constructed, |E*| values at temperatures different
from the reference temperature are calculated using the time-temperature superposition principle.
Level-1 analysis in MEPDG requires |E*| values at minimum 3 temperatures and at minimum 4
frequencies directly measured from laboratory. While inputting |E*| values in MEPDG, the
following guidelines are outlined for selecting the temperatures and frequencies


At least one data at temperature below 30oF.



At least one data at temperature between 40F-100oF.



At least one data at temperature higher than 125oF.



At least one data at frequency less than 1 Hz.



At least two data at frequency between 1-10 Hz.



At least one data at frequency greater than 10 Hz.
Maximum of 8 temperature levels and a maximum of 6 frequency levels are allowed and

inputs at 5 temperature and 4 frequency levels are recommended. The default temperatures and
default frequencies used are: 14, 40, 70, 100 and 130˚F and 0.1, 1, 10 and 25 Hz respectively.
The master curve and shift factor are obtained by fitting the raw |E*| values measured at
different temperatures to a single master curve. Typical steps followed in the development of |E*|
master curve from dynamic modulus test are:
1) Determine dynamic modulus at different temperatures (AASHTO T342 recommends -10,
4, 21, 37 and 54oC). At each temperature, test the specimen at a range of frequencies

27

(AASHTO T342 recommends 25, 10, 5, 1, 0.5, 0.1 Hz). Figure 6 shows an illustrative
diagram of raw |E*| data generated in lab.
2) Select a reference temperature,

for the dynamic modulus master curve.

3) Using time-temperature superposition (TTS) principle, shift the measured raw |E*| data to
a single master curve. The master curve is usually defined as a sigmoid:
|

where,

,

,

and

∗|

	

	

.

are sigmoid coefficients and

is the test frequency and

		

During the shifting process, the shift factor coefficients,
and

defined at the reference

	

	

	

	 13 	

is the shift factor for a given temperature T.

	

10

,

	 12 	

. The equation used in the shift is:
.

where,

	

is reduced frequency. Note that

the shifted master curve is obtained for the reduce frequency
temperature

	

and

	

	

	 14 	

and sigmoid coefficients

,

are numerically optimized until a good sigmoid fit is obtained. Figure 6 shows an

illustrative diagram of the shifting in the optimization process.
However, since the primary input of the layered viscoelastic algorithm utilized in this
research is E(t), it is backcalculated first, allowing then for |E*| master curves to be derived via
interconversion for comparison. In the next subsection, conversion predominantly used in the
study, E(t) to |E*| is explained.
28

10

5

Dynamic Modulus (MPa)

(a)
10

4

10

3

o

-10 C Hz
o

10

10 C Hz

2

o

21 C Hz
o

37 C Hz
o

10

54 C Hz

1

10

-6

10

-4

-2

10

0

2

10

10

4

6

10

10

Frequency (Hz)

5

(b)

(c)

104

o

-10 C

4

o

10

10 C
o

10

21 C

2

o

T

Shift factor, a (T)

Dynamic Modulus (MPa)

10

3

10

o

-10 C Hz

2

10

o

10 C Hz
o

21 C Hz

1

10

o

37 C
o

54 C

100
10-2

37 C Hz
o

54 C Hz

0

10

-7

10

-5

10

-3

10

-1

10

1

10

3

10

5

10

-4

10

7

10

-10

0

10

20

30

40

50

60

Temperature (oC)

o

Reduced Frequency at 21 C (Hz)

Figure 6: Figure illustrating steps in developing |E*| master curve using laboratory
measured |E*| test data :(a) raw |E*| data, (b) master curve development (c) timetemperature shift factors
2.5.1 Interconversion of E(t) to |E*|
Relaxation modulus is the ratio of stress response (with time) to a unit stain input. The
relaxation behavior is typically fitted as Prony series (i.e., generalized Maxwell), which is a
summation of sequence of exponentially decaying functions expressed as (Kim, 2009)

29

N

E (t )  E    E i e  t /  i 	

	

	 15 	

= modulus of each maxwell spring,

i i E is

i 1

where,

= long term relaxation modulus,

	

	

i

relaxation time, and  i is the viscosity of each dashpot element in the generalized Maxwell
model. A numerical optimization is performed to determine the Prony series coefficient

and

in equation. Typically, a series of relaxation time from a range of 10-8 sec to 10+8 seconds are
selected.

and

values in Equation 15 are varied until a good fit to the E(t) is obtained.

Behavior of viscoelastic materials are typically modeled using prudent combination of springs
and dashpot. The mathematical expression in Equation 1 can be shown to have physical
representation of springs and dashpot as shown in Figure 7, often referred to as Generalized
Maxwell model.

σ 
E 1 

E2 

E3 

η1 

η2 

η3 

…...

E ∞ 

En 
ηn 

`

σ 
Figure 7: Generalized Maxwell model
Mathematically the complex modulus can be defined as a complex number as shown in
Equation (10). The real part (E’(f) storage modulus) and the complex part (E”(f) loss modulus)
30

of complex modulus can be obtained for a Generalized Maxwell model (fitting Prony series to
relaxation modulus) shown in Figure 7, using the following equations:

E '  f   E *  f  * cos    E 

2f i 2
  Ei
1  2f  i 2 	
i 1

N

E"  f   E *  f  * sin    Ei
i 1

N

	

	 16 	

	

2f i

1  2f i 2 	

	

	

	

17 	

where  is the phase angle, |E*| is the absolute value of the complex modulus E* function, Ei =
modulus of each Maxwell spring,

i i E is relaxation time, and  i is the viscosity of each
i

dashpot element in the generalized Maxwell model. The dynamic modulus master curves were
calculated from the storage modulus and the loss modulus using Equation (11).
2.5.2 Interconversion D(t) to E(t)
Interconversion from D(t) to E(t) can be performed mathematically owing to the fact that
the product of dynamic creep compliance D* and dynamic modulus E* in frequency domain is
equal to one. An appropriate conversion using this method would require D(t) data over a large
range of time (~10-6 to 10+6 seconds). However, due to testing constraints, creep data are
available only over a smaller range of time and temperature, which is not sufficient for
conversion using the method. Creep compliance data can be converted to E(t) assuming a
classical power-law function for the creep compliance (see Equation (18)). The measured data at
different temperatures can be fitted separately and the associated relaxation modulus is then
obtained using the mathematically exact formula given in Equation (19) (Kim, 2009):

31

D(t)  D1tn 	
E ( t ) D (t ) 

	

	

	

	

	

sin n
n 	

	

	

	

	

	 18 	
	 19 	

where, D1 and n are the power function coefficients of D(t). Once E(t) at each temperature are
obtained, E(t) master curve can be generated using a similar procedure described for |E*| (see
Equation (12) and (13)).
2.6 Motivation and Scope
From the past studies it can be seen that most of the work on backcalculation have been
carried out considering AC as elastic layer. Backcalculation programs which can consider
viscoelastic properties of AC are computationally too expensive. Some of the issues related to
the existing studies can be identified as follows:


The existing backcalculation programs consider AC as elastic layer and cannot be used to
predict viscoelastic properties of HMA in situ.



In the existing backcalculation procedures the FWD temperature corrections are carried
out using empirically developed equation. The correction procedure does not recognize
the actual viscoelastic characteristics of the HMA.



It is well known that FWD deflections and backcalculated properties are influenced by
AC temperature. At present there are no guidelines for FWD testing temperature.



Existing studies have shown that both the viscoelastic properties of AC and stress
sensitivity of unbound layers can play a crucial role in pavement response. None of the

32

existing backcalculation programs consider both the effects simultaneously in
backcalculation.
These shortcomings in the current practice have motivated the author to develop new
forward and backcalculation models for FWD analysis. In the present work, a procedure was
developed for backcalculating these properties from conventional Falling Weight Deflectometer
(FWD) test data in combination with knowledge of the testing temperature profile in the asphalt
layer. The present work has two main objectives: (1) Attempt to simulate more realistic FWD
test conditions with respect to the existence of a non-uniform temperature profile across the
depth of the AC layer, and (2) develop a procedure for viscoelastic backcalculation, taking into
account the field measured temperatures. The models are based on Quasi Linear Viscoelasticity
(QLV) theory developed by Schapery (Schapery, 1965). The current version of the algorithm is
not able to consider the dynamics, i.e., the effects of wave propagation. In pavements where the
bedrock is close to the surface and/or there is a shallow groundwater table, the FWD test
deflection time history may exhibit significant wave propagation effects. In such cases, the
current version of the algorithm should not be used.

33

CHAPTER 3

LAYERED VISCOELASTIC PAVEMENT MODEL
3.1 Introduction
Traditionally, flexible pavements are analyzed using analytic multilayered elastic models
(KENLAYER (Huang, 2004); BISAR (De Jong et al., 1973); CHEVRONX (Warren and
Dieckmann, 1963)), which are based on Burmister’s elastic solution of multilayered structures
(Burmister, 1943; Burmister, 1945a,b,c). These models assume the material in each pavement
layer to be linearly elastic.
In the proposed approach, the (asphalt) pavement system is modeled as a layered halfspace, with top layer as a linear viscoelastic solid. All other layers (base, subbase, subgrade,
bedrock) in the pavement are assumed either linear or nonlinear elastic. Assuming there is fullbonding between the asphalt layer and the underlying base and subgrade layers, the overall
response of the entire pavement system becomes viscoelastic. Therefore, its response under
arbitrary loading can be obtained using Boltzman’s superposition principle (i.e., the convolution
integral) (Kutay et al., 2011; and Levenberg, 2008):

t

R ( x, y, z, t )   RHve ( x, y, z, t  )
ve

 0

34

dI ( )
d
d

		

	

	

	 20 	

ve
where, R ve ( x, y, z , t ) is linear viscoelastic response at coordinate (x,y,z) and time t, RH ( x, y, z, t ) is

the (unit) viscoelastic response of the pavement system to a Heaviside step function input (H(t))
and dI ( ) is the change in input at time τ.
It is worth noting that for a uniaxial viscoelastic system (e.g., a cylindrical asphalt
ve
mixture), if response Rve = ε(t) = strain, then RH = D(t) = creep compliance and I(t) = σ(t) =

stress. Using Schapery’s quasi-elastic theory, the viscoelastic response at time t to a unit input
function, can be efficiently and accurately approximated by elastic response obtained using
relaxation modulus at time t (Schapery, 1965; 1974) as follows:

R Hve ( x , y , z , t )  R He x , y , z , E t  	

	

	

	 21 	

where, RHe x, y, z, E t  is unit elastic response calculated using the elastic modulus equal to relaxed
modulus (E(t)) at time=t. Flexible pavements are exposed to different temperatures over time,
which in turn influence their response. For thermorheologically simple materials, this variation in
response can be predicted by extending Equation (20) and Equation (21) as follows:

R ( x, y , z , t ) 
ve

tR

 R x, y, z, E t

e
H

R

  

0

dI ( )
d
d

		

	

	 22 	

where, t R  t a T T  , where aT T  is shift factor at temperature T defined in Equation (11) and
Tref is reference temperature.





aT T   a1 T 2  Tref2  a 2 T  Tref

35

 		

	

	

	

23 	

where, a1 and a2 are shift factor’s polynomial coefficients. Using Equation (22), formulation for
predicting vertical deflection of a linear viscoelastic asphalt pavement system subjected to an
axisymmetric loading can be expressed as:

ve
vertical

u

(r , z , t ) 

tR

 u

e
H vertical

EtR  , r, z  d ( ) d
d

0

	

	

	

24 	

ve

where, uvertical(r, z, t) is the viscoelastic response of viscoelastic multilayered structure at time t
and coordinate (r,z), u He  vertical

E t R

  , r , z  is the elastic unit response of the

pavement system at reduced time tR due to unit (Heaviside step) contact stress (i.e., σ(t)=1), and
 ( ) is the applied stress at the pavement surface.
ve
e
In this implementation, the vertical surface displacements, i.e., u H ( t )  u H ( t ) values at

the points of interest were computed using CHEVRONX which was later replaced by an inhouse multilayer elastic analysis program called LayerE. Then, the convolution integral in
Equation (12) is used to calculate the viscoelastic deflection u ve (t ) . Detailed stepwise
description of the algorithm is given in the following section.
3.2 Layered Viscoelastic (Forward) Algorithm (LAVA)
As explained earlier, using TTS principle, flexible pavement response at any temperature
and loading frequency can be obtained. An algorithm that numerically calculates the convolution
integral described in Equation 24 has been developed and referred to as LAVA. Algorithm steps
followed in LAVA are as follows (kutay et al. 2011):

36

1. Define the geometric (layer thicknesses, contact radius) and material (E(t), Ebase, Esubgrade,
Poisson’s Ratio) properties of a layered system as shown in Figure 8.
2. Discrete stress time history σ(t) into Ns intervals as shown in Figure 9.
r=5.9”

(t)
h1

AC
Base

h2

1

2

3

4

5

6

E(t)
Ebase

Subgrade Esubgrade, 3

h3=Inf

Figure 8: Typical flexible pavement geometry for analysis
120

Stress (psi)

100

80

60

40

20

0
0

0.02

0.04

0.06

0.08

0.1

Time (sec)

Figure 9: Discretization of stress history in forward analysis
3. As shown in Figure 10, divide the relaxation modulus master curve into NE number of
time steps in log scale. The relaxation modulus E(t) can be approximated with a sigmoid
function as follows:
37

logE(t)  c1 

c2
1  exp(c3  c4 log(tR )) 					 	

	

	 25 	

where tR is reduced time ( t R  t / a T (T ) ) and ci are sigmoid coefficients. The shift factor

Relaxation modulus, E(t) psi

coefficients are computed using the second order polynomial given in Equation (23).

10

7

10

6

10

5

10

4

10

3

logE (t )   c1 

10

c2
1  exp(c3  c4 log(t R ))

Number of time interval N =100
E

-9

10

-7

10

-5

10

-3

10

-1

10

1

10

3

10

5

10

7

10

9

o

Reduced time, t at 19 C (sec)
R

Figure 10: Discretization of the relaxation modulus master curve
4. For a unit stress load, calculate the elastic response at every discretized reduced time tR1,
tR2, tR3 ….tRNE, using the modulus values as E(tRi) obtained from Figure 10 (Schapery
1965, 1974). Figure 11 shows the vertical surface deflection u He values calculated using
LayeredE at various radial distances similar to one shown Figure 8. Since, the computed
u He values are now the fundamental behavior of the viscoelastic multilayer system these

curves are herein called “unit response master curves”.
5. Viscoelastic response is calculated numerically using the discrete form of Equation (24)
given in Equation (26).
38

i

u ve (t i )   u He (t i   j ) d ( j )

where i  1,2,...N s

j 0

		

	

26 	 	

One of the primary reasons for implementing Schapery’s ‘quasi-elastic’ solution above is
due its extreme computational efficiency. Table 3 shows the computation times using a Pentium
2.66 GHz computer with 3.25 GB ram for different number of discrete time steps in the 3 layered
system shown in Figure 8.

rc=0 in
rc=8 in

-2

rc=12 in
rc= 8 in

rc=24 in
rc=36 in

rc=48 in
rc=60 in

-3

10

-4

10

U

ve

H

e

H

=U (t), (inches)

10

-5

10

-9

10

-7

10

-5

-3

10

10

-1

10

1

10

3

10

5

10

7

10

9

10

o

Reduced time, t at 19 C (sec)
R

Figure 11: Deflections calculated under unit stress for points at different distances from the
centerline of the circular load at the surface.
Table 3: LAVA computation times for different numbers of discrete time steps
Ns

NE

50
24
50
100
100
200

50
100
100
100
200
200

Computation
time (sec.)
1.96
2.88
3.03
3.05
5.01
5.13
39

3.3 Implementation of Temperature Profile in LAVA
Temperature in pavements typically varies with depth, which affects the response of the
HMA to the applied load. As shown in Figure 12, the temperature may be increasing with depth
(profile 1-linear, 2-piecewise, and 3-nonlinear) or decreasing with depth (profile 4-linear, 5piecewise, and 6-nonlinear) depending on the time of the day. This variation in temperature with
depth can be approximated with a piecewise continuous temperature profile function as shown in
Figure 12 (profile 2 and 5). The advantage of using a piecewise function is that it may be used to
approximate any arbitrary function.

Depth, z (in)

Temperature (T)

Figure 12: Schematic diagram of temperature profile
3.3.1 Layered Viscoelastic (Forward) Algorithm for Temperature Profile (LAVAP)
An algorithm that considers HMA sublayers with different temperatures within the HMA
layer has been developed, and referred to as “LAVAP” (LAVA profile). The steps involve in the
algorithm were same as for LAVA except that multiple viscoelastic sublayers exist in LAVAP.
Algorithm steps followed in LAVAP are as follows (Chatti et al., 2014):
1. Define the geometric (layer thicknesses, contact radius) and material (E(t), aT(T), Ebase,
Esubgrade, Poisson’s Ratio) properties of a layered system as shown in Figure 8.
2. Define number of sublayers in the AC layer NAC and average temperature in each
sublayer.
40

3. Discretize stress time history σ(t) into Ns intervals as shown in Figure 9.
4. As shown in Figure 10 divide the relaxation modulus master curve into NE number of
time steps in log scale. The relaxation modulus E(t) can be approximated with a sigmoid
function given by Equation (25).
5.

Calculate the shift factors aT T j  for each sublayer j, computed using the second order
polynomial given in the Equation (23).

6. For a unit stress load, calculate the elastic response at every discretized reduced time t1,
t2, t3 ….tNE, using the modulus values for jth AC sublayer as E(tRij) obtained from Figure
10 (Schapery 1965, 1974), where

 .

t Rij  t i a T T j

7. Viscoelastic response is calculated numerically using the discretized Equation (26).
The algorithm has been compared with LAVA as well as ABAQUS. Comparison with LAVA
was made for deflection response at all the sensors for constant temperature throughout all the
sublayers. The pavement section and layer properties used in the forward analysis are shown in
Table 4. Figure 13 shows the response obtained from the LAVAP algorithm at 0oC, 30oC and
50oC match very well with LAVA.
Table 4: Pavement properties used in LAVAP validation with ABAQUS.
Property
AC sublayers
Thickness
Granular layers
Poisson ratio {layer 1,2,3…}
Eunbound {layer 2,3…}
Total number of sensors
Sensor spacing from load (inches)
E(t) sigmoid coefficient {AC}
a(T) shift factor coefficients {AC}

Constant temperature
Temperature profile
6 in
2in, 2in, 2in
20, inf (in)
20, inf (in)
0.35, 0.3, 0.45
0.35, 0.3, 0.45
11450, 15000 (psi)
11450, 15000 (psi)
8
0, 7.99, 12.01, 17.99, 24.02, 35.98, 47.99, 60
0.841, 3.54, 0.86, -0.515
0.841, 3.54, 0.86, -0.515
4.42E-04, -1.32E-01
4.42E-04, -1.32E-01
41

o

o

o

o

Temperature = 0 C, 0 C, 0 C

0.012

o

o

Temperature = 30 C, 30 C, 30 C

0.02
LAVA
LAVAP

LAVA
LAVAP

0.01
0.015

U (t) in

0.006

0.01

VE

VE

U (t) in

0.008

0.004
0.005
0.002

0

0
0

0.005 0.01 0.015

0.02 0.025

0.03 0.035 0.04

0

0.005 0.01 0.015

Time (sec)

0.02 0.025

0.03 0.035 0.04

Time (sec)
o

o

o

Temperature = 50 C, 50 C, 50 C

0.03

LAVA
LAVAP
0.025

VE

U (t) in

0.02

0.015

0.01

0.005

0
0

0.005 0.01 0.015

0.02 0.025

0.03 0.035 0.04

Time (sec)

Figure 13: Comparison of response calculated using LAVAP and original LAVA
In order to qualitatively examine the response of flexible pavement predicted using
LAVAP algorithm, the response obtained under temperature profile was compared with the
response obtained under constant temperatures. As an example, a comparison of the response
under a temperature profile of {40-30-20}oC with that corresponding to a constant temperature
of 40oC, 30oC and 20oC for the entire depth is shown in Figure 14, Figure 15 and Figure 16,
respectively. It can be seen from the figures that the effect of AC temperature is most prominent
in sensors closer to the load center (sensors 1, 2, 3, and 4). For sensors away from the loading
center (sensors 5, 6, 7 and 8), the deflection histories are not influenced by the AC temperature.
42

Figure 17 shows the region of E(t) master curve (at 19oC reference temperature) used by the
LAVAP algorithm in calculating time histories.

o

o

o

Temperature = 40 C, 30 C, 20 C

0.025

o

LAVA 40 C
LAVAP

0.015

VE

U (t) in

0.02

0.01

0.005

0
0

0.005 0.01 0.015

0.02 0.025

0.03 0.035 0.04

Time (sec)

Figure 14: Comparison of responses calculated using LAVAP at temperature profile {4030-20}oC and original LAVA at constant 40oC temperature.
o

o

o

Temperature = 40 C, 30 C, 20 C

0.02

o

LAVA 30 C
LAVAP

VE

U (t) in

0.015

0.01

0.005

0
0

0.005 0.01 0.015

0.02 0.025

0.03 0.035 0.04

Time (sec)

Figure 15: Comparison of responses calculated using LAVAP at temperature profile {4030-20}oC and original LAVA at constant 30oC temperature.

43

o

o

o

Temperature = 40 C, 30 C, 20 C

0.02

o

LAVA 20 C
LAVAP

U (t) in

0.015

VE

0.01

0.005

0
0

0.005 0.01 0.015

0.02 0.025

0.03 0.035 0.04

Time (sec)

Figure 16: Comparison of responses calculated using LAVAP at temperature profile {4030-20}oC and original LAVA at constant 20oC temperature.
10

7

o

o

Relaxation modulus, E(t) psi

20 C

10

6

10

5

10

4

o

30 C

40 C

o

40 C
o

30 C
o

10

20 C

3

10

-5

-3

10

10

-1

10

1

10

3

5

10

Reduced time (sec)

Figure 17: Region of E(t) master curve (at 19oC reference temperature) used by LAVAP
for calculating response at temperature profile {40-30-20}oC

44

As expected, it can be seen that for a condition of higher temperature at the top and lower
temperature at the bottom, the response with a higher constant temperature will always be greater
than the response with a temperature profile. The response with a lower constant temperature
will always be less than the response with a temperature profile. The response with a medium
constant temperature may or may not be less than the profile response depending on the
temperature profile and thickness of the sublayering.
3.3.2 Comparison of LAVAP with FEM Software ABAQUS
Next, the LAVA algorithm was validated against the well-known finite element software,
ABAQUS, where temperature profile in AC layer was simulated as two-sublayers of AC with
different temperatures. For this purpose, two different HMA types were considered, namely;
Terpolymer and SBS 64-40. The viscoelastic properties of these two mixes are shown in Figure
18. As shown in Table 5, for both the mixes, the AC layer was divided into two sublayers, where
temperature in the top and bottom sublayer were assumed to be 19oC and 30oC respectively.

7

3

10

Terpolymer
SBS 64-40

Terpolymer
SBS 64-40

2

10

6

10

1

Shft factor, aT

Relaxation modulus, E(t) (psi)

10

105

10

0

10

10-1
10-2

4

10

10-3

103
10-9

10-7

10-5

10-3

10-1

101

103

105

107

109

o

10-4
-10

0

10

20

30

40
o

Reduced time at 19 C (sec)

Temperature ( C)

Figure 18: Relaxation modulus and shift factor master curves at 19oC reference
temperature
45

50

60

Table 5: Pavement section used in LAVAP validation
Poisson’s
Ratio

AC

Mix1: Terpolymer
(E(t) = Figure 18);
Mix 2: SBS 64-40
(E(t) = Figure 18)

Thickness(in)
(Temperature oC)
Sublayer1 = 3.94”
(19oC)
Sublayer2 = 3.94”
(30oC)

Base

15000 psi (linear elastic)

7.88”

0.35

Subgrade

10000 psi (linear elastic)

Infinity

0.45

Layer

Modulus (E(t) or E)

0.45

Comparison of surface deflection time histories measured at radial distances 0, 8, 12, 18,
24, 36, 48, 60 inches for Mix 1 and 2 are shown in Figure 19 and Figure 20. From the figures it
can be observed that the results obtained from LAVAP and ABAQUS match well.

0.001
LAVAP
ABAQUS

0.0006

VE

U (t) in

0.0008

0.0004

0.0002

0
0

0.005 0.01

0.015 0.02 0.025 0.03 0.035

0.04

Time (sec)

Figure 19: Comparison between LAVAP and ABAQUS at a temperature profile of {1930}oC (Terpolymer)

46

0.0012
LAVAP
ABAQUS
0.001

VE

U (t) in

0.0008

0.0006

0.0004

0.0002

0
0

0.005 0.01

0.015 0.02 0.025 0.03 0.035

0.04

Time (sec)

Figure 20: Comparison between LAVAP and ABAQUS at a temperature profile of {1930}oC (SBS 64-40)
As expected, it can be seen from Table 6 that for both the mixes, surface deflection in the
pavement section at two step AC temperature profile of {19-30}oC lies between the deflections
obtained for constant AC temperatures of 19oC and 30oC.
Table 6: Peak deflections at temperature profile {19-30}oC and at constant 30oC
temperature using LAVA.

SBS 64-40

Terpolymer

Mix

Temp (oC)
Constant:
19oC
Profile:
19oC, 30oC
Constant
30oC
Constant
19oC
Profile:
19oC, 30oC
Constant
30oC

Deflection (mx10-4) (Sensors)
3
4
5
6

1

2

7

8

7.14

6.27

5.72

4.95

4.26

3.14

2.36

1.84

8.39

7.13

6.32

5.27

4.39

3.11

2.30

1.79

9.76

7.84

6.81

5.54

4.52

3.10

2.27

1.76

8.94

7.40

6.53

5.41

4.47

3.13

2.30

1.78

1.04

8.26

7.03

5.59

4.50

3.06

2.24

1.75

1.21

8.92

7.42

5.76

4.55

3.04

2.22

1.73

47

3.4 Layered Viscoelastic Nonlinear (LAVAN) Pavement Model
In this section, a computationally efficient layered viscoelastic nonlinear model, herein
called LAVAN, is presented. The LAVAN can consider linear viscoelasticity of AC layers as
well as stress-dependent modulus of granular layers. The formulation is inspired from QuasiLinear-Viscoelastic (QLV) constitutive modeling (Leaderman, 1943; Schapery, 1969; Fung,
1996; Masad et al., 2008; Yong et al., 2010; and Nekouzadehn and Genin, 2013), which is often
used in analyzing nonlinear viscoelastic materials. In literature, the various forms of the model
are also named as Fung’s model/ Schapery’s nonlinearity model/ modified Boltzmann’s
superposition.
It is well known that the unbound granular materials exhibit nonlinearity, i.e. stress
dependent modulus (Hicks and Monismith, 1971; Witczak and Uzan, 1988; and Ooi et al., 2004).
Modeling response of flexible pavements with stress dependent granular layers requires iterative
analysis, which makes the solution computationally intensive. In order to avoid this complexity,
most of the standards (Shell, 1978; Asphalt Institute, 1999, Theyse et al., 1996; IRC, 2001; and
Austroads, 2004) assume the granular layer to be elastic. Wang and Al-Qadi (2013) developed a
FE based model for asphalt pavements in ABAQUS and considered both viscoelasticity of
asphalt and nonlinearity of granular layers. They found that both the viscoelastic and nonlinear
material properties need to be considered in order to accurately predict flexible pavement
response. They also showed that the nonlinear response of the granular layer is influenced by the
viscoelastic stress caused by the AC layer. Although FE-based modeling is very useful, they are
computationally expensive and can be significantly affected by the boundary conditions.
Furthermore, FE-based models are inefficient to be used as an inverse analysis tool, e.g., in

48

backcalculation algorithms. A computationally efficient semi-analytical model is inevitable in
such algorithms.
In the present study, Quasi-Elastic theory is combined with generalized QLV theory, to
develop an approximate model for predicting response of multilayered viscoelastic nonlinear
flexible pavement structures. Before introducing the generalized QLV model, a brief overview of
granular nonlinear pavement models is also presented. Development of a generalized QLV
model for multilayered system is followed by numerical validation, in which response of flexible
pavements under stationary transient loading has been analyzed. The model is implemented in
general purpose FE-based software, ABAQUS, to validate the model.
3.4.1 Nonlinear Multilayer Elastic Solutions
Under constant cyclic loading, granular unbound materials exhibit plastic deformation
during the initial cycles. As the number of load cycles increase, plastic deformation ceases to
occur and the response becomes elastic in further load cycles (a phenomenon known as
shakedown). Often this elastic response is defined by resilient modulus (MR) at that load level,
which is expressed as:

	

where,

	

	

, is the deviatoric stress in a triaxial test and

	

	

	

	 27 	

is recoverable strain. If the

granular layer reaches this steady state condition under repeated vehicular loading, then the
further response can be considered recoverable and Equation (27) can be used to characterize the
material. However, the MR value shown in Equation (27) is affected by the state of stress (or load
level). Typically, unbound granular materials exhibit stress hardening (Yau and Von Quintus,
2002; and Taylor and Timm, 2009), i.e., the

increases with increasing stress. Hicks and
49

Monismith (1971) related hydrostatic stress and the resilient modulus obtained in Equation (27)
to characterize the stress dependency of the material. The model suggested by Uzan (1985) and
Witczak and Uzan (1988) use deviatory stress as well as octahedral stress to incorporate the
distortional shear effect into the model. The model has been further modified by various
researchers. Yau and Von Quintus (2002) analyzed LTPP

test data using generalized form of

Uzan (1985) model expressed as:

	 	

where,
,

is hydrostatic stress,

,

,

and

is atmospheric pressure,

	

	

	 28 	

is octahedral shear stress,

are regression coefficient. They found that parameter k6 in the equation

regressed to zero for more than half of the tests, and hence the coefficient was set to zero for the
subsequent analysis. The subsequently modified equation is shown in equation (29), which has
been used in the present study to define resilient modulus of unbound granular layer:

1

1

where,
regression constants,

2

,

	

	

	

	

is octahedral shear stress,

is the coefficient of earth pressure at rest and

		 29 	

,

,

are

is atmospheric

pressure.
Although the resilient modulus, MR, is not the Young’s Modulus (E) (Hjelmstad and
Taciroglu, 2000), while formulating granular material constitutive equations, it is often used to
replace E in the following equation:

	 	
50

	

	

	 30 	

where,

is Young’s Modulus,
,

is Poisson’s Ratio,

is stress tensor,

is strain tensor,

is Kroenecker delta. Nonlinear stress dependent

has been

implemented in many FE-based models. These include GTPAVE (Tutumluer, 1995), ILLIPAVE
(Raad and Figueroa, 1980), and MICHPAVE (Harichandran et al., 1989). Typically FE-based
nonlinear pavement analysis is performed by defining a user defined material (UMAT) in FEbased softwares such as ABAQUS and ADINA (Hjelmstad and Taciroglu, 2000; Schwartz,
2002; and Kim et al., 2009). Although FE-based solutions are promising, apart from being
computationally expensive, they may exhibit significant influence of boundary condition if the
semi-infinite geometry of the problem is not implemented adequately in the model.
An approximate nonlinear analysis of pavement can also be performed using Burmister’s
multilayered elastic based solution (Burmister, 1945). For incorporating variation in modulus
with depth, Huang (Huang 2004) suggested dividing the nonlinear layer into multiple sublayers.
Furthermore, he suggested choosing a specific location in the nonlinear layers to evaluate
modulus based on the stress state of the point. Zhou (2000) studied stress dependency of base
layer modulus obtained from base layer mid depth stress state. They analyzed FWD testing at
multiple load levels on two different pavement structures. The study showed that reasonable
nonlinearity parameters can be obtained through regression of backcalculated modulus with
stress state at mid-depth of the base layer.
In the present study, the elastic nonlinearity is solved iteratively assuming an initial set of
elastic modulus. The evaluated stresses obtained using the initial values of modulus are used to
evaluate new set of modulus using Equation (29). The iteration is continued until the computed
modulus from the stresses predicted by the layered solution and input modulus used in the

51

layered solution converges. It is noted that the appropriate stress adjustments were made due to
the fact that unbound granular material can’t take tension. This means that in such a case either
residual stress is generated such that stress state obeys a yield criterion or the tensile stresses are
replaced with zero.
3.4.2 Proposed Constitutive Model for Multilayer Pavement Model under Uniaxial Loading:
Combining Linear Viscoelastic AC and Nonlinear Base
Mechanical behavior of elastic and linear viscoelastic materials are well defined in
literature using springs and dashpots. In the derivation presented, the standard mechanical
responses using springs and dashpots is implied, and has not been stressed upon for brevity. As
shown in Figure 21, for the uniaxial formulation, a two layer system comprising of linear
viscoelastic layer that is in perfect bonding with the underlying nonlinear unbound layer is
considered. A third layer of elastic material is avoided, because it can be easily lumped along the
viscoelastic properties in the system. Under uniaxial loading, formulation for vertical surface
response has been presented.

A
Linear
Viscoelastic
AC

,

Nonlinear base
material

Figure 21: Cross section of multilayer viscoelastic nonlinear system used in uniaxial
analysis.
52

The vertical deflection response of the two layer system can be obtained as an addition of
responses of individual layers. That is, the vertical surface deflection at A in Figure 21 can be
expressed as:

		 	

,
where,

	

is the total vertical deflection at point A,

	

	

		 31 	

is the vertical deformation in

is the vertical deformation in the linear viscoelastic AC layer. Unlike

nonlinear layer and

a system comprising of only linear materials, a nonlinear system is expected to have different
unit response functions at different loads. In this study, unit response function of the system at
any load level has been defined as a secant property (like secant modulus), such that at any load
level

the unit response function

,

	
	

can be defined similar to the definition used for

linear viscoelastic materials. Equation (32) show expressions for unit response functions at
different load levels,

	

according to the definition explained. It is worth noting that in the

uniaxial case, the loading stresses are same as the stress state for calculating the nonlinear
modulus in the nonlinear layer.

,

,

where,

		

	

is the nonlinear unit response of the system at stress equal to

is the deformation in the nonlinear layer at stress equal to
layer at stress equal to
,

,

	 32 	

,

,

is the deformation in AC

is the secant stiffness of the nonlinear layer at stress equal to

is the unit response function of the linear viscoelastic AC material. It is clear from

Equation (32) that, for a constant load a constant modulus exists for the nonlinear layer, which
53

can be used in estimating response of the combined viscoelastic-nonlinear system under any
loading through Boltzmann superposition.
However, since the unit response functions are function of stress, the modified
convolution is as shown below:
,

where,

	

			

	

	
,

is the vertical deflection of the system at time, t,

response function at stress equal to

,

	

	 33 	

is unit nonlinear

	

is arbitrary stress function used as succession of

infinitesimal steps. The modified convolution has been further discussed in detail in the next
section. Similar to the numerical discretization of convolution integral in linear viscoelasticity,
the above equation can be numerically expressed as:
,

	

∆

,

∆

,

∆
,

		

∆
∆

where, ∆

⋯
∆

				

	 34 	

													

	 35 	

∆
⋯

is infinitesimal stress increment at

∆

. As shown in Equation (36) to Equation

(38), Equation (35) can be further expressed as addition of independent response from linear
viscoelastic and nonlinear elastic constituents of the system. The first summation in the equations
represents the well-known convolution integration for linear viscoelastic materials, whereas the
second term represents the secant response from the nonlinear constituent.

54

∑

∆

∑

∑

∆

∑

∆ 													

		

∆
∑

	
∆

	

		 36 	

	

			 37 	

													 	

				 38 	

Although the uniaxial response derivation for multilayer viscoelastic nonlinear system
reduces to a simpler form, it cannot be directly generalized for 2D or 3D (axisymmetric)
conditions because of the following reasons: (1) loading is typically concentrated over specific
loading area, and (2) due to the relaxation of the AC layer, stress state in the nonlinear layer can
change with time. However, the impact of these issues can be minimal if the response of
pavement to the variation in nonlinear modulus value in the radial direction is not significant,
which has been shown by researchers to be an adequate assumption (Huang 2004) in multilayer
nonlinear elastic analysis of flexible pavements. In the next section, this assumption has been
used in deriving response for the viscoelastic nonlinear system, and subsequently, error in the
analysis is discussed.
3.4.3 Proposed Generalized Model for Multilayer Pavement Model: Combining Linear
Viscoelastic AC and Nonlinear Base
3.4.3.1 Applicability of Existing Theories of Nonlinear Viscoelasticity
Mechanistic solutions for NLV (Nonlinear Viscoelastic) materials exhibit variation
depending on the type of nonlinearity that is present. Typical NLV equations involve
convolution integrals that are based on integration kernel which are function of stress or strain.
Equation (39) and (40) show typical form of such expression:
,

55

	

	

	

	

39 	

		
where,

is strain,

is stress,

	

,
,

	

	

	

	 40 	
,

is strain dependent relaxation modulus and

is stress

dependent creep compliance. Typically, in many nonlinear materials, the shape of relaxation
modulus of the material is preserved, even though the material presents stress or strain
dependency (Shames and Cozzarelli, 1997; and Nekouzadeh and Genin, 2013). Such NLV
problems are solved by assuming that time dependence and stress (or strain) dependence can be
decomposed into two functions as follows:

		

,
,
where

is a function of stress,

function of strain,

	

	

	

	 41 	

	 											 	

	

	

	 42 	

is the (only) time dependent creep compliance,

	is a

is the (only) time dependent relaxation modulus. This multiplicative

decomposition facilitates an easy application of superposition principle. For such materials the
following expression has been typically used in NLV formulations to develop the convolution
integral (Nekouzadeh and Genin, 2013):

		
where,

			

	

is a relaxation function that remains unchanged at any strain level, and

function of strain, such that

⁄

	

	 43 	
is a

represents the elastic tangent stiffness at different

strain levels. These models are designated as Fung’s NLV material, which was first proposed by
Leaderman in 1943 (Leaderman, 1943). A generalized form of this nonlinearity model was
presented by Schapery (1969) using thermodynamic principals. Yong et al. (2010) used the

56

model to describe nonlinear viscoelastic-viscoplastic behavior of asphalt sand, whereas, Masad
et al. (2008) used the model to describe nonlinear viscoelastic creep behavior of binders. The
model suggests that nonlinear relaxation function can be expressed as a product of function of
time

⁄

and function of strain
⁄

by the elastic component,

. In Equation (43), nonlinearity is introduced

, and the viscoelasticity comes from the

.

A direct extension of the concepts of QLV model to develop formulations for viscoelastic
nonlinear multilayered system, where the unbound layer is nonlinear and the AC layer is linear
viscoelastic lead to the following:

	

, , ,

where,
strain

, , ,

, ,

is relaxation function at location

	

, and

	

	 44 	

is a function of

. Alternatively, to obtain vertical surface deflection in pavements, Equation (44) can

be expressed in terms of vertical deflection response to Heaviside step loading as follows:
,
,

is surface (NLV) displacement,

where

due to a unit stress and

,

,

	

	 45 	

1 is unit nonlinear elastic response

	 	

	

	

is the nonlinear elastic unit displacement due to a given stress

theory (i.e., Equations (44) through (46)) to hold,
order to investigate this, the

	

is a function of stress, which can be expressed as:

,

where,

	

1

	

	 46 	

. For Fung’s

must be purely a function of stress. In

values were computed using Equation (46) and plotted against

57

surface stress and relaxation modulus (i.e. time). The properties of pavement section and material
properties are shown in Figure 22.

Figure 22: Flexible pavement cross section.
The LAVA algorithm was modified to implement an iterative nonlinear solution for the
,

granular base, which was assumed to follow Equation (29). The
calculated at a range of stress values from 0.1 psi to 140 psi and using

in Equation (36) was
values (for AC) at a

range of times, ranging from 10-8 to 108 seconds. Then,

,

stress. Figure 23 shows the variation of

values decrease with increasing

stress

, where the

1 was calculated for unit

.

1
X: 1360
Y: 0.1
Z: 1.035

1

X: 3.301e+06
Y: 0.1
Z: 0.9985

0.95
0.9

g()

0.8
0.6
0.4
0.2
0

0.85

X: 3.301e+06
Y: 140.1
Z: 0.8887

0.8
0.75

X: 8022
Y: 140.1
Z: 0.6052

30

60

10
90

Stress (psi)

120

10
150

10

5

0

E(t) (psi)

Figure 23: Variation of g(σ) with stress and E(t) of AC layer.

58

10

0.7
0.65

This is expected behavior for a nonlinear material since as the stress increases, the unbound layer
moduli will increase. However, Figure 23 also illustrates that the
. This means that

varies with change in

is not solely based on the stress, as a result, Fung’s model cannot be

used in a layered pavement structure. This is meaningful since the change in the stress
distribution within the pavement layers due to viscoelastic effect (as

varies) will impose

changes in the behavior of stress dependent granular layer. Hence, even though the viscoelastic
layer in a nonlinear multilayered system is linear, it cannot be formulated as a Fung’s QLV
model.
3.5 Proposed Model: LAVAN
QLV model can still be approximated as a convolution integral, provided the stress
dependent relaxation function of the multilayered structure under all the load levels are known.
Using modified superposition, such a generalized QLV model for a multilayered structure can be
expressed as NLV equations involving the convolution integrals of unit response function of the
structure, which is a function of stress or strain as follows:
, , ,

where,
, , ,

, , ,
,

is

, , ,

the

NLV

response

	

,

of

the

layered

	

	

	 47 	

pavement

structure,

is the unit response function that is both function of time and input,

,

which is equal to stress applied at the surface of the pavement. Note that in this formulation,
unlike Fung’s QLV model, time dependence and stress (or strain) dependence are not separated.
Equation (47) can be rewritten in terms of vertical surface deflection under axisymmetric surface
loading as follows:

59

	

	

where,
, ,

, ,

, ,

, ,

	 		

,

is the vertical deflection at time

,

, ,

,

pavement at a loading stress level of

/

, where

	

	 48 	
,

and location

; and

is the nonlinear response of the

. The model in Equation (48) can be expressed in

discretized formulation as follows:
	

0,

where

∑

, ,
,

. The

2D matrix pre-computed for

, ,

,

,

								 	

∆

	 49 	

values are computed via interpolation using the
(which was computed at a range of stress values and

values). The developed model has been referred to LAVAN as an abbreviation for LAVANonlinear.
Step by step procedure to numerically compute response is given below:
1. Define a discrete set of surface stress values:
2. Calculate nonlinear elastic response
for each

	

,

= 0.1 psi to 140 psi.
at a range of

values, by using

value.

3. Recursively compute

until the stress in the middle of the base layer results in the same

as the one used in the layered elastic analysis. At this step, Equation 7 is used in the
nonlinear formulation for the base.
4. Calculate the nonlinear unit elastic response,

, ,

,

as

, ,

,

5. Perform convolution shown in Equation (49) to calculate the NLV response.
60

/

.

3.5.1 Validation of the LAVAN
In order to validate the LAVAN algorithm, a flexible pavement is modeled as a three layered
structure, with viscoelastic AC top layer, followed by stress dependent (nonlinear) granular base
layer on elastic half space (subgrade). Figure 22 shows the geometric properties of the pavement
5.9" and

structure utilized in the validation, where

9.84". The viscoelastic

properties of two HMA mixes, namely CRTB and Control (two materials from FHWA’s ALF
2002 experiment – (Gibson et al., 2012)) were used for the AC layer in the analysis as case 1 and
2. The relaxation modulus master curves of the two mixes are shown in Figure 24.

Relaxation modulus, E(t) (psi)

10

7

Control 70-22
CRTB
10

6

10

5

10

4

10

3

10

-9

10

-7

-5

10

10

-3

10

-1

1

10

10

3

10

5

7

10

10

9

o

Reduced time at 19 C (sec)

Figure 24: Relaxation modulus master curve for LAVAN validation (at 19oC reference
temperature).
These curves were computed from their |E*| master curves by following the Prony series-based
interconversion procedure suggested by Park and Shapery (1999). From the theory of
viscoelasticity, creep compliance, relaxation modulus and dynamic modulus are interconvertible. As a result, the developed viscoelastic nonlinear multilayered model can take any of
three viscoelastic properties. The relaxation modulus
61

can be approximated with a sigmoid

function given by Equation (25). Both the HMA mixes were analyzed with granular nonlinear
model as shown in Equation (29).
3.5.2 Comparison of LAVAN with FEM Software ABAQUS
In ABAQUS, the viscoelastic properties of the HMAs were input in the form of
normalized bulk modulus

and normalized shear modulus

(ABAQUS, 2011). For the

unbound nonlinear layer, a user defined material UMAT was written. ABAQUS requires that
any UMAT should have at least two main components; (i) updating the stiffness Jacobean
Matrix, and (ii) stress increment. Equations (50) and (51) show the mathematical expressions for
these two operations implemented in the UMAT.

	

	
	

where,

is the Jacobian matrix,

	

	

	

	

	 50 	

	

	

	

	

	 51 	

is the updated stress. For nonlinear analysis using

LAVAN, unbound modulus was calculated using stress state at the midpoint of unbound base
layer (vertically). Since LAVAN cannot incorporate nonlinearity along horizontal direction, for
comparison, modulus values were calculated using stress state at
where

3.5 (i.e., r in Figure 22),

is the radius of loading. Estimating base modulus using stress state along the centerline

of the loading would result in stiffer base (Kim et al., 2009); hence it was expected to get closer
results using stress state at some location radially away from the centerline.
In ABAQUS, solution is sensitive to the boundary condition used in the analysis. Kim et
al. (2009) found that for multilayer nonlinear axisymmetric problem with the vertical boundaries
supported using roller supports and the bottom fixed (finite boundaries), the domain size of 140
62

in vertical direction and 20 in the horizontal direction produce results comparable to analytic
solution. Alternatively, researchers have used infinite elements as the boundary to imitate semiinfinite geometry. In the present study, both the boundary conditions were analyzed. It was
observed that, when all the vertical boundaries were supported using roller supports and the
bottom was fixed, the domain size of 133

in vertical direction and 53

in the horizontal

direction was found to produce stable surface deflection (with less than 1% error at the center).
The same domain size is later analyzed with all the vertical as well as the bottom boundaries
supported using infinite elements. For the selected domain size, the FE mesh refinement of
10mm in the AC layer and 25mm in the base layer is used. The model response under transient
loading is compared with results obtained from general purpose FE software ABAQUS. For this
purpose, haversine loading in an axisymmetric setup is used. ABAQUS consumed approximately
17 minutes in analyzing a haversine loading of 138 psi and 35 ms, whereas LAVAN could
generate the results in 3.6 minutes in the same desktop computer.
The surface deflection obtained using FE for the two boundary conditions, (a) finite
boundaries: roller support on the vertical boundaries and fixed support on the bottom (b) infinite
elements at boundaries are compared with LAVAN in Table 7. It should be noted that both the
boundary conditions does not strictly represent the semi-infinite geometry of the problem. In
ABAQUS the solution in the infinite element is considered to be linear, which is assumed to
matches the material properties of the adjacent finite element. Hence the infinite elements
provide stiffness to the boundary assuming the deflection at

∞ to be zero. It can be seen

from Table 7 that, for both the Control mix and CRTB mix, the surface deflection predicted by
FE using finite boundaries were higher compared to results when infinite elements were used at

63

the boundaries. Further for both the mixes the surface deflections predicted by LAVAN were
found to be lying between the FE results predicted by the two boundary conditions.
Table 7: Comparison of ABAQUS and LAVAN: Peak surface deflections for different
boundary conditions.
Control mix (µm)
CRTB mix (µm)
Sensor
radial
Finite
Infinite
Finite
Infinite
distance boundary LAVAN boundary boundary LAVAN boundary
0”
810.9
794.7
747.3
1068.9
1079.5
1005.1
8”
757.4
734.6
693.6
966.3
955.0
902.2
12”
712.1
685.0
648.4
883.9
859.1
819.8
18”
639.5
605.5
575.7
758.6
715.5
694.5
24”
565.9
526.9
501.9
639.4
585.0
575.3
36”
433.9
388.4
369.6
445.5
384.3
381.3
48”
329.1
284.0
264.4
311.6
257.7
246.9
60”
252.5
210.4
187.0
227.2
183.4
161.9
Since, satisfactory performance of the model would require predicting a comparable time
response by the model. The surface deflection history computed by LAVAN and ABAQUS
(analysis with finite boundaries) for the Control mix and CRTB mix are plotted in Figure 25 and
Figure 26 respectively. A comparable response is visible in the figures, which has been further
quantified using Equation (52) and Equation (53). As expected, under the same geometric and
loading conditions the stiffer Control mix generated lower deflections compare to softer CRTB
mix. Figure 25a shows the results when stress at

0 is used in LAVAN for nonlinearity

computation, and provided for comparison purpose Figure 25b shows the results when stress at
3.5

is used in LAVAN. It is noted that S1, S2, S3, S4, S5, S6, S7, S8 in the figures

corresponds to surface deflection at Sensor-1 through Sensor-8 spaced at, 0”, 8”, 12”, 18”, 24”,
36”, 48” and 60” away from the centerline of the load. In general a better match for the
deflection basin between the FE and LAVAN results can be found when stress state at

3.5

is used while incorporating nonlinearity. It is worth noting that, for the structure in Figure 22,
64

Huang’s (2004) procedure for the location of stress used in the nonlinear elastic analysis leads to
2.8 , when trapezoidal stress distribution with 0.5 horizontal slope, 1 vertical slope is
assumed.
The difference between the FE results and LAVAN was quantified using the following
two variables:

100	

	

where

	

	

/
	

/

100	

/

is the percent error in the peaks,

ABAQUS,
	

∑

by LAVAN at time

,

				

	 52 	

	

		 53 	

is peak deflection predicted by

is peak deflection predicted by LAVAN,
is deflection predicted by ABAQUS at time

						

,

is average percent error,
	

is deflection predicted

is total number of time intervals in the deflection time history. Since

the model integrates both viscoelastic and nonlinear material properties, both peak deflection as
well as relaxation of deflection time history should be predicted with accuracy. Therefore,
was used to examine the

normalized average error along the entire deflection history
model performance in creeping. As shown in Table 8, the
mix showed a slight improvement in the results when

and

values for Control

3.5 is used. However, it can be seen

from the table that, error for CRTB mix increased at sensor 1, 2 and 3. As expected, due to the
limitations of LAVAN in incorporating horizontal nonlinearity, both the mixes showed
increasing difference as compared to the FE solutions at sensors away from the loading. In both

65

the cases, however, errors in the first 4-5 sensors were within 6%, indicating good accuracy of
LAVAN.
Table 8: Comparison of ABAQUS and LAVAN: Percent error (PEpeak) calculated using
the peaks of the deflections and average percent error (PEavg) calculated using the entire
time history.
Peak deflection error
Sensor
(%) PEpeak
radial
distance Control mix CRTB mix
(inches)
r=3.5
r=3.5
r=0
a
r=0
a
0”
0.81 1.02 2.93 6.26
8”
1.83 0.07 0.77 4.26
12”
2.64 0.69 0.89 2.71
18”
4.18 2.17 3.83 0.18
24”
5.80 3.76 6.78 3.15
36”
9.52 7.53 12.43 9.32
48”
12.98 11.25 16.46 14.18
60”
16.19 14.81 18.96 17.74

1000

1000

(b) LAVAN using stress at r=3.5a

AS1*
AS2
AS3
AS4
AS5
AS6
AS7
AS8
LS1*
LS2
LS3
LS4
LS5
LS6
LS7
LS8

800

600

400

AS1*
AS2
AS3
AS4
AS5
AS6
AS7
AS8
LS1*
LS2
LS3
LS4
LS5
LS6
LS7
LS8

800
Surface deflection (micro-m)

(a) LAVAN using stress at r=0

Surface deflection (micro-m)

Average deflection
history error (%) PEavg
Control mix CRTB mix
r=3.5
r=3.5
r=0
a
r=0
a
1.61 2.09 1.45 2.06
1.47 1.86 1.29 1.75
1.39 1.67 1.30 1.52
1.31 1.48 1.56 1.30
1.45 1.35 1.99 1.43
1.82 1.53 2.58 2.04
2.20 1.87 2.46 2.25
2.27 2.10 1.42 1.77

600

400

200

200

0

0
0

5

10

15

20

25

30

35

40

0

5

10

15

20

25

30

35

Time (msec)

Time (msec)

Figure 25: Surface deflection comparison of ABAQUS and LAVAN for the Control mix
(*AS1= ABAQUS sensor 1; *LS1=LAVAN sensor 1).

66

40

1200

1200

(a) LAVAN using stress at r=0

AS1*
AS2
AS3
AS4
AS5
AS6
AS7
AS8
LS1*
LS2
LS3
LS4
LS5
LS6
LS7
LS8

800

600

400

AS1*
AS2
AS3
AS4
AS5
AS6
AS7
AS8
LS1*
LS2
LS3
LS4
LS5
LS6
LS7
LS8

1000
Surface deflection (micro-m)

Surface deflection (micro-m)

1000

(b) LAVAN using stress at r=3.5a

200

800

600

400

200

0

0
0

5

10

15
20
25
Time (msec)

30

35

40

0

5

10

15

20

25

30

35

40

Time (msec)

Figure 26: Surface deflection comparison of ABAQUS and LAVAN for the CRTB mix
(*AS1= ABAQUS sensor 1; *LS1=LAVAN sensor 1).
3.6 Chapter Summary
Algorithms have been developed to model multilayer viscoelastic forward solution. The
models presented can consider the unbound granular material as both linear elastic as well as
nonlinear-stress dependent material. Depending on the assumed unbound granular material
property, two generalized viscoelastic flexible pavement models were developed. The developed
forward model for linear viscoelastic AC and elastic unbound layers has been referred to as
LAVA. The developed forward model for linear viscoelastic AC and nonlinear elastic unbound
layers have been termed as LAVAN. LAVA assumes a constant temperature along the depth of
the AC layer. The algorithm was subsequently modified for temperature profile in the AC layer
and has been referred to as LAVAP. The models have been validated using FE based solutions.

67

CHAPTER 4

BACKCALCULATION USING VISCOELASTICITY BASED
ALGORITHM
4.1 Introduction
Typically, a load-displacement history of 60 milliseconds is recorded in an FWD test
(which constitutes of 25-35 msec of applied load pulse). This can give only a limited information
about the time varying E(t) behavior of the AC layer. However, in theory, it should be possible to
obtain the two sought after functions (i.e., E(t) and aT(T)). Two different approaches have been
discussed to obtain comprehensive behavior of asphalt: (i) using series of FWD deflection time
histories at different temperature levels (Varma et al., 2013b) and (ii) using uneven temperature
profile information existing across the thickness of the asphaltic layer during a single or multiple
FWD drops deflection histories (Varma et al., 2013a). Backcalculation of pavement properties
using FWD data involves developing an optimization scheme. The analysis is based on
formulating an objective function, which is minimized by varying the pavement properties.
Response obtained from the forward analysis is compared with response obtained from the FWD
test and the ‘difference’ is minimized by adjusting the layer properties of the system until a best
match is achieved. Typically the existing backcalculation methods either use Root Mean Square
(RMS) or percentage error of peak deflections as the objective function. However, since the
viscoelastic properties are time dependent, the entire deflection history needs to be used. Hence,
the primary component of the proposed backcalculation procedure is a layered viscoelastic
forward solution. Such solution should provide accurate and rapid displacement response
68

histories due to a time-varying (stationary) surface loading. Herein, for linear viscoelastic
pavement model, we used the computationally efficient layered viscoelastic algorithm
LAVA/LAVAP (refer Chapter 3) to support the backcalculation algorithm called
BACKLAVA/BACKLAVAP. Whereas, for viscoelastic nonlinear pavement model, we used the
computationally efficient layered viscoelastic algorithm LAVAN (refer chapter 3) to support the
backcalculation algorithm called BACKLAVAN.

10

5

Relaxation Modulus E(t) (MPa)

(a)

10

4

10

3

o

-10 C
o

10

10 C

2

o

21 C
o

37 C

10

o

54 C

1

10

-6

10

-4

10

-2

0

10

2

4

10

6

10

10

Frequency (Hz)

logE (t)  c1 

(b)

c2
1  exp(c3  c4 log(tR ))

10

4

(c)

log( aT (T))  a1(T 2-T02 )  a2(T-T0 )

4

10

102

T

Shift factor, a (T)

Relaxation Modulus, E(t) (MPa)

5

10

3

10

o

-10 C
o

10 C

2

10

o

21 C

10

0
o

-10 C
o

10

-2

o

21 C

o

o

37 C

37 C

o

54 C

1

10

-7

10

-5

10

10 C

-3

10

-1

10

1

3

10

10

5

10

10

7

10

-4

-10

o

54 C

0

10

20

30

40
o

o

Temperature ( C)

Reduced Time at 19 C (sec)

Figure 27: Figure illustrating steps in E(t) master curve development.

69

50

60

Whenever mechanical properties are derived by way of inverse analysis, it is desirable to
minimize the number of undetermined parameters by using an economical scheme. Such
approach is both advantageous from a computational speed perspective and it also addresses the
non-uniqueness issue, as test data may not be detailed, accurate, or precise enough to allow
calibration of a complicated model. Moreover, it is beneficial to have some inherent ‘protection’
within the formulation, forcing the analysis to a meaningful convergence - fully compliant with
the physics of the problem. Therefore, as discussed before, the E(t) master curve (Figure 27) is
initially assumed to follow a sigmoid shape with the Equation (25). Where, as discussed before,
aT(T) is the shift factor coefficient, which is a second order polynomial function of temperature
(T) and t is time (Equation (23)). Similar to the process used in generating |E*| master curve, E(t)
master curve can be generated using the concept of TTS. As shown by the relaxation modulus
and shift factor equations, total six coefficients are needed to develop the E(t) master curve,
including the temperature dependency (i.e., the shift factor coefficients). In theory, it should be
possible to obtain these six coefficients in two ways: (i) using two or more FWD time history
data at different temperature levels and (ii) using uneven temperature profile information existing
across the thickness of the asphaltic layer during a single FWD drop containing time changing
response data.
The objective function, which is based on deflection differences in the current work, is a
multi-dimensional surface that can include many local minima. In elastic backcalculation
methods, the modulus of the AC layer is defined using a single value. However, in the present
problem the AC properties are represented by a sigmoid containing four parameters for E(t) and
by a polynomial containing 2 parameters for aT(T). Hence, it is naturally expected that the
probability of number of local minimums will increase. In traditional methods, due to the
70

presence of multiple local minima, selection of different initial solution may lead to different
solutions. Reliability and accuracy of the backcalculated results depend on the optimization
technique used. In the present work several optimization techniques were tried to formulate a
procedure to backcalculate these six viscoelastic properties along with unbound material
properties. These optimization techniques can be broadly classified as traditional methods and
non-traditional methods. Traditional optimization methods could be computationally less
expensive and can have better control over solution accuracy. However, typically traditional
optimization methods (such as the “fminsearch” in MATLAB©) have the following issues:
1. Solution (local or global) may depend on initial seed values
2. Convergence time may vary depending on the initial seed values.
To extract benefits from both traditional and non-traditional methods, in this study they
have been hybridized to develop a more effective backcalculation procedure. Simplex-based
traditional optimization method was performed using MATLAB function “fminsearch”, whereas,
genetic algorithm-based evolutionary optimization method was performed using MATLAB
function “ga”. It is important to develop a backcalculation process such that FWD data obtained
at relatively small range of pavement temperatures can be sufficient to derive the viscoelastic
properties of AC. Among various optimization techniques, Genetic Algorithm (GA) was chosen
because of its capability to converge to a unique global minimum solution, irrespective of the
presence of local solutions (Fwa et al., 1997; Alksawneh, 2007; and Park et al., 2009). Through
guided random search from one ‘generation’ to another, GA minimizes the desired objective
function. The detailed discussion on genetic algorithm in FWD backcalculation is presented in
chapter 2 and will not be repeated here for brevity.
71

Formulation of the optimization model using GA is:
Objective:

m

n

k 1

i , o 1

Er  100 

d

k
i

 dok
dik


	 	

	

	

	

	 54 	

Subject to the following constraints:

c1l  c1  c1u ,	 c2l  c2  c2u ,	

c3l c3 c3u ,	 c4l  c4  c4u 	
a1l  a1  a1u , a2l  a2  a2u ,	

Ebl  Eb  Ebu ,	 Esl Es Esu 	
where, m is the number of sensors, di is the input deflection information obtained from field at
sensor k, dok is the output deflection information obtained from forward analysis at sensor k, n is
the total number of deflection data points recorded by a sensor, ci are the sigmoid coefficients, Eb
and Es are base and subgrade modulus and ai are the shift factor polynomial coefficients. The
superscript l represents lower limit and u represents upper limit.
In order to obtain the lower and upper limits of ci and ai, values of sigmoid and shift
factor coefficients of numerous HMA mixtures were calculated. Table 9 shows these limits
which were used in the GA constraints shown in Equation (54). Limits to the elastic modulus
were arbitrarily selected (based on typical values presented in the literature). It should be noted
that sigmoid obtained by using the lower limits or the upper limits of the coefficients give larger
range as compared to the actual range of E(t). This can potentially slow down the
72

backcalculation process. Therefore, as described later in the report, additional constraints were
defined to narrow down the search window. Figure 28 shows the steps in the GA based
backcalculation programs.
Table 9: Upper and lower limit values in backcalculation
Limit
Lower
Upper

c1
0.045
2.155

c2
1.80
4.40

c3
-0.523
1.025

c4
-0.845
-0.380

a1
-5.380E-04
1.136E-03

a2
-1.598E-01
-0.770E-01

E1
10000
13000

E2
22000
28000

Start

Limits on unknowns
Initial pool of solution
[C1,C2,C3,C4,a1,a2,Eb,Es]

Population size

Fitness evaluation & 
selection
Selected Parents

Variation & mutation

New generation
(offsprings)

No

No of generations>15

Yes
End

Figure 28: Backcalculation flow chart in BACKLAVA

73

Pavement info:
Layer thickness, 
Poission’s Ratio, 
AC temperature
Measured FWD:
stress history 
deflection history

4.2 Backcalculation of Relaxation Modulus Master Curve Using Series of FWD Tests Run
at Different Temperatures
The duration of a single pulse of an FWD test is very short, which limits the portion of
the E(t) curve used in the forward calculation using LAVA. As a result, it is not possible to
backcalculate the entire E(t) curve accurately using the deflection data of such short duration.
The longer the duration of the pulse, the larger portion of the E(t) curve is used in LAVA in the
forward calculation process. Therefore, one may conclude that FWD tests need to produce longduration deflection-time history. However, owing to the thermorheologically simple behavior of
AC, time-temperature superposition principle can be used to obtain longer duration data by
simply running the FWD tests at different temperatures and using the concept of reduced time.
This process has been further illustrated in Figure 29. The figure shows a schematic example of
the steps involved in the evaluation of a Candid parent solution (refer to step two “Fitness
evaluation and selection” in Figure 28). The figure shows simulation of two different FWD runs
at temperatures T1 and T2 such that T2 > T1. Since the two tests are at different temperatures, they
correspond to different reduce time in the master curve. For T2 > T1, FWD test at T2 would use
E(t) from higher reduced time compare to FWD test at T1. Before getting into the details of the
required number of FWD test temperatures and magnitudes, first, an analysis on the effects of
different FWD deflection sensor data on the backcalculated E(t) master curve is presented in the
following section.

74

3

10

c2
logE (t )   c1 
1  exp( c3  c 4 log(t R ))

log( aT (T))  a1(T 2 -T02 )  a2(T-T0 )

2

10

4

10

1

Shift factor

Relaxation Modulus (MPa)

5

10

3

10

2

10

1

10
-7
-5
-3
-1
1
10 10 10 10
10

3

10

o

10

5

10

7

10

0

10
10

-1

10

-2

10

-3

10

-4

-10

0

Reduced time at 19 C (sec)

10

20

30

10

4

10

3

10

2

60

FWD at temperature T2
10

Maximum E(t) used at 
temperature T1

1

10
-7
-5
-3
-1
10 10
10 10

Maximum E(t) used at 
temperature T2

4

10

3

10

2

10

1

10

1

10

o

3

10

5

10

10

7

-7

10

-5

10

-3

10

-1

10

1

3

10

10

5

10

7

10

o

Reduced time at 19 C (sec)

Reduced time at 19 C (sec)

0.03

0.012

Calculated 
deflections at 
temperature T1

0.008

Calculated 
deflections at 
temperature T2

0.025

Deflection (in)

0.01

Deflection (in)

50

5

5

Relaxation Modulus (MPa)

Relaxation Modulus (MPa)

FWD at temperature T1
10

40

o

Temperature C

0.006

0.004

0.02
0.015
0.01
0.005

0.002

0

0
0

0.005 0.01 0.015 0.02 0.025

0

0.03 0.035 0.04

0.005

Time (sec)

0.01

0.015

0.02

0.025

Time (sec)

0.03

0.035

0.04

Fitness evaluation=
Error in calculated deflection history at T1 +Error in calculated deflection history at T2 

Figure 29: Schematic of fitness evaluation in BACKLAVA

75

4.2.1 Sensitivity of E(t) Backcalculation to the use of Data from Different FWD Sensors
This section presents an analysis of contribution of individual and group of sensors on the
backcalculation of E(t) master curve. The backcalculation process was run using a populationgeneration of 70 and 15 respectively (selected after trying various combinations), using FWD
time histories obtained at a temperature set of {0, 10, 20, 30, 40, 50, 60, 70 and 80}oC. The
pavement properties used (Table 10) have been kept same throughout the study.
Table 10: Pavement Properties in viscoelastic backcalculation of optimal number of sensors
Property
Thickness (AC followed by granular layers)
Poisson ratio {layer 1,2,3…}
Eunbound {layer 2,3…}
E(t) sigmoid coefficient {layer 1}
aT(T) shift factor polynomial coefficients {layer 1}
Sensor spacing from the center of load (inches)

Case 1
10, 20 (in)
0.35, 0.3, 0.45
11450, 15000 (psi)
0.841, 3.54, 0.86, -0.515
4.42E-04, -1.32E-01
0, 8, 12,18,24, 36,48, 60

Convergence was evaluated based on the backcalculated modulus of base and subgrade layers as
well as the E(t) curve of the AC layer. Average error in the modulus of base and subgrade are
defined as follows:

 E act  E bc
 E act

 unbound  


 x100 	


	

	

	

55 	

where  unbound is the absolute value of the error in the backcalculated unbound layer modulus,
Eact and Ebc are the actual and backcalculated modulus (of the unbound layer), respectively. The
variation of error in the backcalculated E(t) at different reduced times is defined as follows:

 AC (ti ) 

Eact (ti )  Ebc (ti )
x100
Eact (ti )
	 	

76

	

	

	 56 	

where  AC (ti ) is the E(t) error at reduced time ti, where i ranges from 1 to n such that t1 = 10-8
and tn = 108 sec, n is the total number of discrete points on E(t) curve, Eact(ti) is the actual E(t)
value at point i, and Ebc(ti) is the backcalculated E(t) value at i. Finally, average error in E(t) is
defined as follows:



avg
AC

1 n

    AC (ti ) 
n  i 1
	 	

	

	

	 57 	

avg
where AC
is the average error in the E(t) of the AC layer. Figure 30 shows the variation of

 unbound when data from different FWD sensors are used. As shown, the error decreases as data
from farther sensors are incorporated in the backcalculation. This may be because at farther
sensors the deflections are primarily, if not solely, due to deformation in the lower layers.
14
12
Layer 2
Layer 3

unbound (%)

10
8
6
4
2
0
1

1-2

1-3

1-4

1-5

1-7

1-8

1-9

Range of sensors

Figure 30: Error in unbound layer modulus in optimal number of sensor analysis
Figure 31a shows the actual and backcalculated E(t) curve, which is only based on
deflection history obtained from sensor 1 (at the center of load plate). As shown, there is a very
good match between the backcalculated and actual curves. Figure 31b shows the variation of
77

percentage error in E(t) with time, calculated using Equation (57). The magnitude of percent
error ranges from about -9% to 23% and increases with reduced time. This is expected since the
E(t) in longer durations (>10-6 sec) are not used in the forward computations. It is noted that the
result is shown over a time range from 10-8 to 108 sec. However, the forward calculations were
actually made using temperatures ranging from 0 to 80oC, which corresponds to a reduced time
range of approximately from 10-6 to 106 sec.
In order to investigate if using just the farther sensors improves the backcalculated
Eunbound values, backcalculations were performed using data from different combinations of
farther sensors. Figure 32 shows the error in backcalculation of modulus of base (Layer 2) and
subgrade (Layer 3), when data from only further sensors are used. As shown, for Layer 3, error
ranges between 0.27 to 1.43%, with no specific trend. The error in the modulus of base (Layer 2)
is higher, ranging from 1 to 8.96%. However, a clear trend was not observed. Compared to
Figure 30, one can conclude that using all the sensors produces the least error in Eunbound.

78

Relaxation modulus, E(t) (psi)

10

7

Actual E(t)
Backcalculated E(t)
10

6

10

5

10

4

10

3

10

(a)

-9

10

-7

-5

10

10

-3

10

-1

1

10

10

3

10

5

7

10

10

9

o

Reduced time at 19 C (sec)
25
Entire backcalculated E(t)
E(t) used in backcalculation
20

(b)

AC i

 (t )

15

10

5

0

-5

-10
-9
10

10

-7

10

-5

10

-3

10

-1

10

1

10

3

5

10

10

7

10

9

o

Reduced time, 19 C
Figure 31: Backcalculation results using only sensor 1: (a) Backcalculated and actual E(t)
master curve at the reference temperature of 19oC (b) Variation of error,  AC (ti ) .

79

10
Layer 2
Layer 3



unbound

(%)

8

6

4

2

0
5-9

6-9

7-9

8-9

9

Range of sensors

Figure 32: Error in unbound layer modulus using FWD data only from farther sensors
4.2.2 Effect of Temperature Range of Different FWD Tests on Backcalculation
It is typically not feasible to run the FWD test over a wide range of temperatures (e.g.,
from 0oC to 80oC). Depending on the region and the month of the year, the variation of
temperature in a day can be anywhere between 10oC and 30oC during the Fall, Summer and
Spring when most data collection is done. This means that the performance of the
backcalculation algorithm needs to be checked for various narrow temperature ranges. The
purpose of the study explained in this section was to determine the effect of different temperature
ranges on the backcalculated E(t) values. Further it was realized that the results obtained from
GA may not be exact but only an approximation of the global solution. Hence a local search
method was carried out through “fminsearch” using the results obtained from GA as seed. Figure
33 shows the error in the backcalculated elastic modulus values of base and subgrade when
different pairs of temperatures are used. As shown, in most cases, the error was less than 0.1 %.
It is noted that the errors shown in Figure 33 are less than the ones shown in Figure 30 (when all
sensors are used). This is because in Figure 30, only GA was used, whereas in Figure 33,
“fminsearch” is used after the GA, which improved the results.

80

0.3

0.25

unbound (%)

Layer 2
Layer 3
0.2

0.15

0.1

0.05

0
0-10

10-20

10-30

30-40

40-50

50-60

o

Temperature, C

Figure 33: Variation of error in backcalculated unbound layer moduli when FWD data run
at different sets of pavement temperatures are used.
avg
Figure 34 shows the average error in E(t) (i.e., AC
given in Equation (57)), where a

pattern was observed. The error was the least when intermediate temperatures (i.e., {10-30},
{10-20-30}, {20-30-40}, {30-40}, {30-40-50}oC) were used. At low temperatures, it seems the
error increase. This is meaningful because at low temperatures, a small portion (upper left in
Figure 31) is utilized in BACKLAVA. Therefore, the chance of mismatch at the later portions of
the curve (lower right in Figure 31) is high. At high temperatures, error also seems to increase.
Theoretically, the higher the temperature, the larger portion of the E(t) curve is used because of
the nature of the convolution integral, which starts from zero. However, if, only the high
temperatures are used, discrete nature of load and deflection time history leads to a big ‘jump’
from zero to the next time ti, during evaluation of the convolution integral. This is because when
the physical time at high temperatures is converted to reduced time, actual magnitudes become
large and, in a sense, a large portion at the upper left side of the E(t) curve is skipped during the
convolution integral. At intermediate temperatures, however, a more ‘balanced’ use of E(t) curve
in BACKLAVA improves the results.

81

80
70

50

AC

40



avg

(%)

60

30
20
10
0
0-10

10-20

10-30

10-20-30 20-30-40

30-40

30-40-50

40-50

40-50-60

50-60

50-60-70

o

Temperature, C

Figure 34: Error in backcalculated E(t) curve in optimal backcalculation temperature set
analysis minimizing percent error
When results from GA were used as seed values in fminsearch, it was observed that in general
error in E(t) was reduced. Figure 35, Figure 36 and Figure 37 show backcalculated E(t) master
curve using GA and corresponding backcalculated E(t) master curves obtained using GA +
“fminsearch”. As shown, combined use of GA and “fminsearch” results in improved
backcalculation. Table 11 shows the time it takes to run the genetic algorithm for populationgeneration size of 70-15, followed by “fminsearch”. The results are shown for a computer that
has Intel Core 2, 2.40 GHz, 1.98GB RAM.
Table 11: Backcalculation run time for ga-fminsearch seed Runs
Number of
Temperature data
Backcalculation time

Two
(e.g. {10, 30}oC)
30 min

82

Three
(e.g. {10, 20, 30}oC)
40 min

Seed
Run
(“fminsearch”)
15-20 min

107
Actual E(t)
Backcalculated E(t)

Relaxation modulus, E(t) (psi)

Relaxation modulus, E(t) (psi)

107

6

10

105

4

10

103
10-9

10-7

10-5

10-3

10-1

101

103

105

107

Actual E(t)
Backcalculated E(t)
6

10

105

4

10

103
10-9

109

10-7

o

10-5

10-3

10-1

101

103

105

107

109

o

Reduced time at 19 C (sec)

Reduced time at 19 C (sec)

a. Backcalculated E(t) curve using GA only

c. Backcalculated E(t) curve using GA +
fminsearch
5

30
Entire backcalculated E(t)
E(t) used in backcalculation

0

20

-5
-10

AC i

AC i

 (t )

 (t )

10

0

-15
-20
-25

-10

-30
-20
-9
10

10

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

-35
-9
10

9

10

o

10

-7

10

-5

10

-3

10

-1

10

1

10

3

10

5

10

7

10

9

o

Reduced time, 19 C

b. Variation of error,  AC (ti ) : GA only.

Entire backcalculated E(t)
E(t) used in backcalculation

Reduced time, 19 C

d. Variation of error,  AC (ti ) : GA+fminsearch

Figure 35: Results for backcalculation at {10, 30}oC temperature set: (a) and (b) Only GA
is used, (c) and (d): GA+fminsearch used.

83

107
Actual E(t)
Backcalculated E(t)

Relaxation modulus, E(t) (psi)

Relaxation modulus, E(t) (psi)

107

6

10

105

4

10

103
10-9

10-7

10-5

10-3

10-1

101

103

105

107

Actual E(t)
Backcalculated E(t)
6

10

105

4

10

103
10-9

109

10-7

o

10-5

10-3

10-1

101

103

105

107

109

o

Reduced time at 19 C (sec)

Reduced time at 19 C (sec)

a. Backcalculated E(t) curve using GA only

c. Backcalculated E(t) curve using
GA+fminsearch
25

15
Entire backcalculated E(t)
E(t) used in backcalculation

Entire backcalculated E(t)
E(t) used in backcalculation

10

20

5

0

AC i

 (t )

AC i

 (t )

15

-5

10

5
-10
0

-15

-20
-9
10

10

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

-5
-9
10

9

10

o

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

9

10

o

Reduced time, 19 C

b. Variation of error,  AC (ti ) : GA only

10

Reduced time, 19 C

d. Variation of error,  AC (ti ) : GA+fminsearch

Figure 36: Results for backcalculation at {30, 40}oC temperature set

84

107
Actual E(t)
Backcalculated E(t)

6

10

105

4

10

103
10-9

10-7

10-5

10-3

10-1

101

103

105

107

Actual E(t)
Backcalculated E(t)

Relaxation modulus, E(t) (psi)

Relaxation modulus, E(t) (psi)

107

6

10

105

4

10

103
10-9

109

10-7

o

10-5

10-3

10-1

101

103

105

107

109

o

Reduced time at 19 C (sec)

Reduced time at 19 C (sec)

a: Backcalculated E(t) curve using GA only

c: Backcalculated E(t) curve using GA +
fminsearch

40

1.5
Entire backcalculated E(t)
E(t) used in backcalculation

Entire backcalculated E(t)
E(t) used in backcalculation
1

30
0.5

AC i

 (t )

AC i

 (t )

20

10

0

-0.5

-1
0
-1.5

-10
-9
10

10

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

-2
-9
10

9

10

o

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

9

10

o

Reduced time, 19 C

b: Variation of error,  AC (ti ) : GA only

10

Reduced time, 19 C

d: Variation of error,  AC (ti ) : GA+ fminsearch

Figure 37: Results for backcalculation at {30, 40, 50}oC temperature set
4.2.3 Normalization of Error Function (Objective Function) to Evaluate Range of
Temperatures
In the analysis presented in the previous sections, percentage error between the computed
and measured displacement was used as the minimizing error. But deflection curve obtained
from the field often includes noise, especially after the end of load pulse. If percentage error is
used as the minimizing objective, this may lead to over emphasis of lower magnitudes of
deflections at the later portion of the time history, which typically includes noise and integration
85

errors. Hence another fit function was proposed in which, the percentage error was calculated
with respect to the peak of deflection at each sensor. This penalizes the tail data by normalizing it
with respect to the peak.

d  d 
Er   100 
d  	
m

n

k 1

i , o 1

k
i

k
o

k

	

	

	

	 58 	

max

where, d k max is the peak response at sensor k, m is the number of sensors, d ik is the measured
deflection at sensor k, d ok is the output (calculated) deflection from forward analysis at sensor k,
n is the total number of deflection data points recorded by a sensor. The limits considered for
E(t) so far were the limits on the individual parameters of the sigmoid curve (Table 9). The E(t)
curves obtained by considering the upper and lower limits of the parameters represent curves
well beyond the actual data base domain. To curtail this problem, constraints were introduced
putting limit on sum of the sigmoid coefficients c1 and c2 as c1  c2  s1 and c1  c2  s2 ,where s1
and s2 are arbitrary constants, which physically represents the maximum and minimum
instantaneous modulus value of any HMA mix respectively. The arbitrary constants s1 and s2
were obtained by calculating maximum and minimum values of the sum of sigmoid coefficients
c1 and c2 from numerous HMA mixes. Alternatively the problem was reframed by incorporating
the constraints in limit form by redefining the variables as
	1
	2

;	

	
	

;
	3

;	

	

	4

;	

	

86

	

	

	

	 59 	

The problem was then resolved after replacing the inequality constrain with limits on the
variables. Figure 38 and Figure 39 show the search domains in unconstrained and constrained
optimization models. The unconstrained search zone comprise of many unrealistic solutions,
which means the GA would require larger population and generations to search the entire
domain, leading to high computation time. It can be seen from the figures that the search domain
significantly reduces when the constraint is added to the model. The new function gave good
results at temperature sets of {10, 30}oC, {30, 40}oC, {10, 20, 30}oC, {20, 30, 40}oC and {30,
40, 50}oC. The backcalculated E(t) curves were then converted to |E*| using the interconversion
relationship given in Equation (11) and (12). Backcalculated E(t) and |E*| master curves are
compared with the actual curves for temperature sets {10, 30}oC and {10, 20, 30}oC in Figure 40
and Figure 41 respectively.

Relaxation modulus, E(t) (psi)

10

9

Minimum
Maximum

108

Search domain

107
106
105

Search domain

104
10

3

10

2
-9

10

10

-7

10

-5

10

-3

10

-1

10

1

10

3

10

5

Reduced time at 19oC (sec)
Figure 38: Search domain for unconstrained constrain
87

10

7

9

10

Relaxation modulus, E(t) (psi)

10

9

10

8

Minimum
Maximum

Search domain

107
106
105

Search domain

104
10

3

10

2
-9

10

10

-7

10

-5

10

-3

10

-1

10

1

10

3

10

5

10

7

9

10

Reduced time at 19oC (sec)

Dynamic modulus, |E*| (psi)

Figure 39: Search domain for constrained constrain

10

7

10

6

10

5

10

4

Actual |E*|
Backcalculated |E*| 10,30C
10

Backcalculated |E*| 10, 20, 30C

3

10

-9

10

-7

-5

10

10

-3

10

-1

1

10

10

3

10

5

7

10

10

9

o

Reduced frequency, at 19 C (Hz)

Figure 40: Backcalculated lE*l master curve using FWD data at temperature set {10, 30}oC
minimizing normalized error

88

Relaxation modulus, E(t) (psi)

10

7

10

6

10

5

10

4

Actual E(t)
Backcalculated E(t) 10,30C
10

Backcalculated E(t) 10, 20, 30 C

3

10

-9

10

-7

-5

10

10

-3

10

-1

1

10

10

3

10

5

7

10

10

9

o

Reduced time, at 19 C (sec)

Figure 41: Backcalculated E(t) master curve using FWD data at temperature set {10, 20,
30}oC minimizing normalized error
It can be seen from Figure 42 that the results obtained for E(t) errors over temperature sets
distinctly showed a pattern. The E(t) and Eunbound errors with respect to temperature sets showed
trend similar to what was observed in case of percentage error (see Figure 33 and Figure 34).
The error was observed to be high at sets of low {0-10}oC and high {50-60}oC, {40-50-60}oC,
{50-60-70}oC temperatures. This is because the backcalculated E(t) at lower temperatures
represents the left portion of the sigmoidal E(t) curve and higher temperatures represents the
right. As explained earlier both the regions are fairly flat and hence represent constant values of
E(t), which may not optimize to the actual E(t) curve. Better results were obtained for the
temperature range of {10, 20}oC to {30, 40, 50}oC (Figure 42). The backcalculated E(t) master
curves and corresponding error obtained at {10, 30}oC and {20, 30, 40}oC for the proposed
backcalculation model are shown in Figure 44.

89

80
70

50

AC

40



avg

(%)

60

30
20
10
0
0-10

10-20

10-30

10-20-30 20-30-40 30-40-50

30-40

40-50

40-50-60

50-60

50-60-70

o

Temperature, C
avg
Figure 42: Variation of AC
at different FWD temperature sets using modified sigmoidal
variables.

5

unbound (%)

4
Layer 2
Layer 3
3

2

1

0
0-10

10-20

10-30

30-40

40-50

50-60

o

Temperature range C

Figure 43: Variation of  unbound at different FWD temperature sets using modified
sigmoidal variable.

90

107
Actual E(t)
Backcalculated E(t)

6

10

105

4

10

103
10-9

10-7

10-5

10-3

10-1

101

103

105

107

Actual E(t)
Backcalculated E(t)

Relaxation modulus, E(t) (psi)

Relaxation modulus, E(t) (psi)

107

6

10

105

4

10

103
10-9

109

10-7

o

10-5

10-3

10-1

101

103

105

107

109

o

Reduced time at 19 C (sec)

Reduced time at 19 C (sec)

c. Backcalculated E(t) curve using GA (at
temperature {20, 30, 40}oC)

a. Backcalculated E(t) curve using GA (at
temperature {10, 30}oC)

1.5

10
Entire backcalculated E(t)
E(t) used in backcalculation

Entire backcalculated E(t)
E(t) used in backcalculation

1

8

0.5

6

AC i

 (t )

AC i

 (t )

0
4

2

-0.5
-1

0

-1.5

-2

-4
-9
10

-2

10

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

-2.5
-9
10

9

10

o

-7

-5

10

-3

10

10

-1

10

1

3

10

10

5

10

7

9

10

o

Reduced time, 19 C

b. Variation of error,  AC (ti ) : (result at
temperature {10, 30}oC)

10

Reduced time, 19 C

d. Variation of error,  AC (ti ) : (result at
temperature {20, 30, 40}oC)

Figure 44: Backcalculation results obtained using modified sigmoid variables
4.2.4 Backcalculation of Viscoelastic Properties using Various Asphalt Mixtures
In the previous sections, analyses were performed using only a single mix. In order to
verify the conclusions made in the previous sections regarding the optimum range of
temperatures of FWD testing, backcalculations have been performed on 9 typical mixtures.
Actual viscoelastic properties: relaxation modulus and shift factors of the selected mixtures are
shown in Figure 45. Comparison of average error in the backcalculated relaxation modulus
91

function calculated over four time ranges: (a) 10-5 to 1 sec (b) 10-5 to 102 sec and (c) 10-5 to 103
sec are shown in Figure 46. It can be seen from the figure that, for all the mixes, relaxation
modulus curve can be predicted close to less than 15% over a range of relaxation time less than
10+3 seconds. Furthermore, it can be seen that, as suggested, the backcalculated relaxation
modulus prediction provides a good match over an approximate temperature range of 10oC to
30oC.

7

Relaxation modulus, E(t) (psi)

10

106

5

10

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

104

3

10

10-9

10-7

10-5

10-3

10-1

101

103

105

107

109

o

Reduced time at 19 C (sec)
103
ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

2

10

Shift factor, aT

101
100
-1

10

10-2
10-3
10-4
-10

0

10

20

30

40

50

60

o

Temperature ( C)

Figure 45: Viscoelastic properties of field mix in optimal temperature analysis: (a)
Relaxation modulus at 19oC, (b) Time-temperature shift factor.
92

35

25
20

30

Average Error (%)

Average Error (%)

30

35

(a)

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

15
10

25
20

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

15
10

5

5

0

0
10,20

20,30

30,40

40,50

10,20,30

o

20

35

(b)

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

30

Average Error (%)

Average Error (%)

25

15
10

25
20

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

15
10

5

5

0

0
10,20

20,30

30,40

40,50

10,20,30

o

35

(c)

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

30

Average Error (%)

Average Error (%)

20

30,40,50

Temperature set C

35

25

20,30,40
o

Temperature set C

30

30,40,50

Temperature set C

35
30

20,30,40
o

Temperature set C

15
10

25
20

ALF Control 70-22
SBS 64-40
Air Blown
PPA+Elvaloy
CRTB
WAM Foam
SBS LG
Terpolymer
PG6422 12GTR

15
10

5

5

0

0
10,20

20,30

30,40

40,50

10,20,30

o

20,30,40

30,40,50
o

Temperature set C

Temperature set C

Figure 46: Variation of error,  AC (ti ) : (a) ti = 10-5 to ti = 1 sec used in  AC (ti ) , (b) ti = 10-5
to ti = 10+2 sec used in  AC (ti ) and (c) ti =10-5 to ti = 10+3 sec used in  AC (ti ) computation.
93

4.3 Backcalculation of Relaxation Modulus Master Curve using a Single FWD Test and
Known Pavement Temperature Profile
The uneven temperature profile existing across the thickness of the asphaltic layer during
a single FWD drop can theoretically be used to backcalculate the E(t) master curve and the shift
factor coefficients (aT(T)). The AC layer may be divided into several sublayers of same
viscoelastic properties, with different temperature levels. This process has been further illustrated
in Figure 47. The figure shows a schematic example of steps involved in the evaluation of a
candid parent solution (refer to step two “Fitness evaluation and selection” in Figure 28). The
figure shows simulation of a single FWD ran under a temperature profile {T1, T2 and T3} such
that T3 > T2 > T1. The response is expected to be a resultant of AC behavior from all the three
temperatures (which would corresponds to three different reduce time in the master curve). Two
different approaches of backcalculation have been discussed in the present section. In the first
approach all the unknown variables (sigmoid coefficients, shift factor coefficients and unbound
modulus) in the forward algorithm were varied during backcalculation. Whereas in the second
approach, a two-staged backcalculation procedure was adopted. The two-stage method involved
elastic backcalculation in the first stage (unbound modulus assuming elastic AC layer) followed
by viscoelastic backcalculation in the second (sigmoid and shift factor coefficients). Both the
approaches were explored in the present study.

94

3

logE (t )   c1 

10

c2
1  exp( c3  c 4 log(t R ))

log( aT (T))  a1(T 2 -T02 )  a2(T-T0 )

2

10

4

10

1

Shift factor

Relaxation Modulus (MPa)

5

10

3

10

2

10

1

10
-7
-5
-3
-1
1
10 10 10 10
10

3

10

o

10

5

10

7

10

0

10
10

-1

10

-2

10

-3

10

-4

-10

0

10

Reduced time at 19 C (sec)

20

30

40

o

50

60

Temperature C

5

10

4

10

3

10

2

5

5

10

10

Maximum E(t) used at 
temperature T1

1

10
-7
-5
-3
-1
10 10
10 10

Maximum E(t) used at 
temperature T2

4

10

3

10

2

10

Relaxation Modulus (MPa)

10

Relaxation Modulus (MPa)

Relaxation Modulus (MPa)

FWD at temperature profile T1, T2, T3

1

10

o

3

10

5

10

7

10

3

10

2

10

1

1

10

Maximum E(t) used at 
temperature T3

4

10

-7

-5

10

10

-3

10

-1

10

1

3

10

10

5

10

7

10

10

-7

10

Reduced time at 19 C (sec)

-5

10

-3

10

-1

10

1

3

10

10

5

10

7

10

o

o

Reduced time at 19 C (sec)

Reduced time at 19 C (sec)

0.03

Calculated deflections at 
temperature profile T1, T2, T3

Deflection (in)

0.025
0.02
0.015
0.01
0.005
0
0

0.005

0.01

0.015 0.02

0.025 0.03

0.035

0.04

Time (sec)

Fitness evaluation=Error in calculated deflection history

Figure 47: Schematic of fitness evaluation in BACKLAVA
4.3.1 Linear Viscoelastic Backcalculation using Single Stage Method
As discussed earlier, total six coefficients are needed to represent the relaxation
properties of the AC layer, including the temperature dependency. The backcalculation
procedure used was same as used in the previous section (BACKLAVA), except the forward
95

analysis was replaced by LAVAP, which can consider varying temperature along the depth of
the AC layer. Subsequently the backcalculation algorithm was referred as BACKLAVAP. For
executing the GA, lower and upper limits of ci and ai, (sigmoid and shift factor coefficients) and
other specifications were retained the same.
As a first step, the backcalculation algorithm was validated with a synthetic FWD
deflection history, under two different temperature profiles. The data were generated using
LAVAP, and then used in BACKLAVAP for backcalculation of E(t). The AC layer was divided
into three equal sublayers with three different temperatures. Pavement section, properties and
temperature used in the forward analysis were as shown in Table 12.
Table 12: Details of the pavement properties used in single FWD test backcalculation at
known temperature profile
Asphalt Concrete Layer
Subgrade
Sublayer Sublayer Sublayer Granular Base
1
2
3
Thickness
51 mm
51 mm
51 mm
508 mm
Semi-inf
o
o
o
Case 1
20 C
15 C
10 C
N/A
N/A
Temperature
o
o
o
Case 2
30 C
25 C
20 C
N/A
N/A
Poisson’s Ratio
0.35
0.4
0.45
E(t) Coefficients (c1,c2,c3,c4)
Relaxation Modulus
Backcalculated Backcalculated
Backcalculated
Time-Temperature
N/A
N/A
(a1, a2) Backcalculated
Shifting Coefficients
Sensor spacing from the center of load (inches): 0, 8, 12, 18, 24, 36,48, 60
Property

For the case of backcalculation using temperature profile, the GA parameters, namely
population size and generation numbers were again selected after several trials of combinations.
It was observed that at population size of 300, improvement in the best solution was marginal
after 12 to 15 generations, and the population converged to the best solution at around 45
generations. Similarly for population size of 400, improvement in the best solution was marginal
after 10 to 15 generations. Figure 48 shows the backcalculation results at the temperature sets
96

given in Table 12, where a good match is visible. Error in the backcalculated E(t) was quantified
relative to the actual E(t) using  AC (ti ) given in Equation (56). The  AC (ti ) calculation was
avg
performed over a reduced time interval from 10-8 to 10+8 seconds. Then, the average error ( AC
)

was computed using Equation (57). The average error level for the first temperature profile was
found to be 5.2% and for the second: 4.4%.
In order to further investigate the effect of magnitude of the pavement temperature profile
on backcalculation of E(t) master curve, synthetic FWD deflection histories were generated. The
synthetic data was then used in backcalculation. The structure was divided into three layers with
different temperatures, and E(t) was backcalculated using these data. The pavement section
properties used in the study were same as shown in Table 12. Backcalculation was performed
assuming the temperature of the top, middle and bottom sublayers of the asphalt layer as {20, 15,
10}oC, {30, 25, 20}oC, {40, 35, 30}oC and {50, 45, 40}oC. It was again observed that the problem
converged well with 300 GA populations at 45 GA generations. Backcalculated E(t) and
deflection histories for temperature profile {20, 15, 10}oC and {30, 25, 20}oC are shown in
Figure 48. The results shown in Figure 49 did exhibit some trend, suggesting that there is a good
potential for backcalculation of E(t) using a single FWD response for the lower temperature
ranges, assuming that the temperature profile of the pavement is known.

97

o

o

o

o

Temperature = 20 C, 15 C, 10 C

0.02

o

Measured
Backcalculated

Measured
Backcalculated
0.015

Deflection, in

0.015

Deflection, in

o

Temperature = 30 C, 25 C, 20 C

0.02

0.01

0.005

0.01

0.005

0

0
0

0.005

0.01

0.015

0.02

0.025

0.03 0.035

0.04

0

0.005

0.01

Time (sec)

0.015

0.02

0.025

0.03 0.035

Time (sec)

(a) Case 1 deflection history

(b) Case 2 deflection history
7

Relaxation modulus, E(t) (psi)

10

6

10

105

4

10

Actual E(t)
Backcalculated case-1
103
10-9

Backcalculated case-2
-7

10

10-5

10-3

10-1

101

103

105

107

109

o

Reduced frequency, at 19 C (sec)

(c) Relaxation modulus
Figure 48: Comparison of actual and backcalculated values in backcalculation using
temperature profile.

98

0.04

E(t) avergae error (%)

70
60
50

Reduced time 10e-8 to 10e8
Reduced time 10e-5 to 10e5

40
30
20
10
0
20_15_10

30_25_20

40_35_30

50_45_40

o

Temperature set C

Figure 49: Error in backcalculated E(t) curve for a three-temperature profile

4.3.2 Backcalculation of the viscoelastic Properties of the LTPP Sections using a Single FWD
Test with Known Temperature Profile
The BACKLAVAP algorithm was next used with field data to backcalculate the
viscoelastic properties of three LTPP sections. The selection of the sites was done based on the
following rules: (a) Section comprising of three layers with only one AC layer, (b) total number
of constructions of the section to be one, (c) SPS section type (experiment number 1 or 8), (d)
flexible pavement and (e) presence of dynamics.
It was observed that the displacement readings obtained from the LTPP sections at all the
sensors (including the sensor under the load), showed time delay with respect to loading. For an
ideal multilayer elastic pavement, displacement peaks in all the sensors will be perfectly aligned
with load peak. Whereas for a viscoelastic multilayer pavement system, displacement peaks are
expected to be slightly misaligned due to viscoelasticity. However, in the measured FWD data
the displacement peaks are completely staggered. In the absence of measurement error, this
misaligned displacement peak is mainly due to dynamics. Typically this delay in displacement
peaks increases with increase in the sensor location distance. It was assumed that the lag in
99

displacement can be removed by shifting of displacement curves assuming first non-zero value
as the initial starting point of the displacement curves. The sections were checked for presence of
time delay using the following procedure:
1) Perform elastic backcalculation to obtain an approximate estimate of unbound layers
modulus values.
2) Perform forward calculation using creep converted E(t) and unbound modulus values
obtained in step 1.
3) Shift the measured and calculated deflection histories assuming first non-zero value as
the initial starting point.
4) If the peaks of the shifted curves do not reach at the same time then eliminate the section
from the analysis.
Figure 50 shows misaligned displacement peaks after shifting the displacement curves in LTPP
section 06A805. This misalignment after the shifting could be because of the following reasons:


Measurement error.



Error in D(t) measured values.



Presence of significant dynamics.



Erroneous temperature measurement in field.

100

0.04
Measured S1
Measured S3

0.035

Deflection (in)

0.03

Backcalculated S1
Backcalculated S3
Measured peak
for sensor1

Calculated peak
for sensor1

0.025
0.02
0.015
0.01
0.005
0

0

0.005

0.01

0.015

0.02

Time (sec)

0.025

0.03

 

Figure 50: Time lag in LTPP section 06A805 (shifted based on constant displacement) 
Table 13 and Table 14 contain general and structural information of the selected LTPP
site. As shown in Table 14, section 010101 has total four layers which include two AC layers.
However, since the D(t) of the two AC layers reported in the LTPP database are very close, the
section was included in the list, assuming the two layers as a single AC layer in the analysis.
However, it is not clear from the LTPP database whether the D(t) was measured before or after
the constructions were done.
Table 13: List of LTPP sections used in the analysis
State

Section

1
34
46

0101
0802
0804

Year of
Total no. of
Section Experiment
Test date
construction constructions
type
no.
4/30/1991
1
03/11/1993 SPS
1
SPS
1/1/1993
1
02/05/1997
8
SPS
1/1/1992
1
06/18/1993
8

In the LTPP program, each section is tested according to a specific FWD testing plan,
which consists of one or more test passes. Both SPS 1 and SPS 8 are tested along two test passes
(test pass 1 and test pass 3) using test plan 4 in LTPP. Test pass 1 data includes FWD testing
performed along the midlane (ML) whereas test pass 2 data includes FWD testing performed
along outer wheel (OW) path. Since testing along the midlane test pass represent the
101

axisymmetric assumption better, it was used here in the analysis. In LTPP testing protocol,
temperature gradient measurements are taken every 30 min, plus or minus 10 minutes. The
necessary temperature profile data was obtained by interpolating the temperature measured
during the FWD testing. The AC layer was divided into three equal sublayers and a constant
temperature for each sublayer has been estimated.
Table 14: Structural properties of the LTPP sections used in the analysis
State

Section

1
34
46

0101
0802
0804

Total
no. of
layers

No. of
AC
layers
4
3
3

AC layer
thickness
2
1
1

AC1=1.2, AC2=6.2
6.7
6.9

Base layer
thickness
7.9
11.6
12

The FWD deflection data in the selected LTPP sections showed no or minimal waviness
at the end of the load pulse, which indicated that there was no shallow stiff layer. The presence
of a stiff layer was further evaluated using a graphical method suggested by Ullidtz (1988). The
method involves plotting peak deflections obtained from FWD testing versus the reciprocal of
the corresponding sensor location (measured from the center of loading) (Appea, 2003). Depths
of stiff layer in each LTPP section estimated using Ullidtz’s method is shown in Table 15. It
should be noted that negative depth to stiff layer is interpreted as absence of stiff layer in the
method. The results indicate that stiff layers are generally deeper than 18 ft. It was suggested by
Lei (2011) that if the stiff layer is below 18ft, effect of dynamics is not observed on the surface
deflections.
Section properties used for elastic and viscoelastic backcalculations were the same (see
Table 13) except that the modulus of the AC layer in the elastic backcalculation was assumed
constant (modulus unknown). For elastic backcalculation, an in-house genetic algorithm was
102

developed. The Poisson’s ratio for AC, granular base and subgrade layers were assumed to be
0.3, 0.4 and 0.45, respectively.
Table 15: Depths of stiff layer in each LTPP section estimated using Ullidtz’s method
State

Section

1
34
46

0101
0802
0804

Depths of stiff layer from
surface (ft)
25.13 ft
185.61
78.58

As noted above, the results obtained from elastic backcalculation were used to define
bounds for base and subgrade modulus in BACKLAVAP. The elastic backcalculation results are
presented in Table 16. The static backcalculated base modulus values varied between 16157 psi
and 47202 psi, and the subgrade modulus values varied between 15789 psi and 36473 psi. Next,
the viscoelastic backcalculation was performed. The backcalculated unbound layer moduli for
the sections obtained from viscoelastic backcalculation are presented in Table 17. For the
viscoelastic backcalculation, the GA algorithm in BACKLAVAP utilized 300 populations in
each of the 15 generations.
Table 16: Elastic backcalculation results for LTPP unbound layers
State

Section

Test Date

1
34
46

0101
0802
0804

3/11/1993
09/08/1993
6/18/1993

Elastic modulus
(psi)
Ebase
Esubgrade
16157
35127
47202
36473
20877
17547

As shown in Table 17, the viscoelastic backcalculated base modulus values varied between
12117 psi and 39292, and the subgrade modulus values varied between 15299 psi and 34795 psi.

103

Table 17: Viscoelastic backcalculation results for LTPP unbound layers
State

Section

Test Date

1
34
46

0101
0802
0804

3/11/1993
09/08/1993
6/18/1993

Elastic modulus
(psi)
Ebase Esubgrade
12808 34362
39292 34795
15283 16876

In order to validate the backcalculated results, creep compliance data available in the
LTPP database were converted into relaxation modulus E(t). Creep data were available in
tabulated form at three temperatures: -10oC, 5oC, and 25oC, and seven different times: 1s, 2s, 5s,
10s, 20s, 50s and 100s. Assuming classical power law function for the creep compliance
(Equation (18)), the available data was fitted separately to each temperature. The associated
relaxation modulus was then obtained using the mathematically exact formula given in Equation
(19) (Kim, 2009). Finally, for comparison, dynamic modulus master curve was calculated from
the relaxation modulus via interconversion (Kutay et al., 2011). For further verification,
estimated dynamic modulus obtained from Artificial Neural Network (ANN) based model
ANNACAP is also compared. The ANNACAP predictive model developed under FHWA long
term pavement program research (Kim et al., 2011) uses artificial neural network based five
different models to predict dynamic modulus and time-temperature shift factor. Based on the
inputs used each model is a separately trained neural network. In the present study the estimated
dynamic modulus master curve and time-temperature shift factors obtained from ANNACAP are
based on MR model in ANNACAP. From the results it was found that the dynamic modulus
curves estimated using ANNACAP model, especially at higher frequencies, aggress well with
the dynamic modulus curves obtained through interconverted creep data. However, although the
dynamic modulus master curve predicted by ANNACAP matched well at the higher frequencies,
it typically predicted higher values at reduced frequencies less than 10-2 Hz.
104

The section property and interpolated temperature profile for LTPP section 010101 is
shown in Figure 51. The section was first evaluated to assess presence of any time delay. For this
the E(t) obtained from D(t) and the unbound modulus values obtained from elastic
backcalculation (Table 16) were used to calculate deflections. The forward calculated deflection
curves were then shifted assuming first non-zero value as the initial starting point. Assuming that
the E(t) obtained from D(t) and the elastic backcalculated modulus values approximately
represent the AC and the unbound layer properties, the peaks of the measured and the calculated
curves are expected to reach at the same time. The measured and calculated deflection histories
are shown in Figure 52. The curves were found to reach peak at the same time and hence the
section was chosen for further backcalculation analysis. Further it also indicates that the E(t)
obtained from D(t) approximates the viscoelastic properties of the AC layer. As shown in Figure
53a, for section 010101, the backcalculated relaxation modulus master curve and the creep
converted relaxation modulus matched well till 102 seconds. However, the curves showed
disagreement for time greater than 102 seconds. The backcalculated time-temperature shift
factors were compared with creep and ANNACAP computed results in Figure 53b. It can be seen
from the figure that the backcalculated time-temperature shift factor functions showed a good
match over the entire temperature range of 0oC to 50oC.

105

7.4"

7.9"

30.01oC

Asphalt Concrete (Linear Viscoelastic)
ν=0.35
Granular Base (Ebase)
ν=0.4

27.13oC
24.13oC

Subgrade (Esubgrade)
ν=0.45
∞ 
Figure 51: LTPP section 01-0101 cross section and temperature profile in AC layer

0.03

0.03

Measured S1
Measured S3
Calculated peak
for sensor1

0.02

Measured S1
Measured S3

0.025

Measured peak
for sensor1

Deflection in

Deflection in

0.025

Backcalculated S1
Backcalculated S3

0.015
0.01
0.005

Measured peak
for sensor1

Calculated peak
for sensor1

0.02

Backcalculated S1
Backcalculated S3

0.015
0.01
0.005

0
0

0.003

0.006

0.009

0.012

0.015

Time (sec)

0

0

0.003

0.006

0.009

Time (sec)

0.012

0.015

 

Figure 52: Comparison of measured and forward calculated deflection histories for section
010101 (a) Viscoelastic backcalculation (b) Forward calculated using E(t) obtained from
LTPP measured D(t) and elastic backcalculated modulus values.
As shown in Figure 53c, backcalculated, creep converted and ANNACAP dynamic modulus
curves matched well over frequencies greater than 10-3 Hz. A better agreement was observed
between ANNACAP and backcalculated results at frequencies less than 10-3 Hz.

106

7

10

6

10

5

10

4

10

3

10

2

Backcalculation-1
Backcalculated-2
From D(t)
ANNACAP

3

Backcalculated-1

10

Backcalculated-2

Shift factor, a (T)

From D(t)

101

T

Relaxation modulus, E(t) (psi)

10

10-1
-3

10

-5

10

-5

10

-4

10

-3

10

-2

10

-1

10

0

10

1

10

2

10

10

3

5

10

15

20

Dynamic modulus, |E*| (psi)

(a) Relaxation modulus
10

7

10

6

10

5

10

4

10

3

10

2

10

25

30

35

40

45

Temperature (oC)

Reduced time at 19oC (sec)

(b) Shift factor

Backcalculated-1
Backcalculated-2
From D(t)
ANNACAP
-3

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

Reduced frequency at 19 oC (Hz)

(c) Dynamic modulus
Figure 53: Comparison of measured and backcalculated E(t), aT(T) and |E*| at section
010101.
The section property and interpolated temperature profile for LTPP section 340802 is
shown in Figure 54. The section was first evaluated to assess presence of any time delay. For this
the E(t) obtained from D(t) and the unbound modulus values obtained from elastic
backcalculation (Table 16) were used to calculate deflections. The forward calculated deflection
curves were then shifted assuming first non-zero value as the initial starting point. Assuming that
the E(t) obtained from D(t) and the elastic backcalculated modulus values approximately
represent the AC and the unbound layer properties, the peaks of the measured and the calculated
curves are expected to reach at the same time. The measured and calculated deflection histories
are shown in Figure 55. The curves were found to reach peak at the same time and hence the
107

section was chosen for further backcalculation analysis. Further it also indicates that the E(t)
obtained from D(t) approximates the viscoelastic properties of the AC layer. Relaxation
modulus, time-temperature shift factor and dynamic modulus curves for section 340802 are
compared in Figure 56, respectively.

24.06oC

Asphalt Concrete (Linear Viscoelastic)
ν=0.35
Granular Base (Ebase)
11.6" ν=0.4
6.7"

24.66oC
25.45oC

Subgrade (Esubgrade)
ν=0.45
∞ 
Figure 54: LTPP section 340802 cross section and temperature profile in AC layer
0.02

0.02
Measured S1
Measured S3

Measured S1
Measured S3
0.015

Measured peak
for sensor1

Calculated peak
for sensor1

Deflection in

Deflection in

0.015

Backcalculated S1
Backcalculated S3

0.01

0.005

0

Backcalculated S1
Backcalculated S3
Measured peak

Calculated peak

0.01

0.005

0

0.003

0.006

0.009

0.012

0.015

Time (sec)

0

0

0.003

0.006

0.009

Time (sec)

0.012

0.015

 

Figure 55: Comparison of measured and forward calculated deflection histories for section
340802 (a) Viscoelastic backcalculation (b) Forward calculated using E(t) obtained from
LTPP measured D(t) and elastic backcalculated modulus values.
Although, time-temperature shift factor functions showed a good match over the entire
temperature range, the backcalculated relaxation modulus and dynamic modulus curves showed

108

significant deviation at higher reduced time and lower frequencies. Better agreement was found
at time less than 10-2 seconds and frequencies greater than 100 Hz.

7

10

6

10

5

10

4

10

3

10

2

Backcalculated-1
Backcalculated-2
From D(t)

102
101

T

10

Shift factor, a (T)

Relaxation modulus, E(t) (psi)

103

10

-5

10

-4

10

-3

10

-2

10

-1

10

0

10

1

10

2

10

3

0

10

10

-1

10-2
10-3
10

-4

10

-5

Backcalculated-1
Backcalculated-2
From D(t)
ANNACAP
5

10

15

20

Dynamic modulus, |E*| (psi)

(a) Relaxation modulus
10

7

10

6

10

5

10

4

10

3

10

2

10

25

30

35

40

45

Temperature (oC)

Reduced time at 19oC (sec)

(b) Shift factor

Backcalculated-1
Backcalculated-2
From D(t)
ANNACAP
-3

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

o

Reduced frequency at 19 C (Hz)

(c) Dynamic modulus
Figure 56: Comparison of measured and backcalculated E(t), aT(T) and |E*| at section
340802.
The section property and interpolated temperature profile for LTPP section 460804 is
shown in Figure 57. The section was first evaluated to assess presence of any time delay. For this
the E(t) obtained from D(t) and the unbound modulus values obtained from elastic
backcalculation (Table 16) were used to calculate deflections. The forward calculated deflection
curves were then shifted assuming first non-zero value as the initial starting point. Assuming that
the E(t) obtained from D(t) and the elastic backcalculated modulus values approximately
109

represent the unbound layer properties, the peaks of the shifted curves are expected to reach at
the same time. The displaced measured and calculated deflection histories are shown in Figure
58. The curves were found to reach peak at the same time and hence the section was chosen for
further backcalculation analysis. Further it also indicates that the E(t) obtained from D(t)
approximates the viscoelastic properties of the AC layer. Relaxation modulus, time-temperature
shift factor and dynamic modulus curves for section 460804 are compared in Figure 59a, b and c
respectively. The predicted aT(T) curve showed good agreement with ANNACAP at
temperatures less than 25oC. However at temperatures greater than 25oC significant deviation
was observed in shift factor which was also reflected in the relaxation modulus and dynamic
modulus curves at time greater than 100 seconds and less than 101 Hz.

6.9"

12"

Asphalt Concrete (Linear Viscoelastic)
ν=0.35
Granular Base (Ebase)
ν=0.4

31.18oC
28.17oC
25.48oC

Subgrade (Esubgrade)
ν=0.45
∞ 
Figure 57: LTPP section 460804 cross section and temperature profile in AC layer

110

0.04

0.04
Measured S1
Measured S3

0.035

Backcalculated S1
Backcalculated S3

0.03

Deflection (in)

Deflection (in)

0.03
0.025
0.02
0.015

0.015

0.005

0.005
0.003

0.006

0.009

0.012

0

0.015

Time (sec)

Measured peak
for sensor1

0.02

0.01

0

Backcalculated S1
Backcalculated S3

Calculated peak
for sensor1

0.025

0.01

0

Measured S1
Measured S3

0.035

0

0.003

0.006

0.009

0.012

0.015

Time (sec)

 

 

Figure 58: Comparison of measured and forward calculated deflection histories for section
460804 (a) Viscoelastic backcalculation (b) Forward calculated using E(t) obtained from
LTPP measured D(t) and elastic backcalculated modulus values. 

7

10

6

10

5

10

4

10

3

10

2

102
101

T

10

Shift factor, a (T)

Relaxation modulus, E(t) (psi)

103

10

Backcalculated-1
Backcalculated-2
From D(t)
-5

10

-4

10

-3

10

-2

10

-1

10

0

10

1

10

2

10

3

0

10

10

-1

10-2
10-3
10

-4

10

-5

Backcalculated-1
Backcalculated-2
From D(t)
ANNACAP
5

10

15

20

Dynamic modulus, |E*| (psi)

(a) Relaxation modulus
10

7

10

6

10

5

10

4

10

3

10

2

10

25

30

35

40

Temperature (oC)

Reduced time at 19oC (sec)

(b) Shift factor

Backcalculated-1
Backcalculated-2
From D(t)
ANNACAP
-3

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

o

Reduced frequency at 19 C (Hz)

(c) Dynamic modulus
Figure 59: Comparison of measured and backcalculated E(t), aT(T) and |E*| at section
460804.
111

45

4.3.3 Backcalculation of Linear Viscoelastic Pavement Properties using Two-Stage Method
In the previous backcalculation process, viscoelastic and unbound properties were
calculated at the same step; however in this section a two stage linear viscoelastic
backcalculation scheme is presented. The first stage was to perform linear elastic backcalculation
of unbound material properties, which was followed by linear viscoelastic backcalculation (using
BACKLAVA/BACKLAVAP) of AC layer viscoelastic properties (E(t) sigmoid coefficients: c1,
c2, c3 and c4 and shift factor aT(T) coefficients a1 and a2). Details of stage 1 and stage 2 steps are
presented in the following sections.
Stage-1, Elastic backcalculation for Unbound Layer Properties: It is important to verify
that the elastic backcalculation (Stage-1) gives unbound granular modulus values close to the
actual values. If this is verified, the backcalculated Eunbound values can be fixed in viscoelastic
backcalculation (Stage-2) and only the 6 unknowns of AC layer can be backcalculated. Known
and unknown variables in the first and second stages of backcalculation are listed in Table 18. In
the first stage, elastic backcalculation was performed assuming the AC layer as linear elastic. In
the second stage, viscoelastic backcalculation was performed keeping the unbound granular layer
modulus values obtained in the first stage fixed.
In order to perform the verification, first, various synthetic deflection time histories were
obtained by running LAVA on the structure shown in Table 19 at various temperature profiles
(also shown in Table 19). These synthetic deflections were used in “Stage 1 – Elastic
Backcalculation”, which computed Eunbound values. Then these backcalculated Eunbound values
were compared to the original Eunbound values used in the original layered viscoelastic forward
computation.

112

Table 18: Variables in two stage linear viscoelastic backcalculation analysis
Stage 1: Elastic backcalculation
Unknown (backcalculated) parameters

Known parameters

Thickness & Poisson’s ratio of each layer
Eac, elastic modulus of AC layer
FWD parameters (contact radius, pressure, Eunbound (i), unbound layer moduli, i=1…NL ,
locations of the sensors, etc.)
NL = number of unbound layers.
Stage 2: Viscoelastic backcalculation
Known parameters
Unknown (backcalculated) parameters
Thickness & Poisson’s ratio of each layer E(t) sigmoid coefficients: c1, c2, c3 and c4
and FWD parameters
Eunbound (i), unbound layer moduli Shift factor aT(T) coefficients a1 and a2
backcalculated in stage 1
Table 19: Pavement properties in two stage linear viscoelastic backcalculation analysis
Property
Thickness (AC followed by granular layers)
6, 20, inf (in)
Poisson ratio {layer 1,2,3…}
0.35, 0.3, 0.45
Eunbound {layer 2,3…}
25560, 11450 (psi)
E(t) sigmoid coefficient {layer 1}
0.841, 3.54, 0.86, -0.515
aT(T) shift factor polynomial coefficients {layer 1}
4.42E-04, -1.32E-01
Total number of sensors
8
Sensor spacing from the center of load (inches)
0, 8, 12, 18, 24, 36, 48, 60
AC layer temperature profile {T1-T2}: {10-0}, {15-10}, {20-10}, {20-15}, {25-20}, {3020}, {30-25}, {35-30}, {15-10-5}, {20-15-10}, {25-20-15}, {30-25-20}, {35-30-25},
{40-35-30}
4

3 x 10

Average modulus (psi)

Base

Subgrade

Actual E=25560 psi

4

2.5 x 10

4

2 x 10

4

1.5 x 10

Actual E=11450 psi

4

1 x 10

3

5 x 10

0
10_0

15_10

20_10

20_15

25_20

30_20

30_25

35_30

40_30

40_35

Temperature profile oC

Figure 60: Elastic backcalculation of two-step temperature profile FWD data, assuming
AC as a single layer in two stage backcalculation

113

4

3 x 10

Average modulus (psi)

Base

Subgrade

Actual E=25560 psi

4

2.5 x 10

4

2 x 10

4

1.5 x 10

Actual E=11450 psi

4

1 x 10

3

5 x 10

0
15_10_5

20_15_10

25_20_15

30_25_20

35_30_25

40_35_30

o

Temperature profile C

Figure 61: Elastic backcalculation of three-step temperature profile FWD data, assuming
AC` as a single layer in two stage backcalculation
The analysis results shown in Figure 60 and Figure 61 were based on elastic
backcalculations that assumes a single AC layer. However, in the LAVA forward computations,
because of different temperatures with depth, multiple layers of AC (2 layers for Figure 60 and 3
layers for Figure 61 analysis) were used. In order to investigate if selection of number of AC
layers on the elastic backcalculation, the computations were repeated assuming the AC layer to
be consisting of 2 or 3 independent elastic layers. Average backcalculated base and subgrade
modulus values for two-step and three-step temperature profiles are shown in Figure 62 and
Figure 63, respectively. Comparing Figure 60 to Figure 62 and Figure 61 to Figure 63, it can be
seen that, assuming single or multiple AC layers didn’t significantly affect backcalculation of
base and subgrade elastic modulus. From these analyses (Figure 60 through Figure 63), it can be
concluded that, it is possible to first perform elastic backcalculation (Stage-1) for the unbound
layer properties and fix these in Stage-2.

114

4

3 x 10

Average modulus (psi)

Base

Subgrade

Actual E=25560 psi

4

2.5 x 10

4

2 x 10

4

1.5 x 10

Actual E=11450 psi

4

1 x 10

3

5 x 10

0
10_0

15_10

20_10

20_15

25_20

30_20

30_25

35_30

40_30

40_35

o

Temperature profile C

Figure 62: Elastic backcalculation of two-step temperature profile FWD data, assuming
two AC sublayers in two stage backcalculation
4

3 x 10

Average modulus (psi)

Base

Subgrade

Actual E=25560 psi

4

2.5 x 10

4

2 x 10

4

1.5 x 10

Actual E=11450 psi

4

1 x 10

3

5 x 10

0
15-10-5

20-15-10

25-20-15

30-25-20

35-30-25

40-35-30

o

Temperature profile C

Figure 63: Elastic backcalculation of three-step temperature profile FWD data, assuming
three AC sublayers in two stage backcalculation
Stage-2, Viscoelastic Backcalculation for E(t) of AC Layer: After fixing the unbound
layer modulus values, the AC layer properties (E(t) sigmoid coefficients: c1, c2, c3 and c4 and
shift factor aT(T) coefficients a1 and a2) are backcalculated using viscoelastic backcalculation
algorithm (BACKLAVA). It should be noted that for viscoelastic backcalculation, as done
earlier, set of FWD test data at different temperature can be used for backcalculation. This is due
to the fact that even though the temperatures are different, the characteristic properties of AC
layer (E(t) or |E*| master curves) remains the same. In this stage, viscoelastic backcalculation
was performed on a set of temperature profiles keeping the actual unbound modulus values
constant.
115

Average error (over reduced times from 10-8 to 108 seconds) in E(t) master curve,
obtained from set of two two-step and two three-step temperature profiles are shown in Figure 64
and Figure 65 respectively. It can be observed from Figure 64 that, for all the cases of presented
2-step temperature profile sets, average error in backcalculated E(t) is below 10-12% except for
FWD test at {30-20}&{40-30}. It can be observed from Figure 65 that, for all the cases of
presented for 3-step temperature profile sets, average error in backcalculated E(t) is below 5.5%
except for FWD test at {30-25-20}&{25-20-15}. These results indicate that the 2-stage
algorithm works well in backcalculating the E(t) of AC layer. From Figure 64 and Figure 65, E(t)
errors obtained in 2-stage backcalculation are less when compared to single stage
backcalculation (Figure 49). However, it should be noted that results presented in Figure 64 and
Figure 65 are from backcalculation using a set of two FWD test data each obtained at different
temperature profile, whereas results in Figure 49 are from backcalculation using single FWD
data. But, the result does indicate that backcalculation using a set of FWD test data each obtained
at different temperature profile may improve the accuracy.

Average error (%)

20

15

10

5

0

10-0 &
20-10

20-15 &
15-10

20-10 &
30-20

30-25 &
25-20

35-30 &
30-25

30-20 &
40-30

40-35 &
35-30

o

Temperature profile C

Figure 64: Error in backcalculated E(t) curve from two-step temperature profile.

116

Average error (%)

20

15

10

5

0

20-15-10 &
15-10-05

25-20-15 &
20-15-10

30-25-20 &
25-20-15

o

30-25-20 &
35-30-25

40-35-30 &
35-30-25

Temperature profile C

Figure 65: Error in backcalculated E(t) curve from three-step temperature profile.
4.4 Layered Viscoelastic-Nonlinear Backcalculation (BACKLAVAN) Algorithm
4.4.1 Introduction
In the previous section, analyses were performed assuming the unbound layers to be
linear elastic. Due to the linear elastic assumption the forward engine in the backcalculation
algorithms were computationally efficient. Because of inherent iterative nature of the
viscoelastic-nonlinear forward solution (i.e., LAVAN), it can take very long time to
backcalculate all the parameters (i.e.,

,

,

,

,

,

of the AC and

,

,

	of the unbound

layer) if all the parameters are varied simultaneously during the backcalculation. Furthermore,
like many other complex engineering applications involving multi-variable optimization,
backcalculation of pavement properties can also involve multiple local minima solutions. It is
well known that even linear elastic backcalculation (considering all the pavement layers to be
linear elastic) of pavement possess multiple local solutions, and hence, a unique solution using
traditional optimization methods is difficult to obtain.
In this study, optimization tool GA along with a simplex method was utilized together.
The GA-based optimization ensures that the solution approaches to a global minima, provided
that proper GA parameters are used. GA can iteratively sweep search the entire search domain
117

(made by the variables), moving towards the global solution by systematically eliminating the
local optimal solutions. The main advantage of GA algorithm is that it can avoid local minima
and find the global solution; however it can be a relatively slow method. The main advantage of
the simplex method is its speed; however, simplex method is influenced by the seed values and
susceptible to reaching to a local minimum. A combination of GA and simplex may be a good
alternative to reach an accurate solution.
In order to improve the computational efficiency of the backcalculation algorithm, three
basic stages are proposed. These three stages optimize the use of GA and simplex methods to
improve the speed and accuracy.
Stage-1: Nonlinear elastic backcalculation: The first stage involves nonlinear elastic
backcalculation of the unbound granular layer properties (i.e.,

,

,

). At this stage, the AC

layer is assumed to be elastic and the unbound layers are assumed to be nonlinear. Due to the
linear elastic assumption (instead of linear viscoelastic) of AC layer, a quick approximation of
unbound layer properties is obtained. It is noted that this step ignores the effects of stress
redistribution because of the viscoelasticity of AC layer, thus induces error in backcalculated k1,
k2 and k3. However, this error may not significantly affect the backcalculated |E*| master curve
(if the ultimate goal is to obtain the |E*|), as illustrated in later parts of the chapter.
Formulation of the optimization model used in the first stage of genetic algorithm is as
follows:

∑

	

∑,

118

	

	 60

GA
where 1 is the error to be minimized,

is the number of FWD sensors,

input deflection information obtained from field at sensor
(predicted) response at sensor ,

,

is the peak output

is the total number of drops. In this stage, the following

bounds were imposed for each of the unknowns:
,

,

, where

shown in Equation 9,

is the peak

,

,
,

is elastic AC modulus,

is subgrade modulus, and the superscripts and

,

are constants

represent lower and

upper limits. For backcalculation of these six variables, population size of 200 and 15 numbers
of generations were used. The default MATLAB values were used for the rest of the GA
parameters.
Stage-2, Viscoelastic - nonlinear backcalculation with fixed unbound layer properties: In
the second stage, the backcalculated unbound properties (i.e.,	 ,

,

) are fixed, and the

layered viscoelastic-nonlinear (LAVAN) forward model is used to backcalculate the linear
,

viscoelastic properties of AC layer (i.e.,

,

,

,

,

). Fixing the unbound layer properties

in the viscoelastic-nonlinear response calculation significantly reduces the computational effort
required in the backcalculation algorithm. Formulation of the optimization model used in the
second stage of genetic algorithm based inverse analysis is

∑

where,

is the total number of sensors,

field at sensor ,

∑,

	

	 61 	

is the input deflection information obtained from

is the output deflection obtained from forward analysis at sensor ,

is the

total number of deflection data points recorded by a sensor. In this stage, the following bounds
were imposed for each of the unknowns:
119

,		

	

	

	

	

	

,	

,

,	
,

The superscripts and

represent lower and upper limits and

are the sigmoid coefficients and

are the shift factor polynomial coefficients for the AC layer (see Equations (23) and (25)). The
genetic algorithm (GA) results were obtained for population size of 100 and 15 number of
generations.
Stage-3, Viscoelastic - nonlinear backcalculation, all parameters varied: In stage three, all
the viscoelastic parameters (i.e., ,
parameters ( ,

,

,

,

,

,

) of the AC layer as well as the unbound

) of the unbound layer are treated as unknowns. In this stage, a simplex

method is employed. For this purpose, MATLAB’s “fminsearch” algorithm, which is an
unconstrained nonlinear optimization tool, is used (Lagarias 1998). It is well known that
solutions obtained using such traditional methods depend on the seed value used in the
optimization (Chatti et al., 2014). However, using “fminsearch” in the third stage offers three
benefits:
1)

The seed values obtained in stage-1 and stage-2 are generally very close to the

actual values;

120

2)

As GA in stage-2 has already moved toward the global optimal solution, the third

stage improves the precision of the solution, which cannot be achieved using the GA
alone, due to the inherent randomness and discreteness involved in the technique;
3)

It reduces the number of iterations needed in the computationally expensive GA-

based stage-2.
A detailed discussion on various hybrid optimization schemes and benefits can be found
in Grosan and Abraham (2007) and El-Mihoub et al. (2006). It is worth noting that, in this study,
the two-stage optimization may be sufficient, and the stage-3 slightly improves the results.
Details of known and unknown properties used during all the stages are shown in Table 20.
Table 20: Pavement properties and test inputs in staged nonlinear viscoelastic
backcalculation.

Esubgrade (MPa)

Unknown

E(t)AC (MPa)

Unknown
(E(t) = Constant)

Stage-2:
Nonlinear Viscoelastic
Known (AC),
Known (BASE),
inf (SUBGRADE)
Known (AC),
Known (BASE),
Known (SUBGRADE)
Obtained
from Stage-1
Obtained
from Stage-1
Unknown
(sigmoid coefficient)

Genetic Algorithm

Genetic Algorithm

Simplex Algorithm

Known
peak stress

Known
load history

Known
load history

Known
peak deflection

Known
deflection history

Known
deflection history

Property
Thickness (cm)
Poisson ratio
Ebase (MPa)

Backcalculation
algorithm used
Test Inputs
Surface loading
(kPa)
Surface
deflection
(µm)

Stage-1:
Nonlinear elastic
Known (AC),
Known (BASE),
inf (SUBGRADE)
Known (AC),
Known (BASE),
Known (SUBGRADE)
Unknowns
(k1, k2, k3)

121

Stage-3:
Nonlinear Viscoelastic
Known (AC),
Known (BASE),
inf (SUBGRADE)
Known (AC),
Known (BASE),
Known (SUBGRADE)
Unknowns
(k1, k2, k3)
Unknowns
Unknown
(sigmoid coefficient)

4.4.2 Verification of the Nonlinear Backcalculation Procedure
The BACKLAVAN algorithm was verified using synthetic deflection histories for two
HMA mixes, namely, Control and CRTB (for mix properties refer to Table 22). The synthetic
FWD test data (using 35milisecond haversine load) were generated at five test temperatures
(10oC, 20oC, 30oC, 40oC and 50oC) for the pavement section shown in Table 22. Stresses at
distance r=0 were used in calculating unbound base modulus value for both forward calculation
and backcalculation. An FWD test is generally composed of four independent test drops, where
each drop corresponds to a different stress level. Typical ranges of stress levels in each drop are
shown in Table 21. In this study, four FWD drops at peak stress levels of 379, 552, 689 and 950
kPa were simulated to generate synthetic FWD deflection-time history using LAVAN algorithm.
In stage-1 of the BACKLAVAN algorithm, peak stress and deflection values during all the test
drops are used as input in nonlinear elastic backcalculation. The AC layer modulus and unbound
layer properties ( ,

,

) backcalculated from synthetic deflection at different temperatures are

shown in Figure 66 and Figure 67, respectively. As expected, for both Control and CRTB mixes
the backcalculated elastic AC modulus values dropped with increase in temperature. It is noted
that the nonlinear elastic backcalculation was performed using an in-house nonlinear elastic
backcalculation program developed based on the same iterative procedure used in LAVAN.
Table 21: Typical FWD test load levels
Load Level
Drop 1
Drop 2
Drop 3
Drop 4

Allowable range (for 11.81”
diameter plate), psi
49 - 60
74 - 96
99 - 120
132 - 161

122

Used surface load,
psi
55
80
110
137.8

Table 22: Pavement geometric and material properties in two stage nonlinear viscoelastic
backcalculation.
Property
Thickness (cm)
Poisson ratio (ν)
Density (kg/m3)
Nonlinear Ebase (MPa)
Esubgrade(MPa)
AC: E(t) sigmoid
coefficient (ci) (MPa)
Shift factor coefficients (ai)
(35 msec) Haversine peak
stress (kPa)

Value used
15 (AC), 25 (Base), (Subgrade)
0.35 (AC), 0.4 (Base), 0.4
(Subgrade)
2082 (AC), 2082 (Base),
2082 (Subgrade)
Ko=0.6, k1=25, k2=0.5, k3=-0.5
68.95
Control: 1.598, 2.937, -0.272, -0.562
CRTB: 0.895,3.411,0.634,-0.428
Control: 5.74E-04,-1.55E-01
CRTB: 4.42E-04,-1.32E-01
Peak stress = 950

As shown, the backcalculated results are close to the actual values but they are generally under
predicted by the algorithm.
In stage-2, the backcalculated unbound layer properties from stage-1 were used as known
fixed values, and the viscoelastic layer properties of the AC layer were obtained. The forward
algorithm used in stage-2 backcalculation and stage-3 backcalculation are same (LAVAN)
except that, in stage-2 the nonlinear properties are fixed and not varied whereas in stage-3 all the
AC as well as unbound layer properties are varied.

123

Asphalt concrete modulus, Eac (MPa)

1.2 x 104
Control
CRTB

1 x 104

8000

6000

4000

2000

0
10

20

30

40

50

Temperature ( oC)

Figure 66: Nonlinear elastic backcalculated AC modulus (in stage-1) for Control and
CRTB Mixes, using FWD data at different test temperatures.

124

30

Actual value of k = 25 MPa
1

Control
CRTB

Actual value of k = 0.5
2

0.5

0.4

15

0.3

k

2

20

1

k (MPa)

25

0.6

Control
CRTB

10

0.2

5

0.1

0

0

10

20

30

40

o

50

10

20

Temperature ( C)
0

40

o

50

80
Actual value of subgrade modulus
E = 68.95 MPa

-0.1

Control
CRTB

Subgrade modulus, E (MPa)

70

-0.2

k

3

-0.3

-0.4

-0.5

-0.6

30

Temperature ( C)

Actual value of
k = -0.5
3

10

20

30

40

50
40
30
20
10

Control
CRTB

-0.7

60

50

0
10

o

Temperature ( C)

20

30

o

40

50

Temperature ( C)

Figure 67: Nonlinear elastic backcalculated unbound layer properties (in stage-1) for
Control and CRTB Mixes, using FWD data at different test temperature.
It should be noted that the relaxation modulus of linear viscoelastic AC layer is a function
of time, therefore in nonlinear viscoelastic backcalculation, the error function is defined as the
normalized percentage error along the entire deflection history. Although, the entire deflection
history (~35-60msec) is used in the inverse analysis, it may not be sufficient to predict the entire
relaxation modulus curve.

125

In order to obtain both the sigmoid coefficients (Equation (25)) and the shift factor
coefficients (Equation (23)), FWD deflection histories obtained at (at least) two AC temperatures
need to be used simultaneously during the backcalculation. In order to see the effect of AC
temperature on backcalculated E(t), the backcalculation was conducted at the following
temperature pairs: a) 10oC & 20oC, b) 20oC & 30oC, c) 30oC & 40oC, and d) 40oC & 50oC. The
average of backcalculated unbound layer properties obtained in stage-1 at each independent
temperature were used in stage-2.
Figure 68 and Figure 69 show the backcalculation results for Control and CRTB
mixtures, respectively. Figure 68 shows that the backcalculated E(t) master curve for Control
mix was quite good for most of the FWD test temperature pairs. At the temperature pair of 10oC
& 20oC, the backcalculated E(t) deviated from the actual E(t) after reduced time of about 104sec,
which is somewhat expected since the entire E(t) curve is not used during forward analyses at
these temperatures. Figure 69 shows that the backcalculated E(t) master curve for CRTB mix
was also quite good for most of the FWD test temperature pairs. In this mix, it was observed that
the some deviation occurred at low reduced times (<10-4sec) when FWD test temperature pairs of
30oC & 40oC and 40oC & 50oC are used in the backcalculation. This deviation was attributed to
coarse discretization of reduced time evaluated in the convolution integral (Equation 8) at these
high temperatures. Theoretically, the results can be improved by increasing the number of time
steps in the convolution integral given in Equation 8, at the expense of increasing the
computational time. Accuracy of the backcalculation was evaluated based on average error in
backcalculated

master curve using Equation (57) over reduce time of

1

10 sec to

avg
10 	sec. Figure 68a and Figure 69a shows average error ( AC
) calculated for the Control

and CRTB mixtures, respectively.
126

5

10

25

(b)

4

10

Average percentage error (%)

Relaxation modulus, E(t) (MPa)

(a)

Actual

3

10

{10,20} Stage-2
{20,30} Stage-2
{30,40} Stage-2
{40,50} Stage-2

2

10

{10,20} Stage-3
{20,30} Stage-3
{30,40} Stage-3
{40,50} Stage-3

1

10

-7

10

-5

10

-3

-1

1

3

5

10
10
10
10
10
o
Reduced time, at 19 C (Sec)

10

3

10

2

10

1

10

0

20

15

10

7

5

0

10

{10,20}

{30,40}

{40,50}
o

(c)

Shift factor, aT

{20,30}

FWD test Temperatures ( C)

Actual
{10,20}
{20,30}
{30,40}
{40,50}
{10,20}
{20,30}
{30,40}
{40,50}

Stage-2
Stage-2
Stage-2
Stage-2
Stage-3
Stage-3
Stage-3
Stage-3

-1

10

-2

10

-3

10

-4

10

0

10

20

30

40

50

o

Temperature ( C)

Figure 68: (a) Comparison of backcalculated and actual E(t) master curves (b) average
avg
percentage error,	  AC
for the Control mixture and (c) time-temperature shift factor.

127

10

5

10

10

10

10

4

25

Average percentage error (%)

Relaxation modulus, E(t) (MPa)

(a)

3

Actual
{10,20}
{20,30}
{30,40}
{40,50}
{10,20}
{20,30}
{30,40}
{40,50}

2

1

10

-7

10

-5

Stage-2
Stage-2
Stage-2
Stage-2
Stage-3
Stage-3
Stage-3
Stage-3

10

-3

10

-1

10

1

10

3

10

5

10

(b)
20

15

10

5

7

0
{10,20}

o

{20,30}

{30,40}

{40,50}
o

Reduced time, at 19 C (Sec)

FWD test temperatures ( C)

3

10

(c)

Actual
{10,20} Stage-2
{20,30} Stage-2
{30,40} Stage-2
{40,50} Stage-2
{10,20} Stage-3
{20,30} Stage-3
{30,40} Stage-3
{40,50} Stage-3

2

10

1

Shift factor, aT

10

0

10

-1

10

-2

10

-3

10

-4

10

0

10

20

30

40

50

o

Temperature ( C)

Figure 69: (a) Comparison of backcalculated and actual E(t) master curves (b) average
avg
percentage error,	  AC
for the CRTB mixture and (c) time-temperature shift factor.

128

4.4.3 Validation of BACKLAVAN using Field FWD Data
The backcalculation algorithm BACLAVAN was validated using field FWD data
obtained at the LTPP section 0101 from state 01 (Alabama). This section was selected since
FWD tests were conducted in years 2004 and 2005, and the temperature of the AC layer in 2005
FWD test was about 23oC more than the AC layer temperature during the FWD test run in 2004.
This helps in obtaining the shift factor coefficients of the AC layer. In addition, the section was
not modified (i.e., there was no re-construction or rehabilitation) between the two tests.
In LTPP database, viscoelastic properties of the AC layer of the section 0101 were
available from two different sources. The first source was the creep compliance (D(t)) curve of
the AC layer, which was measured using the cores from the field. The second source of
viscoelastic properties in LTPP database is the ANNACAP predictive model developed under
FHWA long term pavement program research (Kim et al. 2011).
In stage-1 of BACKLAVAN algorithm, peak stress and deflection values during all the
test drops shown in Table 23 were used. As shown in Table 24, the separately backcalculated
unbound base properties ( ,

,

) from 2004 and 2005 FWD data were found to be very close.

As expected, the backcalculated elastic AC modulus values dropped with increase in
temperature. However, the effect of the temperature on the unbound layer was not significant.
The nonlinear unbound layer properties obtained in stage-1 were used in stage-2, which
uses the viscoelastic-nonlinear forward algorithm LAVAN during backcalculation of the E(t) of
the AC layer. It should be noted that viscoelastic backcalculation requires the entire time history
for backcalculation. In the FWD test run in 2004, entire deflection history was available only for
drop 1, hence the backcalculation was performed only using drop 1.
129

Table 23: FWD test data from LTPP section 01-0101 for years 2004-2005.

28-Apr-05 23-Feb-04

Test Drop
date Level
1
2
3
4
1
2
3
4

Peak
stress(kPa)
373.0
575.7
781.2
1025.3
359.9
553.0
766.7
959.1

Deflection (in µm) at radial distance r cm
r=0
r=20.3 r=30.5 r=45.7
127.0
108.0 94.0
75.9
62.0
206.0
176.0 154.9 125.0
103.1
295.9
254.0 224.0 181.1
149.1
410.0
353.1 309.1 252.0
208.0
226.1
164.1 128.0 89.9
63.0
355.1
265.9 215.9 148.1
105.9
513.1
395.0 327.9 226.1
162.1
661.9
515.1 430.0 300.0
213.1

r=91.4
39.9
65.0
96.0
133.1
37.1
61.0
95.0
122.9

r=121.9
27.9
46.0
66.0
91.9
23.9
43.9
70.1
88.9

Properties

Table 24: Nonlinear elastic backcalculation results for section 0101.
FWD test year =
Avg. AC temperature =
AC modulus (MPa)
k1 (dimensionless)
k2 (dimensionless)
k3 (dimensionless)
Esubgrade (MPa)

2004
11.9oC
6491.59
1223.79
0.156
-0.594
205.68

2005
35.2oC
1567.50
1086.80
0.165
-0.577
179.93

Average
N/A
N/A
1155.29
0.161
-0.586
192.81

Comparison between predicted and measured deflection history for both the 2004 and
2005 FWD testing are shown in Figure 70. Although the two drops had very close peak loading
stress (Table 23), the peak deflection in 2005 FWD test was 76% higher. Further the test
temperature effect was also clearly visible in the deflection creeping between the two tests.
Therefore, as expected the backcalculated deflection history for FWD test at higher temperature
(i.e., in 2005 test) showed a better match compared to 2004.

130

x 10

-3

10

(a)

Backcalculated

6
4
2
0
0

x 10

-3

Measured 2005
Backcalculated

(b)

Measured 2004

8

Deflection (micro-m)

Deflection (micro-m)

10

8
6
4
2
0

0.01

0.02

0.03

0

Time (seconds)

0.01

0.02

0.03

Time (seconds)

Figure 70: Backcalculated and measured FWD deflection time histories for LTPP section
0101: (a) test run in 2004 and (b) test run in 2005.
Figure 71 shows two relaxation modulus master curves “stage-2 run 1” and “stage-2 run
2” obtained from two independent backcalculation trials in stage-2. Both the stage-2 run 1” and
“stage-2 run 2” results were obtained using the same population size and number of generations.
The only difference between the two runs was in the initial values chosen for the generation of
the first population to start the backcalculation. From the Figure 71, it can be seen that both the
results almost overlap, illustrating a good convergence. It should be noted that E(t) from
ANACAP was not available for comparison. ANACAP cannot predict phase angle and hence
cannot be used to obtain E(t) from |E*|. The dynamic modulus and phase angle master curves
were calculated from the relaxation master curve at 19oC reference temperature via prony seriesbased interconversion procedure given by Park and Schapery (1999). As shown in Table 25,
similar to the relaxation curves, the dynamic results obtained by two independent
backcalculation attempts were identical. As shown in Figure 71b, the time-temperature shift
factor curve computed using ANNACAP was found to be higher as compared to both the
backcalculated and measured curves. It can be seen from the Figure 71b that the backcalculated
time-temperature shift factor curves matched well with the measured for temperature less than
30oC. Further, it can be seen from Table 25 that that the dynamic modulus curves estimated
using ANNACAP model, especially at higher frequencies, agrees well with the dynamic
131

modulus curves obtained through measured as well as the backcalculated results. However,
ANNACAP predicted higher values at reduced frequencies less than 10-2 Hz.

10

5

4

Measured
Stage-2 (Run 1)
Stage-2 (Run 2)
Stage-3
ANNACAP

3

Measured
Stage-2 (Run 1)
Stage-2 (Run 2)
Stage-3

10

1

10

Shift factor, aT

Relaxation modulus, E(t) (MPa)

10

3

10

2

10

1

10

-1

10

-3

10

10

-5

-4

10

-3

10

-2

10

-1

10

0

10

1

10

2

10

3

10

0

1

1 10

1

2 10

1

3 10

4 10

1

1

5 10

1

6 10

o

o

Temperature ( C)

Reduced time, at 19 C (Sec)

Figure 71: Nonlinear viscoelastic backcalculation of section 0101 (a) Relaxation modulus
(b) Shift factor at 19oC.
Table 25: Comparison of viscoelastic nonlinear viscoelastic backcalculated results obtained
at different stages.
Stage 2 (GA)
Run 2

Stage 3
(fminsearch)

ANNACAP

Measured

Stage 2 (GA)
Run 1

Stage 2 (GA)
Run 2

Stage 3
(fminsearch)

f
(Hz)
10
1
0.1
10
1
0.1
10
1
0.1
10
1
0.1
10
1
0.1

Stage 2 (GA)
Run 1

T
(C)
-10
-10
-10
10
10
10
21
21
21
37
37
37
54
54
54

Phase angle (Degrees)

Measured

|E*| (MPa)

21268
18698
15220
13120
8498
4468
6713
3223
1206
1098
357
120
87
39
22

18259
16435
13548
11513
7558
3717
6474
3155
1125
1785
594
177
373
114
42

14230
13502
12119
10492
7718
4582
6501
3531
1379
1750
550
156
250
68
23

13482
11926
10282
8680
6039
3323
5049
2620
1006
1370
468
136
250
82
29

24920
21804
18045
12149
7556
4208
5786
2993
1353
1728
758
346
625
285
150

4.5
6.4
10.9
13.8
21.3
30.1
24.8
33.6
40.9
41.3
43.3
39.3
36.7
27.0
16.3

3.7
5.9
10.2
13.8
21.8
32.3
24.4
34.2
42.9
39.7
45.5
44.7
46.1
42.4
32.1

1.6
2.8
5.8
9.4
16.3
26.2
19.8
30.2
41.4
39.0
47.5
49.7
49.7
46.4
35.9

3.9
3.8
8.5
11.2
18.7
28.3
22.0
31.4
40.5
38.1
44.2
43.4
47.4
40.9
29.2

132

The nonlinear unbound layer properties obtained in stage-1 and linear viscoelastic
properties obtained in stage-2 were used in stage-3 as seed values in MATLAB function
“fminsearch”. It should be noted that similar to stage-2, stage-3 also involves viscoelastic
nonlinear backcalculation, and so requires entire time history for backcalculation. The stage-3
results are labeled as “fminsearch” in Figure 71 and Table 25. As shown in the Table 25, the
backcalculated relaxation modulus, shift factor, dynamic modulus and phase angle does not
change much, as compare to the stage-2 results. Since the unbound layer properties were also
considered as unknown variables in stage-3, new set of values for

,

,

and 	

were

obtained. A marginal increase of 2.5%, 0.11%, 1.5% and 14.0% in

,

,

and 	

were

observed. This agrees with the results obtain in synthetic data analysis, where although the
backcalculated results obtained from stage-2 were close to the actual values, they were in general
under predicted by the algorithm.
4.5 Chapter Summary
A summary of the work presented in chapter 4 is as follows:


Genetic Algorithm based backcalculation procedures (BACKLAVA, BACKLAVAP and
BACKLAVAN) were developed, which are capable of backcalculating relaxation and
dynamic modulus master curves as well as time-temperature shift factor coefficients of
AC pavements, and nonlinear k1, k2, and k3.



For backcalculation using multiple FWD (BACKLAVA), the simulated FWD model
shows that, accuracy of backcalculation depends on the set of temperatures at which
FWD test have been conducted. The present study suggests that FWD data collected at a

133

set of temperatures, between 20-40oC would be useful in estimating E(t) or |E*| master
curve.


In a sensitivity analysis, it is shown that the influence of unbound layer properties
increases with incorporation of data from further sensors and with increase in test
temperature. But using only further sensors is not recommended.



For backcalculation using single FWD at a known temperature profile (BACKLAVAP),
the analyses indicated that, the accuracy of E(t) backcalculation is the highest if FWD is
run at a temperature profile between 30oC and 20oC, with preferably ±5oC or more.



For backcalculation using temperature profile (BACKLAVAP), it is recommended to use
data from two or more FWD tests ran at different temperature profiles.

134

CHAPTER 5

EVOLUTION OF DYNAMIC MODULUS OVER TIME
5.1 Introduction
Chemical and physical changes in HMA, due to construction and in-service use may
significantly affect its mechanical properties. It is well known that one of the significant changes
encountered by asphalt pavements is age hardening. Several studies have shown that factors such
as mean annual air temperature, moisture, binder type, volumetrics of mixture and pavement
structure can significantly influence aging of HMA (Houston et al., 2005; Dickson, 1980; Jung,
2006; Raghavendra et al., 2006l; Mirza and Witczak, 1995).
Aging may lead to higher stiffness and embrittlement in HMA. Embrittlement of HMA
adversely affects its healing potential which leads to an increased rate of cracking. Similarly,
higher stiffness may reduce the fatigue life (Nf at constant strain level) of the pavement. This is
also illustrated by the Asphalt Institute fatigue formulation, were the fatigue life is inversely
related to the stiffness of asphalt (Asphalt Institute (1999)).
18.4
where,

4.32

10

.

is horizantal tensile strain at the bottom of asphalt layer and

.

		

	 62 	

is elastic modulus of

asphalt mix.
Traditionally, in field aging is determined using core samples obtained from the field.
Field coring is resource intensive (time consuming and cumbersome) and the limited amount of
samples yields insufficient information. In the present study, a simple method based on typical
135

FWD backcalculation is proposed to quantify hardening of an in-service pavement. In the next
section popular methods to simulate field aging and their limitations are discussed, which is
followed by the FWD based proposed procedure to quantify age hardening.
5.2 Aging Models
Aging in field is mainly attributed to physical and chemical changes in asphalt binder.
These changes are mainly caused by oxidation and loss of volatile oils from the asphalt binder
(Houston et al., 2005). Aging in asphalt mixtures is typically investigated using two approaches:
(1) interpolating the aging potential of asphalt mixtures based on aging potential of asphalt
binders (2) interpolating the aging potential of asphalt mixtures based on laboratory aging under
different aging conditions. Two popularly used aging models under the above categories are the
Global Aging System (GAS) model and AASHTO R30 (AASHTO R30, 2008). Both the
methods, however, fail to account for the effects of in situ field conditions on age hardening such
as compaction, air voids, moisture, temperature fluctuations and influence of aggregates.
5.2.1 GAS Model
The GAS predictive equation used in Mechanistic Pavement Design Guide (NCHRP 137A) is based on the work by Mirza and Witczak (Mirza and Witczak, 1995). The model predicts
field aged viscosity of binder with time and depth using the following equation

	

	

	 63 	

where,	
0.004166
0.197725

	

1.41213
0.068384

	
136

.

10

.

14.5521

.

	

10.47662

	

1.88161

	

	

	

	

	

	

	

=viscosity at mix/laydown (centipoise) 	

	

	

	

	

	

	

	

	

	

	

	

	

	

	

	

	

	

	

=aged viscosity (centipoise)	

	

=mean annual air temperature (oF)

	

=temperature in Rankine (oR= oF+459.7)

	

,

.

23.82

	

	

	

	

	

	 64 	

	

	 65 	

Equation (64) represents in situ viscosity at a depth of 0.25 inches, which is used in Equation 52
to predict viscosity at other depths. The predicted age hardening of the binder (in situ viscosity in
Equation (65)) is then used to estimate the aged |E*| of HMA as follows:
|

∗

|

1.249937

0.029232

0.802208

.

0.001767
.

.
.

where,
|

∗

|

Dynamic modulus, in 105 psi
Aged bitumen viscosity in 106 poise
Loading frequency in Hz
% air void in HMA
137

0.002841
.

.

.
.

	

0.058097

	 66 	

% effective binder content (by volume)
Cumulative % retained on the ¾ inch sieve (by total aggregate weight)
Cumulative % retained on the No. 4 sieve (by total aggregate weight)
Cumulative % passing the No. 200 sieve (by total aggregate weight)
The GAS model predicts an exponential decrease in hardening with most of the changes
in the viscosity in the top 25 mm. Farrar et al. (2006) conducted verification of the GAS model
using field cores collected four years after construction and found that, the GAS method under
predicts hardening in field, particularly in top 13 mm of pavement. Furthermore due to the
limited data, the existing GAS method incorporated in MEPDG ignores the influence of air voids
in field hardening predictions.
5.2.2 Mixture Conditioning of HMA: AASHTO R30
The AASHTO R30 standard was developed under SHRP studies to simulate field
oxidative age hardening in HMA (Houston et al., 2005). To simulate long term aging, the
compacted samples are required to be prepared from short term aged samples as specified in the
standard (AASHTO R30). The test aims to simulate field aging of five to seven years. The
standard requires conditioning of the compacted test specimen for 5 days in a forced-draft oven
at 85oC. In an NCHRP study to verify the field applicability of the AASHTO R30 standard,
Houston et al. (2005) compared dynamic moduli of field and laboratory aged samples. They
concluded that the existing AASHTO standard is not sufficient to simulate field aging. It was
pointed out that the existing laboratory method ignores the effect of volumetric properties, field

138

temperature (MAAT) and other environmental conditions which limits the applicability of the
method.
Both the GAS and AASHTO methods in there curent form posses serious limitations in
estimating in field age hardening of HMA. Hence more research on in field age hardening of
HMA is needed to:
1. Generate field data in developing/calibrating aging models/tests.
2. Estimate loss in pavement life due to hardening (Roque et al., 2010 and Jung, 2006).
5.3 Effect of Aging on Viscoelastic Properties of HMA
As previously discussed, viscoelastic behavior of HMA requires determining both the
time/frequency and time-temperature functions. Each parameter in the sigmoid function and the
time-temperature shift factor polynomial function influence their shape (property). It is expected
that, with time, the HMA will get stiffer and less temperature dependent. In an NCHRP study
Roque et al. (2010) showed the effect of aging on parameters in the dynamic modulus and timetemperature shift factor. The study was performed on a single mixture aged for three separate
durations in lab. The laboratory aging time was related to field aging time using a lab to field
aging relationship proposed in SHRP (Bell et al., 1994). The effect of aging was shown as the
ratio of aged values to unaged values. They concluded that the viscoelastic behavior in HMA
diminishes with aging. To further illustrate this point, the physical meaning of each parameter
and the hypothesized effect of hardening on these parameters is shown in Figure 72 to Figure 77.
CRTB mix was used as the base mixture in the illustrations. The dynamic modulus and the shift
factor polynomial coefficient of the mix are (c1=0.807, c2=3.519, c3=1.066, c4=0.407) and
(a1=4.42E-04, a2=-1.32E-01) respectively.
139

5.3.1 Effect of Variation of c1
The |E*| plots in Figure 72 show the effect of coefficient c1 on the sigmoid. Three trial
mixes with constant c1, c3, c4 coefficients and varying the c2 coefficient were generated. Trial 1 is
the figure refers to CRTB mix which was used as the base mix in the illustration. It can be seen
from the figure that, an increase or decrease in the coefficient c1 shifts the entire curve up or
down without significantly changing the shape of the function. Further, it should be noted that
the coefficient represents the long term modulus of a mixture. Since HMA hardens with aging, it
is expected that the coefficient should increases with time.

7

Dynamic modulus, |E*| (MPa)

10

6

10

5

10

4

10

C1=0.807
C1=0.600
C1=1.000

3

10

-9

10

-7

10

-5

-3

-1

1

3

5

7

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

9

10

Figure 72: Effect of variation in c1 on |E*| master curve
5.3.2 Effect of Variation of c2
The |E*| plots in Figure 73 show the effect of coefficient c2 on the sigmoid. Three trial
mixes with constant c1, c3, c4 coefficients and varying the c2 coefficient were generated. Trial 1 is
the figure refers to CRTB mix which was used as the base mix in the illustration. It can be seen
from the figure that, an increase or decrease in the coefficient shifts the instantaneous modulus of

140

the mix up or down. It should be noted that in the present study instead of c2, c (defined as sum
of c1 and c2, refer to Equation 47) has been used in the analysis. This is because c is much easier
to physicaly interpretate (the coefficient represents long term modulus of a mixture) and use in
backcalculation. Since the HMA hardens with aging it is expected that the coefficient should
increase with time.

7

Dynamic modulus, |E*| (MPa)

10

6

10

5

10

4

10

C2=1.066
C2=0.600
C2=1.400

3

10

-9

10

-7

10

-5

-3

-1

1

3

5

7

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

9

10

Figure 73: Effect of variation in c2 on |E*| master curve
5.3.2 Effect of Variation of c3
The |E*| plots in Figure 74 show the effect of coefficient c3 on the sigmoid. Three trial
mixes with constant c1, c2, c4 coefficients and varying the c3 coefficient were generated. It can be
seen from the figure that, an increase or decrease in the coefficient shifts the entire curve to the
left or to the right without changing the shape of the function. Since HMA hardens with aging it
is expected that the coefficient should increase with time.

141

7

Dynamic modulus, |E*| (MPa)

10

6

10

5

10

4

10

C3=1.066
C3=0.600
C3=1.400

3

10

-9

10

-7

10

-5

-3

-1

1

3

5

7

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

9

10

Figure 74: Effect of variation in c3 on |E*| master curve
5.3.2 Effect of Variation of c4
The |E*| plots in Figure 75 show the effect of coefficient c4 on the sigmoid. Three trial
mixes with constant c1, c2, c3 coefficients and varying the c4 coefficient were generated. It can be
seen from the figure that, change in the coefficient c4 increases or decreases the slope of the
entire curve. Since HMA hardens with aging, it is expected that the coefficient would decrease
with time.
5.3.2 Effect of Variation of a1
The time-temperature shift factor plots in Figure 76 show the effect of coefficient a1 on
the polynomial function. Three trial mixes with constant a2 coefficients and varying the a1
coefficient were generated. It can be seen from the figure that a change in the coefficient a1
increases or decreases the slope of the entire curve. Since HMA hardens with aging, it is
expected to get less temperature dependent; the coefficient should increase with time.

142

7

Dynamic modulus, |E*| (MPa)

10

6

10

5

10

4

10

C4=0.407
C4=0.300
C4=0.500

3

10

-9

10

-7

10

-5

-3

-1

1

3

5

7

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

9

10

Figure 75: Effect of variation in c4 on |E*| master curve
3

10

2

10

1

0

T

Shift factor, a (T)

10
10

-1

10

-2

10

-3

10

-4

10

a1=4.42E-04
a1=0.42E-04
a1=9.42E-04

-5

10

0

10

20

30

40

50

o

Temperature ( C)

Figure 76: Effect of variation in a1 on aT(T) master curve
5.3.2 Effect of Variation of a2
The time-temperature shift factor plots in Figure 77 show the effect of coefficient a2 on
the polynomial function. Three trial mixes with constant a2 coefficients and varying the a2
coefficient were generated. It can be seen from the figure that a change in the coefficient a2
143

increases or decreases the slope of the entire curve. Since HMA hardens with aging, it is
expected to get less temperature dependent; the coefficient should increase with time.

5

10

3

1

10

T

Shift factor, a (T)

10

-1

10

-3

10

-5

10

a2=-1.32E-01
a2=-2.32E-01
a2=-0.32E-01

-7

10

0

10

20

30

40

50

o

Temperature ( C)

Figure 77: Effect of variation in a2 on aT(T) master curve

5.4 Proposed Solution
In this study a new methodology is suggested that allows investigation of in-service age
hardening in pavement using data obtained from commonly used FWD Test.
Steps in proposed solution:
(1) Collect FWD deflection data and HMA layer temperature profile over the analysis period.
(2) Backcalculate E(t) and shift factor for each FWD test using BACKLAVAP.
(3) Convert E(t) to |E*|.
(4) Calculate normalized |E*| and shift factor coefficients to quantify aging with time.
144

It is well known that, with time, apart from aging, HMA may undergo other various physical
changes such as densification and damage. It should be noted that the above procedure ignores
these factors and lumps everything as age hardening. To recognize when damage could become
dominant in the analysis, backcalculation results are plotted along with fatigue distress data.
Furthermore, a method proposed by MEPDG to analyze damage of an in-service pavement was
used for comparison.
For rehabilitation design in MEPDG the determination of the HMA layer dynamic
modulus follows a modified procedure to account for damage incurred in the HMA layer during
the life of the existing pavement.

The procedure therefore determines a “field damaged”

dynamic modulus master curve as follows:
For Level 1 analysis, the MEPDG calls for the following procedure:
1. Conduct FWD tests on the pavement to be rehabilitated; calculate the mean
backcalculated HMA modulus, Ei.
2. Determine mix volumetric parameters and asphalt viscosity parameters from cores.
Determine binder viscosity-temperature properties through extraction/recovery.
3. Develop an undamaged dynamic modulus master curve using the data from step 2 and the
modified Witczak equation. Calculate |E*| for the same temperature recorded in the field
at the frequency equivalent to the FWD load pulse.
4. Estimate damage, di = 1-Ei/|E*|, with Ei obtained from step 1 and |E*| obtained from step
3.

145

5. Calculate c = (1-di)c2; where c2 is a function of mix gradation parameters and represents
a parameter in the sigmoidal curve definition of |E*| expressed as sigmoid function.
6. Determine field damaged dynamic modulus master curve using c instead of c2 in the
modified Witczak equation.
The above MEPDG method has the following limitation:
1. Inappropriate conversion of time to frequency in step 3. Many time-frequency
relationships have been proposed by various researches. The two most popular
conversions are as follows.

	 67 	

	

	

	

	 68 	

Equation (67) and Equation (68) would work best for harmonic/sinusoidal pulses (Al-Qadi et al.,
2008). For computing a frequency equivalent to the loading pulse Seo et al. (2012) suggested
using Equation (67) where t was assumed to be loading pulse. Al-Qadi et al. 2008 suggested
using Equation (67) where t was assumed to be the compressive stress pulse duration developed
by the surface loading in the AC layer. However, such a conversion does not apply to FWD
pulse loading.
2. Inappropriate calculation of damage. In MEPDG damage in the HMA is defined through
reduction in dynamic modulus value at computed equivalent frequency. As discussed in
the previous paragraph, the computed equivalent frequency does not represent the actual
loading frequency and hence the calculated reduction in dynamic modulus would result in
146

misleading stiffness of the HMA. Further, it was observed that the MEPDG method may
not show any damage up until damage in existing pavement is dominant over the age
hardening effect.
Due to the above limitations the MEPDG method was used only as a reference for comparison.
5.4.1 Results and Discussion of Age Hardening in LTPP Test Section
FWD test data over 5-12 years from three different LTPP sections were used to study infield age hardening. The genetic algorithm based viscoelastic backcalculation model
BACKLAVAP was used to backcalculate both dynamic modulus as well as the time-temperature
shift factors from the FWD test data. The investigation in this chapter focused on quantifying the
effect of aging on dynamic modulus and time-temperature properties of HMA. The LTPP
sections and pavement properties selected in the analysis were the same as listed in Table 13 and
Table 14. Backcalculated |E*| and aT (T) for the sections are shown in Figure 78 to Figure 80.
Dynamic modulus and time-temperature shift factor curves for section 010101 were
backcalculated for 3 FWD tests conducted in year 1993, 1996 and 2000. Fatigue cracking for the
section was observed prior to FWD testing in 2000. Hence it is expected to experience hardening
in year 1996 followed by damage in 2000. As shown in Figure 78 this was observed at higher
frequencies, however, at lower frequencies no trend in the results was found.
Dynamic modulus and time-temperature shift factor curves for section 340802 were
backcalculated for 2 FWD tests conducted in year 1993 and 1998. The backcalculation dynamic
modulus and shift factor from both the results (Figure 79) gave similar results indicating no
change in AC properties.

147

Dynamic modulus and time-temperature shift factor curves for section 460802 were
backcalculated for 5 FWD tests conducted in year 1993, 1995, 1999, 2004 and 2011. Fatigue
cracking for the section was observed prior to FWD testing in 2007. Hence it is expected to
experience hardening in year 1995 followed by 1999 and 2004.

7

10

6

10

5

10

4

10

3

10

Shift factor, aT

Dynamic modulus, |E*| (psi)

10

|E*| from Creep
03.11.93
04.19.96
05.19.00

2

10

-9

10

-7

-5

-3

-1

1

3

5

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

7

10

10

3

10

2

10

1

10

0

10

-1

10

-2

10

-3

10

-4

10

-5

aT from Creep
03.11.93
04.19.96
05.19.00

0

9

10

20

30

40

50

o

Temperature ( C)

Figure 78: Backcalculated |E*| and aT(T) at section 01-0101 for year 1993, 1996 and 2000.
10

3

7

10

2

10

6

10

5

1

10

10

4

10

3

10

Shift factor, aT

Dynamic modulus, |E*| (psi)

10

|E*| from Creep
08.09.93
26.08.98

2

10

-9

10

-7

-5

-3

-1

1

3

5

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

7

10

9

0

10
10

-1

10

-2

10

-3

10

-4

10

-5

aT from Creep
08.09.93
26.08.98

0

10

20

30

40
o

Temperature ( C)

Figure 79: Backcalculated |E*| and aT(T) at section 34-0802 for year 1993 and 1998.

148

50

10

7

10

6

10

5

3

10

2
1

10

10

4

10

3

10

Shift factor, aT

Dynamic modulus, |E*| (psi)

10

|E*| from Creep
06.18.93
06.27.95
11.05.99
08.09.04
05.23.11

2

10

-9

10

-7

-5

-3

-1

1

3

5

10 10 10 10 10 10 10
o
Reduced frequency, at 19 C (Hz)

7

10

9

0

10
10

-1

10

-2

10

-3

10

-4

10

-5

aT from Creep
06.18.93
06.27.95
11.05.99
08.09.04
05.23.11

0

10

20

30

40

50

o

Temperature ( C)

Figure 80: Backcalculated |E*| and aT(T) at section 46-0802 for year 1993, 1995, 1999, 2004
and 2011.
As discussed earlier, MEPDG suggests calculating change in dynamic modulus in an inservice pavement using modular ratio, where the ratio is defined as

	

	

where,
|

∗

| ∗|

	

	

	

	

58 	

is backcalculated AC modulus value (obtained from elastic backcalculation) and

|	is original dynamic modulus value for the same temperature recorded in the field at

frequency equivalent to the FWD load pulse. MEPDG suggest using volumetric and binder
rheological properties in predicting original dynamic modulus using Witczak model (Equation
58). However, the LTPP database does not provide all the necessary inputs to the model and
hence the dynamic modulus obtained from creep compliance data was used. Further, for
temperature correction, the shift factor obtained from the creep data was used.
In field, pavement is subjected to different temperature, vehicle speed and loading.
Hence, the viscoelastic response in HMA is a resultant of both time and temperature. In

149

literature, moving load in field is typically assumed to correspond to 10Hz. Hence the dynamic
modulus value at 10 Hz is assumed to be critical in fatigue cracking. For comparison the
Modular Ratio was plotted along with dynamic modulus at 10 Hz. Measured field fatigue
cracking, backcalculated |E*| at 10 Hz at temperature 10oC, 20oC and 30oC and Modular Ratio
for the sections are shown in Figure 81 to Figure 83.
In general it was observed that, before the initiation of fatigue cracking, both the modular
ratio and backcalculated dynamic modulus (at 10 Hz) shows hardening trend, with exception in
section 46-0804. As oppose to the age hardening trend in section 46-0804 a drop in modulus
value was observed in 1999. It should be noted that Modular Ratio was calculated based on 44
independent FWD test performed at 10 different locations in each test section, whereas,
viscoelastic backcalculation was performed only at a single station using a single FWD.

150

70
Damaged

Fatigue Cracking m

2

60
50
40
30
20
10
0
1/1/92

1/1/94

1/1/96

1/1/98

1/1/00

1/1/02

1/1/04

12/31/01

1/1/04

1/1/02

1/1/04

Year
Dynamic Modulus (psi), at 10 Hz

3.5x106
10 C

3x106
2.5x106

Damaged

20 C
30 C

2x106
1.5x106
1x106
5x105
0
1/1/92

12/31/93

1/1/96

12/31/97

1/1/00

Year
5
Damaged

Modular Ratio

4

3

2

1

0
1/1/92

1/1/94

1/1/96

1/1/98

1/1/00

Year

Figure 81: Comparison of static and viscoelastic backcalculation (|E*| at 10 Hz) for section
01-0101
151

80

(a)

Fatigue Cracking m

2

70
60
50
40
30
20
10
0
9/1/93

9/1/94

9/1/95

9/1/96

9/1/97

9/1/98

9/1/99

9/1/00

Year
Dynamic Modulus (psi), at 10 Hz

3x10

6

(b)

2.5x106
2x106
1.5x10

6

1x106
5x105
0
9/1/93

10 C
20 C
30 C

9/1/94

9/1/95

9/1/96

9/1/97

9/1/98

9/1/99

9/1/00

Year
5

(c)

Modular Ratio

4

3

2

1

0
9/1/93

9/1/94

9/1/95

9/1/96

9/1/97

9/1/98

9/1/99

9/1/00

Year

Figure 82: Comparison of static and viscoelastic backcalculation (|E*| at 10 Hz) for section
34-0802
152

80

(a)

Fatigue Cracking m

2

70
Damaged

60
50
40
30
20
10
0
1/1/93

5/18/95

10/1/97

2/16/00

7/2/02

11/15/04

4/2/07

8/16/09

1/1/12

Year
Dynamic Modulus (psi), at 10 Hz

3x106
2.5x106

(b)

10 C
20 C
30 C

Damaged

2x106
1.5x106
1x106
5x105

0
1/1/93

5/18/95

10/1/97

2/16/00

7/2/02

11/15/04

4/2/07

8/16/09

1/1/12

Year
5

(c)

Modular Ratio

4

3

2
Damaged
1

0
1/1/93

5/18/95

10/1/97

2/16/00

7/2/02

11/15/04

4/2/07

8/16/09

1/1/12

Year

Figure 83: Comparison of static and viscoelastic backcalculation (|E*| at 10 Hz) for section
46-0802
153

Next, the normalized dynamic modulus and shift factor coefficients were plotted over
time. It was observed that none of section gave the expected age hardening trend in coefficient
c1. It should be noted that the backcalculation algorithm, due to the data limitations is not
sensitive to capture modulus value at very low frequencies. Coefficient c1 physically represents
the long term modulus value of mixture, and may not be captured with accuracy using the
existing backcalculation procedure. In two out of three sections the hardening trend in c2 and c3
was captured by the backcalculation and one section demonstrated hardening effect in coefficient
c4. Summary of the aging prediction performance in each LTPP section is presented in Table 26.
Table 26: Backcalculation performance in predicting aging of LTPP section.
Station
01-0101
34-0802
46-0804

C1
X
X
X

C’2

X


C3


X

154

C4


X

a1
X



a2
X
X


1.1

Normalized sigmoid coefficient C' =C + C

1

1

1

Normalized sigmoid coefficient C

1

2

1.04

0.9
Damaged
0.8

0.7

0.6

0.5

0.4

1.03

Damaged
1.02

1.01

1

0.99
0

1

2

3

4

5

6

7

8

0

1

2

3

Years

6

7

8

1.5

4

3

5

4

Normalized sigmoid coefficient C

Normalized sigmoid coefficient C

5

Years

6

Damaged

3

2

1

0

1.4

1.3

Damaged

1.2

1.1

1

0.9
0

1

2

3

4

5

6

7

8

0

1

2

3

Years

4

5

6

7

8

Years
1.5

Normalized shift factor coefficient a

1

2

1.5

Normalized shift factor coefficient a

4

1

0.5

Damaged

0

-0.5

-1

1.4

1.3

Damaged

1.2

1.1

1

0.9

-1.5
0

1

2

3

4

5

6

7

8

0

1

2

3

4

5

6

7

8

Years

Years

Figure 84: Normalized |E*| and aT(T) coefficients with ageing time for section 01-0101

155

1.002

1

1

1

1

Normalized sigmoid coefficient C

1

Normalized sigmoid coefficient C' =C + C

2

1.005

0.995

0.99

0.985

0.98

0.998

0.996

0.994

0.992

0.99

0.975

0.988
0

1

2

3

4

0

5

1

2

Years
1.06

1.03

1.05

Normalized sigmoid coefficient C

4

1.035

3

Normalized sigmoid coefficient C

4

5

3

4

5

3

4

5

Years

1.025
1.02
1.015
1.01
1.005
1

1.04

1.03

1.02

1.01

1

0.995

0.99
0

1

2

3

4

0

5

1

2

Years

Years

1.05

Normalized shift factor coefficient a

2

1.035

1

Normalized shift factor coefficient a

3

1.04

1.03

1.02

1.01

1

0.99

1.03

1.025
1.02
1.015
1.01
1.005
1
0.995

0

1

2

3

4

5

Years

0

1

2

Years

Figure 85: Normalized |E*| and aT(T) coefficients with aging time for section 34-0802

156

1.1

Normalized sigmoid coefficient C'1=C1+ C2

Normalized sigmoid coefficient C1

1.05

1

Damaged

0.95

0.9

0.85

1.08

1.06

Damaged

1.04

1.02

1

0.98
0

5

10

15

20

0

5

Years

20

Normalized sigmoid coefficient C4

1.05

3

Normalized sigmoid coefficient C

15

Years

1.1

1

0.9

0.8
Damaged
0.7

0.6

0.5

1

0.95
Damaged

0.9

0.85
0

5

10

15

20

0

5

Years

10

15

20

Years

10

1.1

Normalized shift factor coefficient a2

Normalized shift factor coefficient a1

10

8

6

4

Damaged

2

0

1

0.9

0.8

Damaged

0.7

0.6
0

5

10

15

20

Years

0

5

10

15

20

Years

Figure 86: Normalized |E*| and aT(T) coefficients with aging time for section 46-0802

157

5.5 Summary
FWD data over 5-12 years were backcalculated using the backcalculation procedure
(BACKLAVAP) developed in Chapter 4. The coefficients of the backcalculated dynamic
modulus and time-temperature shift factor functions were normalized with respect to the
backcalculated results from the first year and plotted with time. The following observations were
made:


Preliminary investigation indicated that the proposed method was able to show the as
expected age stiffening effect for 8 coefficients out of total 18 backcalculated coefficients
(4 sigmoidal coefficients and 2 time-temperature shift factor coefficients for each LTPP
section) (refer Table 26).



The MEPDG based modular ratio method was able to predict the age hardening trend in
the pavement.

158

CHAPTER 6

SENSITIVITY ANALYSIS FOR LAYER THICKNESS
6.1 Introduction
It is well known that pavement response and performance is sensitive to layer thickness.
Therefore thickness is considered to be one of the most important factors both in design and
quality control. Elevation measurement is one of the quickest methods in pavement thickness
control during construction. The method possesses a main disadvantage: with construction of
new layers the underlying layers undergo compaction. Changes in the underlying layer
thicknesses are reflected in the thickness of subsequent layers, producing thinner lower layers
and thicker upper layers compared to the design thicknesses.
In a study to analyze variability in LTPP pavement layer thickness, Selezneva et al.
(2002) found that elevation measured thickness values were statistically different from core
measured values. Important observations made in the study were:


For 84 percent of layers, thickness variation within a section follows a normal
distribution.



For the same layer and material type, the mean constructed layer thickness tends to be
above the design value for thinner layers and below the design values for thicker
pavements.



AC surface and binder layers show the highest number of sections with significant
deviation from the design thickness.
159

Layer thickness variability has prompted many states to conduct inventory analysis using
ground penetration radar (GPR) based non-destructive test method of pavement thickness
evaluation (Noureldin et al., 2005; Uddin, 2006; Marc, 2007; Gucunski et al., 2008; Nazef,
2011). Occasionally, the FWD vehicles are also fitted with GPR to collect pavement surface
thickness along with the FWD measurements. However, thickness measurement using GPR has
the following issues:
1. GPR estimated HMA thickness may deviate from core thickness (in situ) by
approximately 8%-10% (or ±1 in), which could be even higher compared to actual design
thickness. Furthermore GPR requires calibration using field cores to reach this accuracy.
2. The results are influenced by various other factors such as dielectric properties, moisture,
and thickness of the constituent layer material. As an example, two layers with similar
dielectric properties may not be separately recognized.
3. Typically GPR is used only for surface layer thickness evaluation; accuracy of GPR
diminishes significantly with depth of the structure. As an example, occasionally
interface between the unbound layers or surface and the unbound layer is not identifiable.
Because of the inevitable variability in pavement thickness it should be considered as a
variable in FWD analysis. But considering layer thickness as a variable would increase the
number of unknowns and reduce time efficiency of the analysis. Hence in most of the
backcalculation studies, thickness is typically assumed to be a known (constant) parameter. This
has generated interest in researchers in estimating error in backcalculation due to variation in
pavement layer thickness. Irwin et al. (1989) studied the effect of random variation in layer
thickness on elastic backcalculation. Backcalculation program, MODCMP2, was used to
160

backcalculate elastic modulus of pavement structure composed of 3 inches AC, 6 inches base,
and 12 inches subbase layer course. They reported +37% to -21% error in AC modulus and
+19% to -21% errors in base modulus. In another independent study, Rwebangira et al. (1987)
studied the effect of variation in AC and base thickness on backcalculation. Backcalculation
program BISDEF was used to backcalculate elastic modulus of a pavement structure composed
of 2.5 inches of AC and 14.5 inches of base layer. They reported 60% error in AC and base
modulus for a 1 inch change in AC layer thickness and 18% error in base modulus for 2 inch
change in base layer thickness. Both the studies concluded that variation in AC thickness is most
reflected in backcalculation of AC modulus, followed by base and subgrade. Both the studies
reported insignificant error in subgrade modulus backcalculation. In a SHRP study, Briggs et al.
(1992) evaluated the effect of AC and base layer thickness variation in GPS sites on
backcalculation. They reported error up to 100% in AC layer modulus and 80 percent in the base
layer modulus.
Given the focus of the dissertation is to backcalculate E(t), in the present study, effect of
variation in AC and base thickness on viscoelastic backcalculation is studied. Three cases were
considered:
1.

Case-1: Effect of variation in AC thickness (keeping the base thickness constant) on
backcalculation of AC relaxation modulus, base modulus, and subgrade modulus. In
case-1 the AC layer thickness was varied based on equation:
	

	

	
	

161

		

%	
	

	 	

	

	 69 	

	

	 70 	

	

where,

and

	

and

	

	

compared to the

	

are as constructed AC and base thicknesses,

are design AC and base thicknesses, %	
	

is error in

.

2. Case-2: Effect of variation in base thickness (keeping the AC thickness constant) on
backcalculation of AC relaxation modulus, base modulus, and subgrade modulus. In
case-2 the base layer thickness was varied based on equation:
	

	

	

		

%	

	

	

	

	

	 71 	

	

	 72 	

3. Case-3: Effect of variation in AC and base thickness (keeping the sum of the AC and
base layer constant) on backcalculation of AC relaxation modulus, base modulus, and
subgrade modulus.
	
	

	

	
	

		

%	
	

	

	 73 	
	

	 74 	

The three cases were analyzed for two hypothetical pavement structures. The difference
between the two structures was in the AC layer thickness: (i) 6 inch AC layer and (ii) 12 inch AC
layer. All the analysis was performed on five different FHWA field HMA mixes as shown in
Figure 87. Backcalculation was performed using the procedure BACKLAVAP developed in
Chapter 4. To keep all the cases comparable same temperature profile of 30-25-20oC was
assumed in the AC layer for all cases.

162

Steps followed in the sensitivity analysis are as followed:
Step-1: Use forward calculation algorithm LAVAP to calculate deflection histories considering
actual (as built) pavement structure (i.e. using the thicknesses obtained from Equation (69) for
case-1, Equation (71) for case-2 and Equation (73) for case-3).
Step-2: Use BACKLAVAP to backcalculate layer properties (E(t), base and subgrade elastic
modulus) assuming pavement layer thicknesses to be equal to the design thickness.
Step-3: Convert E(t) to |E*|; and calculate error in backcalculated |E*| and base and subgrade
elastic modulus.

7

Relaxation modulus, E(t) (psi)

10

ALF Control
Terpolymer
SBSLG
CRTB
AirBlown

6

10

5

10

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

5

7

10 10 10 10 10 10 10
o
Reduced time, at 19 C (sec)

9

10

Figure 87: Relaxation modulus master curves for the sensivity analysis
6.2 Sensitivity Analysis for Layer Thickness: 6 Inches AC Layer
The pavement properties used in the study are given in Table 27. E(t), base modulus, and
subgrade modulus were backcalculated assuming -10%, -20%, 0%, +10%, +20%, and 30% error
163

in pavement layer thicknesses. Results for the three case studies are presented in the following
sections. The backcalculated E(t) master curves for all the mixes are presented in Appendix A.
Table 27: Pavement properties in sensitivity analysis
Property
Structure-A
Structure-B
Design thickness
6, 12 (in)
12, 12 (in)
Poisson ratio {layer 1,2,3}
0.35, 0.4, 0.45
Eunbound {layer 2,3}
35000, 20000 (psi)
E(t) sigmoid coefficients {layer 1}
5 AC mixtures: Control,
aT(T) shift factor polynomial coefficients
Terpolymer, SBS-LG, CRTB, Air blown
{layer 1}
Sensor spacing from the center of load
0, 8, 12,18,24, 36,48, 60 (in)
6.2.1 Sensitivity Analysis for AC Layer Thickness
Figure 88 and Figure 89 show error in backcalculated relaxation modulus for -10%, 20%, 0%, +10%, +20%, and +30% variation in 6 inches AC layer and 12 inches AC layer
thicknesses. Error in the relaxation modulus was calculated using Equation (56). It should be
noted that error corresponding to 0% thickness variation in AC layer corresponds to the case
when no variation in thickness was considered, and hence represents the inherent error in the
backcaulation procedure itself. It can be seen from both the figures that except for CRTB
structure-B, error in the backcalculated relaxation modulus increased with increased error in AC
layer thickness. Comparison of backcalculated dynamic moduli for all five HMA mixtures, for
structures-A and structure-B are shown in Figure 90 and Figure 91. For comparison the
backcalculated values were normalized with actual modulus values, defined as modular ratio in
the figures. It can be seen from the figure that the dynamic modulus is overestimated with
positive variation in AC layer thickness, were as the relaxation modulus is underestimated with
negative variation in AC layer thickness. This trend was observed for the dynamic modulus at
frequencies 10+8, 10+6, 10+4 and 100 Hz. There was no trend observed at 10-4 Hz and beyond.

164

The modular ratio for the five mixes at -20% and +30% error for frequencies greater than
100 Hz in structure-A ranged from 0.40-0.65 and 1.25-2.18 respectively. The modular ratio for
the five mixes at -20% and +30% error for frequencies greater than 100 Hz in structure-B ranged
from 0.35-0.65 and 0.92-2.34 respectively.

-20%

-10%

0%

10%

20%

30%

Error in relaxation modulus (%)

120

Structure-A
100

80

60

40

20

0
ALF Control

Terpolymer

SBLG

CRTB

Air Blown

AC mixture type

Figure 88: Error in backcalculated relaxation modulus from variation in 6 inch AC
thickness.
-20%

-10%

0%

10%

20%

30%

Error in relaxation modulus (%)

100

Structure-B
80

60

40

20

0
ALF Control

Terpolymer

SBLG

CRTB

Air Blown

AC mixture type

Figure 89: Error in backcalculated relaxation modulus from variation in 12 inch AC
thickness.
165

-20%

-10%

0%

10%

20%

30%

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2

1.5

1

0.5

0

10

8

10

6

10

4

10

0

10

-20%

2

0%

10%

0

10

8

30%

1
0.5

8

10

6

10

4

10

0

10

-20%

2.5

Modular Ratio
(backcalculated |E*|/actual |E*|)

10

4

10

0

10

-4

-10%

0%

10%

20%

30%

1.5
1
0.5
0

-4

10

8

10

6

10

4

10

0

10

-4

Reduced Frequency, at 19 o (Hz)

Reduced Frequency, at 19 (Hz)
-20%

6

2

o

2

10

Reduced Frequency, at 19 (Hz)

20%

1.5

10

30%

o

2

0

20%

0.5

-4

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

-10%

10%

1

Reduced Frequency, at 19 (Hz)
-20%

0%

1.5

o

2.5

-10%

-10%

0%

10%

20%

30%

1.5

1

0.5

0

10

8

10

6

10

4

10

0

10

-4

Reduced Frequency, at 19 o (Hz)

Figure 90: Backcalculated relaxation modulus in 6 inch AC structure: Case-1 sensitivity
analysis.

166

-20%

-10%

0%

10%

20%

30%

2

1.5
1
0.5
0

10

8

10

6

10

4

10

0

10

-20%

2.5

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2.5

0%

10%

0

10

8

10

30%

0.5

10

6

10

4

10

0

10

-20%

2.5

Modular Ratio
(backcalculated |E*|/actual |E*|)

10

0

10

-4

-10%

0%

10%

20%

30%

1.5
1
0.5
0

-4

10

8

10

6

10

4

10

0

10

-4

Reduced Frequency, at 19 o (Hz)

Reduced Frequency, at 19 (Hz)
-20%

4

2

o

2.5

10

Reduced Frequency, at 19 (Hz)

20%

1

8

6

o

1.5

10

30%

0.5

-4

2

0

20%

1

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

-10%

10%

1.5

Reduced Frequency, at 19 (Hz)
-20%

0%

2

o

2.5

-10%

-10%

0%

10%

20%

30%

2
1.5
1
0.5
0

8

10

10

6

10

4

0

10

10

-4

o

Reduced Frequency, at 19 (Hz)

Figure 91: Backcalculated relaxation modulus in 12 inch AC structure: Case-1 sensitivity
analysis.
6.2.2 Sensitivity Analysis for Base Layer Thickness
Figure 92 and Figure 93 show error in backcalculated relaxation modulus for -10%, 20%, 0%, +10%, +20%, and +30% error in base layer thickness for 6 inch AC layer and 12 inch
AC layer thickness structures respectively. It should be noted that error corresponding to 0%
variation in base thickness corresponds to the case when no variation in thickness was

167

considered, and hence represents the inherent error in the backcalculation procedure itself. It can
be seen from both the figures that error in the backcalculated relaxation modulus was not
significantly influenced with error in base layer thickness (compared to case-1). Comparison of
backcalculated dynamic moduli for all five HMA mixtures, for structures-A (Figure 94) and
structure-B (Figure 95) showed no trend with variation in base thickness.

-20%

-10%

0%

10%

20%

30%

Error in relaxation modulus (%)

70

Structure-A
60
50
40
30
20
10
0
ALF Control

Terpolymer

SBLG

CRTB

Air Blown

AC miture type

Figure 92: Error in backcalculated relaxation modulus from variation in base thickness for
6 inch AC structure.
-20%

-10%

0%

10%

20%

30%

Error in relaxation modulus (%)

70

Structure-B
60
50
40
30
20
10
0
ALF Control

Terpolymer

SBLG

CRTB

Air Blown

AC miture type

Figure 93: Error in backcalculated relaxation modulus from variation in base thickness for
12 inch AC structure.
168

-20%

-10%

0%

10%

20%

30%

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

1.5
1.25
1
0.75
0.5
0.25
0

8

10

10

6

10

4

0

10

10

-20%

1.5

10%

0.25

8

10

30%

0.5
0.25

10

6

10

4

0

10

10

-20%

1.5

Modular Ratio
(backcalculated |E*|/actual |E*|)

4

0

10

10

-4

-10%

0%

10%

20%

30%

1
0.75
0.5
0.25
0

-4

8

10

10

6

10

4

0

10

10

-4

o

Reduced Frequency, at 19 (Hz)
2

10

1.25

o

-20%

6

Reduced Frequency, at 19 (Hz)

20%

0.75

8

10

o

1

10

30%

0.5

0

1.25

0

20%

0.75

-4

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

1.5

0%

10%

1

Reduced Frequency, at 19 (Hz)
-10%

0%

1.25

o

-20%

-10%

Reduced Frequency, at 19 (Hz)
-10%

0%

10%

20%

30%

1.5

1

0.5

0

8

10

10

6

10

4

0

10

10

-4

o

Reduced Frequency, at 19 (Hz)

Figure 94: Backcalculated relaxation modulus in 6 inch AC structure, case-2 sensitivity
analysis.

169

-20%

-10%

0%

10%

20%

30%

1.25
1
0.75
0.5
0.25
0
10

8

10

6

10

4

10

0

10

-20%

1.5

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

1.5

10%

0.25
0

10

8

10

10

6

10

4

10

0

30%

10

-20%

2

Modular Ratio
(backcalculated |E*|/actual |E*|)

10

0

10

-4

-10%

0%

10%

20%

30%

1

0.5

0

-4

10

8

10

6

10

4

10

0

10

-4

o

Reduced Frequency, at 19 (Hz)
-20%

4

1.5

o

2

10

Reduced Frequency, at 19 (Hz)

20%

0.5

8

6

o

1

10

30%

0.5

-4

1.5

0

20%

0.75

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2

0%

10%

1

Reduced Frequency, at 19 (Hz)
-10%

0%

1.25

o

-20%

-10%

Reduced Frequency, at 19 (Hz)
-10%

0%

10%

20%

30%

1.5

1

0.5

0

10

8

10

6

10

4

10

0

10

-4

o

Reduced Frequency, at 19 (Hz)

Figure 95: Backcalculated relaxation modulus in 12 inch AC structure, case-2 sensitivity
analysis.

170

6.2.3 Sensitivity Analysis for Total Pavement Thickness
Figure 96 and Figure 97 show error in backcalculated relaxation modulus for -10%, 20%, 0%, +10%, +20%, and +30% error in 6 inch and 12 inch AC layer thickness in case-3
sensitivity analysis. It can be seen from both figures that similar to case-1, error in the
backcalculated relaxation modulus increased with increased variation in AC layer thickness.
Comparison of backcalculated dynamic moduli for all five HMA mixtures, for structures-A and
structure-B are shown in Figure 98 and Figure 99. Similar to case-1, the dynamic modulus is
overestimated with positive variation in AC layer thickness, were as the relaxation modulus is
underestimated with negative variation in AC layer thickness. Again this trend was observed for
the dynamic modulus at frequencies 10+8, 10+6, 10+4, and 100 Hz and there was no trend observed
at 10-4 Hz and beyond.
The modular ratio for the five mixes at -20% and +30% error for frequencies greater than
100 Hz in structure-A ranged from 0.42-0.78 and 1.35-2.44 respectively. The modular ratio for
the five mixes at -20% and +30% error for frequencies greater than 100 Hz in structure-B ranged
from 0.31-0.93 and 0.92-2.29 respectively.

171

-20%

-10%

0%

10%

20%

30%

Error in relaxation modulus (%)

120

Structure-A
100

80

60

40

20

0
ALF Control

Terpolymer

SBLG

CRTB

Air Blown

AC miture type

Figure 96: Error in backcalculated relaxation modulus in 6 inch AC structure, case-3
sensitivity analysis.
-20%

-10%

0%

10%

20%

30%

Error in relaxation modulus (%)

100

Structure-B
80

60

40

20

0
ALF Control

Terpolymer

SBLG

CRTB

Air Blown

AC miture type

Figure 97: Error in backcalculated relaxation modulus in 12 inch AC structure, case-3
sensitivity analysis.

172

-20%

-10%

0%

10%

20%

30%

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2

1.5

1

0.5

0

8

10

10

6

10

4

0

10

10

-20%

2.5

10%

8

10

30%

1
0.5

8

10

6

10

4

0

10

10

-20%

2.5

Modular Ratio
(backcalculated |E*|/actual |E*|)

0

10

10

0%

10%

20%

30%

1
0.5
0

-4

-10%

8

10

10

6

10

4

0

10

10

o

Reduced Frequency, at 19 (Hz)
-10%

0%

10%

20%

30%

2
1.5
1
0.5
0

8

10

-4

1.5

Reduced Frequency, at 19 (Hz)
2.5

4

2

o

-20%

10

Reduced Frequency, at 19 (Hz)

20%

1.5

10

6

o

2

0

30%

0.5

10

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2.5

0%

20%

1

Reduced Frequency, at 19 (Hz)
-10%

10%

1.5

o

-20%

0%

2

0

-4

-10%

10

6

10

4

0

10

10

-4

o

Reduced Frequency, at 19 (Hz)

Figure 98: Backcalculated relaxation modulus in 6 inch AC structure, case-3 sensitivity
analysis.

173

-4

-20%

-10%

0%

10%

20%

30%

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2.5
2
1.5
1
0.5
0

8

10

10

6

10

4

0

10

10

-20%

2.5

10%

8

10

30%

1
0.5

8

10

6

10

4

0

10

10

-20%

2

Modular Ratio
(backcalculated |E*|/actual |E*|)

0

10

10

-4

0%

10%

20%

30%

0.5

0

-4

-10%

1

8

10

10

6

10

4

0

10

10

-4

o

Reduced Frequency, at 19 (Hz)
2.5

4

1.5

o

-20%

10

Reduced Frequency, at 19 (Hz)

20%

1.5

10

6

o

2

0

30%

0.5

10

Modular Ratio
(backcalculated |E*|/actual |E*|)

Modular Ratio
(backcalculated |E*|/actual |E*|)

2.5

0%

20%

1

Reduced Frequency, at 19 (Hz)
-10%

10%

1.5

o

-20%

0%

2

0

-4

-10%

Reduced Frequency, at 19 (Hz)
-10%

0%

10%

20%

30%

2
1.5
1
0.5
0

8

10

10

6

10

4

0

10

10

-4

o

Reduced Frequency, at 19 (Hz)

Figure 99: Backcalculated relaxation modulus in 12 inch AC structure, case-3 sensitivity
analysis.
6.2.4 Summary and Discussion
Average error in backcalculated AC relaxation modulus for 6 inch and 12 inch AC
structures are shown in Table 28 and Table 29 respectively. The observations in the study are
similar to the observations made by Irwin el al. (1989) for elastic backcalculation. It was
observed that AC layer relaxation modulus was most sensitive to variation in AC layer thickness.
174

As shown in Table 28 and Table 29 variation in AC thickness in case-3 was found to be more
critical compared to variation in AC thickness in Case-1. This is due to the fact that, in case-3,
change in AC layer thickness is compensated by thickness in base layer, keeping the overall
structural thickness constant, which further underestimates AC modulus with positive variation
in AC thickness and further overestimates AC modulus with negative variation in AC thickness.
Table 28: Error in backcalculated relaxation modulus in structure-A
Error in Case‐1 
Error in Case‐2 
Error in Case‐3 
%change in 
thickness  Average  Stdev 
Average  Stdev  Average  Stdev 
‐20% 
40.1 
13.8
15.5
8.4
42.4 
5.9 
‐10% 
26.4 
11.5
11.2
9.3
34.2 
6.9 
0% 
10.6 
3.6
10.6
3.6
10.6 
3.6 
10% 
19.0 
8.8
12.8
10.1
17.2 
5.0 
20% 
35.7 
6.4
17.8
15.3
44.0  12.7 
30% 
51.2 
15.5
19.6
10.3
77.8  21.7 
Table 29: Error in backcalculated relaxation modulus in structure-B
Error in Case‐1 
Error in Case‐2 
Error in Case‐3 
%change in 
thickness  Average 
Stdev 
Average  Stdev  Average  Stdev 
‐20% 
41.5
15.5
18.3
11.0
44.8  10.8
‐10% 
21.7
10.7
13.7
11.3
36.2 
8.6
0% 
11.9
5.8
11.9
5.8
11.9 
5.8
10% 
24.2
8.3
14.4
6.0
29.9 
9.6
20% 
38.7
21.4
12.7
7.0
45.3 
5.7
30% 
62.8
39.8
13.8
5.1
68.8  18.1
Average error in backcalculated base modulus for 6 in and 12 in AC structure is shown in
Table 30 and Table 31 respectively. It was observed that base modulus was over predicted for
negative variation in thickness in all the three cases and under predicted for positive variation.
Error in subgrade modulus was insignificant with variation in thickness in all the three cases. The
maximum error in subgrade modulus was found to be 3%.

175

Table 30: Error in backcalculated base modulus in structure-A
%change in 
thickness 
‐20% 
‐10% 
0% 
10% 
20% 
30% 

Error in Case‐1 
Error in Case‐2 
Error in Case‐3 
Average 
Stdev 
Average 
Stdev  Average 
Stdev 
9.3
5.2
15.2
3.9
2.8 
7.6
5.4
4.6
8.2
3.9
2.1 
5.8
3.1
5.4
3.1
5.4
3.1 
5.4
‐11.3
10.1
‐6.0
9.3
‐7.3 
5.8
‐20.6
11.3
‐11.9
1.8
‐14.0 
13.7
‐32.8
9.4
‐19.4
5.8
‐16.7 
8.2

Table 31: Error in backcalculated base modulus in structure-B
%change in 
thickness 
‐20% 
‐10% 
0% 
10% 
20% 
30% 

Error in Case‐1 
Error in Case‐2 
Error in Case‐3 
Average 
Stdev 
Average 
Stdev  Average 
Stdev 
20.9
3.2
13.6
9.1
21.0 
1.8
21.7
3.0
7.4
6.3
16.3 
3.4
‐4.4
15.0
‐4.4
15.0
‐4.4 
15.0
‐25.3
11.1
‐2.6
7.7
‐19.4 
7.6
‐37.8
2.9
‐11.9
6.7
‐28.9 
11.6
‐34.2
11.9
‐10.5
3.8
‐32.7 
12.3

176

CHAPTER 7

SUMMARY AND CONCLUSION
In an effort towards obtaing in-situ |E*| of AC layer and linear and non-linear properties
of base/subgrade, several time domain algorithms were developed. The models developed can
consider temperature in AC layer as both constant (LAVA) as well as varying with depth
(LAVAP). The developed model LAVAN also assumes the AC layer as a linear viscoelastic
material; however, it can consider the nonlinear (stress-dependent) elastic moduli of the unbound
layers. The models have been validated using Finite Element Method (FEM) based solutions.
LAVA was used to develop a genetic algorithm-based backcalculation algorithm
(BACKLAVA) that is capable of backcalculating relaxation (and complex) modulus master
curves as well as time-temperature shift factor coefficients of AC pavements. The algorithm uses
multiple FWD drop time history data and AC layer temperature during each test. Simulated
(theoretical) data showed an excellent match between the backcalculated and actual values of
relaxation moduli. The algorithm was subsequently modified to incorporate FWDs at known
temperature profile in the AC layer and have been referred to as BACKLAVAP, which uses a
single FWD drop time history data and variation of temperature with depth of AC layer. In order
to further validate the algorithm FWD test data from LTPP sections were utilized. To evaluate
the backcalculation performance, laboratory creep compliance data available in the LTPP
database and the |E*| data from ANN-based software ANNACAP were compared.
The LAVAN was used to develop a genetic algorithm-based three stage backcalculation
methodology (called BACKLAVAN). The study showed that FWDs run at two different
177

temperatures can be sufficient to compute |E*| master curve of asphalt pavements as well as
nonlinear properties of unbound layers. The BACKLAVAN algorithm was further validated
using two FWD test runs at an LTPP test section. The backcalculated viscoelastic properties
were compared with measured values converted from laboratory creep tests and the ANN-based
software ANNACAP.
In order to be able to backcalculate the entire |E*| master curve of the asphalt layer,
accurate measurement of deflection time history from FWD is crucial. Measurement errors in the
time history of the deflections may lead to significant errors. The current versions of the
algorithms developed in this study are not able to consider the dynamics, i.e., the effects of wave
propagation. In pavements where the bedrock is close to the surface and/or there is a shallow
groundwater table, the FWD test deflection time history may exhibit significant wave
propagation effects. In such cases the current version of the algorithms should not be used.
Following conclusions can be drawn regarding the FWD data collection:
1. Accurate measurement of FWD deflection data is crucial. As a minimum, highly accurate
deflection time history at least until the end of the load pulse duration is needed for
E(t)/|E*| master curve backcalculation. The longer the duration of the deflection time
history, the better.
2. The temperature of the AC layer needs to be collected during the FWD testing.
Preferably temperatures at every 2” depth of AC layer should be collected.

178

3. Either a single FWD run on an AC with large temperature gradient or multiple FWD runs
at different temperatures can be sufficient to compute E(t)/|E*| master curve of asphalt
pavements.
4. For backcalculation using multiple FWD test data, tests should be conducted at least two
different temperatures, preferably 10oC or more apart. FWD data collected at a set of
temperatures between 20oC and 40oC should maximize the accuracy of backcalculated
E(t)/|E*| master curve with less than 10% error.
5. For backcalculation using single FWD test data at a known AC temperature profile, FWD
testing should be conducted under a temperature variation of preferably ±5oC or more.
6. Study on the effect of FWD sensor data on backcalculation indicates that, the influence of
the unbound layer properties increases with incorporation of data from sensors further
from the load and with an increase in test temperature. Further it could be concluded that
all sensors in the standard FWD configuration are needed for accurate backcalculation of
viscoelastic AC layer and unbound layers.
Following conclusions can be made for the backcalculation procedure:
1. Viscoelastic properties of the AC layer can be obtained using a two-stage scheme. The
first stage being an elastic backcalculation to determine unbound layer properties, which
is followed by viscoelastic backcalculation of E(t) of the AC layer, while keeping the
unbound layer properties fixed.

179

2. In the case of presence of considerable dynamic effects in FWD data (e.g. shallow
bedrock), the algorithms (BACKLAVA/ BACKLAVAP and BACKLAVAN) should not
be used.
3. For the genetic algorithm-based backcalculation procedures, following population and
generation sizes are recommended:
(i)

BACKLAVA model, using a set of FWD tests run at different (but constant)
asphalt layer temperatures: population size = 70 and generations = 15.

(ii)

BACKLAVAP model, using a single FWD test with known AC temperature
profile: population size = 300 and generation 15.

(iii)

BACKLAVAN (nonlinear) model, using FWD tests run at different (but constant)
asphalt layer temperatures: population size = 100 and generations = 15.

Following conclusions can be drawn from field backcalculation study:
1. For in-service pavement, AC relaxation modulus ages with time, which is reflected as
stiffening of the dynamic modulus at high frequencies. However, this was not observed at
lower frequencies.
2. Very limited data was available for analyzing the proposed procedure. More research is
needed for a comprehensive conclusion.
Following conclusions can be drawn from the sensitivity analysis:
1. Backcalculated AC layer relaxation modulus is most influenced by variation in AC layer
thickness.
180

2. Backcalculated relaxation modulus is not significantly influenced with error in base layer
thickness.
3. Subgrade modulus is least affected by variation in AC or base thickness. The maximum
error in subgrade modulus was found to be 3%.

181

APPENDIX

182

APPENDIX
Sensitivity Analysis for AC Layer Thickness
Backcalculated relaxation moduli for all the five HMA mixtures, for the structure-A and

10

7

10

6

10

5

10

4

10

3

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

structure-B are shown in Figure 100 and Figure 101.

Actual
-20%
-10%
0%
10%
20%
30%
-9

-7

10

-5

-3

-1

1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

3

5

10

7

10

6

10

5

10

4

10

3

10

10

Actual
-20%
-10%
0%
10%
20%
30%
-9

-7

10

(a) Control

-1

1

3

10

5

1

3

10

7

10

Relaxatin modulus, E(t) (psi)

10

Relaxatin modulus, E(t) (psi)

-3

(b) Terpolymer

7

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-5

10
10
10
10
10
o
Reduced time, at 19 C (sec)

-9

10

-7

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

10

(c) SBS-LG

-7

10

-5

-3

-1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

5

(d) CRTB

Figure 100: Backcalculated relaxation modulus for structure-A: Case-1 sensitivity analysis.

183

Figure 100 (cont’d)
7

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

(e) Air blown

184

10

5

7

7

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

3

10

6

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

-7

10

10

-5

-3

(a) Control

3

10

5

7

10

Relaxatin modulus, E(t) (psi)

10

Relaxatin modulus, E(t) (psi)

1

(b) Terpolymer

7

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

-9

10

-7

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

-7

10

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

(c) SBS-LG

10

5

(d) CRTB
7

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

(e) Air blown
Figure 101: Backcalculated relaxation modulus for structure-B: Case-1 sensitivity analysis.
185

Sensitivity analysis for Base layer thickness
Backcalculated relaxation moduli for all the five HMA mixtures, for the structure-A and

10

7

10

6

10

5

10

4

10

3

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

structure-B are shown in Figure 102 and Figure 103.

Actual
-20%
-10%
0%
10%
20%
30%
-9

-7

10

-5

-3

-1

1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

3

5

10

7

10

6

10

5

10

4

10

3

10

Actual
-20%
-10%
0%
10%
20%
30%

10

-9

-7

10

-5

(a) Control

1

3

5

10

7

10

Relaxatin modulus, E(t) (psi)

10

Relaxatin modulus, E(t) (psi)

-1

(b) Terpolymer

7

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

10

-9

10

-7

10

(c) SBS-LG

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

(d) CRTB

Figure 102: Backcalculated relaxation modulus for structure-A: case-2 sensitivity analysis.

186

Figure 102 (cont’d)

7

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

(e) Air blown

187

10

5

7

7

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

-7

10

10

-5

-3

(a) Control
7

3

10

5

10

5

7

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

1

(b) Terpolymer

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

-7

10

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

(c) SBS-LG

(d) CRTB
7

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

(e) Air blown
Figure 103: Backcalculated relaxation modulus for structure-B: case-2 sensitivity analysis.
188

Sensitivity analysis for total layer thickness
Backcalculated relaxation moduli for all the five HMA mixtures, for the structure-A and

10

7

10

6

10

5

10

4

10

3

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

structure-B are shown in Figure 104 and Figure 105.

Actual
-20%
-10%
0%
10%
20%
30%
-9

-7

10

-5

-3

-1

1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

3

5

10

7

10

6

10

5

10

4

10

3

10

Actual
-20%
-10%
0%
10%
20%
30%

10

-9

-7

10

-5

(a) Control

1

3

5

10

7

10

Relaxatin modulus, E(t) (psi)

10

Relaxatin modulus, E(t) (psi)

-1

(b) Terpolymer

7

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

10

-9

10

-7

10

(c) SBS-LG

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

(d) CRTB

Figure 104: Backcalculated relaxation modulus for structure-A, case-3 sensitivity analysis.

189

Figure 104 (cont’d)
7

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

(e) Air blown

190

10

5

7

7

10

Relaxatin modulus, E(t) (psi)

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

-7

10

10

(a) Control

-1

1

3

10

5

1

3

10

5

7

10

Relaxatin modulus, E(t) (psi)

10

Relaxatin modulus, E(t) (psi)

-3

(b) Terpolymer

7

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-5

10
10
10
10
10
o
Reduced time, at 19 C (sec)

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

10

5

-9

-7

10

10

-5

-3

-1

10
10
10
10
10
o
Reduced time, at 19 C (sec)

(c) SBS-LG

(d) CRTB
7

Relaxatin modulus, E(t) (psi)

10

6

10

5

10

Actual
-20%
-10%
0%
10%
20%
30%

4

10

3

10

-9

10

-7

10

-5

-3

-1

1

3

10
10
10
10
10
o
Reduced time, at 19 C (sec)

10

5

(e) Air blown
Figure 105: Backcalculated relaxation modulus for structure-B, case-3 sensitivity analysis.
191

BIBLOGRAPHY

192

BIBLIOGRAPHY

AASHTO (1993). “AASHTO Guide for Design of Pavement Structures.” American
Association of State and Highway Transportation Officials. Washington, D. C.
AASHTO R30 – 02. (2008). “Standard Practice for Mixture Conditioning of Hot-Mix
Asphalt (HMA)”. Washington, D.C.
ABAQUS. (2011). “Abaqus analysis user’s Manual.” version 6.11.
Alkasawneh, W., Pan, E., Han, F., Zhu, R., and Green, R. (2007). “Effect of Temperature
Variation on Pavement Response using 3D Multilayered Elastic Analysis.” International
Journal of Pavement Engineering, Vol. 8, No. 3, pp. 203-212.
Al-Khoury, R., Kasbergen, C., Scarpus, A., and Blaauwendraad, J. (2001). “Spectral Element
Technique for Efficient Parameter Identification of Layered Media. Part II: Inverse
calculation.” Int. J. of Solids and Struct., Vol. 38, pp. 8753-8772.
Al-Khoury, R., Scarpus, A., Kasbergen, C., and Blaauwendraad, J. (2002). “Spectral Element
Technique for Efficient Parameter Identification of Layered Media. Part III: Viscoelastic
Aspects.” Int. J. of Solids and Struct., Vol. 39, pp. 2189-2201.
Alksawneh, W. (2007). “Backcalculation of Pavement Moduli Using Genetic Algorithms.”
PhD Thesis, The University of Akron.
Al-Qadi, I., Wei, X., and Mostafa, E. (2008). “Frequency Determination from Vehicular
Loading Time Pulse to Predict Appropriate Complex Modulus in MEPDG.” Journal of the
Association of Asphalt Paving Technologists, Vol.77, pp. 739-772.
Anderson, M. (1989). “A Data Base Method for Backcalculation of Composite Pavement
Layer Moduli.” Nondestructive Testing of Pavements and Backcalculation of Moduli, ASTM
STP 1026, A. J. Bush III and G. Y. Baladi, Eds., American Society for Testing and Materials,
Philadelphia, pp. 201-216.
Appea, A. K. (2003). “Validation of FWD Testing Results At The Virginia Smart Road:
Theoretically and by Instrument Responses.” PhD Thesis, Virginia Polytechnic Institute and
State University.
Asphalt Institute (AI). (1999). “Thickness design-Asphalt pavements for highways and
streets.” 9th Ed., Manual Series No.1, Lexington.
Austroads. (2004). “Pavement design.” Australian Department of Transportation, Sydney,
Australia.
Bell, C. A., Wieder, A. J., and Fellin, M. J. (1994). “Laboratory Aging of Asphalt-Aggregate
Mixtures: Field Validation.” SHRP-A-390, Strategic Highway Research Program, National
Research Council, Washington, D.C.
193

Bremermann, H. J. (1958). “The evolution of intelligence. The nervous system as a model of
its environment.” Technical report, no.1, contract no. 477(17), Dept. Mathematics, Univ.
Washington, Seattle.
Briggs, R., Scullion, T., Maser, K. (1992). “Asphalt thickness variation on Texas Strategic
Highway Research Program Sections and Effect on Backcalculation Moduli.” Transportation
Research Record 1377, National Research Council, Washington, D.C.
Burmister, D. M. (1943). “The theory of stress and displacements in layered systems and
applications to the design of airport runways.” Highway Research Board, Proc. Annual
Meet., Vol. 23, pp. 126–144.
Burmister, D. M. (1945a). “The general theory of stress and displacements in layered soil
systems. I.” Journal of Applied Physics, Vol. 16(2), pp. 89–94.
Burmister, D. M. (1945b). “The general theory of stress and displacements in layered soil
systems. II.” Journal of Applied Physics, Vol. 16(3), pp. 126–127.
Burmister, D. M. (1945c). “The general theory of stress and displacements in layered soil
systems. III.” Journal of Applied Physics, Vol. 16(5), pp. 296–302.
Bush, A. J. (1985). “Computer Program BISDEF.” US Army Corps of Engineer Waterways
Experiment Station, Vicksburg, MS.
Chatti, K. (2004) “Use of Dynamic Analysis for Interpreting Pavement Response in Falling
Weight Deflectometer Testing.” Material Evaluation, Vol. 62, No. 7, pp. 764-774.
Chatti, K., Ji, Y., and Harichandran, R. S. (2006). “Dynamic Backcalculation of Pavement
Layer Parameters Using Frequency and Time Domain Methods.” Proc., 10th International
Conference on Asphalt Pavements, Vol. 3, Canada, pp. 5-14.
Chatti, K., Kutay, E., Lajnef, N., Zaabar, I., Varma, S., and Lee, S. (2014). “Enhanced
Analysis of Falling Weight Deflectometer Data for use with Mechanistic-Empirical Flexible
Pavement Design and Analysis and Recommendations for Improvements to Falling Weight
Deflectometer.” FHWA Report, Federal Highway Administration, Washington, DC.
Chen, D., Bilyeu, J., Lin, H., and Murphy, M. (2000). “Temperature correction on falling
weight deflectometer measurements.” In Transportation Research Record: Journal of the
Transportation Research Board, 1716, Transportation Research Board of the National
Academies, pp. 30-39.
Chou, Y. J., and Lytton, R. L. (1991). “Accuracy and Consistency of Backcalculated
Pavement Layer Moduli.” Transportation Research Record 1022, TRB, National
ResearchCouncil, Washington, D.C., pp. 1-7.
De Jong, D. L., Peutz, M. G. F., and Korswagen, A. R. (1979). “Computer program BISAR,
layered systems under normal and tangential surface loads.” Koninklijke-Shell Laboratorium,
Amsterdam, Netherlands.
194

Dickson, E. J. (1980). “The hardening of Middle East Petroleum Asphalts in Pavement
Surfacing.” Journal of the Association of Asphalt Paving Technologists, Vol. 49, pp. 30-63.
El-Mihoub T., Hopgood A., Nolle, L., and Battersby A. (2006). “Hybrid genetic algorithms:
A review.” Engineering Letters, 13, pp. 124-137.
Farrar, M., Harnsberger, M., Thomas, K., and Wiser, W. (2006). “Evaluation of oxidation in
asphalt pavement test sections after four years of service”, Proceedings of the international
conference on perpetual pavement, Columbus, Ohio.
Fogel L. J. (1962). “Autonomous automata.” Ind Res 4:14–19.
Fogel L. J., Owens A. J., and Walsh M. J. (1966). “Artificial intelligence through simulated
evolution.” Wiley, New York.
Foinquinos, R., Roesset, J. M., and Stokoe, K. H. (1995). “Response of Pavement Systems to
Dynamic Loads Imposed by Nondestructive Tests.” Transportation Research Record 1504,
pp. 57-67.
Fung, Y. C. (1996). “Biomechanics mechanical properties of living tissues.” 2nd edition,
Springer, New York.
Fwa, T. F., Tan C. Y., and Chan, W. T. (1997). “Backcalculation analysis of Pavement-Layer
Moduli Using Genetic Algorithms.” Transportation Research Record, 1570, pp. 134- 142.
Gibson, N., Qi, X., Shenoy, A., Al-Khateeb, G., Kutay, M. E., Andriescu, A., Stuart, K.,
Youtcheff, J., Harman, T. (2012). “Performance Testing for Superpave and Structural
Validation.” Federal Highway Administration Publication, FHWA-HRT-11-045.
Goldberg, D. E. (1989). “Genetic Algorithms in Search, Optimization, and Machine
Learning.” Addison-Wesley.
Gopalakrishnan, K., Thompson, M. R., and Manik, A. (2006). “Rapid Finite-Element Based
Airport Pavement Moduli Solutions using Neural Networks.” International Journal of
Computational Intelligence. Vol. 3, No. 1. World Enformatika Society.
Grosan C., and Abraham A. “Hybrid evolutionary algorithms: Methodologies, Architectures,
and Reviews.” Hybrid Evolutionary Algorithms Studies in Computational Intelligence. Vol.
75, 2007, pp. 1-17.
Gucunski, N., Rascoe, C., and Maher, A. (2008). “New Jersey Department of
Transportation” Report no 166-RU9309.
Harichandran, R. S., Baladi, G. Y., and Yeh, M. S. (1989). “MICH-PAVE user's manual,
Final Report.” FHWA-MI-RD-89-023, Department of Civil and Environmental Engineering,
Michigan State University, East Lansing, MI.

195

Harichandran, S., Mahmood, T., Raab, R., and Baladi, G. (1993). “Modified Newton
Algorithm for Backcalculation of Pavement Layer Properties.” Transportation Research
Record, 1384, pp. 15-22.
Hermansson, A. (2001). “Mathematical model for calculation of pavement temperatures
comparison of calculated and measured temperatures.” Transportation Research Record:
Journal of the Transportation Research Board, 1764, pp. 180-188.
Hicks, R. G., and Monismith, C. L. (1971). “Factors influencing the resilient properties of
granular materials.” Highway Research Record, Vol. 345, pp. 15-31.
Hjelmstad, K. D., and Taciroglu, E. (2000). “Analysis and implementation of resilient
modulus models for granular solids.” Journal of Engineering Mechanics, Vol. 126 (8), pp.
821-830.
Holland, J. H. (1975). “Adaptation in natural and artificial systems.” Ann Arbor: The
University of Michigan Press.
Horak, E. (1987). “The Use of Surface Deflection Basin Measurements in the Mechanistic
Analysis of Flexible Pavements.” Proceedings, Sixth International Conference on Structural
Design of Asphalt Pavements, Volume I. University of Michigan, Ann Arbor, MI.
Houston, W., Mirza, M., Zapata, C., and Raghavendra, S. (2005). “Environmental effects in
pavement mix and structural design systems.” NCHRP web-only document 113, NCHRP
Project 9-23.
Huang, Y. H. (2004). “Pavement analysis and design.” Second Edition, Prentice-Hall.
Indian Roads Congress (IRC). (2001). “Guidelines for the design of flexible pavements.”
IRC: 37-2001, 2nd Revision, New Delhi, India.
Irwin, H., Yang, W., and Stubstad, R. (1989). “Deflection reading accuracy and layer
thickness accuracy in backcalculation of pavement layer moduli” Nondestructive testing of
pavements and backcalculation of pavement layer moduli, American Society for Testing and
Materials, ASTM STP 1026, Philadelphia, pp. 229-224.
Irwin, L. H. (1994). “Instruction Guide for Backcalculation and The Use of MODCOMP.”
CLRP Publication No. 94-10, Cornell Univ., Local Road Program, NY.
Jamrah, A., Kutay, M. E., and Ozturk, H. I. (2014) “Characterization of Asphalt Materials
Common to Michigan in Support of the Implementation of the Mechanistic-Empirical
Pavement Design Guide.” Proceedings of 93rd Annual Transportation Research Board
Conference, Washington, D.C., January 12-16.
Ji, Y. (2005). “Frequency and Time domain Backcalculation of Flexible Pavement Layer
Parameters” Michigan State University, Department of Civil and Environmental
Engineering, Ph.D. thesis.
196

Jung, S. (2006). “The effects of asphalt binder oxidation on hot mix asphalt concrete mixture
rheology and fatigue performance.” Phd Thesis, Chemical Enginering, Texas A&M
University.
Khazanovich, L., and Roesler, J. (1997). “DIPLOBACK: neural-network-based
backcalculation program for composite pavements.” Transportation Research Board, 1570
pp. 143-150.
Kim, M., Tutumluer, E., and Know, J. (2009). “Nonlinear pavement foundation modeling for
three-dimensional finite-element analysis of flexible pavements.” International Journal of
Geomechanics, Vol 9(5), pp. 195-208.
Kim, Y. R. (2009). “Modeling of Asphalt Concrete.” 1st ed., McGraw Hill, ASCE press.
Kim, Y. R., Underwood, B., Sakhaei, F., Jackson,N., and Puccinelli, J. (2011). “LTPP
Computed Parameters: Dynamic Modulus.” FHWA-HRT-10-035, U.S. Department of
Transportation, McLean, VA 22101-2296.
Kutay, M.E., Chatti, K. and Lei, L. (2011) “Backcalculation of Dynamic Modulus from
FWD Deflection Data.” Transportation Research Record: Journal of the Transportation
Research Board, 2227, Vol. 3, pp. 87-96.
Lagarias, J.C., Reeds, J. A., Wright, M. H., and Wright, P. E. (1998). “Convergence
Properties of the Nelder-Mead Simplex Method in Low Dimensions.” SIAM Journal of
Optimization, Vol. 9 Number 1, pp. 112-147.
Leaderman, H. (1943). “Elastic and Creep Properties of Filamentous Materials and other
Higher Polymers.” The Textile Foundation, Washington, DC.
Lee, H. (2014). “Viscowave-a new solution for viscoelastic wave propagation of layered
structures subjected to an impact load.” International Journal of Pavement Engineering, Vol.
15, No. 6, pp. 542-557.
Lei, L. (2011). “Backcalculation of Asphalt Concrete Complex Modulus Curve by Layered
Viscoelastic Solution.” PhD Dissertation, Michigan State University.
Levenberg, E. (2008). “Validation of NCAT Structural Test Track Experiment using INDOT
APT Facility.” Final Report, North Central Superpave Center, Joint Transportation Research
Program Project No. C-36-31R SPR-2813, Purdue University.
Levenberg, E. (2013). “Inverse Analysis of Viscoelastic Pavement Properties using Data
from Embedded Instrumentation.” International Journal for Numerical and Analytical
Methods in Geomechanics, Vol. 37, pp. 1016-1033.
Losa, M. (2002). “The influence of Asphalt Pavement Layer Properties on Vibration
Transmission.” Int. J. of Pavements, Vol. 1, No.1, pp. 67-76.

197

Lukanen E., Stubstad, R., and Briggs, R. (2000). “Temperature predictions and adjustment
factors for asphalt pavement.” Publication No. FHWA-RD-98-085, Federal Highway
Administration (FHWA), McLean, VA.
Marc, L. (2007). “Use of Ground Penetrating Radar to Evaluate Minnesota Roads.”
Minnesota Department of Transportation, Report No MN/RC-2007-01, Maplewood, MN.
Masad, E., Huang, C., Airey, G., and Muliana, A. (2008). “Nonlinear viscoelastic analysis of
unaged and aged asphalt binders.” Construction and Building Materials, Vol. 22, pp. 21702179.
Meier R.W., Rix G.J. (1994). Backcalculation of Flexible Pavement Moduli Using Artificial
Neural Networks. Transportation Research Record 1448. Transportation Research Board.
Washington, D.C.
Meier R.W., Rix G.J. (1995). Backcalculation of Flexible Pavement Moduli From Dynamic
Deflection Basins Using Artificial Neural Networks. Transportation Research Record
1473.Transportation Research Board. Washington, D.C., 1995.
MEPDG (2004). “Guide for Mechanistic-Empirical Design of New and Rehabilitated
Pavement Structures.” National Cooperative Highway Research Program Transportation
Research Board National Research Council.
Michelow, J. (1963). “Analysis of stress and displacements in an N-Layered elastic system
under a load uniformly distributed on a circular area” Chevron Research Corporation,
Richmond, California.
Mirza, M. W., and Witczak. M.W. (1995). “Development of a Global Aging System for
Short and Long Term Aging of Asphalt Cements.” Journal of the Association of the Asphalt
Paving Technologists, Vol. 64, pp. 393-424.
Monismith, C. L. (2004). “Evolution of Long-Lasting Asphalt Pavement Design
Methodology: A Perspective. Distinguished Lecture, International Society for Asphalt
Pavements.” International Symphosium on Design and Construction of Long Lasting Asphalt
Pavements. Auburn University.
Nazef, A. (2011). “US 27/SR 25 Pavement Evaluation.” State of Florida Department of
Transportation, Research Report no. FL/DOT/SMO/11-542.
Nekouzadeh, A., and Genin, G. M. (2013). “Adaptive Quasi-Linear Viscoelastic Modeling.”
Computational Modeling in Tissue Engineering. Studies in Mechanobiology, Tissue
Engineering and Biomaterials, Vol. 10, pp. 47-83.
Newcomb, D. E. (1986). “Development and Evaluation of Regression Methods to Interpret
Dynamic Deflections.” Ph.D. Dissertation. Department of Civil Engineering. University of
Washington, Seattle, WA.

198

Noureldin, S., Zhu, K., Harri, D. and Li, S. (2005). “Non-destructive estimation of pavement
thickness, structural number and subgrade resilience along INDOT highways” Indiana
Department of Transportation, Publication No.: FHWA/IN/JTRP-2004/35.
Ooi, P. S., Archilla, A. A., and Sandefur, K. G. (2004). “Resilient modulus models for
compacted cohesive soils.” Transportation Research Record: Journal of The Transportation
Research Board, Vol. 1874, pp. 115-124.
Park, S. W. and Schapery, R. A. (1999). “Methods of Interconversion between Linear
Viscoelastic Material Functions. Part I – a Numerical Method Based on Prony Series.”
International Journal of Solids and Structures, Vol. 36, pp. 1653-1675.
Park, S. W., Park, H. M., and Hwang, J. J. (2009). “Application of Genetic Algorithm and
Finite Element Method for Backcalculating Layer Moduli of Flexible Pavements.” KSCE
Journal of Civil Engineering, Vol. 14(2), pp. 183-190.
Park, D.Y, Buch, N. and Chatti, K. (2001). “Effective layer temperature prediction model
and temperature correction via falling-weight deflectometer deflections.” In Transportation
Research Record: Journal of the Transportation Research Board, 1764, pp. 97-111.
Park, H. M., Kim, Y. R., and Park, S. (2002). “Temperature correction of multiload-level
falling weight deflectometer deflections.” Transportation Research Record: Journal of the
Transportation Research Board, 1806, pp. 3-8.
Raad, L., and Figueroa, J. L. (1980). “Load response of transportation support systems.”
Transp. Engg. J., 106 (1), pp. 111–128.
Raghavendra, S., Zapata, C., Mirza, M., Houston, W., and Witczak. M. (2006), “Verification
of the Rate of AsphaltMix Aging Simulated by AASHTO PP2-99 Protocol.” In TRB
Meeting, Washington DC.
Rechenberg, I. (1965). “Cybernetic solution path of an experimental problem.” Royal
Aircraft Establishment, Library translation NO. 1122, Farnborough, Hants., UK, August.
Reddy, M. A., Reddy, K. S., and Pandey, B. B. (2004). “Selection of Genetic Algorithm
Parameters for Backcalculation of Pavement Moduli.” The International Journal of Pavement
Engineering, Vol. 5 (2), pp. 81-90.
Roesset, J. M., Stokoe, K. H., and Seng, C. (1995). “Determination of Depth to Bedrock from
Falling Weight Deflectometer Test Data.” Transportation Research Record, 1504, pp. 68-78.
Rohde, G. T. and T. Scullion. (1990). “Modulus 4.0: Expansion and Validation of the
Modulus Backcalculation System.” FHWA/TX-91/1123-3. Texas State Department of
Highways & Public Transportation, Austin, TX.
Roque, R., Zou, J., Kim, R., Baek, C., Thirunavukkarasu, S., Underwood, S. and Guddati, M.
(2010). “Top-down cracking of hot-mix asphalt layers: models for initiation and
propagation.” NCHRP web-only document 162, NCHRP Project 1-42A.
199

Rwebangira, T., Hicks, R., and Truebe, M. (1987). “Sensitivity analysis of selected
backcalculation procedures.” Transportation Research Record 1117, National Research
Council, Washington, D.C.
Sangghaleh, A., Pan, E., Green, R., and Wang, R. (2014). “Backcalculation of pavement
layer elastic modulus and thickness with measurement errors.” International Journal of
Pavement Engineering, Vol. 15 (6), pp. 521-531.
Schapery, R. A. (1969). “On the characterization of nonlinear viscoelastic materials.”
Polymer Engineering and Science, Vol. 9, No 4.
Schapery, R. A. (1965). “Method of Viscoelastic Stress Analysis using Elastic Solutions”.
Journal of the Franklin Institute, Vol. 279(4), pp. 268-289.
Schapery, R. A. (1974). “Viscoelastic Behavior and Analysis of Composite Materials.”
Mechanics of Composite Materials, New York, Academic Press, Inc., pp. 85-168.
Schiffman, R. L. (1962). “General Solution of Stresses and Displacements in Layered Elastic
Systems.” Proceedings of the International Conference on the Structural Design of Asphalt
Pavement, University of Michigan, Ann Arbor, USA.
Schmalzer, P. (2006). “LTPP Manual for Falling Weight Deflectometer Measurements,
Version 4.1.” FHWA-HRT-06-132, U.S. Department of Transportation, McLean, VA 221012296.
Schwartz, C. W. (2002). “Effect of stress-dependent base layer on the superposition of
flexible pavement solutions.” International Journal of Geomechanics, 2(3), pp. 331–352.
Seo, J., Kim, Y., Cho, J., and Jeong, S. (2012). “Estimation of in situ dynamic modulus by
using MEPDG dynamic modulus and FWD data at different temperatures.” International
Journal of Pavement Engineering, Vol. 1, pp. 1-11.
Selezneva O., Jiang, Y., and Mladenovic, G. (2002). “Evaluation and analysis of LTPP
pavement layer thickness data.” Publication no. FHWA-RD-03-041, McLean, Virginia.
Shames, I. H., and Cozzarelli, F. A. (1997). “Elastic and Inelastic Stress Analysis”, revised
printing. Taylor and Francis.
Shell. (1978). “Shell Pavement Design Manual-Asphalt Pavements and Overlay for Road
Traffic.” Shell International Petroleum Company Limited, London.
Shook, J. F., Finn, F. N., Witczak, M. W., and Monismith, C. L. (1982). “Thickness Design
of Asphalt Pavements-The Asphalt Institute Method.” Proceeding of 5th International
Conference on Structural Design of Asphalt Pavements, Vol. I, pp. 17-44.
Taylor, A. J., and Timm, D. H. (2009). “Mechanistic characterization of resilient moduli for
unbound pavement layer materials.” NCAT Report 09-06.

200

Theyse, H. L., Beer, M., and Rust, F. C. (1996). “Overview of the South African mechanistic
pavement design analysis method.” Transportation Research Record 1539, Transportation
Research Board, Washington, DC, pp. 6–17.
Tsai, B., Kannekanti, V., and Harvey, J.T. (2004). Application of Genetic Algorithm in
Asphalt Pavement Design, Transportation Research Board, No. 1891, National Research
Council, Washington, D.C., pp. 112-120.
Tutumluer, E. (1995). “Predicting behavior of flexible pavements with granular bases,” Ph.D.
dissertation, School of Civil and Environmental Engineering, Georgia Institute of
Technology, Atlanta.
Ullidtz, P. (1988). “Pavement Analysis.” Elsevier Science. New York, NY.
Uddin, W. (2006). “Ground penetrating radar study ─ Phase I technology review and
evaluation” FHWA/MS-DOT-RD-06-182, Mississippi Department of Transportation.
Uzan, J. (1994). “Advance Backcalculation Techniques. Nondestructive Testing of
Pavements and Backcalculation of Moduli (Second Volume).” ASTM STP 1198, American
Society for Testing and Materials, Philadelphia.
Uzan, J. (1994). “Dynamic Linear Backcalculation of Pavement Material Parameters.” J. of
Transp. Eng., 120(1), pp. 109-126.
Uzan, J. (1985). “Characterization of granular material.” Transportation Research Record
1022, Transportation Research Board, Washington, D.C. pp. 52–59.
Varma, S. , Kutay, M. E., and Chatti, K. (2013b). “Data Requirements from Falling Weight
Deflectometer Tests for Accurate Backcalculation of Dynamic Modulus Master curve of
Asphalt Pavements,” Airfield & Highway Pavement Conference, Los Angeles, California,
pp. 445-455.
Varma, S., Kutay, M. E. and Levenberg, E. (2013a). “A Viscoelastic Genetic Algorithm for
Inverse Analysis of Asphalt Layer Properties from Falling Weight Deflections,”
Transportation Research Record: Journal of the Transportation Research Board, Vol. 4, pp.
38-46.
Wang, H., and Al-Qadi, I. L. (2013). “Importance of nonlinear anisotropic modeling of
granular base for predicting maximum viscoelastic pavement responses under moving
vehicular loading.” Journal of Engineering Mechanics, Vol. 139, pp. 29-38.
Warren, H., and Dieckmann, W. L. (1963). “Numerical computation of stresses and strains in
a multiple-layered asphalt pavement system,” California Research Corporation, Richmond,
California.
Witczak, M. W., and Uzan, J. (1988). “The universal airport pavement design system, Report
I of IV: Granular material characterization.” University of Maryland, College Park, Md.
201

Yan, A., and Von Quintus, H. L. (2002). “Study of laboratory resilient modulus test data and
response characteristics.” Final Report, Federal Highway Administration. Publication No.
FHWA-RD-02-051.
Yong, Y., Yang, X., and Chen, C. (2010). “Modified Schapery’s model for asphalt sand.”
Journal of Engineering Mechanics, Vol. 136(4), pp. 448-454.
Zaabar, I., Chatti, K., and Lajnef, N. (2014). “Backcalculation of asphalt concrete modulus
master curve from field measured falling weight deflectometer data using a new time domain
viscoelastic dynamic solution and hybrid optimization scheme.” Sustainability, Ecoefficiency and Construction in Transportation Infrastructure Asset Management-Losa and
Papagiannakis (Eds). Taylor & Francis, ISBN 978-1-138-00147-3.
Zhou, H. (2000). “Comparison of backcalculated and laboratory measured moduli on AC and
granular base layer materials.” Nondestructive Testing of Pavements and Backcalculation of
Moduli: Third Vol., ASTM STP 1375, S. D. Tayabji and E. O. Lukamen, Eds., American
Society of Testing and Materials, West Conshohcken, PA.

202