DEVELOPMENT OF A NEW SOLUTION FOR VISCOELASTIC WAVE PROPAGATION OF
PAVEMENT STRUCTURES AND ITS USE IN DYNAMIC BACKCALCULATION
By
Hyung Suk Lee

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

ABSTRACT
DEVELOPMENT OF A NEW SOLUTION FOR VISCOELASTIC WAVE PROPAGATION
OF PAVEMENT STRUCTURES AND ITS USE IN DYNAMIC BACKCALCULATION
By
Hyung Suk Lee
Due to the viscoelastic nature of asphalt materials and the dynamic nature of pavement structures,
it is important to consider both effects simultaneously in modeling of asphalt pavements. In this
study, a new computational algorithm, namely ViscoWave, has been developed and implemented
for modeling the pavement dynamics and viscoelasticity under an impact load generated by a
Falling Weight Deflectometer (FWD). The primary advantage of the proposed solution over
some of the existing solutions is that it uses continuous integral transforms (Laplace and Hankel
transforms) that are more appropriate for the FWD time histories whose signal characteristics are
transient, nonperiodic, and truncated.

Prior to the mathematical formulation of the developed algorithm, the fundamental properties of
a viscoelastic material and the theory of uniaxial viscoelasticity are reviewed. Then, the theory of
linear, uniaxial viscoelasticity is extended to multi-axial viscoelasticity. The multi-axial theory of
viscoelasticity is, in turn, applied to develop a methodology for analyzing the laboratory Indirect
Tensile (IDT) test data.

The theoretical development of ViscoWave follows similar steps to those used for the
development of the spectral element method. However, in place of the discrete transforms
adopted in the spectral element method, ViscoWave utilizes the continuous integral transforms

(namely Laplace and Hankel transforms) that are more appropriate for transient, nonperiodic
signals. The theory behind ViscoWave was verified by comparing the ViscoWave simulation
results to other existing solutions such as the Finite Element Analysis (FEA) and spectral
element method.

To backcalculate the pavement layer parameters, two of the well known unconstrained
optimization algorithms (Gauss-Newton and Levenberg-Marquardt methods) were adopted for
use with ViscoWave. The backcalculation was conducted using both theoretically-generated and
field-obtained FWD time histories. The results indicate that ViscoWave has great potential for
modeling the viscoelastic and dynamic effects of a pavement structure under an impact load.

To my beloved daughter Lena Dah-Ye Lee,
my mother Ok Hee Seo,
my father Doo Seon Lee,
and my sister Seung Eun Lee

iv

TABLE OF CONTENTS

LIST OF TABLES ...................................................................................................................... vii
LIST OF FIGURES ................................................................................................................... viii
CHAPTER 1 - INTRODUCTION ............................................................................................... 1
1.1. General ................................................................................................................................. 1
1.2. Problem Statement ............................................................................................................... 2
1.3. Research Objective .............................................................................................................. 3
1.4. Report Layout ...................................................................................................................... 3
CHAPTER 2 - LITERATURE REVIEW .................................................................................. 6
CHAPTER 3 - REVIEW OF UNIAXIAL THEORY OF LINEAR VISCOELASTICITY 10
3.1. Fundamental Properties of a Viscoelastic Material in Time Domain ................................ 10
3.2. Fundamental Properties of a Viscoelastic Material in Frequency Domain ....................... 11
3.3. Relationship between Viscoelastic Properties ................................................................... 13
3.4. Time-Temperature Superposition Principle....................................................................... 15
3.5. Analytical Models for Viscoelastic Properties .................................................................. 19
3.6. Prony Series ....................................................................................................................... 21
3.7. Elastic-Viscoelastic Correspondence Principle ................................................................. 26
CHAPTER 4 - EXTENSION TO MULTI-AXIAL VISCOELASTICITY AND ANALYSIS
OF INDIRECT TENSILE TEST DATA.............................................................................. 28
4.1. Viscoelastic Poisson’s Ratio .............................................................................................. 28
4.1.1. Viscoelastic Poisson’s Ratio in Time Domain............................................................ 28
4.1.2. Viscoelastic Poisson’s Ratio in Frequency Domain ................................................... 30
4.1.3. Relationship Between Time and Frequency Dependent Poisson’s Ratios ................. 31
4.2. Extension of the Elastic-Viscoelastic Correspondence Principle ...................................... 34
4.3. Determining the Viscoelastic Property from the Indirect Tensile (IDT) Test ................... 36
4.3.1. Stress Distribution within an Indirect Tensile (IDT) Specimen ................................. 38
4.3.2. Indirect Tensile Creep Test ......................................................................................... 42
4.3.3. Indirect Tensile Complex Modulus Test..................................................................... 43
CHAPTER 5 - THEORETICAL DEVELOPMENT OF A tIME DOMAIN FORWARD
SOLUTION ............................................................................................................................. 45
5.1. Governing Equations For Viscoelastic Wave Propagation................................................ 45
5.2. Solutions For The Wave Equations In The Laplace-Hankel Domain ............................... 50
5.3. Formulation of the Stiffness Matrices for the Layer Elements .......................................... 55
5.3.1. Two Noded Element for a Layer with a Finite Thickness .......................................... 55
5.3.2. One Noded Semi-Infinite Element ............................................................................. 59
5.4. Incorporating Elastic And Viscoelastic Layer Properties .................................................. 61
5.5. Construction Of The Global Stiffness Matrix.................................................................... 64
5.6. Boundary Conditions for a Circular Unit Impulse Loading at the Ground Surface .......... 67
v

5.7. Inversion Of Laplace And Hankel Transforms.................................................................. 68
5.8. Numerical Inversion of the Hankel Transform .................................................................. 68
5.9. Numerical Inversion of the Laplace Transform ................................................................. 71
5.10. System Response To Arbitrary Loading .......................................................................... 72
CHAPTER 6 - IMPLEMENTATION AND VALIDATION OF ALGORITHM –
VISCOWAVE ......................................................................................................................... 74
6.1. Simulation of an Elastic Structure using ViscoWave and LAMDA ................................. 74
6.2. Simulation of Viscoelastic Structures using ViscoWave .................................................. 75
CHAPTER 7 - THEORETICAL BACKCALCULATION USING VISCOWAVE............. 81
7.1. Review of Selected Optimization Algorithms ................................................................... 81
7.1.1. Objective Function for Optimization and its Derivatives ........................................... 81
7.1.2. Basics of Iteration ....................................................................................................... 84
7.1.3. Gauss-Newton Method ............................................................................................... 85
7.1.4. Levenberg-Marquardt Method .................................................................................... 86
7.1.5. Remarks on the Optimization Routine ........................................................................ 86
7.2. Theoretical Backcalculation Using ViscoWave ................................................................ 87
7.2.1. Single Temperature Backcalculation .......................................................................... 87
7.2.1.1. Reference Pavement Structures and Seed Values for Backcalculation ............... 87
7.2.1.2. Backcalculation Results for Pavement Structure With a Halfspace .................... 89
7.2.1.3. Backcalculation Results for Pavement Structure on a Bedrock ........................... 99
7.2.2. Multi-Temperature Backcalculation ......................................................................... 109
7.2.2.1. Theoretical Development of Generalized Power Function Relationship at
Multiple Temperatures .................................................................................................... 109
7.2.2.2. Reference Pavement Structures and Seed Values for Backcalculation ............. 110
7.2.2.3. Backcalculation Results for Pavement Structure With a Halfspace .................. 112
7.2.2.4. Backcalculation Results for Pavement Structure on a Bedrock ......................... 117
CHAPTER 8 - FIELD BACKCALCULATION USING VISCOWAVE ............................ 123
8.1. Test Site and Protocol ...................................................................................................... 123
8.2. Laboratory Dynamic Modulus Test ................................................................................. 123
8.3. FWD Backcalculation Using ViscoWave ........................................................................ 126
8.3.1. Backcalculation Results from the 3 Layer System ................................................... 127
8.3.2. Backcalculation Results from the 5 Layer System ................................................... 133
CHAPTER 9 - CONCLUSIONS AND RECOMMENDATIONS ........................................ 140
BIBLIOGRAPHY ..................................................................................................................... 142

vi

LIST OF TABLES

Table 2-1 Existing Dynamic Forward Solution and Backcalculation Pairs ................................... 8
Table 6-1 Layer Properties for Elastic Simulation of LAMDA and ViscoWave ......................... 74
Table 7-1 Reference Pavement Structures .................................................................................... 88
Table 7-2 Seed Values Used for Backcalculation......................................................................... 88
Table 7-3 Backcalculated Values and Percent Error – Pavement on a Halfspace ........................ 93
Table 7-4 Backcalculated Values and Percent Error – Pavement on a Bedrock ........................ 103
Table 7-5 Asphalt Properties Used for Reference Pavement Structures .................................... 111
Table 7-6 Base and Subgrade Properties Used for Reference Pavement Structures ................. 112
Table 7-7 Seed Values Used for Backcalculation....................................................................... 112
Table 7-8 Backcalculated Values and Percent Error – Pavement on a Halfspace ...................... 115
Table 7-9 Backcalculated Values and Percent Error – Pavement on a Halfspace ...................... 120
Table 8-1 Assumed Pavement Properties for Backcalculation ................................................... 127
Table 8-2 Seed Values Used for Backcalculation....................................................................... 127
Table 8-3 Backcalculated Values from 3 Layer System ............................................................. 128
Table 8-4 Backcalculated Values from 5 Layer System ............................................................. 134

vii

LIST OF FIGURES

Figure 3-1. Concept of Time-Temperature Superposition Principle and Mastercurve ................ 18
Figure 4-1. (a) Indirect tensile test setup and (b) specimen dimensions ....................................... 37
Figure 4-2 (a) Schematic illustration of the IDT test and (b) plane stress distributions ............... 39
Figure 5-1 Coordinate System for Axisymmetric Layers on a Halfspace .................................... 48
Figure 5-2 Bessel Functions of the First Kind .............................................................................. 53
Figure 5-3 Construction of the Global Stiffness Matrix for Structures (a) with and (b) without a
halfspace ..................................................................................................................... 65
Figure 6-1 Surface Deflections of a Layered Elastic Structure from (a) ViscoWave and (b)
LAMDA ...................................................................................................................... 75
Figure 6-2 (a) Low Temperature and (b) High Temperature Asphalt Creep Compliance Curves
Used for ViscoWave Simulation ................................................................................ 76
Figure 6-3 (a) Axisymmetric Finite Element Geometry and (b) Finite Element Mesh Used for
Simulation of Pavement Response Under FWD Loading .......................................... 77
Figure 6-4 Surface Deflections of a Layered Viscoelastic Structure with a Halfspace at Low
Temperature Simulated Using (a) ViscoWave and (b) ADINA ................................. 78
Figure 6-5 Surface Deflections of a Layered Viscoelastic Structure with a Bedrock at 3 m below
Surface at Low Temperature Simulated Using (a) ViscoWave and (b) ADINA ....... 79
Figure 6-6 Surface Deflections of a Layered Viscoelastic Structure with a Halfspace at High
Temperature Simulated Using (a) ViscoWave and (b) ADINA ................................. 80
Figure 6-7 Surface Deflections of a Layered Viscoelastic Structure with a Bedrock at 3 m below
Surface at High Temperature Simulated Using (a) ViscoWave and (b) ADINA ....... 80
Figure 7-1 Deflection Time Histories from Pavement Structure on a Halfspace; (a) Reference
Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3..................................... 90
Figure 7-2 History of Root Mean Square Error for (a) Seed Set 1, (b) Seed Set 2, and (c) Seed
Set 3 ............................................................................................................................ 91
Figure 7-3 Backcalculated Deflection Time Histories from Pavement Structure on a Halfspace;
(a) Reference Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3 .............. 92

viii

Figure 7-4 Convergence of Power Function Parameters D0; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Halfspace .................................................. 94
Figure 7-5 Convergence of Power Function Parameters D1; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Halfspace .................................................. 95
Figure 7-6 Convergence of Power Function Parameters m; (a) Seed Set 1, (b) Seed Set 2, and (c)
Seed Set 3 – Pavement Structure on a Halfspace ....................................................... 96
Figure 7-7 Convergence of Base Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set 3 –
Pavement Structure on a Halfspace ............................................................................ 97
Figure 7-8 Convergence of Subgrade Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set 3
– Pavement Structure on a Halfspace ......................................................................... 98
Figure 7-9 Deflection Time Histories from Pavement Structure on a Bedrock; (a) Reference
Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3................................... 100
Figure 7-10 History of Root Mean Square Error for (a) Seed Set 1, (b) Seed Set 2, and (c) Seed
Set 3 .......................................................................................................................... 101
Figure 7-11 Backcalculated Deflection Time Histories from Pavement Structure on a Bedrock;
(a) Reference Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3 ............ 102
Figure 7-12 Convergence of Power Function Parameters D0; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Bedrock .................................................. 104
Figure 7-13 Convergence of Power Function Parameters D1; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Bedrock .................................................. 105
Figure 7-14 Convergence of Power Function Parameters m; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Bedrock .................................................. 106
Figure 7-15 Convergence of Base Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set 3 –
Pavement Structure on a Bedrock ............................................................................. 107
Figure 7-16 Convergence of Subgrade Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set
3 – Pavement Structure on a Bedrock ....................................................................... 108
Figure 7-17 Time-Temperature Shift Factors ............................................................................. 111
Figure 7-18 Creep Compliance Curves at Multiple Temperatures ............................................. 111
Figure 7-19 Simulated FWD Deflection Time Histories for (a) 0⁰C, (b) 10⁰C, (c) 20⁰C, and (d)
40⁰C .......................................................................................................................... 113
Figure 7-20 Seed FWD Deflection Time History for All Temperatures .................................... 114
ix

Figure 7-21 History of Root Mean Square Error ........................................................................ 115
Figure 7-22 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – Pavement
Structure on a Halfspace ........................................................................................... 116
Figure 7-23 Convergence of (a) Base and (b) Subgrade Layer Moduli – Pavement Structure on a
Halfspace................................................................................................................... 117
Figure 7-24 Simulated FWD Deflection Time Histories for (a) 0⁰C, (b) 10⁰C, (c) 20⁰C, and (d)
40⁰C .......................................................................................................................... 118
Figure 7-25 Seed FWD Deflection Time History for All Temperatures .................................... 119
Figure 7-26 History of Root Mean Square Error ........................................................................ 120
Figure 7-27 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – Pavement
Structure on a Bedrock.............................................................................................. 121
Figure 7-28 Convergence of (a) Base and (b) Subgrade Layer Moduli – Pavement Structure on a
Bedrock ..................................................................................................................... 122
Figure 8-1 Laboratory Dynamic Modulus Mastercurves............................................................ 125
Figure 8-2 Laboratory Time-Temperature Shift Factors ............................................................ 125
Figure 8-3 Measured FWD (a) Load and (b) Deflection Time Histories ................................... 126
Figure 8-4 History of Root Mean Square Error for the 3 Layer System .................................... 128
Figure 8-5 Simulated Time Histories for the 3 layer System with (a) Seed and (b) Backcalculated
Values ....................................................................................................................... 129
Figure 8-6 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – 3 Layer
System ....................................................................................................................... 130
Figure 8-7 Convergence of Unbound Pavement Layer Moduli for the 3 Layer System (a) Base
and (b) Subgrade ....................................................................................................... 131
Figure 8-8 (a) Seed and (b) Backcalculated Creep Compliance Functions from 3 Layer System
................................................................................................................................... 132
Figure 8-9 Backcalculated Dynamic Modulus Mastercurves at Reference Temperatures (a) 58⁰C
and (b) 10⁰C – 3 Layer System ................................................................................. 132
Figure 8-10 History of Root Mean Square Error for the 5 Layer System .................................. 133

x

Figure 8-11 Simulated Time Histories for the 5 layer System with (a) Seed and (b)
Backcalculated Values .............................................................................................. 135
Figure 8-12 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – 5 Layer
System ....................................................................................................................... 136
Figure 8-13 Convergence of Unbound Pavement Layer Moduli for the 3 Layer System (a) Base,
(b) Subgrade, (c) Embankment 1, and (d) Embankment 2 ....................................... 137
Figure 8-14 (a) Seed and (b) Backcalculated Creep Compliance Functions from 5 Layer System
................................................................................................................................... 139
Figure 8-15 Backcalculated Dynamic Modulus Mastercurves at Reference Temperatures (a)
58⁰C and (b) 10⁰C – 5 Layer System ....................................................................... 139

xi

CHAPTER 1 - INTRODUCTION

1.1. General

The use of a Falling Weight Deflectometer (FWD) is one of the most frequently employed
nondestructive testing (NDT) methods for evaluating the structural integrity of an existing
pavement. As its full name implies, the FWD is equipped with a falling mass mechanism
capable of inducing an impact load on the pavement surface. Due to the nature of the impact
load generated by a falling mass, the load typically has a short duration (usually 20 ms to 40 ms)
and gives rise to a stress wave that propagates through the pavement structure. The resulting
time dependent response of the pavement structure or more specifically, the vertical deflection at
the pavement surface resulting from the stress wave is measured at various radial distances from
the load and is recorded for the structural analysis of the pavement system.

The process of estimating material parameters from the FWD data can be categorized as an
inverse problem whose objective is to determine the system characteristics (e.g., layer modulus)
from the known input (e.g., applied load) and output (e.g., measured deflection).

In the

pavement engineering community, such inverse problems have typically been solved using a
procedure commonly referred to as backcalculation.

In general, backcalculation of layer parameters is carried out by matching the FWD load and
deflection to those from a theoretical model. Therefore, as is the case for most inverse problems,
the two crucial components of a backcalculation methodology are (1) a forward solution or a
theoretical model capable of simulating the FWD load and deflection, and (2) an iterative or
1

statistical routine capable of determining the optimum layer parameters that minimize the error
between the measured and simulated results.

1.2. Problem Statement

The forward models adopted for use in the dynamic backcalculation methodologies are typically
based on the solutions that stem from the theory of elastic wave propagation. One of the
common characteristics of these available solutions is that the analytical developments for
solving the wave equations were made in the frequency domain. As a result, these solutions
commonly utilized the Discrete Fourier Transform (DFT) algorithm for converting the load and
deflection signals from the time domain to the frequency domain.

However, there has been varying levels of success in dynamic backcalculation using the
theoretical time histories generated by the forward solution. Furthermore, dynamic
backcalculation using the field FWD data has generally shown fair to poor level of success and
still remains a challenge in the field of pavement engineering. Since it is recognized that the
difficulties in dynamic backcalculation arose from the forward solutions derived in the frequency
domain and from the use of the discrete transforms, it is emphasized that there is a need for a
more suitable forward solution for dynamic backcalculation in time domain. It is believed that
such a solution should not only eliminate the difficulties inherent to the solutions derived in the
frequency domain but also allow for modeling the wave propagation phenomenon while
accounting for the fundamental viscoelastic properties in a more appropriate manner.

2

1.3. Research Objective

The primary objective of this research is to develop a new forward solution that could be used
for dynamic backcalculation in time domain. In order to overcome some of the drawbacks
related to the discrete transforms, the new solution will utilize continuous integral transforms that
are more appropriate for transient, nonperiodic time domain signals such as the FWD time
histories.

Then, the resulting algorithm is to be verified against some of the existing solutions that are
capable of modeling the FWD time histories. In addition, sensitivity analyses will be conducted
in order to understand the behavior of various pavement structures subjected to FWD loading.

The research will also look into using the developed algorithm as a forward engine for dynamic
backcalculation.

A feasible backcalculation algorithm will be recommended based on the

lessons learned from the sensitivity analyses.

1.4. Report Layout

The remainder of the dissertation is organized as follows:

Chapter 2 provides a brief review of the existing dynamic solutions typically used for
backcalculation. It also addresses some of the difficulties that are present in the backcalculation
procedures frequently being used currently.

3

Chapter 3 provides a brief review of the uniaxial theory of linear viscoelasticity. The
fundamental properties of a viscoelastic material will be described in both the time and the
frequency domains. Then, the time-temperature superposition of a thermo-rheologically simple
(viscoelastic) material will be reviewed. In addition, a review of the analytical models frequently
used for modeling the viscoelastic constitutive relations is provided with a discussion on the
interconversion of various viscoelastic properties. Finally, the chapter concludes with a brief
discussion on the elastic-viscoelastic correspondence principle in the uniaxial mode.

Chapter 4 builds upon the uniaxial theory of viscoelasticity in Chapter 3 and provides a brief
review of the multi-axial theory of viscoelasticity. The chapter begins with an introduction to the
viscoelastic Poisson’s ratio which is critical in the development of the multi-axial viscoelasticity
and the extended elastic-viscoelastic correspondence principle. Then, a methodology for
applying the multi-axial viscoelasticity in the analysis of the laboratory Indirect Tensile (IDT)
test is presented. This chapter is critical since the viscoelastic properties backcalculated from the
FWD data are compared to those obtained from the IDT tests.

Chapter 5 describes in detail, the mathematical development of the new algorithm. The
theoretical development for the proposed methodology follows similar steps to those used for the
development of LAMDA, the spectral element method which utilized the discrete transforms for
solving the wave equations. However, the proposed solution utilizes the continuous integral
transforms (namely Laplace and Hankel transforms) that are more appropriate for transient,
nonperiodic signals.

4

Chapter 6 presents the verification results for the developed algorithm which is done through
comparing the simulation results from the developed algorithm to some of the other existing
solutions.

Chapter 7 provides a brief mathematical background on the iteration schemes selected for use
with the new forward solution developed in this research. Then the iteration algorithms are tested
using theoretically generated deflection time histories. The results of the theoretical
backcalculation provided in this chapter for both the single and multiple temperature
backcalculation.

Chapter 8 presents the results of the preliminary backcalculation exercise using the fieldmeasured FWD data. Then the backcalculated results are compared to those obtained from the
laboratory.

Chapter 9 summarizes the findings of the current research and some recommendations for future
research.

5

CHAPTER 2 - LITERATURE REVIEW

Evidently, the forward solutions used for the static backcalculation procedures assume that the
pavement structure under the applied load is in static equilibrium.

In other words, these

solutions do not allow for simulating the FWD load and deflection as functions of time.
Correspondingly, the material properties adopted for use in these solutions have been acquired
from the constitutive models that do not depend on time – e.g., linear elastic (Uzan et. al., 1988)
and nonlinear elastic (Irwin, 1977) models. Because the time dependent solution cannot be
achieved, only the peak magnitude of the impact load and the peak deflection measured at
different sensor locations (also known as the deflection basin or deflection bowl) are taken into
account in the backcalculation process. Although the static backcalculation methods are
extremely efficient, these solutions do not properly account for the dynamic (time dependent)
nature of the FWD load and deflection.

As it has been proven by numerous researchers in the past, asphalt mixtures that comprise the top
layer of flexible pavements are viscoelastic in nature. Unlike an elastic material, because their
fundamental properties are time or frequency dependent, viscoelastic materials such as an asphalt
mixture show time dependent response even under a static (or constant) load (A review of the
theory of linear viscoelasticity will be provided in the next chapter and hence it will not be
discussed extensively as part of this chapter). The forward solution such as the ones recently
developed by Kim (2011) and Kutay et. al. (2011) accounts for such material time dependency
through the use of viscoelastic constitutive relations. The primary advantage of such forward
models is that they allow for simulating the FWD time histories without a significant loss of

6

computational efficiency. However, as it was pointed out by Kutay et. al. (2011), because the
viscoelastic solutions do not consider the effect of wave propagation attributed to the impact load,
they are not capable of modeling the different time delays measured at different sensors as well
as the free vibration of the pavement structure which may be significant in the presence of a stiff
layer (e.g., a bedrock) at shallow depth.

Some of the existing dynamic forward and inverse solution pairs are summarized in Table 2-1.
The forward models adopted for use in the dynamic backcalculation methodologies are typically
based on the solutions that stem from the theory of elastic wave propagation. Therefore, the
stress wave propagation is naturally modeled in these solutions. The viscoelasticity of the
material typically has been taken into consideration through the use of a damping ratio which is a
concept derived from the theory of vibrations (Chatti and Yun, 1996, Chatti et. al., 2004, Ji, 2005,
Matsui et. al., 2011). Excluding the solutions based on the Finite Element Analysis (FEA) that
are generally more time consuming and inefficient for backcalculation (Hadidi and Gucunski,
2007), the dynamic forward solutions can be categorized into the finite layer type (Chatti and
Yun, 1996, Chatti et. al., 2004, Ji, 2005, Matsui et. al., 2011) or the spectral element type
solutions (Al-Khoury et. al., 2001a, 2001b, & 2002, Grenier and Konrad, 2007). However, one
of the common characteristics of both solution types is that the analytical developments for
solving the wave equations were made in the frequency domain. As a result, these solutions
commonly utilized the Discrete Fourier Transform (DFT) algorithm for converting the load and
deflection signals from the time domain to the frequency domain.

7

Table 2-1 Existing Dynamic Forward Solution and Backcalculation Pairs
Forward
Program Name

Calculation

Backcalculation Method

Domain

Developed By

Program
PAVE-SID

SCALPOT

System Identification (SID)

Frequency

FEDPAN

SAP IV

Linear Least-Squares

Time

No Name

UTFWIBM

Newton’s Method

Frequency/Time

BKGREEN

GREEN

Nonlinear Least-Squares

Frequency

No Name

FEA

Gauss-Newton

Time

No Name

SAPSI

Levenberg-Marquardt

Frequency

Secant Update,
No Name

LAMDA

Levenberg-Marquardt,

Frequency

Powell Hybrid
DYNABACK-

SAPSI

Newton Raphson Method

Frequency/Time

No Name

FEA

Neighborhood Algorithm

Time

UCODE

DYNAPAV-UL

Levenberg-Marquardt

Time

F/T

Magnuson, et. al.
(1991)
Ong, et. al.
(1991)
Ujan, et. al.
(1994)
Kang, Y. V.
(1998)
Matsui, et. al.
(1998)
Losa
(2002)
Al-Khoury, et. al.
(2002)
Ji
(2005)
Hadidi and
Gucunski (2007)
Grenier and
Konrad (2007)

As it was pointed out by Bendat and Piersol (2010), the aforementioned discrete transform is not
appropriate for transient nonperiodic signals such as those generated by the FWD. Some of the
drawbacks and difficulties arising from the use of the discrete transform that are documented in
8

the dynamic backcalculation literature can be summarized as: (1) the truncation in the FWD time
histories – i.e., the recording of the signal being terminated before the pavement system has
come to a rest, which acts like a box filter applied to the signal, (2) the periodicity of the signal
assumed in the discrete DFT algorithm not being able to accurately disclose the frequency
content of the transient FWD time histories with short duration, (3) the DFT being very sensitive
to noise that is always present in the FWD data (Chatti et. al., 2004, Ji, 2005, Matsui et. al.,
2011), and (4) the DFT being impractical for representing the fundamental properties of a
viscoelastic material such as creep compliance or dynamic modulus, due to a large number of
harmonics necessary for modeling them (Zhang et. al., 1997a, 1997b) – this is also the reason
behind most of the frequency domain solutions utilizing the damping ratio concept and hence
failing to model and/or backcalculate the fundamental properties of a viscoelastic material .

9

CHAPTER 3 - REVIEW OF UNIAXIAL THEORY OF LINEAR VISCOELASTICITY

Prior to the presentation of the mathematical formulation of the developed algorithm, it is
necessary to review the constitutive relation of a viscoelastic material. Therefore, this chapter is
dedicated to providing a theoretical review on the fundamental properties of a viscoelastic
material and the theory of linear viscoelasticity.

3.1. Fundamental Properties of a Viscoelastic Material in Time Domain

According to the theory of linear viscoelasticity, the time-dependent stress-strain relationship of
a viscoelastic material subjected to a constant uniaxial stress, , can be expressed as the creep
compliance, D(t):

Dt  

 t 
for t  0
0

(3-1)

where (t) is the resulting time dependent strain. The strain resulting from any stress history can
be expressed in the form of a convolution integral by means of the Boltzman’s superposition
principle as:

t

 t    Dt   
0

where  is the integral variable.
10

  
d


(3-2)

On the other hand, the relaxation modulus, E(t), is defined as the time dependent stress resulting
from an applied step of constant unit strain, that is:

E t  

 t 
for t  0
0

(3-3)

where (t) is the resulting time dependent stress and  is the magnitude of the constant strain.
With the relaxation modulus defined as above, the stress history resulting from any strain history
can be expressed as:

t

 t    E t   
0

 t 
d


(3-4)

3.2. Fundamental Properties of a Viscoelastic Material in Frequency Domain

In the frequency domain, the sinusoidal oscillations of stress and strain are commonly
represented by complex variables through Euler’s formula:

 t    0 eit   0 cos t  i sin t 

(3-5)

where 0 and  are the amplitude and the angular frequency of the axial stress oscillation. The
resulting strain, (t) , at a steady state oscillates at the same frequency as the stress oscillation but
with a phase angle:

11

 t    0  eit     * eit

(3-6)

where 0 is the strain amplitude,  is the phase angle by which the strain lags behind the stress,
and * is a complex strain amplitude defined as:

 *   0  e i   0 cos   i sin  

(3-7)

The complex compliance defined as the strain over stress can be obtained from Equations (3-5)
and (3-7), analogous to the creep compliance as follows:

D    D   D  

 *  0  i

e
0 0

(3-8)

where D  and D  are the real and imaginary parts of the complex compliance. Similarly,
if the viscoelastic material is subjected to a sinusoidal strain oscillation, that is:

 t    0 eit   0 cos t  i sin t 

(3-9)

Then, the resulting stress is:

 t    0  eit     * eit

where * is the complex stress amplitude defined as:
12

(3-10)

 *   0  ei   0 cos   i sin  

(3-11)

The complex modulus is defined as the ratio between the complex stress amplitude, *, over the
input strain amplitude 0:

E     E    E   

 *  0 i

e
0 0

(3-12)

where E   and E   are the real and imaginary parts of the complex modulus.

It is also worthwhile to note that within the field of pavement engineering, the magnitude of the
complex modulus given as the following equation has frequently been referred to using the term
“Dynamic Modulus” (Kim, 2009).

E    

E  2  E  2

(3-13)

3.3. Relationship between Viscoelastic Properties

Since creep and stress relaxation phenomena are two aspects of the same viscoelastic behavior of
a given material, they are obviously related. In other words, the viscoelastic properties presented
above are not independent of each other. However, determining the relationship between the
fundamental properties in time domain becomes quite tedious (although possible) due to the
convolution integrals shown in Equations (3-2) and (3-4). Instead, it is more convenient to derive
13

the relationships in the Laplace domain. Taking the Laplace transform on Equations (3-2) and (34) yields the following:













 s   sDs   s 

(3-14)

 s   sEs   s 

(3-15)



where s is the Laplace variable and L[ f (t )]  f s    e  st f (t )dt is the Laplace transform of
0

the function f t  . By combining the two equations above, the relationship between creep
compliance and relaxation modulus can be expressed as an algebraic equation in the Laplace
domain as follows:

ˆ
ˆ
D( s ) E ( s )  1 / s 2

(3-16)

The inverse Laplace transform of the above equation reveals the relationship between the creep
compliance and relaxation modulus in time domain:

t

t

 E t   D d   E  Dt   d  t

0

(3-17)

0

In the frequency domain, the relationship between the complex compliance and the complex
modulus can easily be derived by combining Equations (3-8) and (3-12) as:

14

D    E     1

(3-18)

In addition, the relationship between the complex compliance and the creep compliance can be
obtained by substituting Equations (3-5) and (3-7) into Equation (3-2). After rearranging the
variables, the final relationship can be written as (Findley et al, 1976):


D    i  Ds  s  i

(3-19)

The relationship between the complex modulus and the relaxation modulus can be obtained in a
similar manner. Substituting Equations (3-9) and (3-11) into Equation (3-4) results in the
following upon summarizing:


E     i  E s  s  i

(3-20)

Although the detailed derivation for Equations (3-19) and (3-20) is omitted here, similar
derivation will be shown in the following chapter for the viscoelastic Poisson’s ratio (Section
4.1.3).

3.4. Time-Temperature Superposition Principle

As it has been reported by several researchers in the past, the behavior of a viscoelastic material
is strongly affected by temperature (Findley et. al., 1976, Tschoegl, 1989, Wineman and
Rajagopal, 2000). However, the constitutive relations provided so far did not account for the
15

effect of temperature (since the fundamental properties were assumed to be function of time only)
and hence, they should be used only when the material is subjected to a constant temperature. In
order to account for the effect of temperature, it is necessary to generalize the fundamental
properties described above to include temperature as an independent variable. For example, the
relaxation modulus can now be defined as:

E  E t , T 

(3-21)

where the uppercase T has been used for temperature to distinguish from the time variable, t.

Fortunately, there has been sufficient theoretical and experimental evidence showing that most
linearly viscoelastic materials obey the Time-Temperature Superposition Principle (TTSP) which
allows for combining the effect of time and temperature (Findley et. al., 1976, Tschoegl, 1989,
Wineman and Rajagopal, 2000). A viscoelastic material obeying the TTSP is also frequently
referred to as a “Thermorheologically Simple” material. According to the TTSP, the above
relaxation modulus can be written as the following:

E  E t , T   E  , T0 

(3-22)

where T0 is the reference temperature andis called the “Reduced Time” which is related to the
physical time in the following manner:

  t aT T 
16

(3-23)

where aT is the shift factor. The above equation can be also be written as:

log   logt aT T   logt   logaT T 

(3-24)

The above equations indicate that fundamental property of a viscoelastic material at the time t
and at temperature T is equal to that property at the reduced time  and at the reference
temperature T0. This in turn, indicates that the fundamental properties at different temperatures
can be shifted by the amount log(aT) when plotted in log-scale (or stretched/shrunk in arithmeticscale), as shown conceptually in Figure 3-1.

As shown in the figure, once the fundamental properties obtained from multiple temperatures
have been shifted to their equivalent at the reference temperature T0, it results in a curve
representing the viscoelastic material property for a wide range of frequencies (or time) and is
called the mastercurve of a viscoelastic material.

In reality, testing of a viscoelastic material for its characterization is typically conducted at a
limited range of frequencies or for a short duration of time. This is due to (1) the limitations of
the testing equipment (e.g., the duration of the FWD load is typically between 20 ms and 40 ms)
and (2) the specimen being damaged when tested for a long period of time (or a wide range of
frequencies). In that sense, the TTSP provides a great deal of advantage for characterization of
viscoelastic material properties because it provides a methodology for overcoming this problem.
In other words, instead of testing the material at a wide range of frequencies or time, the TTSP

17

allows for testing the material at multiple temperatures (and at a limited range of frequencies or
time) and providing the same mastercurve.

Range of Testing
Frequency

Dynamic Modulus

Horizontal Shift,
Log(aT)

Mastercurve

Log(Reduced Frequency )
Figure 3-1. Concept of Time-Temperature Superposition Principle and Mastercurve (For
interpretation of the references to color in this and all other figures, the reader is referred to the
electronic version of this dissertation).

18

3.5. Analytical Models for Viscoelastic Properties

Numerous closed form equations have been proposed and used to model the fundamental
properties of a viscoelastic material and they will be briefly reviewed in this section. Although
the equations will be presented in the form of creep compliance or relaxation modulus, they can
generally be used to model both (Kim, 2009).

One of the simplest equations for modeling the viscoelastic property is the pure power function
shown below:

Dt   D1  t m

(3-25)

where D1 and m are the power function parameters. The Laplace transform of the above equation
is given as:


D (m  1)
Ds   1
s m 1

(3-26)

where mis the gamma function.

Although the above equation is very simple, the critical disadvantage of the pure power function
is that it does not account for the instantaneous (or elastic) response of the viscoelastic material
which may become important at low temperatures.

19

Alternatively, the generalized power function still maintains a simple form while representing
the viscoelastic property effectively. The function is expressed as:

Dt   D0  D1  t m

(3-27)

where D0, D1, and m are the power function parameters. The disadvantage of the above equation
is that it fails to model the entire creep compliance (or dynamic modulus) mastercurve which
spans over a wide range of temperature and frequency. Nevertheless, it has been shown that the
generalized power function is an excellent analytical representation of the uniaxial viscoelastic
creep compliance at a single temperature (Roque et al. 1997, Roque et al. 1998, Kim et al. 2005,
and Kim et al. 2008). In addition, it has been shown by some of the previous research studies
(Lee et. al., 2012, and Kim et. al. 2010) that the generalized power function can be used
successfully to fit the laboratory test data for each temperature prior to fitting the master curve
using a sigmoidal function proposed by the latest Mechanistic-Empirical Pavement Design
Guide (MEPDG) (ARA, 2004). Therefore, it was decided to use the generalized power function
for modeling the viscoelastic material behavior in this research. Taking the Laplace transform on
the above equation results in:


D s m  D1(m  1)
Ds   0
s m 1

(3-28)

On the other extreme, the modified power function shown below is one of the most complicated
closed form equations that have been used for modeling the viscoelastic creep compliance.
20

M

Di

Dt   D0  

i 1 1   i





t 


n

(3-29)

where D0, Di, i, and n are the parameters for the modified power function. Although the above
function is generalized enough so that the entire creep compliance mastercurve can be fitted,
taking the Laplace transform of the above equation is close to impossible due to its complicated
form and the number of parameters involved.

The recent Mechanistic-Empirical Pavement Design Guide (MEPDG) has adopted the sigmoidal
function for representing the mastercurve which is given as the following (ARA, 2004):

log  E *    



1  e   log 

(3-30)

where , , , and  are the parameters for the sigmoidal function. The advantage of the above
equation is that the mastercurve can be fitted accurately with less number of parameters when
compared to the modified power function. Nonetheless, as it was the case for the modified
power function, the Laplace transform of the above cannot be obtained in closed form.

3.6. Prony Series

Although the analytical functions introduced above have been successfully used as fitting
functions that represent the viscoelastic behavior of materials, mathematical difficulties arise

21

when attempting the strict analytical interconversion of the fundamental properties necessary for
the stress or strain development under various loading conditions in viscoelastic media. For
example, the Laplace transforms of the modified power function and the sigmoidal function
cannot be obtained in closed form and hence the interconversion equations presented earlier
cannot be used directly.

Although the Laplace transforms of the pure power function and the generalized power function
can be obtained analytically, they still fail to provide an analytical interconversion of the
viscoelastic properties. For example, substituting Equation (3-27) into Equation (3-16) results in
the following for the relaxation modulus in the Laplace domain whose inverse transform cannot
be carried out analytically:


E s  

s m 1
D0 s m  D1(m  1)

(3-31)

Furthermore, substituting Equation (3-27) into Equation (3-19) results in a non-integer power of
the (i) term (since 0 < m < 1), which in turn yields a non-unique solution for the complex
compliance (Findley et. al., 1976).

Instead, Prony series (generalized model), which has one Maxwell model (a spring element and
dashpot element connected in series) and several Kelvin elements (a spring element and dashpot
element connected in parallel), has been widely used for analytical representation of viscoelastic
materials due to its remarkable computational efficiency. For creep compliance, the Prony series
representation is of the following form:
22

D(t )  D0 

1



N



t

t   Di (1  e  i )

(3-32)

i 1

where D0, and Di are Prony series parameters,  is the dashpot constant, and i are retardation
times.

Prony series is convenient for viscoelastic analysis in cases where the stress history is prescribed,
whereas the generalized Maxwell model, which is several Maxwell elements connected in
parallel, is rather convenient for predicting stress associated with a prescribed strain. Therefore,
converting creep compliance to relaxation modulus is the process of converting the Prony series
to the generalized Maxwell model. Taking the Laplace transform on the above equation results
in:

N

D
Di
1
D( s )  0 

s   s 2 i 1 s   i  s  1

(3-33)

Substituting the above into Equation (3-16) yields the following equation for the relaxation
modulus in the Laplace domain:


E ( s) 

1
N


Di
1

s D0 


  s i 1  i  s  1 



The above equation can be rearranged as a ratio of two polynomials of s:
23

(3-34)

M 1
 aM  2  s M  2  a1 )
ˆ ( s)  (a M 1  s
E
(bM  s M  bM 1s M 1  b1 )

(3-35)

where, a and b are the coefficients of polynomial functions. Expanding Equation (3-35) by
partial fractions yields the following form:

ˆ
E (s) 

E1
E2
EM


1
1
1
s
s
s

1

2

(3-36)

M

where, 1/1, 1/2, etc., are the roots of the denominator in Equation (3-35), and E1, E2, etc., are
numerators that satisfy the partial fractions (the terms i are frequently called the relaxation times
in the literature). Performing inverse Laplace transformation finally yields the generalized
Maxwell model in parallel (Prony series for relaxation modulus):

E (t )  E1  e



t

1

 E2  e



t

2

  EM  e



t

M

M

  Ei  e



t

i

(3-37)

i 1

Frequently, the Prony series of the following form is used in which the term including the
dashpot constant is eliminated by setting   0 in Equation (3-32), that is:

N



t

D(t )  D0   Di (1  e  i )
i 1

24

(3-38)

Following the same procedure above, it can be shown that the relaxation modulus corresponding
to Equation (3-38) is obtained as:



M

t

E (t )  E0   Ei  e  i

(3-39)

i 1

Using the Prony series presented above, their frequency domain counterparts can also be derived
easily. Substituting Equation (3-33) into Equation (3-19) results in the Prony representation of
the complex compliance:

N

Di

D( )  D0  

2 2
i 1  i  1

D( ) 

1



N



 i Di

2 2
i 1  i  1

(3-40)

(3-41)

Similarly, the Prony representation for the complex modulus can be obtained by substituting
Equation (3-36) into Equation (3-20). Upon summarizing, the outcome is given as:

M

E ( )  

 2  i2 Ei

2 2
i 1  i  1

M

E ( )  

i Ei

i 1

2 2
i

25

1

(3-42)

(3-43)

3.7. Elastic-Viscoelastic Correspondence Principle

The uniaxial constitutive equation for a linear elastic material is given as the well known
Hooke’s law:

  E

(3-44)

where , , and E are the stress, strain, and the Young’s modulus, respectively. Alternatively, the
above equation can also be written in terms of the elastic compliance D=1/E as:

  D

(3-45)

By comparing Equation (3-44) to Equations (3-15) and (3-12), it can be seen that replacing the


Young’s modulus, E, in Equation (3-44) by sE s  and E* results in the uniaxial viscoelastic
constitutive relation in the Laplace domain as shown in Equation (3-15) and in the frequency
domain as shown in Equation (3-12), respectively. Similarly, replacing the elastic compliance, D,


in Equation (3-45) with sDs  yields the viscoelastic equation shown in Equation (3-14) while
replacing it with the complex compliance, D*, results in the constitutive equation shown in
Equation (3-8). This indicates that the viscoelastic constitutive relation can be obtained easily in
the Laplace domain or in the frequency domain if the elastic constitutive relation is determined
and the elastic constants are replaced with their appropriate viscoelastic counterparts. This is
known as the elastic-viscoelastic correspondence principle. Although it seems obvious and trivial

26

in the uniaxial case shown here, the elastic-viscoelastic correspondence principle will prevail a
great advantage in the multi-axial viscoelasticity which will be discussed later.

27

CHAPTER 4 - EXTENSION TO MULTI-AXIAL VISCOELASTICITY AND ANALYSIS
OF INDIRECT TENSILE TEST DATA

In this chapter, the theory of viscoelasticity introduced in the previous chapter will be extended
to multi-axial viscoelasticity. After a brief theoretical development, the methodology for
analyzing the laboratory Indirect Tensile (IDT) test data will be presented.

4.1. Viscoelastic Poisson’s Ratio

4.1.1. Viscoelastic Poisson’s Ratio in Time Domain

Tschoegl (1989), Wineman and Rajagopal (2000), and Tschoegl et al (2002) emphasized that the
time-dependent Poisson’s ratio,  t  , of a linearly isotropic viscoelastic material should be
0
defined as the ratio of the time dependent lateral strain,  2 t   to a constant axial strain,  1 . By

the definition, Poisson’s ratio is expressed as:

 t   

 2 t 
0
1

for t  0

(4-1)

The above definition implies that the theoretical viscoelastic Poisson’s ratio can be obtained
from a uniaxial relaxation test in which the strain in the axial direction is held constant. The
definition also reveals that Poisson’s ratio cannot be obtained directly by taking the ratio of the
time-dependent strains resulting from a uniaxial creep test because in this testing mode, the axial

28

stress is held constant rather than the axial strain. Hence, the ratio of 1 t  and  2 t  directly
measured from a creep test cannot be termed Poisson’s ratio (Tschoegl et al, 2002).

From the definition of the Poisson’s ratio shown in Equation (4-1), one can observe that the
Poisson’s ratio is the lateral strain as a function of time resulting from a unit step strain applied in
the axial direction. Therefore, for any given time-dependent axial strain function, 1 t  , the
resulting lateral strain,  2 t  , is given as a convolution integral (observe the similarities between
Equations (3-3) , (3-4), (4-1), and(4-2)):

t

 2 t    t   
0

1  
d


(4-2)

Taking Laplace transform on the above equation yields:







 2 s   s s 1 s 

(4-3)


Equation (11) can also be expressed in terms of the axial stress,  1 s  , by means of Equation (3):









 2 s   s 2 s  Ds  1 s 

29

(4-4)

4.1.2. Viscoelastic Poisson’s Ratio in Frequency Domain

Similar to Poisson’s ratio in the time domain that was defined under the uniaxial relaxation
testing mode, Poisson’s ratio in the frequency domain can be defined under a uniaxial, straincontrolled complex modulus testing mode as (Tschoegl, 1989 and Tschoegl et al., 2002):

     

0
 2 ei (t  )
0
1 eit

0 i 

 e
 2

0
1

(4-5)

0
where  2   is the amplitude of the lateral strain and   is the phase angle between 1 t  and

 2 t  .

Under a uniaxial, stress-controlled complex modulus test, the axial strain, 1 t  , lags behind the
axial stress,  1 t  , with a phase lag . At the same time, the lateral strain,  2 t  , lags behind the
axial strain, 1 t  , which already has a phage lag. Therefore, the lateral strain can be expressed as:


0

 2 t    2  ei t      2, p  ei t     2,t  ei t

(4-6)

, p
In the above equation, two complex amplitudes are introduced. First,  2 is the complex

amplitude that contains only the phase lag between 1 t  and  2 t  :


0
 2, p   2 e i 

30

(4-7)


whereas  2, t is the complex amplitude that includes both phase angles:



0
 2,t   2, p e i   2 e i  

(4-8)

Then, the definition of the complex Poisson’s ratio,  , can be rewritten in terms of the complex
strains as:

    



 2, t

1



0
 2 e  i  
0
1 e  i



0
 2 e  i 
0
1




 2, p
0
1

(4-9)

It should be noted that the last term shown in Equation (4-9) is the same as the form of that

shown in Equation (4-1). After substituting Equation (3-8) into the above equation for  1 , the

complex lateral strain can be described as a function of the complex compliance and complex
Poisson’s ratio as follows:



0
 2,t       1       D   1

(4-10)

4.1.3. Relationship Between Time and Frequency Dependent Poisson’s Ratios

In order to find a relationship between  t  and     , a new integration variable    t   is
substituted into Equation (4-2):

31

t

 2 t     '
0

1 t   '
d 
 '

(4-11)

Since a stress term is not included in the above equation, it is more convenient to use the
0
complex strain function, 1 e it , than Equation (4-8) because the phase angle between the axial

stress and strain is not included in this equation. According to Euler’s formula, the complex
strain function can be expressed as:

0
0
1 t   1 eit  1 cos t  i sin t 

(4-12)

or

0
0
0
1 t   '  1 ei (t  ')  1 eit e i'  1 eit cos  'i sin  '

(4-13)

Substituting the above strain into Equation (4-11) yields:

t

 2 t     
0





 0 i t
 e cos    i sin   d 
  1

(4-14)

The lateral strain,  2 t  , in the above equation is a result of an axial strain oscillation,
0
1 t   1 eit , and can be expressed as the lateral complex strain function with only the phase

angle between the two strains, that is:
32

*,
 2 t    2 p ei t

(4-15)

Therefore, from Equations (4-14) and (4-15),

t

*,
 2 p e i t    
0





 0 i t
 e cos    i sin   d 
  1

(4-16)

0
After canceling e i t and dividing both sides of the above equation by  1 , the left-hand side of

the equation is left with the complex Poisson’s ratio due to Equation (4-9).

t

      
0
t


cos    i sin  d 
 
t

    sin  d   i    cos  d 
0

0
t

t

 i     cos  d   i   sin  d 
0

0



After changing the upper limit to infinity:

33

(4-17)







0


0




   i     cos  d   i   sin  d 


 i   cos    i sin  d 

(4-18)

0


 i     e  i  d 
0

The last equality in the above equation indicates that the relationship between the Poisson’s ratio
functions in the two domains is obtained as:



     i  s  s  i

(4-19)

4.2. Extension of the Elastic-Viscoelastic Correspondence Principle

In the previous chapter, it was shown that for the uniaxial case, the viscoelastic constitutive
equations can be obtained by replacing the elastic constant by its appropriate viscoelastic
counterpart.

Equation (4-4) indicates that in order to apply the elastic-viscoelastic

correspondence principle for the multiaxial case, the elastic compliance, D, and the Poisson’s



ratio, , should be replaced by sDs  and s s  , respectively. On the other hand, Equation (410) shows that the above elastic constants should be replaced by D* and * in order to yield the
viscoelastic constitutive relations in the frequency domain. For example, starting with the threedimensional constitutive equation for an elastic material given as:

34













 x  D   x   y   z
 y  D   y   z   x
 z  D   z   x   y

(4-20)

The viscoelastic constitutive relation in the Laplace domain can be obtained as:

















































 x s   sDs    x s   s s  y s   s s  z s 

 y s   sDs    y s   s s  z s   s s  x s 

(4-21)

 z s   sDs    z s   s s  x s   s s  y s 

Similarly, in the frequency domain, the following is obtained:













0
0
 x *    D *     x   *   0   *   z
y
0
0
 y *    D *     0   *   z   *   x
y
0
0
 z *    D *     z   *   x   *   0
y

35

(4-22)

4.3. Determining the Viscoelastic Property from the Indirect Tensile (IDT) Test

Although the uniaxial tension or compression test is desirable for testing of the materials, the
problem associated with size requirements makes it difficult to perform the uniaxial test,
especially on cored mixtures from the thin asphalt layers.

As an alternative to the uniaxial tests, the indirect tension test (Buttlar and Roque 1994, Zhang et
al., 1997a, 1997b, and Kim et al., 2005), which is often called a Brazilian test, has been widely
used for testing both laboratory-made and field-cored mixtures because the test requires
relatively thin specimens. Therefore, the remainder of this chapter will be dedicated to
developing a methodology for obtaining the viscoelastic properties from the IDT which is a
biaxial test. This is critical because the viscoelastic property backcalculated from the FWD data
will be compared to the IDT result from the pavement cores which will serve as a ground truth.
A picture of the IDT set up and the IDT specimen geometry are shown in Figure 4-1.

36

(a)
1.0 in.

1.5 in.

6.0 in.

1.5 in.
LVDTs
Gauge Points

(b)
Figure 4-1. (a) Indirect tensile test setup and (b) specimen dimensions

37

4.3.1. Stress Distribution within an Indirect Tensile (IDT) Specimen

For an IDT specimen with radius R and thickness d subjected to a strip load of width a =
2R· ) and magnitude P as shown in Figure 4-2 (Zhang et al., 1997a, 1997b, and Kim et al.,
sin(
2005), the distributions of the tensile stress along the horizontal axis and the compressive stress
along the vertical axis were given by Hondros (1959) under the assumption of plane stress
conditions as:

 1  x R 2

1  x R 2 sin 2
2P 

 tan 1 
tan  
 1  x R 2

ad 1  2x R 2 cos 2  x R 4



2P

m x 
ad

 x  x, 0  

2P
 y 0, y   
ad


2



1   y R 2 sin 2
1  1   y R 

 tan
tan  
 1   y R 2

1  2 y R 2 cos 2   y R 4




2P
n y 
ad

38

(4-23)

(4-24)

y
a


x

R

P
(a)

y

x(0,y)

x(x,0)
x

(b)
Figure 4-2 (a) Schematic illustration of the IDT test and (b) plane stress distributions

It would be ideal if the strain values could be measured at a point located at the center of the
specimen, i.e., x = y = 0. The corresponding stresses then could be obtained easily as  x 0, 0
and  y 0, 0 from the above equations; however, the strain values are obtained based on the
displacements measured by the strain gauges with a length of l and are better represented by the
39

average or the mean of the strain distributions between the gauge points. Therefore, the average
stresses calculated over the gauge length in conjunction with the measured strains were used.
The average stresses are calculated as:

2P 1
 x, avg 

ad l

l 2

 mx dx

(4-25)

l 2

2P 1
 y, avg  

ad l

l 2

 n y dy

(4-26)

l 2

For an IDT creep test, the load P in the above equations can be substituted with
P  Pt   P 0 H t  , where P 0 is the magnitude of the creep load, and H(t) is the Heaviside step

function. Substituting this into Equations (4-25) and (4-26), and taking the Laplace transform
yields:



 x, avg s  



0
x

 y, avg s   

0
y

0
where the magnitudes of the stresses  x and  0 are:
y

40

(4-27)

s

s

(4-28)

0
x

0
y

2P 0 1


ad l

2P 0 1


ad l

l 2

 mx dx

(4-29)

l 2

l 2

 n y dy

(4-30)

l 2

Similarly, if the load is a sinusoidal oscillation with magnitude P 0 and angular frequency ,
then P  Pt   P 0 eit and Equations (4-25) and (4-26) become:

0
 x, avg t    x  eit

(4-31)

 y, avg t    0  eit
y

(4-32)

Another constant that is derived in this section is the ratio of the stress magnitudes. This ratio, ,
will be used later to simplify the equations for the viscleoastic stress-strain relation significantly:

l 2

 n y dy



 y, avg  0  l 2
y



0
l 2
 x, avg  x

(4-33)

 mx dx

l 2

With the above definition for , the relationship between the two stresses can be written as:
41

 y, avg     x, avg

(4-34)

For the remainder of the chapter, the subscript “avg” will be dropped from the above equation
and, x and y will be used to represent the average stresses over the gauge length unless stated
otherwise.

4.3.2. Indirect Tensile Creep Test

From Equation (4-21), the plane stress constitutive equations for a viscoelastic material with
time-dependent Poisson’s ratio can be written in the Laplace domain as:













(4-35)













(4-36)

 x (s)  s Ds    x s   s s  y s 

 y (s)  s Ds    y s   s s  x s 


After eliminating the term s s  from the above equations and reorganizing the variables, one
can solve for an expression for the creep compliance in the Laplace domain. Upon simplifying
the result using Equation (4-34), one arrives at:




1  x    y
D( s )   

s  x    y

42

(4-37)

For an IDT creep test, the stresses in the Laplace domain were given in Equations (4-27) and (428). After substitution of these stresses, the denominator in the above equation simply reduces to
a constant, that is:



 x    y

D( s ) 
0
 x    0
y

(4-38)

Taking the inverse Laplace transform on the above equation results in the creep compliance in
time domain:

D(t ) 

1
0
 x    0
y





  x t      y t 

(4-39)

4.3.3. Indirect Tensile Complex Modulus Test

Equation (4-22) indicates that the plane stress constitutive equations for a viscoelastic material in
the frequency domain can be written as:





(4-40)





(4-41)

0
   D   x  * 0
x
y

0
   D   0  * x
y
y

Combining Equations (4-40) and (4-41) by eliminating * yields:

43

*

D ( ) 

 *    *
x
y

(4-42)

0
 x    0
y

The complex strain amplitudes in the numerator of the above equation can be expressed in a
form shown in Equation (3-7) as:

0
0
  ( )   x  e i x   x  cos x   i sin x 
x

  ( )   0  e
y
y

i y

  

 

  0  cos  y  i sin  y
y

(4-43)

(4-44)

0
where  x and  0 are the horizontal and vertical strain amplitudes and  x and  y are the
y

horizontal and vertical phase lags measured from the IDT test. However, through IDT
experiments, Lee and Kim (2009) showed that the difference between the two phase lags  x and

 y is negligible, that is,    x   y . Substituting this into Equations (4-43) and (4-44), and the
resulting equations into Equation (4-42) results in:

*

D ( ) 

0
 x    0
y
0
 x    0
y

44

 e  i

(4-45)

CHAPTER 5 - THEORETICAL DEVELOPMENT OF A TIME DOMAIN FORWARD
SOLUTION

5.1. Governing Equations For Viscoelastic Wave Propagation

Similar to any other wave propagation problems, the proposed solution begins with the classical
equation of motion for a continuous medium given as the following (Malvern, 1969):


  σ  b  ρu

(5-1)

where  is the stress tensor, b is the vector of body forces per unit volume,  is the mass density
of the material, and u is the displacement vector. According to the theory of linear elasticity, the
stress-strain relationship for a linear, homogenous, and isotropic material is obtained from the
generalized Hooke’s law:

σ   tr ε I  2 ε

(5-2)

where  is the strain tensor,  and  are the lamé constants, and  is the identity tensor. The
strain tensor in the above equation is related to the displacement vector according to the
following:

ε



1
u  u T
2

45



(5-3)

For a viscoelastic material such as an asphalt concrete mixture, the fundamental materials
properties – in this case, the lamé constants – as well as the stresses and strains are timedependent and hence, their relationship can be written as the following in reference to the theory
of linear viscoelasticity (Schapery, 1974, Wineman,and Rajagopal, 2000, Christensen, 2003):

σ    tr ε I  2  ε

(5-4)

where the function  represents the well known Stieltjes convolution integral defined as:

t

      t   
0

  
d


(5-5)

It should be noted that the kinematic strain-displacement relationship shown in Equation (5-3)
also applies to linear viscoelastic materials. The only difference from an elastic material is that
the displacement and hence the strain are functions of not only the material (or spatial)
coordinates but also time. Substituting Equations (5-3) and (5-4) into the equation of motion
shown in Equation (5-1) and ignoring the body forces result in the following equation in terms of
displacements:


       u     2u   u

46

(5-6)

By means of the Helmholtz decomposition, the displacement vector in the above equation can be
expressed in terms of potentials as follows:

u      H

(5-7)

where  represent a scalar potential and H is a vector potential whose divergence vanishes (i.e.,

  H  0 ).

Similar to the spectral element solution provided by Al-Khoury et. al. (2001a), a cylindrical
axisymmetric coordinate system shown in Figure 5-1 shall be employed herein. Then, the
following equations are obtained for the potentials by substituting Equation (5-7) into Equation
(5-6):

 2

1   2  
 2



 2
 r 2 r r z 2 
t



  2    2     2    


(5-8)

  2 H 1 H  2 H H 
H 
 2 H
 2


   
    H 






r r
 r 2
r2 
z 2
r2 
t 2




(5-9)

where His the tangential and also the only component of H that does not vanish. By defining
H as:

H  


r

47

(5-10)

it can be shown that the scalar potential  satisfies the following wave equation. The proof can
be obtained immediately if Equation (5-11) is differentiated with respect to r (Ewing et. al.,
1957).

 2


  2    

 r 2




1   2  
 2


 2
r r
z 2 
t


(5-11)

y
θ

x

h1

ur1

r

uz1

h2

ur2
uz2

··
·

ur3
uz3

h(n-1)

ur(n-1)
uz(n-1)

hn

urn
uzn
z
Figure 5-1 Coordinate System for Axisymmetric Layers on a Halfspace
48

Equations (5-8) and (5-11) are the wave equations that govern the axisymmetric wave motion in
a continuous, linear viscoelastic medium. It is also worthwhile to note that if the lamé constants,

 and , were independent of time (i.e., the material is linear elastic), then the convolution
integral in Equations (5-8) and (5-11) reduce to an arithmetic multiplication and these 2
equations become the well known axisymmetric wave equations for a linear elastic material
(Ewing et. al., 1957, Graff, 1991).

Another immediate consequence of adopting the axisymmetric coordinate system is that the
displacement component in the tangential direction, u, vanishes (Graff, 1991). The remaining
deflections can be written as the following in terms of the scalar potentials and :

ur 

  2 

 r rz

(5-12)

uz 

  2  1 


z
r 2 r r

(5-13)

And the stresses can be written as:

 rz   

    2  1   2  





2
2
2 
r  z
r r
r
z 



 2
 

z   

 r 2







1   2  
    2  1  

 2  




r r
z  z
z 2 
r 2 r r 




49

(5-14)

(5-15)

5.2. Solutions For The Wave Equations In The Laplace-Hankel Domain

The solution to the wave equations presented above can be worked more conveniently by
utilizing the integral transforms. Taking the Laplace transform of Equation (5-8) results in:


2
2
    1    


2
s   2  


  s 
 r 2 r r z 2 









(5-16)



where s is the Laplace variable and f s    e  st f (t )dt is the Laplace transform of a function
0

f(t). Then, taking the Hankel Transform (also known as the Fourier-Bessel transform) of order
zero defined as f k   

 
0

f (r ) J 0 kr rdr on Equation (5-16), one finds:





s   2

 2



 k 2   s 2
 z 2




  


(5-17)

After a simple rearrangement of the terms, the above equation can be written as:


s 
  k 2     0
2
z 2 
c1 



2

where,

50

(5-18)



 2   2
c1 

(5-19)



From Equation (5-18), the solution for the Laplace-Hankel transformed potential function,  , is
obtained as the following, after dropping the term that develops an unbounded result, that is the
wave that propagates in the negative z direction (Ewing et. al., 1957, Graff, 1991):

z k2 

  Ae

s
2
c
1  Ae  zf k , s 

(5-20)

where A is an arbitrary constant.

By following the same mathematical steps shown above, Equation (5-11) can be rewritten as:


s 
  k 2     0
2
z 2 
c2 



2

(5-21)

where,


2 
c2 



(5-22)

Again, after dropping the term leading to unbounded results, the solution for the transformed
potential,  , is obtained as:

51

z k2 

  Ce

s
ˆ
c2
2  Ce  zg k , s 

(5-23)

where C is also an arbitrary constant.

In order to make use of the solutions obtained above for the transformed potentials, it is also
necessary to acquire the equations for the displacements and the stresses in the transformed
domain. While taking the Laplace transform on the displacements and the stresses is straight
forward, additional attention is needed in taking the Hankel transform due to the spatial
symmetry supplied by the cylindrical coordinate system adopted for the solution. Referring back
to Figure 1, one finds that the displacement at any point on the z-axis (i.e., when r = 0) is only
allowed to occur in the z-direction (i.e., uz ≠ 0 when r = 0) but is confined in the r-direction (i.e.,
ur = 0 when r = 0), unless the axisymmetric assumption is to be violated. Due to these physical
characteristics of the axisymmetric displacements, Hankel transforms of different orders need to
be applied to ur and uz.

Figure 5-2 shows the first few cycles of the Bessel functions of the first kind, and of orders zero
(J0) and one (J1) that make up the kernels of the Hankel transform. The primary difference
between the two Bessel functions shown in the figure is that while the Bessel function of order
zero (J0) has a nonzero value at r = 0, the Bessel function of order one (J1) is equal to zero
when r = 0. This implies that the Hankel transform of order zero whose kernel is composed of
J0 is appropriate for transforming the functions that exhibit nonzero values at the origin, whereas

52

the the Hankel transform of order one whose kernel is made up of J1 is more appropriate for
transforming the functions that have zero values at r = 0. Therefore, the appropriate Hankel
transformed that should be applied to ur and uz are of orders one and zero, respectively. Taking
the Laplace and the respective Hankel Transforms on the displacements, ur and uz, shown in
Equations (5-12) and (5-13) results in:

u r  k  k

J0 (r) and J1 (r)

uz 

1.2
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6


z

(5-24)


 k 2
z

(5-25)

J0
J1

0

5

10
r

15

Figure 5-2 Bessel Functions of the First Kind

53

20

It should be noted that although the Hankel transform of order one was used to transform ur
shown in Equation (5-12), the Hankel transform of the potentials shown in Equation (5-24) is
still of order zero. This is a consequence of the partial derivative with respect to r that is present
in both terms of the right hand side of Equation (5-12) and the following property of the Hankel
transform which associates the first order transform of a function’s derivative to the zero order
transform of the original function (Graff, 1991, Sneddon, 1995):



 f r 

0

r

J1 krrdr   k 


0


f (r ) J 0 krrdr  k  f k 

(5-26)

Subsequently, the Laplace-Hankel transforms need to be carried out on the relevant stresses.
Due to the same mathematical arguments presented above for the displacements, the Hankel
transform of orders one and zero should be applied respectively on rz and z to allow for a
solution that is compatible with the axisymmetric coordinate system chosen for the solution.
After simplifying, these equations are obtained as:

 
2 s 


2
 
2
2
ˆ
 z
z
c2 


ˆ
 rz   sk 2

 
s 


ˆ
  sk 2
 2k 2  

2 
ˆ
 z
c2 



54

(5-27)

 


 z  s   k 2  




2 
  

ˆ
 k 2
  2s 
z  z

z 2 


2

s
 

ˆ

     2sk 2    

2
z 
c2



(5-28)

5.3. Formulation of the Stiffness Matrices for the Layer Elements

The solutions presented above for the scalar potentials in the transformed domain are not readily
applicable for a multi layered system such as the one shown in Figure 5-1. In order to allow for
the analysis of a layered system such as an asphalt pavement, it is necessary to develop the
formulations for the layer elements of which the underlying concept origins from the method of
Finite Element Analysis (FEA). In this section, two types of layer elements are developed – a 2
noded element for a layer with a finite thickness (e.g., the top layer inFigure 5-1) and a 1 noded
element for simulation of a semi-infinite halfspace (e.g, the bottom layer in Figure 5-1).

5.3.1. Two Noded Element for a Layer with a Finite Thickness

The solutions for the scalar potentials shown in Equations (5-20) and (5-23) only account for the
incident waves that propagate from the upper boundary of a layer in the direction of the positive
z-axis, i.e., downward direction in Figure1.

However, a layer with a finite thickness also

encompasses the waves that reflect from the lower boundary and propagate in the direction of the
negative z-axis. To account for these reflected waves, an additional term must be added to each
of the potentials, which results in the following equations:
55

  Ae  zf  Be  h  z  f

(5-29)

  Ce  zg  De h  z g

(5-30)

where B and D are arbitrary constants and h is the layer thickness. Substituting the above
equations into Equations (5-24) and (5-25) results in the following equations for the
displacements within a two noded element:

ur   Ake zf  Bke h  z  f  Ckge zg  Dkge h  z g

(5-31)

u z   Afe  zf  Bfe  h  z  f  Ck 2e  zg  Dk 2e  h  z g

(5-32)

For the formulation of a layer element, the displacements at the upper and lower boundaries need
to be extracted from the above equations. The radial and the vertical displacements at the upper
boundary, denoted respectively as u r1 and u z1 , can be obtained by substituting z = 0 in the
above equations.

Similarly, the displacements at the lower boundary ( u r 2 and u z 2 ) are

acquired by substituting z = h. In matrix form, the resulting equations for the displacements can
be written as:

 u r1    k
u  
 z1    f
 
 hf
u r 2    ke
u z 2   fe  hf
  

 ke hf
fe  hf
k
f

kg
k2

kge hg

k 2 e  hg

 A
 kge hg   A 
 
B
 
k 2 e  hg   B 
 S1   
 C 
 kg  
C 

D
k 2  D
 
 

56

(5-33)

It is also necessary to obtain the equations for the stresses. By substituting Equations (5-29) and
(5-30) into Equations (5-27) and (5-28), one arrives at the following equations:



 rz  sk  2 Afe  zf  2Bfe  h  z  f  CKe  zg  DKe h  z g






 z  s AKe  zf  BKe  h  z  f  2Ck 2 ge  zg  2Dk 2 ge  h  z g




(5-34)

(5-35)

where,

s
K  2k 2  
2
c2

(5-36)

Again, the stresses at the upper boundary (  rz1 and  z1 ) are obtained by substituting z = 0 in
Equations (5-34) and (5-35), while those at the lower boundary (  rz 2 and  z 2 ) are found by
substituting z = h, all of which can be summarized in a matrix form as follows:

 2f
 rz1 

 
 z1 
  K

  s    hf
2 fe
 rz 2 
  hf
 z2 
 Ke




 2 fe  hf
Ke  hf
2f
K

K
 2k 2 g

 Ke  hg

 2k 2 ge  hg

 A
 Ke  hg   A 
 
B

 
2k 2 ge  hg   B 
 s  S 2   
 C 
K
C 
 
D
2k 2 g   D 
 
 

(5-37)

Combining Equations (5-33) and (5-37) by eliminating the vector of arbitrary constants, the
stresses can be expressed in terms of the displacements as:

57

 rz1 
 u r1 
 


 z1 

1  u z1 

  s  S 2  S1  

 rz 2 
u r 2 
 z2 
u z 2 





(5-38)

where S1 and S2 are the 4 by 4 matrices defined in Equations (5-33) and (5-37), respectively.
According to the concepts of FEA, the stiffness matrix of an element defines the relationship
between the displacement vector and the boundary traction vector. Owing to the Cauchy stress
principle, the boundary tractions are obtained by taking the dot product between the stress tensor
and a unit vector directed along the outward normal of the boundary. Calculating these tractions
at the upper and the lower boundaries of the element and reorganizing them in a vector form
results in the following relationship between the tractions, stresses, and displacements:

Tr1    rz1 
 u r1 
  

u 
Tz1     z1 
 z1 
 
  S 2  noded  

Tr 2    rz 2 
u r 2 
Tz 2    z 2 
u z 2 



  

(5-39)

From the above equation, it is seen that the 4 by 4 matrix S2-noded is the stiffness matrix of the 2
noded layer element which is calculated as:



S 2 noded  s  N  S 2  S1 1

where

58

(5-40)

 1 0
 0 1
Ν
0 0

0 0

0 0
0 0

1 0

0 1

(5-41)

5.3.2. One Noded Semi-Infinite Element

The axisymmetric one noded element was schematically shown as the bottom layer in Figure 5-1.
As shown in the figure and as its name implies, the one noded element only has a single
boundary at the top of the layer and extends infinitely in all other directions. As a consequence,
the waves in this element are only allowed to propagate away from the upper boundary (which is
also the only boundary) without any waves reflecting back. Therefore, the solutions for the
scalar potentials shown in Equations (5-20) and (5-23) can be used without any modifications.
Substituting these two equations into Equations (5-24) and (5-25) results in the following for the
displacements:

ur   Ake zf  Ckge zg

(5-42)

u z   Afe  zf  Ck 2e  zg

(5-43)

The displacements at the boundary are obtained by substituting z = 0 in the above equations and
can be written as the following in matrix form:

59

u r1    k
 
u z1   f

kg   A
 A
2    S 3   
k  C 
C 

(5-44)

Again, the equations for the shear and normal stresses are obtained by substituting the potentials
(Equations (5-20) and (5-23)) into Equations (5-27) and (5-28), respectively:







 

(5-45)



(5-46)

 rz  s  2 Akfe  zf  Ck k 2  g 2 e  zg







 z  s A k 2  g 2 e  zf  2Ck 2 ge  zg

Substituting z = 0 in Equations (5-45) and (5-46) results in the stresses at the boundary:

 2kf
 r1 
ˆ
   s   2
2
k g
 z1 










 A
 k k 2  g 2   A
ˆ
    s  S 4  
 2k 2 g  C 
C 


(5-47)

From Equations (5-44) and (5-47), the following relationship is attained between the stresses
and the displacements:

 rz1 

1 u r1 

  s  S 4  S 3   
  z1 
u z1 

(5-48)

where S3 and S4 were defined in Equations (5-44) and (5-47), respectively. By applying the
Cauchy stress principle, the following relationship is achieved between the tractions, stresses,
and displacements:
60

Tr1    rz1 
u r1 
 
  S1 noded   
u z1 
Tz1     z1 

(5-49)

where the stiffness matrix for the one noded element can be written in terms of the previously
defined variables as:



S1noded  s  S 3  S 4 1

(5-50)

5.4. Incorporating Elastic And Viscoelastic Layer Properties

For a homogenous, isotropic, elastic material whose properties are independent of time, the
relationship between the elastic modulus, E, and the lamé constant, , is given by the theory of
linear elasticity as:



E
2(1   )

(5-51)

Because the parameters in the above equation are not functions of time, the Laplace transform of
the above equation is simply obtained as:



 s  


s



61

E
1

2(1   ) s

(5-52)

However, as it was noted by the pioneer of the spectral element method for layered media (Rizzi,
1989), it is advantageous to add a small amount of damping to the lamé constant, , as no
realistic material is purely elastic. Following Rizzi (1989), this artificial damping can be added to
the above lamé constant as:





 s    s   1    s 

(5-53)

where  is a damping constant. To simulate the wave propagation through an elastic layer, the
above simple equation can be substituted into the equations for the layer elements presented
earlier.

To incorporate the viscoelastic material effects into the solution derived in the previous sections,
it is necessary to adopt a simple function that is capable of representing the fundamental property
of a viscoelastic material analytically. In addition, because all of the time-dependent variables
including stresses, displacements, and material properties (i.e., lamé constants) were transformed
into the Laplace domain, it is preferable to choose a function that is easily transformable into the
Laplace domain. Among the analytical functions described in Chapter 3, the generalized power
function has been selected because it still maintains a simple form while representing the
viscoelastic property effectively. Again, this function is expressed as:

Dt   D0  D1  t m

62

(5-54)

where D0, D1, and m are the power function parameters. Taking the Laplace transform on the
above equation results in:


D0 s m  D1(m  1)
Ds  
s m 1

(5-55)

where mis the gamma function. By substituting Equation (5-55) into Equation (3-16) and
rearranging, one obtains the following equation for the uniaxial relaxation modulus in the
Laplace domain:

ˆ
E s  

1
s m 1

ˆ
s 2 Ds  D0 s m  D1(m  1)

(5-56)

As it was mentioned in Chapter 4, the Poisson’s ratio of a viscoelastic material is also time
dependent. However, the Poisson’s ratio has typically been assumed to be a time-independent
constant in past literatures (Huang, 2004).

In addition, Lee and Kim (2009) showed that

assuming a reasonable constant value for the Poisson’s ratio still results in accurate viscoelastic
responses.

The assumption of a constant Poisson’s ratio implies that the time dependent

behavior of a viscoelastic material in shear or bulk is identical to the behavior in uniaxial mode,
and simplifies the solution significantly (Kim et. al., 2010).

With this assumption, the

relationship between the viscoelastic lamé constant and the uniaxial relaxation modulus shown in
Equation (5-56) is found to be the following:

63


E s 
s m 1
 s  

2(1   ) 2  1     D0 s m  D1(m  1)






(5-57)

5.5. Construction Of The Global Stiffness Matrix

After the stiffness matrices have been obtained for all the layers that make up the structure, the
global stiffness matrix may be constructed in the same way as the traditional FEA methods
(Cook et. al., 2001). In-depth explanation on the concept of the FEA as well as the relationship
between the element and the global stiffness matrices is beyond the scope of this paper. Hence,
it will not be explained herein and interested readers are referred to a variety of textbooks
available on this subject. In this paper, only the generic conceptual schematics will be outlined
and the discussion will be kept to a minimal for conciseness of the manuscript. Figure 5-3 shows
the schematics of the global stiffness matrices for the two types of layered structures that are
most widely adopted for modeling a pavement system.

Figure 5-3(a) shows how the global stiffness matrix is constructed for a layered system resting
on a halfspace.

As mentioned, this pavement model is capable of dissipating the energy

geometrically through the one noded halfspace and is generally used for simulating the FWD
time histories that do not show free vibration at the end of the loading. On the other hand, Figure
5-3(b) shows the global stiffness matrix that can be used for a layered system sitting on a stiff
bedrock at shallow depth. In general, the bedrock in such a model is generally assumed to have
an infinite stiffness, not allowing any displacements to occur within the layer. This implies that

64

the wave energy is trapped between the pavement surface and the top of bedrock, which may
induce free vibration in the layered system at the end of the loading.

S1 noded
2
S 2 noded
2

S Global 





S n -1noded
2
n
S1 noded

(a)
S1 noded
2
S 2 noded
2



S Global 



S n -1noded
2

The shaded portion is not used

S n noded
2

(b)
Figure 5-3 Construction of the Global Stiffness Matrix for Structures (a) with and (b) without a
halfspace
65

Upon constructing the global stiffness matrix, the displacements at the system nodes can be
found from the following equation:

U  S Global 1  P

(5-58)

where,

U   r1 U z1  U ri
U

U zi


 U rn U zn 

(5-59)

is a vector of system displacements to be calculated in global coordinates with U ri and U zi
th

being the radial and vertical displacements at the i node from the top, respectively. Similarly,

P  Pr1

Pz1  Pri

Pzi

 Prn


Pzn 

is a nodal force vector in global coordinates, with the radial and vertical forces at the i

(5-60)

th

node

denoted as Pri and Pzi , respectively. The nodal forces in this vector should be obtained from the
boundary conditions as will be presented in the next section.

66

5.6. Boundary Conditions for a Circular Unit Impulse Loading at the Ground Surface

For the problem in hand where the loading is induced by an impact of a falling weight at the
ground surface, all components of P in Equation (5-60) vanish except for Pz1 . In other words,
the only external load applied to the system is in the vertical direction at the top node (node 1).
In this paper, this surface force will also be in the form of a unit impulse load acting over a
circular area, for the reasons to be explained in subsequent sections of the paper. In the physical
time and spatial domain, this boundary condition is mathematically expressed as the following:

Pz1 (r, t )  R(r )   (t )

(5-61)

where (t) is the dirac delta function for the impulse loading and,

1 , 0  r  a
R(r )  
0 , r  a

(5-62)

where a is the radius of the circular loaded area. However, it should be noted that the stiffness
matrices were previously derived in the Laplace-Hankel domain rather than the physical domain.
As such, it is also necessary to convert the above boundary condition into the one in the
transformed domain. Because the Laplace transform of (t) is equal to 1, taking the LaplaceHankel transforms on Equation (5-61) simply results in the following equation:

Pz1 (k , s) 

67

a
J1 ka
k

(5-63)

5.7. Inversion Of Laplace And Hankel Transforms

As mentioned, the displacements at all nodes of the system can be obtained through Equation (558) from the global stiffness matrix and the force boundary condition described in the previous
sections. It is noted again that the displacements obtained in this manner are in the LaplaceHankel domain and need to be inverse transformed back to the physical domain. However, it has
been shown that even for an elastic halfspace (which simply has a single boundary) subjected to
a point load, the closed form inversion of the Laplace-Hankel transformed displacement is rather
complicated and is close to impossible for a generalized problem (Graff, 1991). Therefore, the
closed form inversion of the displacements obtained from Equation (5-58) is not even attempted
due to the mathematical complexity arising from the viscoelastic material behavior and the wave
propagation phenomenon. Instead, the inversion will be carried out numerically for both the
Laplace and Hankel transforms.

5.8. Numerical Inversion of the Hankel Transform

As it was mentioned earlier, Hankel transforms of orders zero and one were used to transform
the vertical and radial displacements, respectively. Therefore, the inverse Hankel transform of
respective orders must be carried out for the two displacements. In this paper, the numerical
integration scheme will be outlined for the vertical displacement (i.e, the inverse Hankel
transform of order zero). The inverse transform of the radial displacement can also be evaluated
in a similar manner.

68

The closed form equation for the inverse Hankel transform of the vertical displacement at node i
is given as:



U zi (r )   U zi k J 0 krkdk

(5-64)

0

The above integral can also be written as a series of integrals:

b3
bn 1
b2

U zi (r )   U zi k J 0 krkdk   U zi k J 0 krkdk     U zi k J 0 krkdk   (5-65)
b1

b2

bn

Then, each integral in the right hand side of Equation (5-65) needs to be evaluated numerically.
Upon selecting the 6-point Gaussian quadrature as the numerical scheme to be used, the integral
in the above equation can be evaluated as (Abramowitz and Stegun, 1972):

bn 1

b
b 6
U zi k J 0 krkdk  n 1 n  w pU zi  p J 0  p  p

2
p 1

   

(5-66)

bn

where

b 
 bn 1  bn 
b
 x p   n 1 n 
2
2





p 

and xp and wp are the Gaussian nodes and their corresponding weights, respectively.
69

(5-67)

The parameter bn defines the limits of each integration which can be chosen arbitrarily.
However, Cornille (1972) indicated that the convergence of the Gaussian quadrature is greatly
improved if the limits are selected to be the successive roots of the derivative of the Bessel
function that comprise the kernel of the inverse transform. Based on a sensitivity analysis
conducted by the author, subdividing the region between the successive roots of the Bessel
function of order one (that is, the derivative of the Bessel function of order zero) into ten smaller
regions of equal intervals provided satisfactory results for the numerical integration.

It is also noted that the upper bound of the integral shown in Equation (5-64) is equal to infinity.
This indicates that the summation of integrals shown in Equation (5-65) should also span over an
infinite range. However, as it was indicated by Kim (2011), the numerical integration converges
very rapidly even after the first few cycles of the Bessel function comprising the kernel of the
inverse Hankel transform. Therefore, in his static solution for a viscoelastic layered system, the
first five cycles of the Bessel function were used to invert the Hankel transform near the loaded
area, and less number of cycles in the region far from the loading (Kim, 2011). The developers
of the axisymmetric spectral element method used the Fourier-Bessel series (which is the
discrete version of the Hankel transform) in their solution and the summation was also carried
out approximately for the first five cycles of the Bessel function (Al-Khoury et. al., 2001a).
Although the details will be omitted for the compactness of this paper, the sensitivity analysis
performed for the proposed algorithm also showed that the numerical integration over the first
five cycles of the Bessel function is adequate for the solution.

70

5.9. Numerical Inversion of the Laplace Transform

For the inverse Laplace transform, a multi precision numerical scheme known as the Fixed
Talbot Algorithm is adopted in this paper due to its efficiency, accuracy, and ease for
implementation (Abate and Valko, 2004).

The Bromwich integral which is the standard

equation for the inverse Laplace transform is given as:

U zi t  

1
ˆ
e tsU zi s  ds
2j 

(5-68)

B

where j   1 . The contour, B, chosen for the above integral is along the following path:

s    cot   j ,      

(5-69)

where  is a fixed value calculated as:



2M
5t

(5-70)

In Equation (5-70), M is the number of precision decimal digits to be used for the numerical
analysis. For the sake of accuracy, this value is specified to be equal to the machine precision.
By replacing the contour path in Equation (5-68) with the one shown in Equation (5-69), one
finds:

71









U zi t    Re e ts U zi s 1  j   d


(5-71)

0

where,

      cot   1 cot 

(5-72)

Finally, the inverse Laplace transform is obtained by approximating the integral shown in
Equation (5-71) through the trapezoidal rule:

U zi t  

M 1

ts  

U zi  et   Re e q U zi s  q 1  j  q  



M 2


q 1



 1 


  

 

(5-73)

where,

q 

q
M

(5-74)

5.10. System Response To Arbitrary Loading

As it was described in Equations (5-61) and (5-62), the boundary condition considered in the
previous sections was for a unit impulse load distributed over a circular area. As such, the
vertical displacement, U zi , obtained from Equation (5-73) represents the unit impulse response
72

of the layered system in time domain. The primary advantage of the time domain unit impulse
response is that the system response to any arbitrary loading can be obtained through the
convolution integral (Santamarina and Fratta, 1998, Bendat and Piersol, 2010). Theoretically,
this convolution integral for a continuous function is given as:

t

y zi t   U zi t   T t    U zi t   T  d

(5-75)

0

where T(t) could be any arbitrary time dependent loading function and yzi(t) is the corresponding
vertical displacement at node i. For a discrete signal such a FWD time history, the above
equation needs to be evaluated numerically as (Santamarina and Fratta, 1998, Bendat and Piersol,
2010):

y zi t n  

tn

U zi t n  t p T t p t

t p 1

where t is time interval of the discrete signal and tn = nt for an integer n.

73

(5-76)

CHAPTER 6 - IMPLEMENTATION AND VALIDATION OF ALGORITHM –
VISCOWAVE

The solution presented in the preceeding chapter has been implemented into a computer
algorithm named ViscoWave (for Viscoelastic Wave analysis). The algorithm was used to
simulate the behavior of elastic and viscoelastic structures subjected to a FWD loading. In
addition, other available solutions were also used to simulate the response of the same pavement
structures for validation of the ViscoWave algorithm. The results of these numerical simulations
are presented in this chapter.

6.1. Simulation of an Elastic Structure using ViscoWave and LAMDA

The properties of the pavement layers used for the elastic analysis are shown in Table 6-1. The
FWD loading was idealized to be a half-sine load distributed over a circular area of radius 6 in., a
peak magnitude of 9.0 kips, and a duration of 26 ms. The surface deflections were calculated at
radial distances of 0 in., 8 in., 12 in., 18 in., 24 in., 36 in., and 60 in.from the center of the
loading plate.

Table 6-1 Layer Properties for Elastic Simulation of LAMDA and ViscoWave
Layer
Asphalt
Base
Subgrade

Elastic Modulus
(ksi)
145
30
15

Poisson’s Ratio
0.35
0.40
0.45

74

Thickness
(inch)
6
10
∞

Mass Density
(pcf)
4.503
3.882
3.106

In order to verify the results from ViscoWave, the elastic simulation was also conducted using
the axisymmetric spectral element algorithm, namely LAMDA, which has already been verified
through a comparison with 3-dimensional FEA solution (Al-Khoury et. al., 2001a, Lee, 2011).

The time histories for the resulting surface deflections are shown in Figure 6-1. The figure
indicates that ViscoWave and LAMDA showed almost identical results, which validates the
algorithm behind ViscoWave.

25

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

20
15
10
5

Deflection (mils)

Deflection (mils)

25

15
10
5

0
-5
0

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 n.
r = 36 in.
r = 60 in.

20

0

0.01

0.02 0.03
Time (sec)

0.04

-5
0

0.05

0.01

(a)

0.02 0.03
Time (sec)

0.04

0.05

(b)

Figure 6-1 Surface Deflections of a Layered Elastic Structure from (a) ViscoWave and (b)
LAMDA

6.2. Simulation of Viscoelastic Structures using ViscoWave

The viscoelastic simulation was carried out using the same pavement structure that was used for
the previous elastic simulation (Table 1) with a couple of exceptions. The viscoelasticity of the
asphalt concrete was modeled using two different creep compliance functions: one that
75

represents a low temperature behavior (Figure 6-2a) and the other representing a high
temperature behavior in which the viscoelastic effects are more pronounced (Figure 6-2b). In
addition, for each of the creep compliance functions shown in Figure 6-2, the subgrade layer was
first modeled to be a halfspace (infinite thickness) and then with a shallow bedrock (infinite

1.2E-06

Creep Compliance (1/psi)

Creep Compliance (1/psi)

sitffness) located 9.6 ft. below the pavement surface.

1.0E-06
8.0E-07
6.0E-07
4.0E-07

2.0E-07
0.0E+00
0

0.2

0.4
0.6
Time (sec)

0.8

3.0E-05
2.5E-05
2.0E-05
1.5E-05
1.0E-05

5.0E-06
0.0E+00
0

(a)

0.2

0.4
0.6
Time (sec)

0.8

(b)

Figure 6-2 (a) Low Temperature and (b) High Temperature Asphalt Creep Compliance Curves
Used for ViscoWave Simulation

To verify the results of the viscoelastic simulation from ViscoWave, a commercially available
FEA package, ADINA, was used to simulate the dynamic response of the viscoelastic pavement
subjected to the FWD loading. Figure 6-3 shows the geometry and the FEA mesh that was used
for the analysis. Although the elements in ViscoWave assume that the elements extend to infinity
in the horizontal direction and also in the vertical direction for the one noded element, the FEA
simulation was inevitably conducted with a finite geometry. More specifically, the FEA model
only extended to 20 ft. in the horizontal direction and 41.3 ft. in the vertical direction for the
76

simulation of the halfspace (Figure 6-3a). The FEA mesh was generated in such a way that finer
meshes were used near the loaded area and coarser meshes were used near the geometric
boundaries. A total of approximately 8,600 axisymmetric elements, each consisting of 9 nodes,
were consistently used for all FEA simulations.

Axisymmetric
Asphalt: 6 inch
Base: 10 inch
Subgrade:
480 inch for halfspace
simulation and 100 inch
for bedrock simulation

20 ft.

(a)
Axisymmetric
6 inch

(b)
Figure 6-3 (a) Axisymmetric Finite Element Geometry and (b) Finite Element Mesh Used for
Simulation of Pavement Response Under FWD Loading
77

The simulated deflection time histories generated for the low temperature asphalt pavement with
a halfspace and with a shallow bedrock are shown in Figure 6-4 and Figure 6-5, respectively.
Both figures indicate that the results from ViscoWave and ADINA are in excellent agreement.
In Figure 6-4, the effect of phase characteristics in the viscoelastic material can be clearly
identified as shown by the increased duration (or delayed recovery) of the deflections when
compared to the elastic simulation results shown in Figure 6-1.

However, the deflections

calculated from ADINA at 60 in. from the center of the load plate showed some unexpected
outcome towards the end of the time history, especially for the case with a shallow bedrock
(Figure 6-5b). It is believed that this is primarily due to the horizontal boundary at which the
stress wave is reflected back towards the loaded area and secondarily due to the increased error
coming from the increased size of the mesh used for the FEA analysis (Lee, 2013).

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

Deflection (mils)

10
8
6
4
2

12

0
-2
0

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

10

Deflection (mils)

12

8
6
4
2
0

0.01

0.02 0.03
Time (sec)

0.04

-2
0

0.05

(a)

0.01

0.02 0.03
Time (sec)

0.04

0.05

(b)

Figure 6-4 Surface Deflections of a Layered Viscoelastic Structure with a Halfspace at Low
Temperature Simulated Using (a) ViscoWave and (b) ADINA

78

Deflection (mils)

10
8
6
4
2
0

10

-2
-4
0

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

12

Deflection (mils)

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

12

8
6
4
2
0
-2

0.02

0.04
Time (sec)

-4
0

0.06

(a)

0.02

0.04
Time (sec)

0.06

(b)

Figure 6-5 Surface Deflections of a Layered Viscoelastic Structure with a Bedrock at 3 m below
Surface at Low Temperature Simulated Using (a) ViscoWave and (b) ADINA

The simulated deflection time histories for the high temperature asphalt pavement are shown
inFigure 6-6 for the case with a halfspace and in Figure 6-7 for the case with a shallow bedrock.
As expected, the viscoelastic effect is even more pronounced when compared to the low
temperature results, as evidenced in Figure 6-6 by the delayed recovery and in Figure 6-7 by the
noticeable phase difference seen in the free vibration response. Although it is not as significant
as the one shown in Figure 6-5b, the deflection simulated for the outmost sensor (60 in. from
center of load plate) shown in Figure 6-7b again shows the effect of the horizontal boundary
placed at 20 ft. Nonetheless, both Figure 6-6 and Figure 6-7 again show that the ViscoWave
results are in excellent agreement with the FEA results. This concludes the validation of the
viscoelastic simulation using ViscoWave.

79

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

0.01

0.02 0.03
Time (sec)

0.04

Deflection (mils)

Deflection (mils)

16
14
12
10
8
6
4
2
0
-2
0

0.05

16
14
12
10
8
6
4
2
0
-2
0

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

0.01

(a)

0.02 0.03
Time (sec)

0.04

0.05

(b)

Figure 6-6 Surface Deflections of a Layered Viscoelastic Structure with a Halfspace at High
Temperature Simulated Using (a) ViscoWave and (b) ADINA

Deflection (mils)

15
10
5
0
-5
0

0.02

0.04
Time (sec)

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in
r = 36 in.
r = 60 in.

15

Deflection (mils)

r = 0 in.
r = 8 in.
r = 12 in.
r = 18 in.
r = 24 in.
r = 36 in.
r = 60 in.

10
5
0
-5
0

0.06

(a)

0.02

0.04
Time (sec)

0.06

(b)

Figure 6-7 Surface Deflections of a Layered Viscoelastic Structure with a Bedrock at 3 m below
Surface at High Temperature Simulated Using (a) ViscoWave and (b) ADINA

80

CHAPTER 7 - THEORETICAL BACKCALCULATION USING VISCOWAVE

In this chapter, the ViscoWave algorithm developed in the previous chapters will be used along
with available non-linear optimization algorithms for backcalculation of layer properties from
theoretically generated FWD time histories.

7.1. Review of Selected Optimization Algorithms

The non-linear optimization algorithms that were selected for backcalculation were the GaussNewton and Levenberg-Marquardt methods. As these are classical optimization methods and the
detailed description of these algorithms are available in most optimization textbooks, only a brief
overview will be provided in this report (Scales 1985, Rao 2002, Venkataraman 2002).

7.1.1. Objective Function for Optimization and its Derivatives

If the optimization problem includes n independent scalar variables denoted as x1 through xn, it is
more convenient to gather all the scalar variables together and represent them in a vector form.
That is:

x  x1

x2  xn T

Similarly, the scalar valued functions of x can also be gathered into a vector form:

81

(7-1)

f x   f1 x

f 2 x 

f m xT

(7-2)

In backcalculation, the individual functions, fi, in the above equation is frequently taken as the
scalar difference between the measured and simulated deflections, whereas the xi is referred to
the pavement layer parameter that needs to be backcalculated. Then, the objective of the
backcalculation would be to minimize the error between the measured and simulated deflections.
Mathematically, this can be achieved by finding the vector x that minimizes the following
objective function which is a simple sum of squares of the error:

m

F x    f i2 x   f T x   f x 

(7-3)

i 1

It should be noted that in practice, the differences between the predicted and the measured values
are more frequently presented in terms of the Root Mean Square Error (RMSE) which has the
same units as the measured values. The objective function shown above is related to the RMSE
in the following manner:

RMSE 

F
m

(7-4)

For the purpose of optimization, taking the first partial derivative on Equation (7-3) with respect
to xj results in a gradient vector of F which can be written as:

82

m
f
F
gj 
 2 f i i
x j
x j
i 1

(7-5)

or, equivalently in vector form as the following:

gx  F  2J T xf x

(7-6)

where the Jacobian matrix, J, is defined as the following.

 f1
 x
 1
J 
 f m
 x
 1

f1 
xn 

 
f m 

xn 



(7-7)

Differentiating Equation (7-5) or Equation (7-6) with respect to xk results in the Hessian matrix
of F. That is:

Gkj 

m 
 2 fi 
 f f i

 2  i
 fi

xk
xk x j 
i 1  k
 x x j


g j

(7-8)

which can also be written in matrix form as:

Gx  g  2J T xJx  2Sx

83

(7-9)

where the matrix S is given as:

m

Sx    f i x Ti x 

(7-10)

Ti x   2 f i x

(7-11)

i 1

and Ti is the Hessian matrix of fi:

7.1.2. Basics of Iteration

In order for the objective function given in Equation (7-3) to be minimized, it is evident that its
corresponding gradient vector in Equation (7-6) must be a vector of zeros. Denoting the
th

independent variable, objective function value, and the gradient after the p iteration as xp, Fp,
and gp, respectively, the goal of the following iteration is to find xp+1 that satisfies the zerogradient condition. Mathematically, this can be stated as the following:



 



g p 1  g x p 1  g x p  r p  0

(7-12)

where rp is the search vector (or direction vector). Expanding the above expression by Taylor
series and eliminating the higher order terms results in the following expression:

G p r p  g p
84

(7-13)

After the search vector has been found from the above equation, the vector of independent
variable, x, is updated in the following manner and the iteration is continued until a satisfactory
point has been reached.

x p 1  x p  r p

(7-14)

7.1.3. Gauss-Newton Method

By substituting Equation (7-6) and Equation (7-9) into Equation (7-13), the following
relationship is obtained for the search vector:

J pT J p  S p  rp  J pT f p

(7-15)

The above equation, along with Equation (7-14), defines the Newton’s method. However, the
major problem with the above formulation is that the computational effort needed for calculating
the Sp is extremely expensive. In the Gauss-Newton method, the computationally expensive
matrix Sp is completely neglected and the search vector is found from the following relationship.

J p T J p r p  J p T f p

85

(7-16)

7.1.4. Levenberg-Marquardt Method

Although the Gauss-Newton method is known to be very efficient, the problem that is frequently
encountered is that the matrix J pT J p being singular or ill-conditioned (Scales 1985). The
Levenberg-Marquardt method incorporates a technique for handling the issue of J pT J p being
singular by adding an identity matrix to J pT J p . In other words, the search vector is found from
a modified version of Equation (7-16) which can be written as:

J pT J p   pI p rp  J pT f p

(7-17)

It should be noted that as  p  0 , the above equation becomes identical to the Gauss-Newton
method. On the other extreme, as  p   , the matrix J pT J p becomes negligible and the
search vector r p becomes an infinitesimal step in the direction of steepest descent (Scales 1985).
However, the actual magnitude of  p should be found such that the objective function, F, is
minimized at x p  r p .

7.1.5. Remarks on the Optimization Routine

Both the Gauss-Newton and Levenberg-Marquardt optimization routines used for this study were
readily implemented in the “lsqnonlin” function of Matlab’s Optimization Toolbox (Mathworks
2008). The backcalculation was initiated using the Gauss-Newton method which is known to be
86

more efficient. However, when the problem of J pT J p matrix being singular or ill-conditioned
for inversion, the Gauss-Newton was terminated and the backcalculation was continued using the
Levenberg-Marquardt method.

7.2. Theoretical Backcalculation Using ViscoWave

The theoretical backcalculation was carried out in two folds: (1) using a single set of FWD time
histories representing a single temperature and multiple sets of FWD time histories at various
temperatures. The results and findings of the theoretical backcalculation exercise are provided in
the remainder of the chapter.

7.2.1. Single Temperature Backcalculation

7.2.1.1. Reference Pavement Structures and Seed Values for Backcalculation

The properties of the reference pavement structures sitting on a halfspace and on a bedrock are
summarized in Table 7-1. The generalized power function was used for modeling the viscoelastic
behavior of the asphalt surface with its coefficients shown in the table. The FWD loading was
again idealized to be a half-sine load with a peak magnitude of 9.0 kips and a duration of 26 ms,
distributed over a circular area of radius 6 in. The time histories were modeled with a discrete
time interval of 0.2 ms and a duration of 50 ms for pavement with a halfspace and 70 ms for the
pavement on a bedrock. The surface deflections were calculated at radial distances of 0 in., 8 in.,
12 in., 18 in., 24 in., 36 in., and 60 in. from the center of the loading plate. Table 7-2 shows the

87

three sets of seed values that were assigned to the variables prior to backcalculation. The first set
of seed values set to ±50 percent increase from the true values. In order to study the effect of
seed values, the second and the third sets of seed values were set to 100 percent and -80 percent
increase from the true values, respectively, for all variables. These seed values were consistently
used for both the pavement structures with and without a halfspace.

Table 7-1 Reference Pavement Structures
Layer

D0
(1/psi)

D1
(1/psi)

m

Asphalt 3.00E-07 5.00E-07
0.3
Base
N/A
Subgrade
Note 1: For pavement with a halfspace
Note 2: For pavement with a bedrock

Elastic
Modulus
(psi)
N/A
30,000
15,000

Poisson’s
Ratio

Density
(pcf)

Thickness
(in.)

0.35
0.4
0.45

145
125
100

6
12
1
2
∞ /60

Table 7-2 Seed Values Used for Backcalculation
Variable

Seed Value (Percent Increase from True Value)
Seed Set 1
Seed Set 2
Seed Set 3

D0 (1/psi)

4.50E-07 (50%)

6.00E-07 (100%)

6.00E-08 (-80%)

D1 (1/psi)
m

2.50E-07 (-50%)

1.00E-06 (100%)

1.00E-07 (-80%)

0.45 (50%)

0.60 (100%)

0.06 (-80%)

EBase (psi)

15,000 (-50%)

60,000 (100%)

6,000 (-80%)

ESubgrade (psi)

22,500 (50%)

30,000 (100%)

3,000 (-80%)

88

7.2.1.2. Backcalculation Results for Pavement Structure With a Halfspace

The deflection time histories generated from ViscoWave for the reference pavement with a
halfspace (Table 7-1) are shown in Figure 7-1a, whereas those generated for the pavement
structures with the three sets of seed values (Table 7-2) are shown in Figure 7-1b through Figure
7-1d. As shown by the deflection time histories in Figure 7-1, the first set of seed values resulted
in the deflection time histories that are very similar to those from the reference pavement. The
second set of seed values resulted in deflection time histories that are quite different from those
of the reference pavement in terms of magnitude but the shape of the time histories still resemble
those of the reference pavement. The last set of seed values resulted in the poorest deflection
time histories both in terms of shape and magnitude. Based on these observations, it is expected
that this set of seed values has the greatest potential for converging to the true values with the
fewest number of iterations, followed by the second and third set of seeds.

The backcalculation was conducted using the Gauss-Newton and Levenberg-Marquardt
optimization routines previously described in this chapter. The backcalculation was initiated
using the Gauss-Newton method and was switched over to the Levenberg-Marquardt method
after 2 iterations for seed sets 1 and 2, and after the first iteration for seed set 3. The three
solutions converged after 24 (seed set 1), 60 (seed set 2), and 42 (seed set 3) iterations. The
RMSE histories and the deflection time histories obtained after backcalculation are shown in
Figure 7-2 and Figure 7-3, respectively. As expected, Figure 7-2a and Figure 7-2b show that the
first set of seed values converged faster than the second set, and the backcalculated deflection
time histories shown in Figure 7-3b and Figure 7-3c indicate that the backcalculated solution

89

may be reasonably close to the true values. On the other hand, the terminal RMSE shown in
Figure 7-2c indicates immediately that the backcalculation with the third set of seed values
converged to a wrong solution and this is confirmed by the time histories shown in Figure 7-3d.

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

9
7
5
3

1
-1 0

9
7
5
3
1
-1 0

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(a)

7
5
3
1
-1 0

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

11

Deflection (mils)

Deflection (mils)

9

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(b)

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

11

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

11

Deflection (mils)

Deflection (mils)

11

9
7
5
3

1

-1 0 0.01 0.02 0.03 0.04 0.05
Time (Sec)

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(c)

(d)

Figure 7-1 Deflection Time Histories from Pavement Structure on a Halfspace; (a) Reference
Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3

90

RMSE (mils)

RMSE (mils)

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
0

5
10
15
20
Number of Iterations

25

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0

10 20 30 40 50
Number of Iterations

(b)

RMSE (mils)

(a)

60

3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

10
20
30
Number of Iterations

40

(c)
Figure 7-2 History of Root Mean Square Error for (a) Seed Set 1, (b) Seed Set 2, and (c) Seed
Set 3

91

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

9
7
5
3

1
-1 0

9
7
5
3
1
-1 0

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(a)

7
5
3

1
-1 0

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

11

Deflection (mils)

Deflection (mils)

9

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(b)

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

11

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

11
Deflection (mils)

Deflection (mils)

11

9
7
5
3

1

-1 0 0.01 0.02 0.03 0.04 0.05
Time (Sec)

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(c)

(d)

Figure 7-3 Backcalculated Deflection Time Histories from Pavement Structure on a Halfspace;
(a) Reference Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3

92

The backcalculated results are summarized in Table 7-3, for all sets of seed values. Again, the
table confirms that the first set of seed values resulted in the best outcome in terms of the
backcalculated values and their percent errors. It also confirms that the third set of seed values
converged to a wrong solution, with the backcalculated D1 parameter being negative and the
backcalculated base and subgrade moduli being very close to the initial seed values.

Table 7-3 Backcalculated Values and Percent Error – Pavement on a Halfspace
Variable

True Value

D0 (1/psi)

3.00E-07

D1 (1/psi)

5.00E-07

m

0.30

EBase (psi)

30,000

ESubgrade (psi)

15,000

Backcalculated Value (Absolute Percent Error)
Seed Set 1
Seed Set 2
Seed Set 3
3.00E-07
3.14E-07
2.14E-07
(<0.001%)
(4.7%)
(28.8%)
5.00E-07
5.27E-07
-2.24E-07
(<0.001%)
(5.5%)
(144.8%)
0.30
0.34
0.05
(0.002%)
(12.9%)
(84.7%)
30,000
30,229
5,997
(<0.001%)
(0.8%)
(80.0%)
15,000
14,988
3,003
(<0.001%)
(0.1%)
(80.0%)

The above observations confirm the well-accepted fact that both the Gauss-Newton and
Levenberg-Marquardt methods may only find the local minimum and not necessarily the global
minimum (Scales, 1985). In other words, with the Gauss-Newton and Levenberg-Marquardt
methods, it is important that reasonable seed values are chosen in order for the backcalculation to
converge to the correct solution. Figure 7-4 through Figure 7-8 show the convergence histories
of the backcalculated parameters for all three cases and confirm the observations already made in
the above.

93

5.0E-07
4.0E-07

Backcalculated
Value
True Value

6.0E-07

D0 (1/psi)

4.5E-07

D0 (1/psi)

7.0E-07

Backcalculated
Value
True Value

3.5E-07
3.0E-07

5.0E-07
4.0E-07
3.0E-07

2.5E-07
2.0E-07

2.0E-07
0

5 10 15 20 25
Number of Iterations

0

(a)

10 20 30 40 50 60
Number of Iterations
(b)

5.0E-07

Backcalculated
Value
True Value

D0 (1/psi)

4.0E-07
3.0E-07
2.0E-07
1.0E-07
0.0E+00
0

10 20 30 40
Number of Iterations
(c)

Figure 7-4 Convergence of Power Function Parameters D0; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Halfspace

94

1.2E-06

Backcalculated
Value
True Value

5.2E-06

D1 (1/psi)

1.0E-06

D1 (1/psi)

6.2E-06

Backcalculated
Value
True Value

8.0E-07
6.0E-07
4.0E-07

4.2E-06
3.2E-06

2.2E-06
1.2E-06

2.0E-07

2.0E-07
0

5 10 15 20 25
Number of Iterations

0

(a)

10 20 30 40 50 60
Number of Iterations
(b)

9.0E-07

Backcalculated
Value
True Value

D1 (1/psi)

7.0E-07
5.0E-07
3.0E-07

1.0E-07
-1.0E-07
-3.0E-07
0

10
20
30
40
Number of Iterations
(c)

Figure 7-5 Convergence of Power Function Parameters D1; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Halfspace

95

Backcalculated
Value
True Value

m

m

1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0

5
10
15
20
Number of Iterations

1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

Backcalculated
Value
True Value

0

25

(a)

10 20 30 40 50
Number of Iterations

60

(b)

0.6

Backcalculated
Value
True Value

0.5

m

0.4
0.3
0.2
0.1
0.0
0

10
20
30
Number of Iterations

40

(c)
Figure 7-6 Convergence of Power Function Parameters m; (a) Seed Set 1, (b) Seed Set 2, and (c)
Seed Set 3 – Pavement Structure on a Halfspace

96

130,000

Backcalculated
Value
True Value

Backcalculated
Value

110,000

EBase (psi)

EBase (psi)

50,000
45,000
40,000
35,000
30,000
25,000
20,000
15,000
10,000

True Value

90,000
70,000
50,000
30,000
10,000

0

5
10 15 20 25
Number of Iterations

0

(a)

10 20 30 40 50 60
Number of Iterations
(b)

50,000

Backcalculated
Value

EBase (psi)

40,000

True Value

30,000
20,000
10,000
0
0

10
20
30
40
Number of Iterations
(c)

Figure 7-7 Convergence of Base Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set 3 –
Pavement Structure on a Halfspace

97

ESubgrade (psi)

ESubgrade (psi)

35,000

Backcalculated
Value
True Value

24,000
22,000
20,000
18,000
16,000
14,000
12,000
10,000

Backcalculated
Value
True Value

30,000
25,000
20,000
15,000
10,000

0

5
10 15 20 25
Number of Iterations

0

10 20 30 40 50 60
Number of Iterations

(a)

(b)

ESubgrade (psi)

25,000

Backcalculated
Value
True Value

20,000
15,000
10,000
5,000
0
0

10
20
30
40
Number of Iterations
(c)

Figure 7-8 Convergence of Subgrade Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set 3
– Pavement Structure on a Halfspace

98

7.2.1.3. Backcalculation Results for Pavement Structure on a Bedrock

The deflection time histories generated for the reference pavement and the seed structures on a
bedrock are shown in Figure 7-9. As was the case with the pavement on a halfspace presented
above, the first set of seed values resulted in the deflection time histories that resemble the true
histories both in shape and magnitude. The second set of seed values resulted in deflection time
histories that are similar in shape but reduced in magnitude when compared to the true ones. The
last set of seed values resulted in the poorest deflection time histories both in terms of shape and
magnitude.

The backcalculation was initiated using the Gauss-Newton method. As shown by the RMSE
histories in Figure 7-10, the first, second, and third sets of seed values converged after 17, 42,
and 39 iterations, respectively. The respective backcalculation was switched over to the
Levenberg-Marquardt method after 13, 1, and 2 iterations. The deflection time histories obtained
after backcalculation are shown in Figure 7-11. Both Figure 7-10 and Figure 7-11 indicate that
the first set of seed values were most successful followed by the second set, and that the third set
converged to a wrong solution.

99

0

0.02

0.04 0.06
Time (Sec)

Deflection (mils)

Deflection (mils)

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

12
10
8
6
4
2
0
-2
-4

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

12
10
8
6
4
2
0
-2
-4

0.08

0

0.02

12
10
8
6
4
2
0
-2
-4

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0

0.02

0.04 0.06
Time (Sec)

0.08

(b)

Deflection (mils)

Deflection (mils)

(a)

0.04 0.06
Time (Sec)

0.08

12
10
8
6
4
2
0
-2
-4

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0

(c)

0.02

0.04 0.06
Time (Sec)

0.08

(d)

Figure 7-9 Deflection Time Histories from Pavement Structure on a Bedrock; (a) Reference
Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3

100

1.2

0.4

1.0

RMSE (mils)

RMSE (mils)

0.5

0.3
0.2
0.1

0.8

0.6
0.4
0.2
0.0

0.0
0

5
10
15
Number of Iterations

0

20

(b)

RMSE (mils)

(a)

10
20
30
40
Number of Iterations

4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

10
20
30
Number of Iterations

40

(c)
Figure 7-10 History of Root Mean Square Error for (a) Seed Set 1, (b) Seed Set 2, and (c) Seed
Set 3

101

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0

0.02

0.04 0.06
Time (Sec)

Deflection (mils)

Deflection (mils)

12
10
8
6
4
2
0
-2
-4

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

12
10
8
6
4
2
0
-2
-4
0

0.08

0.02

12
10
8
6
4
2
0
-2
-4

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0

0.02

0.04 0.06
Time (Sec)

0.08

(b)

Deflection (mils)

Deflection (mils)

(a)

0.04 0.06
Time (Sec)

0.08

12
10
8
6
4
2
0
-2
-4

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0

(c)

0.02

0.04 0.06
Time (Sec)

0.08

(d)

Figure 7-11 Backcalculated Deflection Time Histories from Pavement Structure on a Bedrock;
(a) Reference Pavement, (b) Seed Set 1, (c) Seed Set 2, and (d) Seed Set 3

102

The backcalculated results for all sets of seed values are summarized in Table 7-4 which shows
that the first and second sets of seeds were able to backcalculate the true values with reasonable
errors. The third set of seed values converged to a wrong solution, as clearly seen by the negative
power function parameters and the minimal improvement in the unbound layer modulus from the
seed values. Again, these observations confirm that the Gauss-Newton method that was used for
backcalculation is sensitive to the seed values. The convergence histories shown in Figure 7-12
through Figure 7-24 confirm the above observations.

Table 7-4 Backcalculated Values and Percent Error – Pavement on a Bedrock
Variable

True Value

D0 (1/psi)

3.00E-07

D1 (1/psi)

5.00E-07

m

0.3

EBase (psi)

30,000

ESubgrade (psi)

15,000

Backcalculated Value (Absolute Percent Error)
Seed Set 1
Seed Set 2
Seed Set 3
2.96E-07
3.01E-07
5.93E-07
(1.2%)
(0.4%)
(-97.6%)
4.99E-07
5.01E-07
-5.11E-07
(0.2%)
(0.2%)
(202.3%)
0.29
0.30
-0.01
(2.5%)
(0.9%)
(102.9%)
29,936
30,049
5,989
(0.2%)
(0.2%)
(80.0%)
15,007
14,998
3,000
(0.1%)
(0.01%)
(80.0%)

103

Backcalculated
Value
True Value

D0 (1/psi)

4.5E-07
4.0E-07
3.5E-07
3.0E-07

3.0E-06

Backcalculated
Value
True Value

2.5E-06
D0 (1/psi)

5.0E-07

2.0E-06
1.5E-06
1.0E-06

2.5E-07

5.0E-07

2.0E-07

0.0E+00
0

0

5
10
15
20
Number of Iterations
(a)

10 20 30 40
Number of Iterations
(b)

1.0E-06

Backcalculated
Value
True Value

D 0 (1/psi)

8.0E-07

6.0E-07
4.0E-07
2.0E-07
0.0E+00
0

10
20
30
40
Number of Iterations
(c)

Figure 7-12 Convergence of Power Function Parameters D0; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Bedrock

104

1.2E-06

Backcalculated
Value
True Value

D1 (1/psi)

D1 (1/psi)

1.0E-06

Backcalculated
Value
True Value

4.0E-06

8.0E-07
6.0E-07

3.0E-06
2.0E-06
1.0E-06

4.0E-07
2.0E-07

0.0E+00

0

5
10
15
20
Number of Iterations

0

10 20 30 40
Number of Iterations

(a)

(b)

Backcalculated
Value
True Value

D1 (1/psi)

9.0E-07
4.0E-07
-1.0E-07
-6.0E-07
0

10
20
30
40
Number of Iterations
(c)

Figure 7-13 Convergence of Power Function Parameters D1; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Bedrock

105

0.5

1.0

Backcalculated
Value
True Value

0.4

Backcalculated
Value
True Value

0.8

0.3

m

m

0.6

0.2

0.4

0.1

0.2

0.0

0.0
0

5
10
15
Number of Iterations

20

0

10
20
30
40
Number of Iterations

(a)

(b)

0.5

Backcalculated
Value
True Value

0.4

m

0.3
0.2
0.1
0.0
-0.1
0

10
20
30
Number of Iterations

40

(c)
Figure 7-14 Convergence of Power Function Parameters m; (a) Seed Set 1, (b) Seed Set 2, and
(c) Seed Set 3 – Pavement Structure on a Bedrock

106

Backcalculated
Value
True Value

0

EBase (psi)

EBase (psi)

50,000
45,000
40,000
35,000
30,000
25,000
20,000
15,000
10,000

5
10
15
20
Number of Iterations

80,000
70,000
60,000
50,000
40,000
30,000
20,000
10,000
0

Backcalculated
Value
True Value

0

(a)

10
20
30
40
Number of Iterations
(b)

50,000

Backcalculated
Value

EBase (psi)

40,000

True Value

30,000
20,000
10,000
0
0

10
20
30
40
Number of Iterations
(c)

Figure 7-15 Convergence of Base Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set 3 –
Pavement Structure on a Bedrock

107

ESubgrade (psi)

ESubgrade (psi)

50,000

Backcalculated
Value
True Value

24,000
22,000
20,000
18,000
16,000
14,000
12,000
10,000

Backcalculated
Value
True Value

40,000
30,000
20,000
10,000
0

0

5
10
15
20
Number of Iterations

0

10
20
30
40
Number of Iterations

(a)

(b)

ESubgrade (psi)

30,000

Backcalculated
Value
True Value

25,000
20,000
15,000
10,000
5,000

0
0

10
20
30
40
Number of Iterations
(c)

Figure 7-16 Convergence of Subgrade Modulus; (a) Seed Set 1, (b) Seed Set 2, and (c) Seed Set
3 – Pavement Structure on a Bedrock

108

7.2.2. Multi-Temperature Backcalculation

7.2.2.1. Theoretical Development of Generalized Power Function Relationship at Multiple
Temperatures

The purpose of the theoretical backcalculation at multiple temperatures is to study if the dynamic
modulus mastercurve can be obtained for a broader range of time (or frequency) and more
accurately from the FWD time histories. As such, it is important to understand how the
generalized power function used in ViscoWave behaves theoretically under the time-temperature
superposition principle. The theoretical development is provided in this section of the report in a
brief manner.

Rewriting Equation (3-27) given for the generalized power function in terms of the reduced time,

, results in:

D   D0  D1   m

(7-18)

Upon substituting the Equation (3-23) into above, the following is obtained.

 t
D   D0  D1  
a
 T

where,

109






m

 D0  D1T  t m

(7-19)

D1T T  

D1

aT T m

(7-20)

The above equations reveal that D1 is the only parameter affected by the change in material
temperature and the other parameters D0 and m remain constant irrespective of temperature. In
addition, the above equations also indicate that if the creep compliances of a viscoelastic material
at multiple temperatures are given in terms of generalized power functions, the shift factors can
easily be obtained through Equation (7-20), since the parameter m does not vary with
temperature.

7.2.2.2. Reference Pavement Structures and Seed Values for Backcalculation

For the multiple temperature backcalculation exercise, the FWD time histories were first
obtained for a couple of pavement structures (with and without halfspace) exposed to four
different temperatures: 0⁰C, 10⁰C, 20⁰C, and 40⁰C. The D1 values for the viscoelastic asphalt
material corresponding to these temperatures are summarized in Table 7-5 along with the other
properties used for generating the time histories. Figure 7-17 and Figure 7-18 show the shift
factors that were used to obtain the D1 values and the creep compliance curves at the above
mentioned temperatures, respectively. The properties of the base and subgrade layers with and
without a halfspace are summarized in Table 7-6. The pavement structures were imposed to a
half-sine load with the same characteristics provided in the previous section for the single
temperature backcalculation. The surface deflection time histories were again obtained at radial
distances of 0 in., 8 in., 12 in., 18 in., 24 in., 36 in., and 60 in. from the center of the loading plate.
110

Table 7-5 Asphalt Properties Used for Reference Pavement Structures
D1 (1/psi)

m

Poisson’s
Ratio

Density
(pcf)

Thickness
(in.)

1.01E-07

7.52E-06
1.74E-05
3.82E-05
1.57E-04

0.27

0.35

145

6

2.0
Log(Shift Factor)

0
10
20
40

D0 (1/psi)

y = 0.0004x2 - 0.1407x + 1.367

1.0
0.0
-1.0

-2.0
-3.0
-4.0
0

10
20
30
Temperature (⁰C)

40

Figure 7-17 Time-Temperature Shift Factors

Creep Compliance (1/psi)

Temp(⁰C)

6.0E-04
5.0E-04
0C
10C
20C
40C

4.0E-04
3.0E-04
2.0E-04

1.0E-04
0.0E+00
0

20

40 60
Time (sec)

80

100

Figure 7-18 Creep Compliance Curves at Multiple Temperatures

111

Table 7-6 Base and Subgrade Properties Used for Reference Pavement Structures
Layer

Elastic Modulus
(psi)

Poisson’s Ratio

Density (pcf)

Thickness (in.)

0.4
0.45

125
100

12
2
∞ /60

Base
30,000
Subgrade
15,000
Note 1: For pavement with a halfspace
Note 2: For pavement with a bedrock

1

Table 7-7 summarizes the seed values used for backcalculation. Although the D1 parameter was
dependent on temperature, a single seed value was assigned regardless of temperature.

Table 7-7 Seed Values Used for Backcalculation
Variable

Temperature (⁰C)

Seed Value

Percent Increase from
True Value

D0 (1/psi)

N/A

1.52E-06

50.0%

m
EBase (psi)
ESubgrade (psi)

15.9%

0
10
20
40
N/A

D1 (1/psi)

8.71E-06

-50.0%
-77.2%
-94.5%

0.4

50.0%

N/A

15,000

-50.0%

N/A

22,500

50.0%

7.2.2.3. Backcalculation Results for Pavement Structure With a Halfspace

Figure 7-19 shows the deflection time histories of the reference pavement for all temperatures.
The seed deflection history that was used regardless of temperature is shown in Figure 7-20. It is

112

also noted that the seed deflection time histories resemble the shape of the true time histories

40
35
30
25
20
15
10
5
0
-5 0

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

Deflection (mils)

Deflection (mils)

reasonably well.

0.01 0.02 0.03 0.04 0.05
Time (Sec)

40
35
30
25
20
15
10
5
0
-5 0

40
35
30
25
20
15
10
5
0
-5 0

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(b)

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

Deflection (mils)

Deflection (mils)

(a)

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.01 0.02 0.03 0.04 0.05
Time (Sec)
(c)

40
35
30
25
20
15
10
5
0
-5 0

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0.01 0.02 0.03 0.04 0.05
Time (Sec)
(d)

Figure 7-19 Simulated FWD Deflection Time Histories for (a) 0⁰C, (b) 10⁰C, (c) 20⁰C, and (d)
40⁰C

113

40

Sensor
Offset

35

0 in.

Deflection (mils)

30

8 in.

25

12 in.

20

18 in.

15

24 in.
36 in.

10

60 in.

5
0

-5

0

0.01

0.02
0.03
Time (Sec)

0.04

0.05

Figure 7-20 Seed FWD Deflection Time History for All Temperatures

The backcalculation was terminated after 18 iterations for which the Gauss-Newton method was
used throughout. Figure 7-21 shows the history of RMSE and Table 7-8 summarizes the
backcalculated values and their respective errors. The table shows that the backcalculation was
successful, with the percent errors for all variables being less than 0.001 percent expect for the m
value which showed an error of 1.1 percent.

114

2.5

RMSE (mils)

2.0
1.5

1.0
0.5

0.0
0

5

10
15
Number of Iterations

20

Figure 7-21 History of Root Mean Square Error

Table 7-8 Backcalculated Values and Percent Error – Pavement on a Halfspace
Variable

Temperature
(⁰C)

True Value

Backcalculated
Value

Absolute Percent
Error

D0 (1/psi)

N/A

1.01E-07

1.01E-07

<0.001%

0

7.52E-06

7.52E-06

10

1.74E-05

1.74E-05

20

3.82E-05

3.82E-05

40

1.57E-04

1.57E-04

m

N/A

0.27

0.27

<0.001%
<0.001%
<0.001%
<0.001%
1.1%

EBase (psi)

N/A

30,000

30,000

<0.001%

ESubgrade (psi)

N/A

15,000

15,000

<0.001%

D1 (1/psi)

The convergence history of the viscoelastic parameters are shown in Figure 7-22 while those for
the base and subgrade layers are shown in Figure 7-23.

115

Backcalculated
Value
True Value

1.8E-04
1.3E-04

D1 (1/psi)

D0 (1/psi)

3.5E-06
3.0E-06
2.5E-06
2.0E-06
1.5E-06
1.0E-06
5.0E-07
0.0E+00
-5.0E-07

0C
20C

8.0E-05

10C
40C

3.0E-05
-2.0E-05
0

5
10
15
20
Number of Iterations

0

(a)

5
10
15
20
Number of Iterations
(b)

0.5

Backcalculated
Value
True Value

0.4

m

0.3
0.2
0.1
0.0
0

5
10
15
Number of Iterations

20

(c)
Figure 7-22 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – Pavement
Structure on a Halfspace

116

Backcalculated
Value

ESubgrade (psi)

EBase (psi)

70,000
60,000
50,000
40,000
30,000
20,000
10,000
0

True Value

0

Backcalculated
Value
True Value

24,000
22,000
20,000
18,000
16,000
14,000
12,000
10,000
0

5
10
15
20
Number of Iterations
(a)

5
10
15
20
Number of Iterations
(b)

Figure 7-23 Convergence of (a) Base and (b) Subgrade Layer Moduli – Pavement Structure on a
Halfspace

7.2.2.4. Backcalculation Results for Pavement Structure on a Bedrock

Figure 7-24 shows the deflection time histories using ViscoWave for all temperatures. The seed
deflection history that was used regardless of temperature is shown in Figure 7-25. Note that the
free vibration is present in all deflection time histories shown in Figure 7-24 including those
generated for 40⁰C. However, the deflection histories generated with the seed values did not
show any free vibration. Such discrepancies in the shape of the deflection time histories may
cause the solution some additional iterations to converge.

117

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0.02
0.04
Time (Sec)

Deflection (mils)

Deflection (mils)

40
35
30
25
20
15
10
5
0
-5 0

0.06

40
35
30
25
20
15
10
5
0
-5 0

40
35
30
25
20
15
10
5
0
-5 0

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0.02
0.04
Time (Sec)

0.02
0.04
Time (Sec)

0.06

(b)

Deflection (mils)

Deflection (mils)

(a)

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.06

(c)

40
35
30
25
20
15
10
5
0
-5 0

Sensor
Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.
0.02
0.04
Time (Sec)

0.06

(d)

Figure 7-24 Simulated FWD Deflection Time Histories for (a) 0⁰C, (b) 10⁰C, (c) 20⁰C, and (d)
40⁰C

118

40

Sensor
Offset

35

0 in.

Deflection (mils)

30

8 in.

25

12 in.

20

18 in.

15

24 in.
36 in.

10

60 in.

5
0

-5

0

0.02

0.04
Time (Sec)

0.06

Figure 7-25 Seed FWD Deflection Time History for All Temperatures

The backcalculation converged and was terminated after 26 iterations. As expected, the 7
additional iterations when compared to the case with a halfspace, are believed to be due to the
seed deflection time histories not showing any free vibration (i.e., the shape of the deflection
time histories not resembling the true ones). As was the case with the pavement on a halfspace,
the Gauss-Newton method was used throughout the backcalculation process, without switching
over to the Levenberg-Marquardt method. Figure 7-26 shows the history of RMSE and Table 7-9
summarizes the backcalculated values and their respective errors Again, the percent errors for all
variables were less than 0.001 percent except for the m value. The convergence histories for the
asphalt and the unbound layers are shown in Figure 7-27 and Figure 7-28, respectively.

119

2.5

RMSE (mils)

2.0
1.5

1.0
0.5

0.0
0

5

10
15
20
Number of Iterations

25

30

Figure 7-26 History of Root Mean Square Error

Table 7-9 Backcalculated Values and Percent Error – Pavement on a Halfspace
Variable

Temperature
(⁰C)

True Value

Backcalculated
Value

Absolute Percent
Error

D0 (1/psi)

N/A

1.01E-07

1.01E-07

<0.001%

0

7.52E-06

7.52E-06

10

1.74E-05

1.74E-05

20

3.82E-05

3.82E-05

40

1.57E-04

1.57E-04

m

N/A

0.27

0.27

<0.001%
<0.001%
<0.001%
<0.001%
1.1%

EBase (psi)

N/A

30,000

30,000

<0.001%

ESubgrade (psi)

N/A

15,000

15,000

<0.001%

D1 (1/psi)

120

Backcalculated
Value
True Value

D0 (1/psi)

3.0E-06
2.0E-06
1.0E-06
0.0E+00

1.8E-04
1.3E-04

D1 (1/psi)

4.0E-06

0C
20C

8.0E-05

10C
40C

3.0E-05

-1.0E-06
-2.0E-06

-2.0E-05

0

5 10 15 20 25 30
Number of Iterations

0

10
20
30
Number of Iterations

(a)

(b)

0.5

Backcalculated
Value
True Value

0.4

m

0.3
0.2
0.1
0.0
0

5 10 15 20 25
Number of Iterations

30

(c)
Figure 7-27 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – Pavement
Structure on a Bedrock

121

Backcalculated
Value
ESubgrade (psi)

EBase (psi)

50,000
45,000
40,000
35,000
30,000
25,000
20,000
15,000
10,000

True Value

0

Backcalculated
Value
True Value

24,000
22,000
20,000
18,000
16,000
14,000
12,000
10,000
0

5 10 15 20 25 30
Number of Iterations
(a)

5 10 15 20 25 30
Number of Iterations
(b)

Figure 7-28 Convergence of (a) Base and (b) Subgrade Layer Moduli – Pavement Structure on a
Bedrock

122

CHAPTER 8 - FIELD BACKCALCULATION USING VISCOWAVE

In this chapter, the ViscoWave algorithm and the optimization routines described in previous
chapters will be used for backcalculating the pavement properties using the FWD data collected
in the field and the results will be compared to those obtained in the laboratory.

8.1. Test Site and Protocol

The tested pavement was part of the Accelerated Pavement Test (APT) tracks in the State
Materials Office (SMO) facility of the Florida Department of Transportation (FDOT). The
pavement consisted of a 6.5 in. asphalt surface layer on top of a 10 in. limerock base. The asphalt
mixture included a PG-76-22 polymer modified binder which is typically used in regular FDOT
projects. The FWD testing was conducted in June of 2010 and the pavement cores were sampled
immediately after the FWD tests. According to the measurements made at several monitoring
wells installed around the APT facility, it was found that the ground water table at the time of
FWD testing was located approximately 11.0 ft to 11.5 ft. below the surface.

8.2. Laboratory Dynamic Modulus Test

As was mentioned in Chapter 4 of the report, although the uniaxial tension or compression test is
desirable for testing of the viscoelastic materials, the problem associated with size requirements
(4 in. diameter and 6 in. thickness) frequently makes it impossible, especially on cored mixtures
from thin asphalt layers. As an alternative to the uniaxial test, the dynamic modulus tests were

123

conducted in the IDT mode (as described in Chapter 4), with the adoption of the dynamic
modulus testing procedure used in AASHTO TP 62.

A total of four field cores obtained from the test track and were delivered to the FDOT’s asphalt
laboratory. Each core was carefully cut to yield two IDT specimens (one of each from the top
and the bottom lifts) with a diameter of 6.0 in. and a thickness of 1.5 in. To allow for generating
the viscoelastic mastercurve, the dynamic modulus test was conducted at three distinct
temperatures: 0⁰C, 10⁰C, and 20⁰C. For each temperature, the specimens were placed in the
target temperature for at least 12 hours prior to the IDT tests. For each temperature, the dynamic
modulus tests were then conducted at five frequencies: 0.1 Hz, 0.5 Hz, 1.0 Hz, 5.0 Hz, and 10.0
Hz. The resulting IDT dynamic modulus data was analyzed as described in Chapter 4 and the
mastercurve was constructed according to the procedures shown in Chapter 3. Figure 8-1 shows
the dynamic modulus mastercurve for both the top and the bottom lift at a reference temperature
of 10⁰C while the corresponding shift factors are shown in Figure 8-2.

124

4,500

4,000
3,500

Top

Bottom

|E*| (ksi)

3,000

2,500
2,000

1,500
1,000
500

0
1.0E-06

1.0E-03 1.0E+00 1.0E+03
Frequency (Hz)

1.0E+06

Figure 8-1 Laboratory Dynamic Modulus Mastercurves

Top

1.5

Log(Shift Factor)

2.0

Bottom

1.0

0.5
0.0
-0.5
-1.0
-1.5
-2.0
0

5

10
15
Temperature (⁰C)

20

Figure 8-2 Laboratory Time-Temperature Shift Factors

125

25

8.3. FWD Backcalculation Using ViscoWave

The load and deflection time histories measured from the FWD is shown in Figure 8-3. Similar
to the deflection time histories shown in Figure 7-19d, it is noted that due to the increased
temperature of the asphalt at the time of FWD testing (58⁰C), the deflection measured at the
center of the load plate did not return to zero.

10,000

Deflection (mils)

Load (lbs)

8,000
6,000
4,000
2,000
0
-2,000
0

0.02
0.04
Time (Sec)

0.06

14
12
10
8
6
4
2
0
-2 0

Sensor Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.01 0.02 0.03 0.04 0.05 0.06
Time (Sec)

(a)

(b)

Figure 8-3 Measured FWD (a) Load and (b) Deflection Time Histories

The backcalculation was carried out in two folds: (1) with a 3 layer pavement system and (2)
with a 5 layer pavement system to incorporate a ground water table at 11.3 ft. below the
pavement surface. Table 8-1 shows the Poisson’s ratio, material density and the thickness that
was assumed for the backcalculation. The seed values for the generalized power function and the
elastic moduli for the unbound layers are summarized in Table 8-2.

126

Table 8-1 Assumed Pavement Properties for Backcalculation
Layer

Poisson’s Ratio

Density (pcf)

Asphalt
0.35
Base
0.35
Subgrade
0.40
1
Emb. 1
0.45
1
Emb. 2
0.45
Note 1: Only applicable to the 5 layer system.

145
125
115
110
110

Thickness (in.)
3 Layer System
5 Layer System
6.5
10
∞
12
108
N/A
∞
N/A

Table 8-2 Seed Values Used for Backcalculation
Seed Value

Variable

3 Layer System

5 Layer System

D0 (1/psi)

5.00E-06

D1 (1/psi)
m

2.00E-05

EBase (psi)

80,000

0.60
32,000

ESubgrade (psi)
1

EEmb-1 (psi)

1

EEmb-2 (psi)
Note 1: Only applicable to the 5 layer system.

N/A

30,000
30,000

8.3.1. Backcalculation Results from the 3 Layer System

The backcalculation using the 3 layer system converged after 31 iterations with a terminal RMSE
value of 0.23 mils. The history of the RMSE is shown in Figure 8-4 and the backcalculated
values are summarized in Table 8-3.

127

0.8
0.7

RMSE (mils)

0.6
0.5
0.4
0.3
0.2

0.1
0.0
0

10
20
Number of Iterations

30

Figure 8-4 History of Root Mean Square Error for the 3 Layer System

Table 8-3 Backcalculated Values from 3 Layer System
Variable

Backcalculated Value

D0 (1/psi)

1.90E-06

D1 (1/psi)
m

1.02E-04

EBase (psi)

32,732

ESubgrade (psi)

62,225

0.64

The time histories using the seed values as well as the backcalculated values are shown in Figure
8-5 which clearly shows that the backcalculated time histories are more comparable to the
measured histories (Figure 8-3b) than the seed histories. Figure 8-6 and Figure 8-7 show the
convergence histories for the viscoelastic parameters of the asphalt and the unbound layer
moduli, respectively.
128

Deflection (mils)

14
12
10
8
6
4
2
0
-2 0

Sensor Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.01 0.02 0.03 0.04 0.05 0.06
Time (Sec)

Deflection (mils)

(a)

14
12
10
8
6
4
2
0
-2 0

Sensor Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.01 0.02 0.03 0.04 0.05 0.06
Time (Sec)
(b)

Figure 8-5 Simulated Time Histories for the 3 layer System with (a) Seed and (b) Backcalculated
Values

129

1.0E-04

4.0E-06

8.0E-05

D1 (1/psi)

1.2E-04

5.0E-06

D0 (1/psi)

6.0E-06

3.0E-06
2.0E-06

6.0E-05
4.0E-05

1.0E-06

2.0E-05

0.0E+00

0.0E+00
0

10
20
30
Number of Iterations

0

(a)

10
20
30
Number of Iterations
(b)

0.7
0.7

m

0.6
0.6
0.5
0.5
0.4
0

10
20
30
Number of Iterations
(c)

Figure 8-6 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – 3 Layer
System

130

100,000

ESubgrade (psi)

EBase (psi)

80,000
60,000
40,000
20,000
0
0

70,000
60,000
50,000
40,000
30,000
20,000
10,000
0

10
20
30
Number of Iterations

0

(a)

10
20
30
Number of Iterations
(b)

Figure 8-7 Convergence of Unbound Pavement Layer Moduli for the 3 Layer System (a) Base
and (b) Subgrade

Figure 8-8 shows the seed and backcalculated creep compliance functions graphically. As can be
seen from the figure, the backcalculated function shows the increased creep compliance which is
most likely responsible for the viscoelastic response as shown by the delayed recovery of the
deflection time histories in Figure 8-5b.

In order to compare the backcalculated viscoelastic parameters to those obtained from the
laboratory, Prony series was fitted to the backcalculated creep compliance function and then
converted to the dynamic modulus values following the procedures outlined in Chapter 3. Figure
8-9a shows the dynamic modulus mastercurve at 58⁰C with the laboratory mastercurve shifted to
this temperature by extrapolating the shift factors shown in Figure 8-2. Figure 8-9b shows the
same dynamic modulus mastercurve but at a reference temperature of 10⁰C with the
backcalculated mastercurve shifted using the extrapolated laboratory shift factors. Both Figure

131

8-9a and Figure 8-9b show that the backcalculated dynamic modulus master curve is in good
agreement with those obtained from the laboratory.

Creep Compliance (1/psi)

2.5E-03

Seed

2.0E-03

Backcalculated

1.5E-03
1.0E-03
5.0E-04

0.0E+00
0

20

40
60
Time (sec)

80

100

Figure 8-8 (a) Seed and (b) Backcalculated Creep Compliance Functions from 3 Layer System

|E*| (ksi)

4000

Top
Bottom
Backcalculated

2000
5000
4000
1500
|E*| (ksi)
|E*| (ksi)

5000

3000
2000
1000
0
1.0E-06 1.0E-02 1.0E+02 1.0E+06
1.0E-03.0E+00
1
1.0E+03
Frequency (Hz)

Top
Bottom
Backcalculated

3000
1000
2000
500
1000
0
0
1.0E-06 1.0E-02 1.0E+02 1.0E+06
1.0E-06
1.0E-03.0E+00
1
1.0E+03
1.0E+06
Frequency (Hz)
Frequency (Hz)

(a)

(b)

Figure 8-9 Backcalculated Dynamic Modulus Mastercurves at Reference Temperatures (a) 58⁰C
and (b) 10⁰C – 3 Layer System
132

8.3.2. Backcalculation Results from the 5 Layer System

The backcalculation using the 5 layer system converged after 20 iterations with a terminal RMSE
value of 0.18 mils. Compared to the 3 layer backcalculation results, this means that the 5 layer
backcalculation resulted in a decrease of 0.05 mils in RMSE even with 11 less iterations. The
history of the RMSE is shown in Figure 8-10 and the backcalculated values are summarized in
Table 8-4.

0.8
0.7

RMSE (mils)

0.6
0.5

0.4
0.3
0.2
0.1
0.0
0

5

10
15
Number of Iterations

20

Figure 8-10 History of Root Mean Square Error for the 5 Layer System

133

Table 8-4 Backcalculated Values from 5 Layer System
Variable

Backcalculated Value

D0 (1/psi)

2.34E-06

D1 (1/psi)
m

1.20E-04

EBase (psi)

29,199

ESubgrade (psi)

204,522

EEmbankment_1 (psi)

44,390

EEmbankment_2 (psi)

155,565

0.69

The time histories using the seed and the backcalculated values are shown in Figure 8-11. Again
the backcalculated time histories clearly show the delayed recovery of the pavement system
which is also seen in the measured time histories (Figure 8-3b). Figure 8-12 and Figure 8-13
show the convergence histories for the viscoelastic parameters of the asphalt and the unbound
layer moduli, respectively.

134

Deflection (mils)

14
12
10
8
6
4
2
0
-2 0

Sensor Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.01 0.02 0.03 0.04 0.05 0.06
Time (Sec)

Deflection (mils)

(a)

14
12
10
8
6
4
2
0
-2 0

Sensor Offset
0 in.
8 in.
12 in.
18 in.
24 in.
36 in.
60 in.

0.01 0.02 0.03 0.04 0.05 0.06
Time (Sec)
(b)

Figure 8-11 Simulated Time Histories for the 5 layer System with (a) Seed and (b)
Backcalculated Values

135

6.0E-06
4.0E-06

D1 (1/psi)

D0 (1/psi)

5.0E-06
3.0E-06

2.0E-06
1.0E-06
0.0E+00
0

1.4E-04
1.2E-04
1.0E-04
8.0E-05
6.0E-05
4.0E-05
2.0E-05
0.0E+00

5
10
15
20
Number of Iterations

0

(b)

m

(a)

5
10
15
20
Number of Iterations

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
0

5
10
15
Number of Iterations

20

(c)
Figure 8-12 Convergence of Power Function Parameters (a) D0, (b) D1, and (c) m – 5 Layer
System

136

200,000

ESubgrade (psi)

250,000

80,000

EBase (psi)

100,000

60,000
40,000
20,000
0

150,000
100,000
50,000
0

0

5
10
15
20
Number of Iterations

0

(a)

(b)

200,000

EEmbankment_2 (psi)

50,000

EEmbankment_1 (psi)

5
10
15
20
Number of Iterations

40,000
30,000
20,000
10,000

150,000
100,000
50,000
0

0
0

0

5
10
15
20
Number of Iterations
(c)

5
10
15
20
Number of Iterations
(d)

Figure 8-13 Convergence of Unbound Pavement Layer Moduli for the 3 Layer System (a) Base,
(b) Subgrade, (c) Embankment 1, and (d) Embankment 2

137

Figure 8-14 shows the creep compliance functions plotted using the seed and the backcalculated
parameters. It is noted that although the backcalculated creep compliance shows increased
viscoelastic behaviour when compared to the seed creep compliance, the increased viscoelastic
creep behaviour was slightly higher than the one backcalculated using the 3 layer system.

The backcalculated creep compliance was again converted to the dynamic modulus through the
use of Prony series. The dynamic modulus mastercurves at 58⁰C and 10⁰C are shown in Figure
8-15. Although the mastercurve plotted for 58⁰C (Figure 8-15a) may seem fairly reasonable, the
mastercurve at 10⁰C (Figure 8-15b) clearly shows that the backcalculated dynamic modulus
underestimates the laboratory dynamic modulus.

Although the backcalculation with a 5 layer system resulted in a lower RMSE than the 3 layer
system, the 3 layer backcalculation showed better agreement with the laboratory results for the
asphalt modulus.

138

Creep Compliance (1/psi)

3.0E-03

Seed

2.5E-03

Backcalculated

2.0E-03

1.5E-03
1.0E-03
5.0E-04

0.0E+00
0

20

40
60
Time (sec)

80

100

Figure 8-14 (a) Seed and (b) Backcalculated Creep Compliance Functions from 5 Layer System

|E*| (ksi)
|E*| (ksi)

1500
4000

5000
5000

Top
Bottom
Backcalculated

4000
4000

Top
Bottom
Backcalculated

|E*| (ksi)
|E*| (ksi)

2000
5000

3000
3000

3000
1000

2000
2000

2000
500
1000
0
0
1.0E-06 1.0E-02 1.0E+02 1.0E+06
1.0E-06
1.0E-03.0E+00
1
1.0E+03
1.0E+06
Frequency (Hz)
Frequency (Hz)

1000
1000
0
0
1.0E-06 1.0E-02 1.0E+02 1.0E+06
1.0E-06
1.0E-03.0E+00
1
1.0E+03
1.0E+06
Frequency (Hz)
Frequency (Hz)

(a)

(b)

Figure 8-15 Backcalculated Dynamic Modulus Mastercurves at Reference Temperatures (a)
58⁰C and (b) 10⁰C – 5 Layer System

139

CHAPTER 9 - CONCLUSIONS AND RECOMMENDATIONS

Due to the viscoelastic nature of asphalt materials and the dynamic nature of pavement structures,
it is important to consider both effects simultaneously in modeling of asphalt pavements. In this
study, a new computational algorithm, namely ViscoWave, has been developed and implemented
for modeling the pavement dynamics and viscoelasticity under an impact load generated by a
Falling Weight Deflectometer (FWD). The primary advantage of the proposed solution over
some of the existing solutions is that it uses continuous integral transforms (Laplace and Hankel
transforms) that are more appropriate for the FWD time histories whose signal characteristics are
transient, nonperiodic, and truncated.

The theoretical development of ViscoWave follows similar steps to those used for the
development of the spectral element method but in place of the discrete transforms adopted in
the spectral element method, ViscoWave utilizes the continuous integral transforms (namely
Laplace and Hankel transforms) that are more appropriate for transient, nonperiodic signals. A
comparison of the ViscoWave results with other existing solutions such as the Finite Element
Analysis (FEA) and spectral element method indicated that ViscoWave is capable of simulating
the viscoelastic and dynamic effects of asphalt pavements.

Preliminary backcalculation efforts were conducted by adopting the Gauss-Newton and the
Levenberg-Marquardt methods as the optimization routines. Both the theoretical and field
backcalculation indicate that the new solution (ViscoWave) has great potential for
backcalculation, which has remained a major challenge in pavement engineering. Nonetheless, it

140

is recognized that the field data presented in this report only include the FWD data from a single
temperature.

Additional effort is needed to look into the FWD data collected at multiple

temperatures.

It is also recommended that future work be conducted to improve the dynamic backcalculation
using ViscoWave with different optimization routines that may show improved efficiency. It is
also noted that the backcalculation routines used for this study fall into the category of
unconfined optimization, meaning that the backcalculated variables were allowed to take on any
value (even negative values). Future work should include looking into the use of confined
optimization routines which may increase the efficiency and reliability of the backcalculation.

141

BIBLIOGRAPHY

142

BIBLIOGRAPHY
1. Abate, J., and Valko, P. P. 2004. Multi-Precision Laplace Transform Inversion. International
Journal for Numerical Methods in Engineering. 60: 979-993.
2. Abramowitz, M., and Stegun, I. A. 1972. Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables. Dover Publications, Inc. NY.
3. Al-Khoury, R., Scarpas, A., Kasbergen, C., and Blaauwendraad, J. 2001a. Spectral Element
Technique for Efficient Parameter Identification of Layered Media. I. Forward Calculation.
International Journal of Solids and Structures. 38: 1605-1623.
4. Al-Khoury, Kasbergen, C., R., Scarpas, A., and Blaauwendraad, J. 2001b. Spectral Element
Technique for Efficient Parameter Identification of Layered Media Part II. Inverse
Calculation. International Journal of Solids and Structures. 38: 8753-8772.
5. Al-Khoury, R., Scarpas, A., Kasbergen, C., and Blaauwendraad, J. 2002. Spectral Element
Technique for Efficient Parameter Identification of Layered Media Part III. Viscoelastic
Aspects. International Journal of Solids and Structures. 39: 2189-2201.
6. ARA., ERES Division. 2004. Guide for Mechanistic-Empirical Design of New and
Rehabilitated Pavement Structures. National Cooperative Highway Research Program1-37A,
Transportation Research Board, National Research Council.
7. Bendat, J. S., and Piersol, A. G. 2010. Random Data Analysis and Measurement Procedures.
4th Ed. John Wiley & Sons, Inc., NJ.
8. Buttlar, W. G., and Roque, R. 1994. Development and evaluation of the strategic highway
research program measurement and analysis system for indirect tensile testing at low
temperatures. Transportation Research Record, 1454: 163-171.
9. Chatti, K., and Yun, K. K. 1996. SAPSI-M: Computer Program for Analyzing Asphalt
Concrete Pavements Under Moving Arbitrary Loads. Transportation Research Record, 1539:
88-95.
10. Chatti, K., Ji, Y., and Harichandran, R. 2004. Dynamic Time Domain Backcalculation of
Layer Moduli, Damping and Thicknesses in Flexible Pavements. Transportation Research
Record, 1869: 106-116.
11. Christensen, R. M. 2003. Theory of Viscoelasticity. 2nd Ed. Dover Publications, Inc. NY.
12. Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J. 2001. Concepts and Applications
of Finite Element Analysis. 4th Ed. John Wiley & Sons, Inc. NY.
13. Cornille, P. 1972. Computation of Hankel Transforms. SIAM Review. 14(2): 278-285.
143

14. Ewing, W. M., Jardetzky, W. S., and Press, F. 1957. Elastic Waves in Layered Media.
McGraw-Hill Book Company. NY.
15. Findley, W. N., Lai, J. S., and Onaran, K. 1976. Creep and Relaxation of Nonlinear
Viscoelastic Materials. Dover Publications, Inc. NY.
16. Graff, K. 1991. Wave Motion in Elastic Solids. Dover Publications, Inc. NY.
17. Grenier, S., and Konrad, J. 2007. Backcalculation of Asphalt Concrete layer Thickness Using
Static and Dynamic Interpretations of FWD Tests. Transportation Research Board Annual
Meeting CD-ROM, 07-0573, Transportation Research Board, Washington D. C.
18. Hadidi, R., and Gucunski, N. 2007. A Probabilistic Approach to Falling Weight
Deflectometer Backcalculation. Transportation Research Board Annual Meeting CD-ROM,
07-3187, Transportation Research Board, Washington D. C.
19. Hondros, G. 1959. The evaluation of Poisson’s ratio and the modulus of materials of a low
tensile resistance by the brazilian (indirect tensile) test with particular reference to concrete.
Australian Journal of Applied Science, 10(3): 243-268.
20. Huang, Y. H. 2004. Pavement Analysis and Design, 2nd Ed. Pearson Education, Inc. NJ.
21. Irwin, L.H. 1977. Determination of pavement layer moduli from surface deflection data for
pavement performance evaluation. Proc. Fourth Intl. Conf. on Structural Design of Asphalt
Pavements. Ann Arbor, MI, University of Michigan, 831-840.
22. Ji, Y. 2005. Frequency and time domain backcalculation of flexible pavement layer
parameters. Ph. D. Dissertation. Michigan State University.
23. Kang, Y. V. 1990. Effect of Finite Width on Dynamic Deflections of Pavement. Ph. D.
Dissertation. University of Texas, Austin.
24. Kim, J., Roque, R., and Birgisson, B. 2005. Obtaining creep compliance parameters
accurately from static or cyclic creep tests., ASTM International, ASTM STP 1469: 177-197.
25. Kim, J., Sholar, G., and Kim, S. 2008. ‘Determination of accurate creep compliance and
relaxation modulus at a single temperature for viscoelastic solids. Journal of Materials in
Civil Engineering, ASCE, 20(2): 147-156.
26. Kim, J., Lee, H. S., and Kim, N. 2010. Determination of Shear and Bulk Moduli of
Viscoelastic Solids from the Indirect Tension Creep Test. Journal of Engineering Mehcanics.
ASCE. 136:1067-1076.
27. Kim, J. 2011. General Viscoelastic Solutions for the Multilayered Systems Subjected to
Static and Moving Loads. Journal of Materials in Civil Engineering. ASCE. 23:1007-1016.
144

28. Kim, Y. R. 2009. Modeling of Asphalt Concrete. ASCE Press. McGraw Hill.
29. Kutay, E., Chatti, K., and Lei, L. 2011. Backcalculation of Dynamic Modulus (|E*|)
Mastercurve from FWD Surface Deflections. Transportation Research Board Annual
Meeting DVD, 11-3434, Transportation Research Board, Washington D. C.
30. Lee, H. S., and Kim, J. 2009, Determination of viscoelastic Poisson’s ratio and creep
compliance from the indirect tension test. Journal of Materials in Civil Engineering ASCE,
21(8): 416-425.
31. Lee, H. S. 2011, ViscoWave – A New Forward Solution for Dynamic Backcalculation. A
Presentation Made at the 20th Annual FWD Users Group Meeting.
http://pms.nevadadot.com/2011.asp. Last Accessed January 2012.
32. Lee, H. S. 2013. Viscowave – A New Solution for Viscoelastic Wave Propagation of
Layered Structures Subjected to an Impact Load, International Journal of Pavement
Engineering, DOI:10.1080/10298436.2013.782401.
33. Losa, M. 2002. The Influence of Asphalt Pavement Layer Properties on Vibration
Transmission. International Journal of Pavements. 1(1): 67-76.
34. Magnuson, A. H., Lytton, R. L., and Briggs, R. 1991.Comparison of Computer Predictions
and Field Data for Dynamic Analysis of Falling-Weight Deflectometer Data. Transportation
Research Record, 1293: 61-71.
35. Malvern, L. E. 1969. Introduction to the Mechanics of a Continous Medium. Prentice-Hall,
Inc. NJ.
36. Mathworks. 2008. Matlab Optimization Toolbox User’s Guide. The Mathworks Inc., MA.
37. Matsui, K., Maina, J. W., Hachiya, Y., and Ozawa, Y. 2011. Wave Propagation Analysis for
Flexible Pavements and Its Application to Backcalculation Analysis. Transportation
Research Board Annual Meeting DVD, 11-2971, Transportation Research Board,
Washington D. C.
38. Matsui, K., Nishizawa, T., and Kikuta, Y. 1998. Time Domain Backcalculation of Pavements.
Structural Materials Technology III: An NDT Conference. Medlock, R. D., and Laffrey, D.
C. Ed. San Antonio, TX, 3400: 410-419.
39. Ong, C. L., Newcomb, D. E., and Siddharthan, R. 1991. Comparison of Dynamic and Static
Backcalculation Moduli for Three-Layer Pavements. Transportation Research Record, 1293:
86-92.
40. Rao, S. S. 1996. Engineering Optimization: Theory and Practice. 3rd Ed. John Wiley & Sons,
Inc. NY.
145

41. Rizzi, S. 1989. A Spectral Analysis Approach to Wave Propagation in Layered Solids. Ph. D.
Dissertation. Purdue University
42. Roque, R., Buttlar, W. G., Ruth, B. E., Tia, M., Dickison, S. W., and Reid, B. 1997.
Evaluation of SHRP indirect tension tester to mitigate cracking in asphalt concrete
pavements and overlays. Final Report, FDOT B-9885, University of Florida, Gainesville, FL.
43. Roque, R., Buttlar, W. G., Ruth, B. E., and Dickison, S. W. 1998. Short-Loading-Time
Stiffness from Creep, Resilient Modulus, and Strength Tests Using Superpave Indirect
Tension Test. Transportation Research Record, 1630: 10-20.
44. Santamarina, J. C., and Fratta, D. 1998. Introduction to Discrete Signals and Inverse
Problems in Civil Engineering. ASCE Press. VA.
45. Scales, L. E. 1985. Introduction to Non-Linear Optimization. Springer-Verlag, Inc. NY.
46. Schapery, R. A. 1974. Viscoelastic Behavior and Analysis of Composite Materials, in
Mechanics of Composite Materials. Sendeckyj, G. P. Ed. Academic Press. NY.
47. Sneddon, I. N. 1995. Fourier Transforms. Dover Publications, Inc. NY.
48. Tschoegl, N. W., Knauss, W. G., and Emri, I. 2002. Poisson’s ratio in viscoelasticity: a
critical review. Mechanics of Time-Dependent Materials, 6: 3-51.
49. Tschoegl, N. W. 1989. The Phenomenological Theory of Linear Viscoelastic Behavior: An
Introduction. Springer-Verlag, NY.
50. Uzan, J. 1994. Dynamic Linear Backcalculation of Pavement Material Parameters. Journal
of Transportation Engineering ASCE, 120(1): 109-126.
51. Uzan, J., Scullion, T., Parades, M., and Lytton, R. L. 1988. A Microcomputer Based
Procedure for Backcalculating Layer Moduli from FWD Data. Research Report 1123-1,
Texas Transportation Institute.
52. Venkataraman, P. 2002. Applied Optimization with MATLAB Programming. John Wiley &
Sons, Inc. NY.
53. Wineman, A. S., and Rajagopal, K. R. 2000. Mechanical response of polymers: an
introduction. Cambridge University Press, NY.
54. Zhang, W., Drescher, A., and Newcomb, D. E. 1997a. Viscoelastic analysis of diametral
compression of asphalt concrete. Journal of Engineering Mechanics, ASCE, 123(6), 596-603.
55. Zhang, W., Drescher, A., and Newcomb, D. E. 1997b. Viscoelastic behavior of asphalt
concrete in diametral compression. Journal of Transportation Engineering, ASCE, 123(6),
495-502.
146