Plasticity theory for granular media


Material Information

Plasticity theory for granular media
Physical Description:
xxi, 325 leaves : ill. ; 28 cm.
Seereeram, Devo, 1957-
Publication Date:


Subjects / Keywords:
Sand   ( lcsh )
Plasticity   ( lcsh )
Particles   ( lcsh )
Soils -- Testing   ( lcsh )
Civil Engineering thesis Ph. D
Dissertations, Academic -- Civil Engineering -- UF
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1986.
Bibliography: leaves 312-324.
Statement of Responsibility:
by Devo Seereeram.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 030455089
oclc - 17394139
System ID:

This item is only available as the following downloads:

Full Text








I would first like to acknowledge Dr. Frank C. Townsend, the

chairman of my supervisory committee, for his tremendous support and

encouragement during this study. Dr. Townsend was instrumental in

acquiring the hollow cylinder test data and in locating key references.

Professor Daniel C. Drucker made fundamental and frequent

contributions to the basic ideas and their connections and to the

overall intellectual structure of this work. I am deeply grateful for

his vigorous critical readings of early drafts of my dissertation, his

constructive and creative suggestions for revision, and his major

contributions which in many ways influenced the content of this study.

The delight I found in our many discussions is one of my chief rewards

from this project.

I am deeply indebted to the other members of my doctoral committee:

Professors John L. Davidson, Martin A. Eisenberg, William Goldhurst,

Lawrence E. Malvern, and Michael C. McVay for their helpful discussions

and criticism of my work. I would also like to thank Paul Linton for

carefully carrying out the series of in-house triaxial tests.

This acknowledgement would not be complete without mention of the

love and support I received from my family and friends. Particularly, I

would like to thank Charmaine, my fiancee, for her understanding and

patience during the long hours spent at work, and my mother, who put my

education above everything else, and my father, who gave me the

financial freedom and the motivation to seek knowledge.

Finally, I would like to acknowledge the financial support of the

United States Air Force Office of Scientific Research, under Grant No.

AFOSR-84-0108 (M.C. McVay, Principal Investigator), which made this

study possible.



ACKNOWLEDGEMENTS ......................................................i

LIST OF TABLES ....................................................... vii

LIST OF FIGURES. ................ .................................... viii

KEY TO SYMBOLS .......... ..............................................xvi

ABT A T .. .. ......... ... .... .................. ........... .. ...... xx


1 .IITRDODi.i :TI:Ni

1.1 The Role and Nature of Theory ..............................1
1.2 Statement of the Problem*....*................**.............2
1.3 Approach*****......*.. ........*.......................... 4
1.4 Scope* ..... ................... ... ........................ 6


2.1 Introduction ....***...... ... ..... ...................... 10
2.2 Tensors ................................. ..... ..... ...... 11
2.2.1 Background.......***.....*..................... 11
2.2.2 The Stress Tensor ......... .................. 16
2.2.3 The Strain Tensor.............. ............ .. ......28
2.3 Stress-Strain Equations and Constitutive Theory. ....... 33
2.4 A Note on Stress and Strain in Granular Media-........... .38
2.5 Anisotropic Fabric in Granular Material................... 46
2.5.1 Introduction...*........... .................. ...... 46
2.5.2 Common Symmetry Patterns ......................... 47
2.5.3 Fabric Measures......* ..*.. ................. ...... 49
2.6 Elasticity ........ ................... ................. 52
2.6.1 Cauchy Type Elasticity *............................ 53
2.6.2 Hyperelasticity or Green Type Elasticity .......... 57
2.6.3 Hypoelasticity or Incremental Type Elasticity**.....58
2.7 Plasticity ......................... ...................... 61
2.7.1 Yield Surface ...................................... 62
2.7.2 Failure Criteria.................................. ..69
2.7.3 Incremental Plastic Stress-Strain Relation, and
Prager's Theory....*................................76

2.7.4 Drucker's Stability Postulate.......................84
2.7.5 Applicability of the Normality Rule to
Soil Mechanics.*....................................87
2.7.6 Isotropic Hardening.................................91
2.7.7 Anisotropic Hardening...............................99
2.7.8 Incremental Elasto-Plastic Stress-Strain Relation*.102


3.1 Introduction ..............................................106
3.2 Material Behavior Perceived as Most Essential
and Relevant .**........*............................... 113
3.3 Details of the Yield Function And Its Evolution ..........122
3.3.1 Isotropy .................................. ........ 124
3.3.2 Zero Dilation Line ................................126
3.3.3 Consolidation Portion of Yield Surface............ 130
3.3.4 Dilation Portion of Yield Surface*................ 136
3.3.5 Evolutionary Rule for the Yield Surface........... 138
3.4 Choice of the Field of Plastic Moduli.................... 139
3.5 Elastic Characterization ................................142
3.6 Parameter Evaluation Scheme...............................143
3.6.1 Elastic Constants..................................145
3.6.2 Field of Plastic Moduli Parameters ............... 145
3.6.3 Yield Surface or Plastic Flow Parameters ..........147
3.6.4 Interpretation of Model Parameters*...............*148
3.7 Comparison of Measured and Calculated Results Using
the Simple Model.......................................... 148
3.7.1 Simulation of Saada's Hollow Cylinder Tests*.......151
3.7.2 Simulation of Hettler's Triaxial Tests.............168
3.7.3 Simulation of Tatsuoka and Ishihara's
Stress Paths*............ ..........................173
3.8 Modifications to the Simple Theory to Include Hardening...191
3.8.1 Conventional Bounding Surface Adaptation-..........191
3.8.2 Prediction of Cavity Expansion Tests.....*** *......196
3.8.3 Proposed Hardening Modification..........*.........210
3.9 Limitations and Advantages *.. ............................225


4.1 Introduction ............................................230
4.2 Field of Work Hardening Moduli Concept ...................231
4.3 Model Characteristics*...................................237
4.4 Yield Function ..........................................237
4.5 Flow Rule ..................................... ............ 238
4.6 Hardening Rule .......................................... 240
4.7 Initialization of Model Parameters....................... 247
4.8 Verification ....................................... 253

5 CONCLUSIONS AND RECOMMENDATIONS ................................267





OF YIELD SURFACE ..............................................275
PARAMETER B ....................................................301

LIST OF REFERENCES ..*................................ .............312

BIOGRAPHICAL SKETCH ................................................ 325



3.1 Comparison of the Characteristic State and Critical
State Concepts.............................................. 131

3.2 Simple Interpretation of Model Constants....................149

3.3 Expected Trends in the Magnitude of Key Parameters With
Relative Density........................................... 150

3.4 Model Constants for Reid-Bedford Sand at 75% Relative
Density.........* ............................................154

3.5 Computed Isotropic Strength Constants for Saada's Series
of Hollow Cylinder Tests ...................................167

3.6 Model Parameters for Karlsruhe Sand and Dutch Dune Sand.....172

3.7 Model Parameters for Loose Fuji River Sand ................. 186

3.8 Summary of Pressuremeter Tests in Dense Reid-Bedford Sand*..198

3.9 Model Constants Used to Simulate Pressuremeter Tests*.......199

4.1 Prevost Model Parameters for Reid-Bedford Sand..*...........254

5.1 Typical Variation of the Magnitude of n:do Along Axial
Extension and Compression Paths.............................272

A.1 Formulas for Use in Inspecting the Nature of the
Quadratic Function Describing the Dilation Portion
of the Yield Surface ........................................278



2.1 Representation of plane stress state at a "point"............20

2.2 Typical stress-strain response of soil for a conventional
'triaxial' compression test (left) and a hydrostatic
compression test (right)................................ ....40

2.3 Typical stress paths used to investigate the stress-strain
behavior of soil specimens in the triaxial environment.......42

2.4 Components of strain: elastic, irreversible plastic,
and reversible plastic.......................................44

2.5 Common fabric symmetry types.................................48

2.6 Rate-independent idealizations of stress-strain response.....63

2.7 Two dimensional picture of Mohr-Coulomb failure criterion**..66

2.8 Commonly adopted techniques for locating the yield stress**..68

2.9 Yield surface representation in Haigh-Westergaard stress
space ......... .............................................71

2.10 Diagrams illustrating the modifying effects of the
coefficients Ai and A2: (a) A, = A2 = 1; (b) Ai 4 Az;
(c) Ai A2 Z = A- ****. .. .............. ...... .. ..... ......90

2.11 Schematic illustration of isotropic and kinematic
hardening *................................................ 92

2.12 Two dimensional view of an isotropically hardening
yield sphere for hydrostatic loading ............. ........97

3.1 In conventional plasticity (a) path CAC' is purely
elastic; in the proposed formulation (b) path CB'A is
elastic but AB"C' is elastic-plastico..................... .107


3.2 The current yield surface passes through the current stress
point and locally separates the domain of purely elastic
response from the domain of elastic-plastic response*.......108

3.3 Pictorial representation for sand of the nested set of
yield surfaces, the limit line, and the field of plastic
moduli, shown by the dep associated with a constant
value of n do ...........................................110
pq pq
3.4 Path independent limit surface as seen in q-p
stress space ................................................ 115

3.5 Axial compression stress-strain data for Karlsruhe sand
over a range of porosities and at a constant confinement
pressure of 50 kN/m2 ........................................ 116

3.6 Stress-strain response for a cyclic axial compression test
on loose Fuji River sand....................................117

3.7 Medium amplitude axial compression-extension test on loose
Fuji River sand *.... .......... ..............................119

3.8 Plastic strain path obtained from an anisotropic
consolidation test.*.................. .......................120

3.9 Plastic strain direction at common stress point.............121

3.10 Successive stress-strain curves for uniaxial stress or
shear are the initial curve' translated along the strain
axis in simplest model......................................123

3.11 Constant q/p ratio (as given by constant o,/o, ratio) at
zero dilation as observed from axial compression stress-
strain curves on dense Fountainbleau sand. Note that the
peak stress ratio decreases with increasing pressure*.......128

3.12 Characteristic state friction angles in compression and
extension are different, suggesting that the Mohr-Coulomb
criterion is an inappropriate choice to model the zero
dilation locus *........................... ...................129

3.13 Establishment of the yield surfaces from the inclination
of the plastic strain increment observed along axial
compression paths on Ottawa sand at relative densities of
(a) 39% (e=0.665), (b) 70% (e=0.555), and (c) 94% (e=0.465).132

3.14 Typical meridional (q -p) and octahedral sections
(inset) of the yield surface................................134



3.15 Trace of the yield surface on the triaxial q-p plane*.......135

3.16 Stress state in "thin" hollow cylinder......................152

3.17 Saada's hollow cylinder stress paths in Mohr's stress space-155

3.18 Measured vs. fitted response for hydrostatic compression
(HC) test using proposed model (po = 10 psi)*...............156

3.19 Measured vs. fitted response for axial compression test
(DC 0 or CTC of Figure 2.3) @30 psi using proposed model*..*157

3.20 Measured vs. predicted response for axial compression test
(DC 0 or CTC of Figure 2.3) @35 psi using proposed model....158

3.21 Measured vs. predicted response for axial compression test
(DC 0 or CTC of Figure 2.3) @45 psi using proposed model.** 159

3.22 Measured vs. predicted response for constant mean pressure
compression shear test (GC 0 or TC of Figure 2.3) using
proposed model ............................................ 161

3.23 Measured vs. predicted response for reduced triaxial
compression test (RTC of Figure 2.3) using proposed model.**162

3.24 Measured vs. predicted response for axial extension test
(DT 90 or RTE of Figure 2.3) using proposed model...........163

3.25 Volume change comparison for axial extension test*..........165

3.26 Results of axial compression tests on Karlsruhe sand at
various confining pressures and at a relative density
of 99% ....... .......... .......................... 169

3.27 Results of axial compression tests on Dutch dune sand
at various confining pressures and at a relative density
of 60.9% .............................. ...............171

3.28 Measured and predicted response for hydrostatic
compression test on Karlsruhe sand at 99% relative density..174

3.29 Measured and predicted response for axial compression test
(03 = 50 kN/m2) on Karlsruhe sand at 62.5% relative density-175

3.30 Measured and predicted response for axial compression test
(o3 = 50 kN/m2) on Karlsruhe sand at 92.3% relative density-176

3.31 Measured and predicted response for axial compression test
(03 = 50 kN/m2) on Karlsruhe sand at 99.0% relative density-177



3.32 Measured and predicted response for axial compression test
(o3 = 50 kN/m2) on Karlsruhe sand at 106.6% relative
density .................................................... 178

3.33 Measured and predicted response for axial compression test
(03 = 50 kN/m2) on Dutch dune sand at 60.9% relative
density .................................................... 179

3.34 Measured and predicted response for axial compression test
(03 = 200 kN/m2) on Dutch dune sand at 60.9% relative
density .....................................................180

3.35 Measured and predicted response for axial compression test
(o3 = 400 kN/m2) on Dutch dune sand at 60.9% relative
density .................................................. 181

3.36 Type "A" (top) and type "B" (bottom) stress paths of
Tatsuoka and Ishihara (1974b)...*...........................182

3.37 Observed stress-strain response for type "A" loading path
on loose Fuji River sand....................................183

3.38 Observed stress-strain response for type "B" loading path
on loose Fuji River sand .................................. 184

3.39 Simulation of type "A" loading path on loose Fuji River
sand using the simple representation........................187

3.40 Simulation of type "B" loading path on loose Fuji River
sand using the simple representation *......................189

3.41 Simulation of compression-extension cycle on loose Fuji
river sand using the simple representation................. 190

3.42 Conventional bounding surface adaptation with radial
mapping rule ............* ..................... ........ 193

3.43 Finite element mesh used in pressuremeter analysis..........201

3.44 Measured vs. predicted response for pressuremeter
test #1 ............................. ... ....... ..... ........ 202

3.45 Measured vs. predicted response for pressuremeter
test #2 ****............................................ ... 203

3.46 Measured vs. predicted response for pressuremeter
test #3 ................................................... 204




3.47 Measured vs. predicted response for pressuremeter
test #4..................................................... 205

3.48 Measured vs. predicted response for pressuremeter
test #5 ................................................... 206

3.49 Variation of principal stresses and Lode angle with
cavity pressure for element #1 and pressuremeter
test #2 ................................................... 208
3.50 Variation of plastic modulus with cavity pressure for
pressuremeter test #2.*..**........*........................209

3.51 Meridional projection of stress path for element #1,
pressuremeter test #2 *.....****.*** ........................211

3.52 Principal stresses as a function of radial distance
from axis of cavity at end of pressuremeter test #2 ........212

3.53 Experimental stress probes of Tatsuoka and
Ishihara (1974b)* ................. ..........................214

3.54 Shapes of the hardening control surfaces as
evidenced by the study of Tatsuoka and Ishihara (1974b)
on Fuji River sand ...........o...................... ...... 215

3.55 Illustration of proposed hardening control surface
and interpolation rule for reload modulus-*................ 217

3.56 Illustration of the role of the largest yield surface
(established by the prior loading) in determining the
reload plastic modulus on the hydrostatic axis .............218

3.57 Influence of isotropic preloading on an axial compression
test (03 = 200 kN/m2) on Karlsruhe sand at 99% relative
density***** ............. .****...................... 220

3.58 Predicted vs. measured results for hydrostatic
preconsolidation followed by axial shear................... 222

3.59 Shear stress vs. axial strain data for a cyclic axial
compression test on Reid-Bedford sand at 75% relative
density. Nominal stress amplitude q = 70 psi, and
confining pressure 03 = 30 psi******........................223

3.60 Prediction of the buildup of the axial strain data of
Figure 3.59 using proposed cyclic hardening representation-.226


3.61 Any loading starting in the region A and moving to region
B can go beyond the limit line as an elastic unloading or
a neutral loading path......................................227

4.1 Initial (a) and subsequent (b) configurations of the
deviatoric sections of the field of yield surfaces.........*234

4.2 Field of nesting surfaces in p-q (top) and Cp-q subspaces
(bottom) ............................................... .....242

4.3 Measured vs. fitted stress-strain response for axial
compression path using Prevost's model......................257

4.4 Initial and final configurations of yield surfaces for
CTC path (see Fig. 2.3) simulation..........................258

4.5 Measured vs. fitted stress-strain response for axial
extension path using Prevost's model........................259

4.6 Initial and final configurations of yield surfaces for
axial extension simulation*................................ 260

4.7 Measured vs. predicted stress-strain response for constant
mean pressure compression (or TC of Fig. 2.3) path using
Prevost's Model .............................................261

4.8 Initial and final configurations of yield surfaces for
TC simulation *...............................................262

4.9 Measured vs. predicted stress-strain response for reduced
triaxial compression (or RTC of Fig. 2.3) path using
Prevost's model............................................. 263

4.10 Initial and final configurations of yield surfaces for
RTC simulation*.............**...............................264

4.11 Measured vs. predicted stress-strain response for constant
pressure extension (or TE of Fig. 2.3) path using
Prevost's model.......*......................................265

4.12 Initial and final configurations of Yield Surfaces for
TE simulation...............................................266

D.1 Measured vs. predicted stress-strain response for DCR 15
stress path using proposed model............................286

D.2 Measured vs. predicted stress-strain response for DCR 32
stress path using proposed model ...........................287





D.3 Measured vs. predicted stress-strain response for DTR 58
stress path using proposed model............................288

D.4 Measured vs. predicted stress-strain response for DTR 75
stress path using proposed model ******.....................289

D.5 Measured vs. predicted stress-strain response for GCR 15
stress path using proposed model..........................**290

D.6 Measured vs. predicted stress-strain response for GCR 32
stress path using proposed model........................... 291

D.7 Measured vs. predicted stress-strain response for R 45
(or pure torsion) stress path using proposed model.........*292

D.8 Measured vs. predicted stress-strain response for GTR 58
stress path using proposed model ...........................293

D.9 Measured vs. predicted stress-strain response for GTR 75
stress path using proposed model ........................... 294

D.10 Measured vs. predicted stress-strain response for GT 90
stress path using proposed model*.............*.............295

E.1 Measured and predicted response for axial compression test
(03 = 400 kN/m2) on Karlsruhe sand at 92.3% relative
density******............ ........... .......................297

E.2 Measured and predicted response for axial compression test
(03 = 80 kN/m2) on Karlsruhe sand at 99.0% relative
density.................................... ...............298

E.3 Measured and predicted response for axial compression test
(o3 = 200 kN/m2) on Karlsruhe sand at 99.0% relative
density* ....******...*.............. ..................299

E.4 Measured and predicted response for axial compression test
(03 = 300 kN/m2) on Karlsruhe sand at 99.0% relative
density...o.......... ...........................*......... 300

G.1 Measured vs. predicted stress-strain response for DCR 15
stress path using Prevost's model .....*................... 303

G.2 Measured vs. predicted stress-strain response for DCR 32
stress path using Prevost's model ................*** *...... 304

G.3 Measured vs. predicted stress-strain response for DTR 58
stress path using Prevost's model*..........................305




G.4 Measured vs. predicted stress-strain response for DTR 75
stress path using Prevost's model...........................306

G.5 Measured vs. predicted stress-strain response for GCR 15
stress path using Prevost's model*..........................307

G.6 Measured vs. predicted stress-strain response for GCR 32
stress path using Prevost's model...........................308

G.7 Measured vs. predicted stress-strain response for R 45
stress path using Prevost's model ......................... 309

G.8 Measured vs. predicted stress-strain response for GTR 58
stress path using Prevost's model ...........................310

G.9 Measured vs. predicted stress-strain response for GTR 75
stress path using Prevost's model- ..........................311



dge, dop
e p
de de

dee, dE

e p
dekk, de
kk' kk







F (a)



II, 12, I,


parameter controlling shape of dilation portion of
yield surface

compression and swell indices

total, elastic, and plastic (small) strain increments

deviatoric components of dc, dee & de respectively

equal to /(3 de:de), /(3 dee:dee) & /(3 deP:dep)
2 2 2

incremental volumetric strain

incremental elastic and plastic volumetric strains

deviatoric components of do

stress increment

relative density in %

deviatoric components of strain e

initial voids ratio

elastic Young's modulus

failure or limit surface in stress space

yield surface in stress space

bounding surface in stress space

elastic shear modulus

function of Lode angle 6 used to normalize /J2

first, second & third invariants of the stress tensor o

initial magnitude of Ii for virgin hydrostatic loading

Io intersection of yield surface with hydrostatic axis
(the variable used to monitor its size)

(I0) magnitude of 10 for the largest yield surface
established by the prior loading

/Jz square root of second invariant of s
/J2 equivalent octahedral shear stress = /J2/g(6)

k parameter controlling size of limit or failure surface

k maximum magnitude of k established by the prior

kmob current mobilized stress ratio computed by inserting
the current stress state in the function f(o)

K elastic bulk modulus

K dimensionless elastic modulus number
K plastic modulus
K plastic modulus at conjugate point o
(K )o plastic modulus at the origin of mapping

m exponent to model curvature of failure meridion

n unit normal gradient tensor to yield surface

n exponent to control field of plastic moduli
interpolation function

n magnitude of n applicable to compression stress space

N slope of zero dilatancy line in /J2-I, stress space

NREP number of load repetitions

p mean normal pressure (=Ii/3)

Pa atmospheric pressure

po or pO initial mean pressure

q shear stress invariant, = /(3J2) = /(3 s..s..)
q equivalent shear stress invariant, = (3J g(
q equivalent shear stress invariant, = /(3J2)/g(9)










Y1, Y2, Y3





- -e -p
, E

e p
Ekk' kk


parameter controlling shape of consolidation portion of
yield surface
parameter to model the influence of 03 on E

parameter to model deviatoric variation of strength

deviatoric components of a

slope of dilation portion of *
yield surface at the origin of /J2-I. stress space

slope of radial line in /JJ-II stress space (below the
zero dilation line of slope N) beyond which the effects
of preconsolidation are neglected (0 < X < 1)

stress obliquity /J2/I,

scalar mapping parameter linking current stress state a
to image stress state o on hardening control surface

modified magnitude of B in proposed hardening option to
account for preconsolidation effects

reload modulus parameter for bounding surface hardening

reload modulus parameters for proposed cyclic hardening

Lame's elastic constant

distance from current stress state to conjugate stress

distance from origin of mapping to conjugate or image
stress state

Kronecker delta

components of small strain tensor

total, elastic, and plastic shear strain invariants,
/(3 e..e..), etc.
1J 13

total volumetric strain

elastic and plastic volumetric strains

Lode's parameter


A plastic stiffness parameter for hydrostatic compression

W Lame's elastic constant

v Poisson's Ratio

a components of Cauchy stress tensor

a stress tensor at conjugate point on bounding surface

a, 02, 03 major, intermediate, and minor principal stresses

r' 7az a radial, axial, and hoop stress components in
cylindrical coordinates

SMohr-Coulomb friction angle or stress obiliquity

^c Mohr-Coulomb friction angle observed in a compression
test (i.e., one in which a0 = 03)

Oe Mohr-Coulomb friction angle observed in an extension
test (i.e., one in which a0 = a2)

Ocv friction angle at constant volume or zero dilatancy

X ratio of the incremental plastic volumetric to shear
strain (= V3 deP /deP)

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



Devo Seereeram

May 1986

Chairman: Dr. Frank C. Townsend
Major Department: Civil Engineering

A special time-independent or elastic-plastic formulation is

developed through qualitative and quantitative comparisons with

experimental results reported in the literature. It does appear to

provide a simple yet adequate model for a number of key aspects of the

inelastic response of sands over a wide variety of loading paths. In

its simplest form, the material model is purely stress dependent and

exhibits no memory at all of prior inelastic deformation. Elementary

procedures are presented for matching the limit or failure surface, the

yield surface which passes through the current stress point for

unloading as well as loading, and the associated scalar field of purely

stress-dependent plastic moduli. Specific choices are presented for

several sands of different origin and initial density.

Based on well-known experimental investigations, a hardening

modification to the simple theory is proposed. The versatility of this

novel proposal is demonstrated by predicting the cyclic hardening

phenomenon typically observed in a standard resilient modulus test and

the influence of isotropic preconsolidation on a conventional triaxial

test. Another more conventional bounding-surface hardening option is

also described, and it is implemented to predict the results of a series

of cyclic cavity expansion tests.

For comparative evaluation with the proposed theory, a study of the

Prevost effective stress model is also undertaken. This multi-surface

representation was chosen because it is thought of as one of the most

fully developed of the existing soil constitutive theories.


1.1 The Role and Nature of Theory

In most fields of knowledge, from physics to political science, it

is essential to construct a theory or hypothesis to make sense of a

complex reality. The complex reality scrutinized in this dissertation

is the load-deformation behavior of a statistically homogenous

assemblage of unbound particles. More specifically, the mathematical

theory of plasticity is used as the basis for developing a constitutive

model for granular material. Such constitutive relations are of

fundamental importance in a number of areas of science and technology

including soil mechanics, foundation engineering, geophysics, powder

processing, and the handling of bulk materials.

The mathematical theories of plasticity of this study should be

clearly distinguished from the physical or microstructural plasticity

theories which attempt to model the local interaction of the granules.

A mathematical (or phenomenological) theory is only a formalization of

known experimental results and does not inquire very deeply into their

physical basis. It is essential, however, to the solution of problems

in stress analysis and also for the correlation of experimental data

(Drucker, 1950b).


To explain or model the complex phenomenon of particles crushing,

distorting, sliding, and rolling past each other under load, a theory

must simplify and abstract from reality. However, these simplification

and idealizations must lie within the realm of physically and

mathematically permissible stress-strain relations. The test of any

scientific theory is whether it explains or predicts what it is designed

to explain or predict, and not whether it exactly mirrors reality. The

most useful theory is the simplest one which will work for the problem

at hand. A theory can consider only a few of the many factors that

influence real events; the aim is to incorporate the most important

factors into the theory and ignore the rest.

1.2 Statement of the Problem

The characterization of the complex stress-strain response of

granular media is a subject which has generated much interest and

research effort in recent years, as evidenced by the symposia organized

by Cowin and Satake (1978), Yong and Ko (1980a), Pande and Zienkiewicz

(1980), Vermeer and Luger (1982), Gudehus and Darve (1984), and Desai

and Gallagher (1984), among others. This focusing of attention on

constitutive models is a direct consequence of the increasing use of the

finite element method to solve previously intractable boundary value

problems. Solutions obtained from this powerful computer-based method

are often precise to several significant digits, but this impressive

degree of precision loses its significance if the governing equations,

coupled with the constitutive assumptions or the imposed boundary

conditions, are inappropriate idealizations of the physical problem.

Progress in the area of theoretical modelling of soil response has

lagged conspicuously behind the state-of-the-art numerical solution

techniques. An all-encompassing stress-strain model for soil media, or

for that matter any other material, has yet to be formulated and

opinions differ as to whether such a task is even remotely possible. An

apparent drawback of all presently available constitutive relations is

that each has been founded on data gathered from standard laboratory

tests, and as Yong and Ko (1980b, p. 55) succinctly state, "the

relationships developed therefrom have been obviously conditioned to

respond to the soils tested as well as for the particular test system

constraints, and therefore the parameters used and material properties

sensed have been chosen to fit the test circumstance. Extension and

projection into a more general framework for wider use do not appear to

be sufficiently well founded."

Although the evolution of a fundamental set of constitutive

equations will benefit foundation engineering as a science, this

particular research effort was stimulated by the problem of rutting in

pavement base courses--in particular, the existing U.S. Air Force runway

system which is soon expected to be overloaded by a new generation of

heavier aircraft. Dr. Salkind (1984), the director of the Air Force

Office of Scientific Research (AFOSR), elucidates:

The relevance is extraordinarily high for this
nation. There is the obvious deterioration of our
highway system including potholes. The Air Force
has 3700 miles of runways around the world designed
for a 20 year life. Ninety-two percent are more
than 20 years old and 25 percent are significantly
deteriorated. The anticipated replacement cost
with today's technology is $1.9 billions. .The
underlying methodology is empirical and should be
put on a sound analytical basis. .The pavement
system, consisting of supporting soil,

underpavement, and paving material should be
analyzed for loads and moments (and loading
spectrum) recognizing the differing response of the
various layers with different material properties.
A basic science need is the lack of measuring
techniques for fundamental soil properties and
descriptions of soil constitutive properties.
Design is based on empirical values such as the
penetration of a standard cone. As soil is a
multi-phase mixture of solid particles, water, and
air, the challenge is to define what are the basic
fundamental properties (eg. soil "fabric" or
spatial arrangement of particles) and how such
properties change with loading (Personal
communication, October 12).

Ever since the pioneering work of Drucker and Prager (1952),

phenomenological plasticity theory has been developed and applied

extensively to model the mechanical behavior of soil. Constitutive

relations have grown increasingly complex as engineering mechanicians

have attempted to include the details of response for a broader spectrum

of loading paths. However, it is not clear that some of these more

sophisticated idealizations are better approximations of reality, or

whether they do capture the key aspects of soil behavior. The present

situation is complicated further by another problem: practicing

geotechnical engineers, the group most qualified to evaluate the

usefulness of these models, do not, for the most part, have a full and

working knowledge of tensor calculus and basic plasticity precepts.

They therefore tend to shun these potentially useful stress-strain

relations in favor of the simpler elastic and quasi-linear theories.

1.3 Approach

Using concepts recently advanced by Drucker and Seereeram (1986), a

new stress-strain model for granular material is introduced. This

representation incorporates those key aspects of sand behavior

considered most important and relevant, while also attempting to

overcome the conceptual difficulties associated with existing theories.

Many aspects of conventional soil plasticity theory are abandoned in

this novel approach:

1. The material is assumed to remain at yield during unloading in

order to simulate inelastic response (either "virgin" or

partially hardened) on reloading.

2. Plastic deformation is assumed possible at all stress levels

(i.e., there is a vanishing region of elastic response for

loading or reloading). The yield surface is not given by the

traditional permanent strain offset or tangent modulus

definitions, but by its tangent plane normal to the observed

plastic strain increment vector.

3. The consistency condition does not play a central role in the

determination of the plastic modulus. Instead, a scalar field

of moduli in stress space is selected to give the plastic

stiffness desired.

4. The limit surface is not an asymptote of or a member of the

family of yield surfaces. These distinct surfaces intersect at

an appreciable angle.

5. Hardening is controlled solely by changes in the plastic

modulus. Therefore, the surface enclosing the partially or

completely hardened region can be selected independently of the

size and shape of the current yield surface.

In its most elementary form, the model ignores changes in state

caused by the inelastic strain history. The field of plastic moduli

remains fixed and the yield surface expands and contracts isotropically

to stay with the stress point. Supplementary features, including

conventional work-hardening, bounding surface hardening, and cyclic

hardening or softening, can be added as special cases by some simple and

straightforward modifications to the basic hypotheses.

For comparative evaluation, a study of the Prevost (1978, 1980)

pressure-sensitive isotropic/kinematic hardening theory is also

undertaken. This model was chosen because it is thought of as the most

complete analytical statement on elasto-plastic anisotropic 'Ird-r.,in

theories in soil mechanics (Ko and Sture, 1980).

1.4 Scope

Chapter 2 attempts to elucidate the fundamentals of plasticity

theory from the perspective of a geotechnical engineer. It is hoped

that this discussion will help the reader, particularly one who is

unfamiliar with tensors and conventional soil plasticity concepts and

terminology, to understand the fundamentals of plasticity theory and

thus better appreciate the salient features of the new proposal.

Based on well-known observations on the behavior of sand, details

of the new theory are outlined in Chapter 3. Specific choices are

tendered for the analytical representations of 1) the yield surface, 2)

the scalar field of plastic moduli (which implicitly specifies a limit

or failure surface), and 3) the evolution of the yield surface. Several

novel proposals are also embedded in these selections.

A procedure is outlined for computing the model constants from two

standard experiments: a hydrostatic compression test and an axial

compression test. Each parameter is calculated directly from the

stress-strain data, and the initialization procedure involves no trial

and error or curve fitting techniques. Each parameter depends only on

the initial porosity of the sand. What is particularly appealing is

that all model constants can be correlated directly or conceptually to

stress-strain or strength constants, such as friction angle and angle of

dilation (Rowe, 1962), considered fundamental by most geotechnical


A number of hollow cylinder and solid cylinder test paths are used

to demonstrate the predictive capacity of the simple "non-hardening"

version of the theory. These tests include one series with a wide

variety of linear monotonic paths, another consisting of axial

compression paths on specimens prepared over an extended range of

initial densities and tested under different levels of confining

pressure, and still another sequence with more general load-unload-

reload stress paths, including one test in which the direction of the

shear stress is completely reversed. The range of the data permitted

examination of the influence of density, if any, on the magnitudes of

the model constants.

Although most of the predictions appeared satisfactory, many

questions are raised concerning the reliability of the data and the

probable limitations of the mathematical forms chosen for the yield

surface and the field of plastic moduli.

Two hardening modifications to the simple theory are described.

Unfortunately, both options sacrifice one important characteristic of

the simple model: the ability to model "virgin" response in extension

after a prior loading in compression, or vice-versa. The first, less

realistic option is an adaptation of Dafalias and Herrmann's (1980)

bounding surface theory for clay, which is itself an outgrowth of the

nonlinearly hardening model proposed by Dafalias and Popov (1975). Two

modifications to the simple theory transform it to the first hardening

option: 1) the largest yield surface established by the loading history

is prescribed as a locus of "virgin" or prime loading plastic moduli

(i.e., a bounding surface), and 2) for points interior to the bounding

surface, an image point is defined as the point at which a radial line

passing through the current stress state intersects the bounding

surface. Then the plastic modulus at an interior stress state is

rendered a function of the plastic modulus at the image point and the

Euclidean distance between the current stress state and the image point.

These constitutive equations are implemented in a finite element

computer code to predict the results of a series of cyclic cylindrical

cavity expansion tests.

Based on the observations of Poorooshasb et al. (1967) and Tatsuoka

and Ishihara (1974b), a second, more realistic hardening option is

proposed. It differs from the bounding surface formulation in that 1)

the shape of the surface which encloses the "hardened" region differs

from the shape of the yield surface, and 2) a special mapping rule for

locating the conjugate or image point is introduced. The versatility of

this proposed (cyclic) hardening option is demonstrated by predicting a)

the influence of isotropic preconsolidation on an axial compression

test, and b) the buildup of axial strain in a uniaxial cyclic

compression test.

In Chapter 4 the Prevost (1978, 1980) model is described. Althch-i

this theory has been the focus of many studies, the writer believes that

certain computational aspects of the hardening rule may have until now

been overlooked. These equations, appearing here for the first time in

published work, were gleaned from a computer program written by the

progenitors of the model (Hughes and Prevost, 1979).

Three experiments specify the Prevost model parameters: i) an axial

compression test, ii) an axial extension test, and iii) a one-

dimensional consolidation test, and although the initialization

procedure was followed with great care, this model seemed incapable of

realistically simulating stress paths which diverge appreciably from its

calibration paths. Because of this serious limitation, no effort was

expended beyond predicting one of the series of experiments used for

verifying the proposed model.


2.1 Introduction

It is the primary objective of this chapter to present and to

discuss in a methodical fashion the key concepts which form the

foundation of this dissertation. At the risk of composing this section

in a format which is perhaps unduly elementary and prolix to the

mechanicist, the author strives herein to fill what he considers a

conspicuous void in the soil mechanics literature: a discussion of

plasticity theory which is comprehensible to the vast majority of

geotechnical engineers who do not have a full and working knowledge of

classical plasticity or tensor analysis.

The sequence in which the relevant concepts are introduced is

motivated by the writer's background as a geotechnical

engineer--accustomed to the many empirical correlations and conventional

plane strain, limit equilibrium methods of analysis--venturing into the

field of generalized, elasto-plastic stress-strain relations. The terms

"generalized" and "elasto-plastic" will be clarified in the sequel. At

the beginning, it should also be mentioned that, although an attempt

will be made to include as many of the basic precepts of soil plasticity

as possible, this chapter will give only a very condensed and selected

treatment of what is an extensive and complex body of knowledge. In a

less formal setting, this chapter might have been titled "Plain Talk

About Plasticity For The Soils Engineer."

2.2 Tensors

2.2.1 Background

Lack of an intuitive grasp of tensors and tensor notation is

perhaps the foremost reason that many geotechnical engineering

practitioners and students shun the theoretical aspects of work-

hardening plasticity, and its potentially diverse computer-based

applications in geomechanics.

In this chapter, the following terms and elementary operations are

used without definition: scalar, vector, linear functions, rectangular

Cartesian coordinates, orthogonality, components (or coordinates), base

vectors (or basis), domain of definition, and the rules of a vector

space such as the axioms of addition, scalar multiple axioms and scalar

product axioms. Except where noted, rectangular Cartesian coordinates

are used exclusively in this dissertation. This particular set of base

vectors forms an orthonormal basis, which simply means that the vectors

of unit length comprising the basis are mutually orthogonal (i.e.,

mutually perpendicular).

Quoting from Malvern (1969, p.7),

Physical laws, if they really describe the physical
world, should be independent of the position and
orientation of the observer. That is, if two
scientists using different coordinate systems
observe the same physical event, it should be
possible to state a physical law governing the
event in such a way that if the law is true for one
observer, it is also true for the other.

Assume, for instance, that the physical event recorded is a spatial

vector t acting at some point P in a mass of sand, which is in

equilibrium under a system of boundary forces. This vector represents

some geometrical or physical object acting at P, and we can

instinctively reason that this "tangible" entity, t, does not depend on

the coordinate system in which it is viewed. Furthermore, we can

presume that any operations or calculations involving this vector must

always have a physical interpretation. This statement should not be

surprising since many of the early workers in vector analysis, Hamilton

for example, actually sought these tools to describe mathematically real

events. An excellent historical summary of the development of vector

analysis can be found in the book published by Wrede (1972).

Having established that the entities typically observed, such as

the familiar stress and strain vectors, are immutable with changes in

perspective of the viewer, we must now ask: How does one formulate

propositions involving geometrical and physical objects in a way free

from the influence of the underlying arbitrarily chosen coordinate

system? The manner in which this invariance requirement is

automatically fulfilled rests on the representation of physical objects

by tensors. To avoid any loss of clarity from using the word "tensor"

prior to its definition, one should note that a vector is a special case

of a tensor. There are several excellent references which deal with the

subject of vector and tensor analysis in considerably more detail than

the brief overview presented in the following. These include the books

by Akivis and Goldberg (1972), Hay (1953), Jaunzemis (1967), Malvern

(1969), Synge and Schild (1949) and Wrede (1972).

Although the necessity to free our physical law from the

arbitrariness implicit in the selection of a coordinate system has been

set forth, it is important to realize that this assertion is meaningless

without the existence of such coordinate systems and transformation

equations relating them. The transformation idea plays a major role in

the present-day study of physical laws. In fact, the use of tensor

analysis as a descriptive language for theoretical physics is largely

based on the invariant properties of tensor relations under certain

types of transformations. For example, we can imagine that the vector t

was viewed by two observers, each using a different rectangular

Cartesian coordinate system (say rotated about the origin with respect

to each other). As a result, an alternative set of vector components

was recorded by each scientist. Nonetheless, we should expect the

length of the vector--a frame indifferent quantity--computed by both

observers to be identical.

The transformation rules, which guarantee the invariant properties

of vectors and tensors, are actually quite simple, but they are very

important in deciding whether or not a quantity does indeed possess

tensorial characteristics. To illustrate how a vector is converted from

one rectangular Cartesian coordinate system to another, consider the

following example in which the "new" coordinate components and base

vectors are primed (') for distinction. The transformation from the old

basis (11,i2,i3) to the new basis (i:1,i,i1) can be written in the

matrix form
[cos(i1,il) cos(iz,i) cos(i3,ii
[Cl ,l,i, ] = [i i,2,i3] cos(i1 ,i ) cos(i2,i ) cos(i3,i )
cos(ii3) cos(i2,i) COS(i ('2.2.

(2.2.1 .1)

where cos(i,1,), for example, represents the cosine of the angle

between the base vectors i, and il. This is an ideal juncture to

digress in order to introduce two notational conventions which save an

enormous amount of equation writing.

The range convention states that when a small Latin suffix occurs

unrepeated in a term, it is understood to take all the values 1,2,3.

The summation convention specifies that when a small Latin suffix is

repeated in a term, summation with respect to that term is understood,

the range of summation being 1,2,3. To see the economy of this

notation, observe that equation is completely expressed as

i' Q k i (
-m mk k'
where Qm is equal to cos(iki). The index "m" in this equation is

known as the free index since it appears only once on each side. The

index "k" is designated the dummy index because it appears twice in the

summand and implies summation over its admissible values (i.e., 1,2,3).

The corresponding transformation formulas for the vector components

(t to t') can now be derived from the information contained in equation and the condition of invariance, which requires the vector

representations in the two systems to be equivalent. That is,

t = t i = t' = t' i (
k k m -m
Substituting the inverse of equation (i.e., i =Q i')
k kr ~r
into equation leads to

t Q i' = t' i'
k kr-r r -r

(t' tk Qkr) i' = 0,
r k kr '
from which we see

S= t kr (

With the invariance discussion and the vector transformation

example as background information, the following question can now be

asked: What actually is a tensor? It is best perhaps to bypass the

involved mathematical definition of a tensor and to proceed with a

heuristic introduction (modified from Malvern, 1969, and Jaunzemis,

1967). The discussion will focus on the particular type of tensor in

which we are most interested: second order (or second rank), orthogonal


Scalars and vectors are fitted into the hierarchy of tensors by

identifying scalars with tensors of rank (or order) zero and vectors of

rank (or order) one. With reference to indicial notation, we can say

that the rank of a tensor corresponds to the number of indices appearing

in the variable; scalar quantities possess no indices, vectors have one

index, second order tensors have two indices, and higher rank tensors

possess three or more indices. Every variable that can be written in

index notation is not a tensor, however. Remember that a vector has to

obey certain rules of addition, etc. or, equivalently, transform

according to equation These requirements for first order

tensors (or vectors) can be-generalized and extended for higher order


To introduce the tensor concept, let us characterize the state at

the point P (of, say, the representative sand mass) in terms of the

nature of the variable under scrutiny. If the variable can be described

by a scalar point function, it is a scalar quantity which in no way

depends on the orientation of the observer. Mass, density, temperature,

and work are examples of this type of variable.

Suppose now that there exists a scalar v(n) (such as speed)

associated with each direction at the point P, the directions being

described by the variable unit vector n. This multiplicity of scalars

depicts a scalar state, and if we identify this scalar with speed, for

instance, we can write
v = v [n] = v.n. (

where v(n) is the component of speed in the nth direction, and the

square brackets are used to emphasize that v, the velocity vector, is a

linear operator on n. Deferring a more general proof until later, it

can be said that the totality of scalars v(n) at a point is fully known

if the components of v are known for any three mutually orthogonal

directions. At the point P, therefore, the scalar state is completely

represented by a first order tensor, otherwise known as a vector.

The arguments for a second order tensor suggest themselves if one

considers the existence of a vector state at P; that is, a different

vector, t is associated with each direction n. Two important

examples of this type of tensor--the stress tensor and the strain

tensor--are discussed in some detail in the following.

2.2.2 The Stress Tensor

An example of second order tensors in solid mechanics is the stress

tensor. It is the complete set of data needed to predict the totality

of stress (or load intensity) vectors for all planes passing through

point P.

Recalling the routinely used Mohr circle stress representation, we

generally expect different magnitudes of shear stress and normal stress

to act on an arbitrary plane through a point P. The resultant stress

vector (or traction) t(n) is unique on each of these planes and is a

function of n at the point P, where n is the unit vector normal to a

specified plane. In order to describe fully the state of stress at P, a

relationship between the vectors t and n must be established; in

other words, we seek a vector function of a single vector argument n.

It turns out that we are in fact seeking a linear vector function, say

a, which is a rule associating the vector t with each vector n in the

domain of definition. A linear vector function is also called a linear

transformation of the domain or a linear operator acting in the domain

of definition of the function a.

A second order extension of equation is

t = ; [n], (

where again the square brackets imply a linear operation. The linearity

assumption of the function a implies the following relationships:

[(n, + 2n)/nit + n |.] = o[n ] + aCn2] (

for arbitrary unit vectors n, and n2, and

a[an] = a G[n] (

for arbitrary unit vector n and real number a.

Geometrically, equation means that the operator a carries

the diagonal of the parallelogram constructed on the vectors nj and n2

into the diagonal of the parallelogram constructed on the vectors ti =

[oCn] and t2 = 2En2]. Equation means that if the length of the

vector n is multiplied by a factor a, then so is the length of the
vector t = o[n].

Using a rectangular Cartesian coordinate system, the traction

vector t(n) and the unit normal vector n can each be resolved into their

(n) (n) (n)
components ti t2 t3 and n,, n2, n3 respectively. The linear

relationship between t and n can be expressed in the matrix form

(n) (n) (n) '111 012 (J13
[ti ,t2 ,t3 ] = [nln2,n3] 21 022 23 (
31 032 033

or alternatively, in the indicial notation,
t = o. n. (
i 31 J
where the components of the 3x3 matrix a are defined as the stress

tensor acting at point P. Note that the wavy underscore under symbols

such as "o" is used to denote tensorial quantities; however, in cases

where indices are used, the wavy underscore is omitted.

In general, tensors can vary from point to point within the

illustrative sand sample, representing a tensor field or a tensor

function of position. If the components of the stress tensor are

identical at all points in the granular mass, a homogenous state of

stress is said to exist. The implication of homogeneity of stress (and

likewise, strain) is particularly important in laboratory soil tests

where such an assumption is of fundamental (but controversial)

importance in interpreting test data (Saada and Townsend, 1980).

Second order tensors undergo coordinate transformations in an

equivalent manner to vectors (see equation For a pure

rotation of the basis, the transformation formula is derived by employing

a sequence of previous equations. Recall from equation that

t' = t Q
r k kr'
and by combining this equation with equation, we find that

tr = jk nj Qkr (
r 3k jkr*

Furthermore, n in this equation can be transformed to n' resulting in

t' = jk Qjs n' Qr (
r jk s s kr
The left hand side of equation can also be replaced by the

linear transformation so that

o' n' = o Q. n' Q
pr p jk s kr
which when rearranged yields

o' n' Q. n' Q = 0. (
pr p jk js s kr
All the indices in equation are dummy indices except "r"--

the free index. A step that frequently occurs in derivations is the

interchange of summation indices. The set of equations is unc:hanrd if

the dummy index "p" is replaced by the dummy index "s." This

manipulation allows us to rewrite equation in the form

o' n' o Q. n' Q = 0,
sr s jk js s kr '
and by factoring out the common term n' we obtain
(a' a Q Q ) n' = 0.
sr jk js kr ns
From this equation, the tensor transformation rule is seen to be

's = Ojk Qs Qr' (
sr jk is kr'
or in tensor notation,

o' = Q 9 9 (

It was previously stated (without verification) that a vector is

completely defined once its components for any three mutually orthogonal

directions are known. The reciprocal declaration for a second order

tensor will therefore be that the components of a second order tensor

are determined once the vectors acting on three mutually orthogonal

planes are given. For the particular case of the stress tensor, this

statement can be substantiated by inspecting the free body diagram of

Figure 2.1 (note that this is not a general proof). Here, a soil prism



21 (Ty x)

'22 (y)

Representation of plane stress state at a "point"

0*22 (y )

2, (xy)



a0 (0-x)

012 ( Txy)

Lj- '3 --

0/ '21 ( yx )

Figure 2.1


is subject to a plane stress state, plane stress simply meaning there is

no resultant stress vector on one of the three orthogonal planes;

therefore, the non-zero stress components occupy a 2x2 matrix instead of

the generalized 3x3 matrix. Generalized, in this context, refers to a

situation where the full array of the stress tensor is considered in the

problem, and when used as an adjective to describe a stress-strain

relationship, the word tacitly relates all components of strain (or

strain increment) to each stress (or stress increment) component for

arbitrary loading programs.

Figure 2.1 shows the two-dimensional free body diagram of a

material prism with a uniform distribution of stress vectors acting on

each plane; note that the planes AB and BC are perpendicular. By taking

moments about the point D, it can be shown that Txy = Ty and this is

known as the theorem of conjugate shear stresses, a relationship which

is valid whenever no distributed body or surface couple acts on the

element. This two dimensional observation can be generalized to three

dimensions, where as a consequence, the 3x3 stress tensor matrix is

symmetric. Symmetry implies that only six of the nine elements of the

3x3 matrix are independent.

By invoking force equilibrium in the x- and y-directions of Figure

2.1, the two resulting equations can be solved simultaneously for the

unknowns Te and oe, thus verifying that the shear and the normal stress

(or the stress vector in this case) on an arbitrary plane can be

computed when the stress vectors on perpendicular planes are known.

Extension of this two-dimensional result to three dimensions reveals

that the components of three mutually perpendicular traction vectors,

acting on planes whose normals are the reference axes, comprise the rows

of the stress tensor matrix.

Most geotechnical engineers are familiar with the Mohr-Coulomb

strength theory for granular soils. This criterion specifies a limit

state (or a locus in stress space where failure occurs with "infinite"

deformations) based on a combination of principal stresses (oa, 02, and

03). As will be described in a later section on plasticity, even the

more recently proposed failure criteria for soils are also only

functions of the principal stresses. This is the motivation for

presenting the following procedure for computing the principal stresses

from the frame-dependent components of a.

A principal plane is a plane on which there are no shear stresses.

This implies that the normal stress is the sole component of the

traction vector acting on such a plane, and the geometrical

interpretation is that the traction vector and the unit normal vector

(n) to the plane at a point both have the same line of action.

Mathematically, the principal plane requirement can be expressed as
( = A n, (

or in indicial notation,
ti = A n (

where A is the numerical value sought. Remember that there are, in

general, three principal planes and therefore three principal values

(A,, A2, and A3).

Substituting equation into equation and

rearranging leads to

aji nj A ni = O. (

As an aid to solving this equation for A, an extremely useful

algebraic device, known as the Kronecker delta 6, is now introduced. It

is a second order tensor defined as

6. = if i = (
ij 0 if i j
By writing out the terms in long form, one may easily verify that

n. = 6.. n.. (

Equation can now be substituted into equation to


.ji nj A 6ij nj = 0,
31 *J 13 3

(o.. A 6..) n. = 0. (

For clarity, equation is expanded out to

(ol1 A) nl + o12 n2 + 013 n3 = 0

021 n" + (022 A) n. + 023 n3 = 0, (

a31 n, + 032 nz + (o33 A) n3 = 0

which may be organized in the matrix form
o 1-A 012 013 n1 0
021 o22-A 023 n2 = 0 (
031 032 o33-A n 01

and where it is seen to represent a homogenous system of three linear

equations in three unknowns (nl, n2, and n3) and contains the unknown

parameter A. The fourth equation for solving this system is provided by

the knowledge that

n*n = ni n. = 1, (

since n is a unit vector.

Equation has a nontrivial solution if and only if the

determinant of the coefficient matrix in equation is equal to

zero (see, for example, Wylie and Barrett, 1982, p.188). That is,

o11-A o12 013
021 o22-A 23 = 0 (2.2.2.
031 032 033-A

must be true for non-trivial answers.

This determinant can be written out term by term to give a cubic

equation in A,

A3 II A2 12 A 13 = 0, (2.2.2.;

where the coefficients

I1 = 011 + 022 + 033 = okk, (

12 = -(011o22 + 022033 + 033011) + 0 + 31 + 02

= ( o ao I2 ) + 2, (
13 ij






011 012 013
13 = 021 022 023 3 (
031 032 033

Since this cubic expression must give the same roots (principal

stresses) regardless of the imposed reference frame, its coefficients--

the numbers Ii, I, and 13--must also be independent of the coordinate

system. These are therefore invariant with respect to changes in the

perspective of the observer and are the so-called invariants of the

stress tensor a. The notation Ii, 12, and 13 are used for the first,

second, and third invariants (respectively) of the stress tensor a.

When provided with a stress tensor that includes off-diagonal terms

(i.e., shear stress components), it is much simpler to compute the

invariants as an intermediate step in the calculation of the principal

stresses. Of course, writing the failure criterion directly in terms of

the invariants is, from a computational standpoint, the most convenient

approach. In any event, one should bear in mind that the stress

invariants and the principal stresses can be used interchangeably in the

formulation of a failure criterion. The following discussion centers on

a typical methodology for computing the principal stresses from the

stress invariants.

Start by additively decomposing the stress tensor into two

components: 1) a spherical or hydrostatic part (p 6, ), and 2) its

deviatoric components (s i). The first of these tensors represents the

average pressure or "bulk" stress (p) which causes a pure volumetric

strain in an isotropic continuum. The second tensor, s, is associated

with the components of stress which bring about shape changes in an

ideal isotropic continuum. The spherical stress tensor is defined as p

6. where p is the mean normal pressure ( kk/3 or Ii/3) and 6ij is the
1J kk 1J
Kronecker delta. Since, by definition, we know the spherical and

deviatoric stress tensors combine additively to give the stress tensor,

the components of the stress deviator (or deviatoric stress tensor) are

sj = ij p 6ij (2..2.2.5)

where compression is taken as positive. This particular sign convention

applies throughout this dissertation.

The development of the equations for computing the principal values

and the invariants of a apply equally well to the stress deviator s,

with two items of note: a) the principal directions of the stress

deviator are the same as those of the stress tensor since both represent

directions perpendicular to planes having no shear stress (see, for

example, Malvern, 1969, p.91), and b) the first invariant of the stress

deviator (denoted by J1) is equal to zero. The proof of the latter


J1 = s1 + 322 + s33

= 1 22 1 + 33 I,
3 3 3
and by recalling equation, it is clear that

J, = 0. (

From the last equation and equation, observe that the

second invariant of the stress deviation (denoted by J2) is simply

J2 = (s..s. ) + 2. (

Denoting the third invariant of the stress deviation by J3, the

cubic expression for the stress deviator s, in analogy to equation for the stress tensor o, becomes

A3 J2 A J3 = 0, (

where the roots of A are now the principal values (or more formally, the

eigenvalues) sj, s,, and s, of the stress deviator s. Since the

coefficient (i.e., J1) of the quadratic term (A2) is zero, the solution

of equation is considerably easier than that of equation It is therefore more convenient to solve for the principal

values of s and then compute the principal values of a using the


oi = si + p, 02 = s2 + p, and 03 = S3 + p. (

The direct evaluation of the roots, A, of equation is not

obvious until one observes the similarity of this equation to the

trigonometric identity

sin 36 = 3 sine 4 sin3e.

Dividing through by four and rearranging shows the relevancy of this


sin3e 3 sine + 1 sin 36 = 0. (

Replacing A with r sine in equation gives

r3 sin3e J, r sine J3 = 0,

which when divided through by r3 gives

sin3e J_ sine J_ = 0. (
r2 r3
A direct correlation of this equation with equation shows that

o r

r = 2 /J2,

J3 = 1 sin 30,

sin 3 = 4

sin 30 = 4 J3.



Substitution of the negative root of equation into

equation leads to

sin 36 = [3/3 (J3//J23)], (
from which we find that

e = 1 sin-1 [3/3 (J3/JJz3)], (
3 2
where 6 is known as the Lode angle or Lode parameter (Lode, 1926). As

will be described in a later section on plasticity, the Lode angle is an

attractive alternative to the J3 invariant because of its insightful

geometric interpretation in principal stress space. Physically, the

Lode angle is a quantitative indicator of the relative magnitude of the

intermediate principal stress o2 with respect to oi and a3.

Owing to the periodic nature of the sine function, the angles 30,

30 + 2i, and 30 + 47 all give the same sine in terms of the calculated

invariants of the deviator in equation If we further restrict

30 to the range + (i.e., 5< 0 +r), the three independent roots of

the stress deviator are furnished by the equations (after Nayak and

Zienkiewicz, 1972)

s, = 2 /J2 sin(6 + 4 t), (

s, = 2 /J2 sin(e), (

S3 = 2 /J2 sin(0 + 2 1). (2.2.2.
7-3 7
Finally, these relations can be combined with those of equation to give the principal values of the stress tensor a,
0ol sin (e + 4/3 Ti)
oz = 2 /J sin 0 + 1 I, (2.2.2.
[oz- (sin (6 + 2/3 T) I

To gain a clearer understanding of how the Lode angle 0 accounts

for the influence of the intermediate principal stress, observe from

this equation that
6 = sin1 [oi + a3 2 02], 300 < 0 < 300. (2.2.2.
2 /(3 J2)




2.2.3. The Strain Tensor

The mathematical description of strain is considerably more

difficult than the development just presented for stress. Nevertheless,

a brief'introduction to the small strain tensor is attempted herein,

while the interested reader should refer to a continuum mechanics

textbook to understand better the concept and implications of finite

deformation. This presentation has been modified from Synge and Schild


Most soils engineers are familiar with the geometrical measure of

unit extension, e, which is defined as the change in distance between

two points divided by the distance prior to straining or

e = (Li Lo) + Lo, (

where Lo and Li are respectively the distances between say particles P

and Q before and after the deformation. If the coordinates of P and Q

are denoted by x (P) and x (Q) respectively, we know that

L2 = [x (P) x (Q)] [x (P) x (Q)] (

from the geometry of distances.

Further, if the particles P and Q receive displacements u (P) and

ur(Q) respectively, the updated positions (using primed coordinates for

distinction) are

x'(P) = x (P) + ur(P), (


x'(Q) = x (Q) + u (Q). (

The notation u r(P) and u (Q) indicates that the particles undergo

displacments which are dependent on their position. Note that if the

displacement vector, u, is exactly the same for each and every particle

in the medium, the whole body translates without deforming--a rigid body


From equations and, we find that

L2 = [x'(P) x'(Q)] [x'(P) x'(Q)],
r r r r
= [Xr(P) + ur(P) xr(Q) Ur(Q)] x

[x (P) + ur(P) Xr(Q) u (Q)], (

and subtracting equation from this equation leads to

L' L' = [x (P) + ur(P) xr(Q) u (Q)][x (P) + u (P) -

xr(Q) Ur(Q)] [xr(P) (Q(P) (Q)] [(P) (Q)

which when reordered gives

L1 L2 = [ur(Q) u (P)][Cu(Q) u (P)] +

2 [x (Q) xr(P)][Cu(Q) ur(P)]. (

If attention is fixed on point P and an infinitesimally close

particle Q, the description of the state of strain at P can be put in a

more general form than the uniaxial unit extension measure. Since the

distance between P and Q is assumed small, the term

[xr(Q) xr(P)] [Cx(Q) Xr(P)]

and its higher orders are negligible; a Taylor expansion about P is

therefore approximately equal to

u (Q) u (P) = au r/3xsl x [x(Q) x (P)]. (

Substitution of this equation into equation gives

L- L = 3u /xsP [x (Q) x(P)] u / Ext, [xt(Q) xt(P)] +
r SP s s r P t
2 Ex (Q) x (P)] au r/xIp Cxm(Q) xm(P)]. (

Furthermore, we know approximately that

[xr(Q) Xr(P)] = Lo nr, (

where n are the components of the unit vector directed from P to Q;

substitution of this relation into equation gives

L2 L2 = au /ax p Lo ns ur/xt LO nt +
1 s0 s prtP t

2 Lo n 3u u/a3x I Lo

= L [aur /x lp ns Du /xtlp nt +

2 n 3ur /ax p nm],
r r m P m


L1 Lo = [aur/3xslp ns u r/9xtP n +

2 n au r/ax nm]. (

If an assumption is made that the strain is small, aur/axtlp is

small and hence the product

ur /ax s I aur/ax t

in the last equation is negligible. Therefore, for small strain

L L = 2 nr aur /a3x n. (


Lj L2 = L, Lo L, + Lo
LL Lo Lo

and with the


= L Lo L, Lo + 2 L9
Lo Lo

= Li Lo [L1 Lo + 2]
Lo Lo

e ( e + 2 ), (

assumption of small strain, e2 is negligible, which implies

L L = 2 e.


By equating the previous equation with equation, one finds


e = nr u r/ax I n (
r r mP m
If the components of the small strain tensor at point P are now

defined as

Es =1 [au /3x + au /ax ], (
rs r s s r

then the unit extension of every infinitesimal line emanating from P in

the arbitrary direction n is given by

e = E n n (
rs r s
Soil engineers may wonder how the traditional shear strain concept

enters this definition of strain. It can be shown (see, for example,

Malvern, 1969, p.121) that the off-diagonal terms of the tensor e are

approximately equal to half the decrease, Yrs' in the right angle

initially formed by the sides of an element initially parallel to the

directions specified by the indices r and s. This only holds for small

strains where the angle changes are small compared to one radian.

Another important geometrical measure in studying soil deformation

is the volume change or dilatation. The reader can easily verify that

the volume strain is equal to the first invariant (or trace) of the

strain tensor E (or in indicial notation, mm).

In analogy to the stress deviator, the strain deviator (denoted by

e) is given by

e.. = .. 1 mm .j' (

and since, like stress, strain is a symmetric second order tensor, the

corresponding discussion for principal strains and invariants parallels

the previous development for the stress tensor. In analogy to the

stress invariant /(3J2), the shear strain intensity is given by

E = /(3 e..e .). (
1J J1

2.3 Stress-Strain Equations and Constitutive Theory

To solve statically indeterminate problems, the engineer utilizes

the equations of equilibrium, the kinematic compatibility conditions,

and a knowledge of the load-deformation response (or stress-strain

constitution) of the engineering material under consideration. As an

aside, it is useful to remind the soils engineer of two elementary

definitions which are not part of the everyday soil mechanics

vocabulary. Kinematics is the study of the motion of a system of

material particles without reference to the forces which act on the

system. Dynamics is that branch of mechanics which deals with the

motion of a system of material particles under the influence of forces,

especially those which originate outside the system under consideration.

For general applicability, the load-deformation characterization of

the solid media is usually expressed in the form of a constitutive law

relating the force-type measure (stress) to the measure of change in

shape and/or volume (strain) of the medium. A constitutive law

therefore expresses an exact correspondence between an action (force)

and an effect (deformation). The correspondence is functional--it is a

mathematical representation of the physical processes which take place

in a material as it passes from one state to another. This is an

appropriate point to interject and to briefly clarify the meaning of

another word not commonly encountered by the soils engineer: functional.

Let us return to the sand mass which contains particle P and extend

the discussion to include M discrete granules (Pi, i= 1,2,. .,M). Say

the body of sand was subjected to a system of boundary loads which

induced a motion of the granular assembly, while an observer, using a

spatial reference frame x, painstakingly recorded at N prescribed time

intervals the location of each of the M particles. His data log

therefore consists of the location of each particle M (xM) and the time

at which each measurement was made (t'). At the current time t (a t'),

we are interested in formulating a constitutive relationship which gives

us the stress at point P, and in our attempt to construct a model of

reality, we propose that such a relation be based on the MN discrete

vector variables we have observed; i.e., the M locations x (in the

locality of point P) at N different times t' (< t). In other words,

stress at P is a function of these MN variables. This function

converges to the definition of a functional as the number of particles M

and the discrete events in the time set t' approach infinity.

For our simplest idealization, we can neglect both history and time

dependence, and postulate that each component of current stress o

depends on every component of the current strain tensor E and tender a

stress-strain relationship of the form

a ij Cijkl kl' (2.3.1)

or inversely,

kl = Dklij ij' (2.3.2)
where the fourth order tensors Cikl and Dkli (each with 81 components)

are called the stiffness and compliance tensors respectively. Both of

these constitutive tensors are discussed in detail later in this

section. Note that the number of components necessary to define a

tensor of arbitrary order "n" is equal to 3"

Because the behavior of geologic media is strongly non-linear and

stress path dependent, the most useful constitutive equations for this

type of material are formulated in incremental form,

ij = Cijki l' (2.3.3)

or conversely,

kl =D klij ij, (2.3.4)
k1 klij ij'
where the superposed dot above the stress and strain tensors denote a

differentiation with respect to time. In these equations, C and D are

now tangent constitutive tensors. The terms a and Z are the stress rate

and strain rate respectively.

If the "step by step" stress-strain model is further idealized to

be insensitive to the rate of loading, the incremental relationship may

be written in the form

doij = Cijkl dEkl' (2.3.5)

or inversely,

dEkl = Dklij doij (2.3.6)

where do and dE are the stress increment and strain increment

respectively, and C and D are independent of the rate of loading. Only

rate-independent constitutive equations are considered in this thesis.

The formulation, determination, and implementation of these

constitutive C and D tensors are the primary concern of this research.

In the formulation of generalized, rate independent, incremental

stress-strain models, the objective is one of identifying the variables

that influence the instantaneous magnitudes of the components of the

stiffness (C) or compliance (D) tensors. Such a study bears resemblance

to many other specialized disciplines of civil engineering. The

econometrician, for instance, may determine by a selective process that

the following variables influence the price of highway construction in a

state for any given year: cost of labor, cost of equipment, material

costs, business climate, and a host of other tangible and intangible

factors. The soils engineer, perhaps using the econometrician's

techniques of regression analysis and his personal experience, can

easily identify several factors which influence soil behavior. From our

basic knowledge of soil mechanics, we might make the following

preliminary list: 1) the void ratio or dry unit weight--perhaps the most

important measure of overall stiffness and strength of the material; 2)

the composition of the grains, which includes information on the mineral

type (soft or hard), particle shape, angularity of particles, surface

texture of particles, grain size distribution, etc.; 3) the orientation

fabric description or anisotropy of the microstructure; 4) the stress
history a or stress path, which may be used, for example, to indicate

how close the current stress state is to the failure line, the number of

cycles of loading, the degree of overconsolidation, etc.; 5) the

magnitude and direction of the stress increment do; 6) the rate of

application of the stress increment; and 7) the history of the strain
S, from which one may compute, for example, the current void ratio and

magnitude of the cumulative permanent distortion.

In writing general mathematical formulations, it is convenient to

t t
lump all variables--except for a e and do--as a group known as the

set of n internal variables 9i (i = 1,2,3,. .,n). These internal

variables represent the microstructural properties of the material. A

generalized, rate-independent, incremental stress-strain functional dE

can therefore be put in the form

^ t t
dg = dE (o do, ). (2.3.7)

This means that the components of the compliance (or stiffness) tensor

t t
depends on o c do (and its higher orders), and g .

One basic difference between the econometrician's model and the

mechanician's load-deformation model must be emphasized: the mechanician

is dealing with dependent and independent variables which are physically

significant, but the econometrician uses variables which may frequently

be intangible. Therefore, in the selection of constitutive variables

(such as stress and strain) and in the actual formulation of the stress-

strain equations, certain physical notions (leading to mathematical

constraints) must be satisfied. These conditions are embodied in the

so-called axioms or principles of constitutive theory. An axiom is a

well-established basis for theoretical development. Since geotechnical

engineers are, for the most part, interested in isothermal processes,

the principles linked to thermomechanical behavior are suppressed in the


The Axiom of Causality states that the motion of the material

points of a body is to be considered a self-evident, observable effect

in the mechanical behavior of the body. Any remaining quantities (such

as the stress) that enter the entropy production and the balance

equations--i.e., the equations of conservation of mass, balance of

momentum, and conservation of energy--are the causes or dependent

variables. In other words, there can occur no deformation (effect)

without an external force (cause).

The Principle of Determinism is that the stress in a body is

determined by the history of the motion of that body. This axiom

excludes the dependence of the stress at a point P on any point outside

the body and on any future events. This phenomenon is sometimes

referred to as the Principle of Heredity.

In the purely mechanical sense, the Axiom of Neighborhood or Local

Action rules out any appreciable effects on the stress at P that may be

caused by the motion of points distant from P; "actions at a distance"

are excluded from constitutive equations.

During the discussion of stress and strain, it was made quite clear

that the tensor measures should be independent of the perspective of the

observer. It is therefore natural to suggest a similar constraint for

the constitutive equations: C and D must be form-invariant with respect

to rigid motions (rotation and/or translation) of the spatial frame of

reference. This is termed the Principle of Material Frame Indifference

or Objectivity.

Finally, the Axiom of Admissibility states that all constitutive

equations must be consistent with the basic principles of continuum

mechanics; i.e., they are subject to the principles of conservation of

mass, balance of moment, conservation of energy, and the entropy


2.4 A Note on Stress and Strain in Granular Media

The concepts of stress and strain discussed in the previous

sections are closely associated to the concept of a continuum, which

effectively disregards the molecular structure of matter and treats the

medium as if there were no holes or gaps. The following quotation from

Lambe and Whitman (1969, p.98) succintly summarizes the applicability of

the continuum stress measure to granular materials:

when we speak of the stress acting at a
point, we envision the forces against the sides of
an infinitesimally small cube which is composed of
some homogenous material. At first sight we may
therefore wonder whether it makes sense to apply
the concept of stress to a particulate system such
as soil. However, the concept of stress as applied
to soil is no more abstract than the same concept
applied to metals. A metal is actually composed of
many small crystals, and on the submicroscopic
scale the magnitude of the forces vary randomly
from crystal to crystal. For any material, the
inside of the infinitesimally small cube is thus
only statistically homogenous. In a sense all
matter is particulate, and it is meaningful to talk
about macroscopic stress only if this stress varies
little over distances which are of the order of
magnitude of the size of the largest particle.
When we talk about about stresses at a "point"
within a soil, we often must envision a rather
large "point."

Local strains within a statistically homogenous mass of sand are

the result of distortion and crushing of individual particles, and the

relative sliding and rolling velocities between particles. These local

strains are much larger than the overall (continuum) strain described in

section 2.2.3. The magnitude of the generated strain will, as mentioned

before, depend on the composition, void ratio, anisotropic fabric, past

stress history, and the stress increment. Composition is a term used in

soil mechanics to refer to the average particle size, the surface

texture and angularity of the typical grain, the grain size

distribution, and the mineral type.

Figure 2.2 illustrates typical qualitative load-deformation

response of loose and dense soil media subject to two conventional

laboratory stress paths: hydrostatic compression, and conventional

w I


u u
Z u.






1 ^ -


o 2


Figure 2.2


'6k V0

Typical stress-strain response of soil for a conventional 'triaxial' compression test (Left)
and a hydrostatic compression test (right)

triaxial compression. Figure 2.3 shows these paths together with an

assortment of other "triaxial" stress paths used for research as well as

routine purposes. In this context, note that the adjective "triaxial"

is somewhat ambiguous since this particular test scenario dictates that

the circumferential stress always be equal to the radial stress. The

stress state is therefore not truly triaxial, but biaxial. As we can

gather from Figure 2.2, the stress-strain behavior of soil is quite

complicated, and in order to model approximately real behavior, drastic

idealizations and simplifications are necessary. More complex details

of soil response are mentioned in Chapter 3.

The major assumptions in most present idealizations are that: a)

soil response is independent of the rate of loading, b) behavior may be

interpreted in terms of effective stresses, c) the interaction between

the mechanical and thermal processes is negligible, and d) the strain

tensor can be decomposed into an elastic part (Ee) and a plastic

conjugate (P ) without any interaction between the two simultaneously

occurring strain types,

e p
= + (2.4.1)

or in incremental form,
e p
de = de + dE (2.4.2)

The elastic behavior (e or dee) is modeled within the broad

framework of elasticity theory, while the plastic part (E or de ) is

computed from plasticity theory. Both these theories will be elaborated

later in this chapter.

With the introduction of the strain decomposition into elastic and

plastic components, it is now important to emphasize the difference

between irreversible strains and plastic strains for cyclic loading on

q-y- Ox


Figure 2.3 Typical stress paths used to investigate the stress-
strain behavior of soil specimens in the triaxial


Conventional Triaxial CTC Acrx = AOz = 0 A cy > 0

Hydrostatic Compression HC Acrx AOz AO, > 0

Conventional Triaxial
Extension CTE d =Acr>0; ACT, -0

Mean Normal Pressure TC crx + Aocz + Aoy = 0;
Triaxial Compression ATy>Ao-x (=Aaz )

Mean Normal Pressure TE aCT +* Cz + AOy =0;
Triaxial Extension A cTx =Aoz>Aaoy

Reduced Triaxial RTC Ox=Az<; A 0
Compression TC A = <; r =

Reduced Triaxial RTE A Extension


soils. Consider a uniaxial cyclic test consisting of a virgin loading,

an unloading back to the initial hydrostatic state of stress, and a

final reloading to the previous maximum deviatoric stress level. During

the first virgin loading both elastic and plastic strains are generated,

and these components may be calculated using an elastic and a plastic

theory respectively. If at the end of this segment of the stress path

we terminate the simulation and output the total, elastic, and plastic

axial strains, one may be tempted to think that the plastic component

represents the irrecoverable portion of the strain. However, when the

stress path returns to the hydrostatic state, the hysteresis loop in

Figure 2.4 indicates that reverse plastic strains are actually generated

on the unload and a (small) portion of the plastic strain at the end of

the virgin loading cycle is, in fact, recovered. This is an

illustration of the Bauschinger effect (Bauschinger, 1887). Therefore,

for such a closed stress cycle, the total strain can more generally be

broken down into the three components:

P P e
E E e + E ,
i~ rrev -rev
where is the irreversible plastic strain e is the reverse
irrev -rev
plastic strain, and as before, e denotes the elastic strain, which is

by definition recoverable. Some complicated models of soil behavior,

such as the one described in Chapter Four, allow for reverse plastic

strains on such "unloading" paths. However, ignoring this aspect of

reality, as is done in Chapter Three, can lead to very rewarding


Three broad classes of continuum theories have evolved in the

development and advancement of soil stress-strain models (Cowin, 1978):



irrev rev



SE rev E-

Figure 2.4 Components of strain: elastic, irreversible
elastic, and reversible plastic

1) the kinematically ambiguous theories, 2) the phenomenological

theories, and 3) the microstructural theories.

The kinematically ambiguous hypotheses employ the stress equations

of equilibrium in conjunction with a failure criterion to form a system

of equations relating the components of the stress tensor. This

category is referred to as kinematically ambiguous because displacements

and strains do not appear in and are therefore not computed from the

basic equations of the theory. They assume the entire medium to be at a

state of incipient yielding. A modern example of this type of

formulation can be found in Cambou (1982).

A phenomenological continuum theory endeavors to devise

constitutive relations based on experimentally observed stress-strain

curves. It is presently the most popular class of the theories and it

concentrates on the macroscopically discernible stress and strain

measures. This theory does not inquire very deeply into the mechanisms

which control the process of deformation. A controversial assumption of

these phenomenological continuum theories, as applied to granular media,

is that the laboratory tests, such as the standard triaxial test,

achieve homogenous states of strain and stress. Many researchers are

now seeking the answer to the question of when bifurcation of the

deformation mode becomes acute enough to render interpretation of the

supposedly "homogenous state" data troublesome (see, for example, Lade,

1982, and Hettler et al., 1984).

Microstructural theories attempt to incorporate geometric measures

of local granular structure into the continuum theory. Local granular

structure is also called fabric, which is defined as the spatial

arrangement and contact areas of the solid granular particles and

associated voids. For clarity, fabric is subdivided into isotropic

fabric measures (such as porosity, density, etc.) and anisotropic fabric

measures (which are mentioned in the next section). In this

dissertation, unless otherwise stated, the word fabric refers to

anisotropic fabric. Perhaps the best known microstructural formulation

is that proposed by Nemat-Nasser and Mehrabadi (1984).

2.5 Anisotropic Fabric in Granular Material

2.5.1 Introduction

The fabric of earthen materials is intimately related to the

mechanical processes occurring during natural formation (or test sample

preparation) and the subsequent application of boundary forces and/or

displacements. Fabric evolution can be examined in terms of the

deformations that occur as a result of applied tractions (strain-induced

anisotropy), or the stresses which cause rearrangement of the

microstructure (stress-induced anisotropy). Strains are influenced to

some extent by the relative symmetry of the applied stress with respect

to the anisotropic fabric symmetry (or directional stiffness). If

straining continues to a relatively high level, it seems logical to

expect that the initial fabric will be wiped out and the intensity and

pattern of the induced fabric will align itself with the symmetry (or

principal) axes of stress. Before introducing and discussing a select

group of microscopic fabric measures, some of the commonly encountered

symmetry patterns, caused by combined kinematic/dynamic boundary

conditions, will be reviewed.

2.5.2 Common Symmetry Patterns

Triclinic symmetry implies that the medium possesses no plane or

axis of symmetry. This fabric pattern is produced by complex

deformations. Gerrard (1977) presents a simple example of how this most

general and least symmetric system may arise. Consider the sketch in

the upper left hand corner of Figure 2.5: triclinic symmetry develops as

a result of the simultaneous application of compression in direction 1,

differential restraint in directions 2 and 3, and shear stress

components acting in directions 2 and 3 on the plane having axis 1 as

its normal.

Monoclinic symmetry is characterized by a single plane of symmetry.

Any two directions symmetric with respect to this plane are equivalent.

An example of this symmetry group is shown in the lower left of Figure

2.5. The concurrent events leading to it are compression in direction

1, no deformation in the 2 and 3 directions, and a shear stress

component acting in the 2-direction and on the plane with axis 1 as its


A slight modification of the previous example permits a

demonstration of a case of n-fold axis symmetry or cross-anisotropy.

Exclusion of the shear stress component causes an axis of fabric

symmetry to develop such that all directions normal to this axis are

equivalent, bottom right of Figure 2.5.

The orthorhombic symmetry group can best be described by bringing

to mind the true triaxial device. Here for example (top right of Figure

2.5), three mutually perpendicular planes of symmetry are produced by

normal stresses of different magnitudes on the faces of the cubical sand



2 02

Materials = Materials
properties properties
(3) (2)

Figure 2.5 Common fabric symmetry types (after Gerrard, 1977)



n'- fold ox

Lastly, the rarest natural case is spherical symmetry or material

isotropy which implies that all directions in the material are

equivalent. However, because of its simplicity, isotropy is a major and

a very common simplifying assumption in many of the current

representations of soil behavior.

2.5.3 Fabric Measures

The selection of the internal variables, gn, to characterize the

mechanical state of a sand medium (see equation 2.3.7) has been a

provocative subject in recent times (Cowin and Satake, 1978; and Vermeer

and Luger, 1982). There is no doubt that the initial void ratio is the

most dominant geometric measure, but as Cowin (1978) poses: "Given that

porosity is the first measure of local granular structure or isotropicc]

fabric, what is the best second measure of local granular structure or

[anisotropic] fabric?" Trends suggest that the next generation of

constitutive models will include this second measure. It is therefore

worthwhile to review some of these variables.

An anthropomorphic approach is perhaps most congenial for

introducing the reader to the concept of anisotropic fabric in granular

material. Let us assume for illustrative purposes that, through a

detailed experimental investigation, we have identified a microscopic

geometric or physical measure (say variable X), which serves as the

secondary controlling factor to the void ratio in interpreting the

stress-strain response of sand. Some of the suggestions offered for the

variable X are 1) the spatial gradient of the void ratio ae (Goodman and

Cowin, 1972); 2) the orientation of the long axes of the grains (Parkin

et al., 1968); 3) the distribution of the magnitude and orientation of

the inter-particle contact forces (Cambou, 1982); 4) the distribution of

the inter-particle contact normals (see, for example, Oda, 1982); 5) the

distribution of branches [note: a branch is defined as the vector

connecting the centroids of neighboring particles, and it is thus

possible to replace a granular mass by a system of lines or branches

(Satake, 1978)]; 6) the mean projected solid path (Horne, 1964); and 7)

mathematical representations in the form of second order tensors

(Gudehus, 1968).

A commander (mother nature) of an army (the set representing the

internal variable of the sand medium) stations her troops (variable X)

in a configuration which provides maximum repulsive effort to an

invading force (boundary tractions). The highest concentration of

variable X will therefore tend to point in the direction of the imposed

major principal stress. If the invading army (boundary tractions)

withdraws (unloading), we should expect the general (mother nature) to

keep her distribution of soldiers (X) practically unaltered. It is an

experimental fact that there is always some strain recovery upon

unloading, and this rebound is caused partly by elastic energy stored

within individual particles as the soil was loaded and partly by

inelastic reverse sliding between particles (Figure 2.4).

Traditionally, it has been convenient to regard this unloading strain as

purely elastic, but in reality, it stems from microstructural changes

due to changes of the fabric and should be considered a dissipative

thermodynamically irreversible process (Nemat-Nasser, 1982). Returning

to our anthropomorphic description, we can therefore say that the

general (mother nature) has an intrinsic command to modify slightly the

arrangement of her troops (X) once the offensive army (boundary

tractions) decamps. The configuration of the defensive forces

(distribution of X) after complete or partial withdrawal of the

aggressor (complete or partial removal of the boundary loads) still,

however, reflects the intensity and direction of the earlier attack

(prior application of the system of boundary loads). This represents an

induced fabric or stress-induced anisotropy in the granular material.

We can create additional scenarios with our anthropomorphic model

to illustrate other features of fabric anisotropy. During the initial

placement of the forces (initial distribution of the variable X during

sample preparation or during natural formation of the soil deposit)

under the general's command, there is a bias in this arrangement which

is directly related to the general's personality (gravity as a law of

nature). This is the so-called inherent anisotropy (Casagrande and

Carillo, 1944) of soil which differs from the stress-induced anisotropy

mentioned previously. Say the invading army (boundary tractions)

attacks the defensive fortress (sand mass) with a uniform distribution

of troops (uniform distribution of stress vectors), we will expect

maximum penetration (strain) at the weakest defensive locations

(smallest concentration of X), but our rational general (mother nature)

should take corrective measures to prevent intrusion by the enemy forces

(boundary tractions) through the inherently vulnerable sites (points of

initially low X concentration). We can relate this situation to the

effect of increasing hydrostatic pressure on an inherently cross-

anisotropic sand specimen; the results of such a test carried out by

Parkin et al. (1968) indicate that the ratio of the incremental

horizontal strain to incremental vertical strain decreases from about 6

to 2.5. Increasing the hydrostatic pressure decreases the degree of

anisotropy, but it does not completely wipe out the inherent fabric. We

may infer that the general (mother nature) cannot reorient her forces at

will since she is faced by the annoying internal constraints (particles

obstructing each other) which plague most large and complex

organizations (the microscopic world of particles sliding and rolling

over each other).

It may seem logical to assume that if the demise of anisotropy is

inhibited in somd way, then so is its induction, but experimental

evidence reported by Oda et al. (1980) indicates that the principal

directions of fabric (i.e., principal directions of the distribution of

X or the second order tensor representation) match the principal

directions of the applied stress tensor during a virgin or prime

loading, even with continuous rotation of the principal stress axes.

There appears to be no lag effect. Data presented by Oda (1972)

describing the evolution of the contact normal distribution suggests

that fabric induction practically ceases once the material starts to

dilate. However, no firm conclusions can be drawn until many tests have

been repeated and verified by the soil mechanics community as a whole.

2.6 Elasticity

We now turn our attention to the mathematical models used to

simulate the stress-strain response of soil. In this section, the

essential features of the three types of elasticity-based stress-strain

relations are summarized (Eringen, 1962): 1) the Cauchy type, 2) the

Hyperelastic (or Green) type, and 3) the incremental (or Hypoelastic)

type. Although, in the strict sense, elastic implies fully recoverable

response, it is sometimes convenient to pretend that total deformations

are "elastic" and to disregard the elastic-plastic decomposition set

forth in equations 2.4.1 and 2.4.2. This approach has some practical

applications to generally monotonic outward loading paths. However, for

unload-reload paths, this class of formulation will fail to predict the

irrecoverable component of strain. Furthermore, one should not be

misled into believing that elasticity theory should be used exclusively

for predicting one-way loading paths because even in its most

complicated forms, elasticity theory may fail to predict critical

aspects of stress-strain behavior, many of which can be captured

elegantly in plasticity theory.

2.6.1 Cauchy Type Elasticity

A Cauchy elastic material is one in which the current state of

stress depends only on the current state of strain. Each stress

component is a single-valued function of the strain tensor,

ij fij ( kl) (

where f are nine elastic response functions of the material. Since

the stress tensor is symmetric, fkl fk and the number of these

independent functions reduces from nine to six. The choice of the

functions fj must also satisfy the Principle of Material Frame

Indifference previously mentioned in section 2.3; such functions are

called hemitropic functions of their arguments. The stress o is an

analytic isotropic function of E if and only if it can be expressed as

ai = mo 6i + m1 Eij + 02 im Emj' (

where (,, 1, and 02 are functions only of the three strain invariants

(see, for example, Eringen, 1962; p. 158).

For a first order Cauchy elastic model, the second order strain

terms vanish (2 0 = 0) and $( is a linear function of the first strain

invariant emm'

ij = (ao + a Cmm) 6j + a2 Eij. (

where a,, al, and a2 are response coefficients. At zero strain, ao 6ij

is the initial spherical stress. Higher order Cauchy elastic models can

be formulated by letting the response functions 00, 41, and (2 depend on

strain invariant polynomials of corresponding order. For example, the

second order Cauchy elastic material is constructed by selecting as the

response functions

40 = a, E + a2 (E )2 + a3 (1 E.
mm mm i3 ij

41 = a4 + as emm'


02 = a6,

where a,, a2,. .., a6 are material constants (Desai and Siriwardane,


An alternative interpretation of the first order Cauchy model is

presented in order to show the link between the elastic bulk and shear

moduli (K and G respectively) and Lame's constants (r and U). For this

material classification,

oij = Cijkl kl'
where the components of Cijkl are each a function of the strain

components, or if isotropy is assumed, the strain invariants. Since

both a.. and e are symmetric, the matrix Cijl is also symmetric in

"ij" and in "kl." A generalization of the second order tensor

transformation formula (equation to its fourth order analogue


C' =Q. Q. Q Q C (
ijkl = p jq kr is pqrs 6
as the transformation rule for the "elastic" stiffness tensor C. With

the isotropy assumption, the material response must be indifferent to

the orientation of the observer, and hence we must also insist that C be

equal to C'. A fourth order isotropic tensor which obeys this

transformation rule can be constructed from Kronecker deltas 6 (see, for

example, Synge and Schild, 1949, p.211); the most general of these is

ijkl = r 6ij 6kl + i 6ik + jk (

where r, p, and v are invariants. From the symmetry requirement,

Cijkl Cijlk (


r 6 j 6kl + 6ik 6jl + 6i 6jk

S6 ij 61k + 6il jk + v 6ik 6jl (

and collecting terms,

(P v) (6ik 6j 6 6j ) = 0, (

which implies that v v. With this equality, equation

simplifies to

Cjkl = r 6ij 61 + (6ik 6jl + 6il 6 ), (

where r and i are Lame's elastic constants.

The incremental form of the first-order, isotropic, elastic stress-

strain relation is therefore

doij r 6ij 6k1 + v (6ik 6j + 6l 6jk) ] dkl

= 6j demm + 2 deij. (

Multiplication of both sides of this equation by the Kronecker delta 6..
results in
results in

dkk = 3 r dm + 2 p dEmm,
kk mm mm


do kk3 dE = K = r + 2 p,
kk mm

where K is the elastic bulk modulus.

Substituting the identities

do = dsj + 1 dokk 6
ij ij kk ij


de.. = de + 1 dE, 6,
di deij kk ij
into equation results in

ds.. + 1 do 6 j = r 6 de + 2 P (de.. + 1 de 6ij),
i kk j ij mm 3 kk 1
3 3
and using equation in this expression shows that

dsij/2 deij = G = j, (

where G is the elastic shear modulus.

Combining equations and gives a more familiar

form of the isotropic, elastic stiffness tensor, namely

Cjkl = (K 2 G) 6 j 6 + G (6 6k + 6 6k ) (
ijkl ij 1 kl ik jl il jk
Many researchers have adapted this equation to simulate, on an

incremental basis, the non-linear response of soil; they have all

essentially made K and G functions of the stress or strain level. Some

of the better-known applications can be found in Clough and Woodward,

1967; Girijavallabhan and Reese, 1968; Kulhawy et al., 1969; and Duncan

and Chang, 1970.



2.6.2 Hyperelasticity or Green Type Elasticity

Green defined an elastic material as one for which a strain energy

function, W (or a complementary energy function, Q) exists (quoted from

Malvern, 1969, p. 282). The development of this theory was motivated by

a need to satisfy thermodynamic admissibility, a major drawback of the

Cauchy elastic formulation. Stresses or strains are computed from the

energy functions as follows:

ai. = W (

and conversely,

E.. a 82 (

For an initially isotropic material, the strain energy function, W,

can be written out in the form (see, for example, Eringen, 1962)

W = W(I~, 12, I) = Ao + A1 I + A2 2 +

A, 11 ~2 + As I3 + A7 I~ + A, IY 12 +

A9 Ii 13 + A10 Ij, (

where Ii, I2, and 1, are invariants of E,

I Ekk I = ij E ij i' 3 km Ekn mn'
2 3
and Ak (k =0,2,..,10) are material constants determined from curve

fitting. The stress components are obtained by partial differentiation,

oij = 3W a + + Wa a3 + W aI (
31, aj 3e ei 31, a

= 1 6ij + 2 ij + 03 Eim Emj' (

where Di (i = 1,2,3) are the response functions which must satisfy the

condition 3D /3I = 3$ /3.i in order to guarantee symmetry of the

predicted stress tensor.

Different orders of hyperelastic models can be devised based on the

powers of the independent variables retained in equation If,

for instance, we keep terms up to the third power, we obtain a second-

order hyperelastic law. These different orders can account for various

aspects of soil behavior; dilatancy, for instance, can be realistically

simulated by including the third term of equation Green's

method and Cauchy's method lead to the same form of the stress-strain

relationship if the material is assumed to be isotropic and the strains

are small, but the existence of the strain energy function in

hyperelasticity imposes certain restrictions on the choice of the

constitutive parameters. These are not pursued here, but the interested

reader can find an in-depth discussion of these constraints in Eringen

(1962). Also, detailed descriptions--including initialization

procedures--for various orders of hyperelastic models can be found in

Saleeb and Chen (1980), and Desai and Siriwardane (1984).

2.6.3 Hypoelasticity or Incremental Type Elasticity

This constitutive relation was introduced by Truesdell (1955) to

describe a class of materials for which the current state of stress

depends on the current state of strain and the history of the stress ot

(or the stress path). The incremental stress-strain relationship is

usually written in the form

do = f(o de), (

where f is a tensor valued function of the current stress a, and the

strain increment de. The principle of material frame indifference (or

objectivity) imposes a restriction on f: it must obey the transformation

Q f(o, de) Q = f(Q d Q Q o QT) (

for any rotation Q of the spatial reference frame. When f satisfies

this stipulation, it is, as mentioned in the previous section, a

hemitropic function of a and de. A hemitropic polynomial representation

of f is

do' = f(o, dE) = ao tr(de) 6 + al de + a2 tr(dE) a' +

a3 tr(o' de) 6 + 1 a, (de o' + a' de) + as tr(de) a'2 +

as tr(a' de) a' + a7 tr(a'2 de) 6 +

1 as (de o'2 + 012 dE) + aS tr(o' dE) o'2 +
a, tr(a'2 de) a' + a,1 tr(a'2 dE) ',2, (

where a' is the nondimensional stress o/2V (p being the Lame shear

modulus of equation, ak (k = 0,2,..,11) are the constitutive

constants (see, for example, Eringen, 1962, p.256), and "tr" denotes the

trace operator of a matrix (i.e., the sum of the diagonal terms). The

constants ak are usually dimensionless analytic functions of the three

invariants of a', and these are determined by fitting curves to.

experimental results.

Various grades of hypoelastic idealizations can be extracted from

equation This is accomplished by retaining up to and including

certain powers of the dimensionless stress tensor a'. A hypoelastic

body of grade zero is independent of a', and in this case, the general

form simplifies to

do' = f(g, de) = ao tr(de) 6 + a, de. (

Comparing this equation with the first order Cauchy elastic model

(equation shows that

ao = F and a, = 1.
2 u

Similarly, a hypoelastic constitutive equation of grade one can be

elicited from the general equation by keeping only the terms up to and

including the first power of o',

do' = f(a, de) = ao tr(de) 6 + ai de + a2 tr(dE) a' +

a3 tr(a' de) 6 + 1 a4 (de a' + a' de).

By a similar procedure, the description can be extended up to grade two,

with the penalty being the task of fitting a larger number of parameters

to the experimental data. These parameters must be determined from

representative laboratory tests using curve fitting and optimization

techniques, which often leads to uniqueness questions since it may be

possible to fit more than one set of parameters to a given data set.

Romano (1974) proposed the following special form of the general

hypoelastic equation to model the behavior of granular media:

doij = [a, d m+ a3 pq d pq] 6i + a dEij +
l odmm pq pq ij 13
[a2 de + a ar dE rsij. (
mr ds rs i
This particular choice ensures that the predicted stress increment is a

linear function of the strain increment; in other words, if the input

strain increment is doubled, then so is the output stress increment.

Imposing linearity of the incremental stress-strain relation is one way

of compelling the stress-strain relation to be rate-independent; a more

general procedure for specifying rate independence will be described in

the section on plasticity theory.

Davis and Mullenger (1978), working from Romano's equation, have

developed a model which can simulate many aspects of real soil behavior.

Essentially, they have used well-established empirical stress-strain

relations and merged them with concepts from plasticity to arrive at

restrictions on and the interdependency of the constitutive parameters.

2.7 Plasticity

Having outlined the theories used to compute the elastic, or

sometimes pseudo-elastic component dE of the total strain increment de,

the next topic deals with the computation of its plastic conjugate dE .

This section prefaces the mathematical theory of plasticity, a framework

for constitutive laws, which until 1952 (Drucker and Prager, 1952)

remained strictly in the domain of metals. Over the past three decades,

the role of elastic-plastic constitutive equations in soil mechanics has

grown in importance with the development of sophisticated computers and

computer-based numerical techniques. These tools have significantly

increased the geotechnical engineer's capacity to solve complicated

boundary value problems. The three main ingredients for these modern

solution techniques are computer hardware, numerical schemes, and

stress-strain equations, and, of these, the development of constitutive

laws for soils has lagged frustratingly behind.

The fundamentals of plasticity theory still remain a mystery to

many geotechnical engineers. It is very likely that a newcomer to this

field will find considerable difficulty in understanding the literature,

usually written in highly abstruse language. The chief objective of

this section is to provide some insight into plasticity theory by

highlighting the basic postulates, with special emphasis on their

applicability and applications to soil mechanics.

In brief, plasticity theory answers these questions:

a) When does a material plastically flow or yield? Or more directly,

how do we specify all possible stress states where plastic

deformation starts? The answer to this question lies in the

representation of these stress states by yield surfaces. Also

underlying this discussion are the definitions of and the possible

interpretations of yield.

b) Once the material reaches a yield stress state, how are the plastic

strains computed? And, if the stress path goes beyond the initial

yield surface (if an initial one is postulated), what happens to

the original yield surface (if anything)? The first question is

addressed in the discussion on the flow rule (or the incremental

plastic stress-strain relation), while the second is treated in the

discussion on hardening rules.

2.7.1 Yield Surface

Perhaps the best starting point for a discussion of plasticity is

to introduce, or rather draw attention to, the concept of a yield

surface in stress space. At the outset, it should be noted that yield

is a matter of definition, and only the conventional interpretations

will be mentioned in this chapter. The reader is, however, urged to

keep an open mind on this subject since a different perspective, within

the framework of a new theory for sands, will be proposed in Chapter 3.

Since strength of materials is a concept that is familiar to

geotechnical engineers, it is used here as the stimulus for the

introduction to yield surfaces. Figure 2.6 shows a variety of uniaxial

rate-insensitive stress-strain idealizations. In particular, Figures






(d) RIGID,

a a


Figure 2.6 Rate-independent idealizations of stress-strain

2.6 (d) and (e) show examples of perfectly plastic response, and one may

infer from this that, for homogenous stress fields, yield and failure

are equivalent concepts for this simplest idealization of plastic


In the calculation of the stability of earth structures, the Mohr-

Coulomb failure criterion is typically used to estimate the maximum

loads a structure can support. That is, when this load is reached, the

shear stress to normal stress ratio is assumed to be at its peak value

at all points within certain zones of failure. This method of analysis

is known as the limit equilibrium method. Using the classification set

forth in section 2.4, it is a kinematically ambiguous theory in that no

strains are predicted. Another common method of analysis is the wedge

analysis method. This is a trial and error procedure to find the

critical failure plane, a failure plane being a plane on which the full

strength of the material is mobilized and the critical plane being the

one that minimizes the magnitude of the imposed load.

A feature common to both the limiting equilibrium and the wedge

analysis methods is the need to provide a link between the shear and

normal stress at failure. A constitutive law, which is a manifestation

of the internal constitution of the material, provides this information.

More generally, the kinematically ambiguous theories for a perfectly

plastic solid must specify the coordinates of all possible failure

points in a nine dimensional stress space. Mathematically, this is

accomplished by writing a failure function or criterion in the form

F(oij) = 0; many well-established forms of the yield function are

previewed in the following.

The Mohr-Coulomb frictional failure criterion states that shear

strength increases linearly with increasing normal stress, Figure 2.7.

For states of stress below the failure or limit or yield line, the

material may be considered rigid [Fig. 2.6 (d)] or elastic [Fig.

2.6 (e)]. For a more general description, it is necessary to extend the

two-dimensional yield curve of Figure 2.7 to a nine-dimensional stress

space. Although such a space need not be regarded as having an actual

physical existence, it is an extremely valuable concept because the

language of geometry may be applied with reference to it (Synge and

Schild, 1949). The set of values Oal, 012, 013, 021, 022, 023, 031, 032

and 033 is called a point, and the variables oi. are the coordinates.

The totality of points corresponding to all values of say N coordinates

within certain ranges constitute a space of N dimensions denoted by VN.

Other terms commonly used for VN are hyperspace, manifold, or variety.

Inspection of, say, the equation of a sphere in rectangular

cartesian coordinates (x,y,z),

F(x,y,z) = (x a)2 + (y b)2 + (z c)2 k2 = 0

where a, b, and c are the center coordinates and k is the radius, is a

simple way of showing that the nine-dimensional equivalent of a

stationary surface in stress space may be expressed as

F(oij) = 0. (

A surface in four or more dimensions is called a hypersurface. The

theoretician must therefore postulate a mechanism of yield which leads

directly to the formulation of a yield surface in stress space or he

must fit a surface through observed yield points.

Rigorously speaking, a yield stress (or point) is a stress state

which marks the onset of plastic or irrecoverable strain and which may





(b) o0

(C) 0,+ 03

Figure 2.7 Two dimensional picture of Mohr-Coulomb failure

lie within the failure surface. Yield surfaces specify the coordinates

of the entirety of yield stress states. These (not necessarily closed)

surfaces bound a region in stress space where the material behavior is

elastic. But an all-important practical question still looms: How can

we tell exactly where plastic deformation begins? Is the transition

from elastic to elastic-plastic response distinct? At least for soils,

it is not that simple a task. The stress-strain curves continuously

turn, and plastic deformation probably occurs to some extent at all

stress states for outward loading paths. However, for the perfectly

plastic idealization, there should be no major difficulty since the

limit states are usually easy to identify.

Among the techniques used to locate the inception of yield are:

a) for materials like steel with a sharp yield point, the yield

stress is usually taken as the plateau in stress that occurs just

after the yield point;

b) for soft metals like aluminium, the yield stress is defined as the

stress corresponding to a small value of permanent strain (usually


c) a large offset definition may be chosen which more or less gives

the failure stress;

d) a tangent modulus definition may be used, but it must be

normalized if mean stress influences response; and

e) for materials like sand which apparently yield even at low stress

levels, a Taylor-Quinney (1931) definition is used. This and some

of the alternative definitions are illustrated in Figure 2.8.

q /




Figure 2.8 Commonly adopted techniques for locating the yield



( i
i /



Soil mechanicians will identify the Taylor-Quinney definition with

the Casagrande procedure (Casagrande, 1936) for estimating the

preconsolidation pressure of clays.

Defining a yield surface using the methods outlined above usually

leads to one with a shape similar to that of the failure or limit

surface. However, in Chapter 3, an alternative approach will be

suggested for determining the shape of the yield surface based on the

observed trajectory of the plastic strain increment--for sands, these

surfaces have shapes much different from the typical failure or yield


2.7.2 Failure Criteria

If an existing testing device had the capability to apply

simultaneously the six independent components of stress to a specimen,

the yield function F(aij) = 0 could be fitted to a comprehensive data

set. Unfortunately, such equipment is not available at present, and

most researchers still rely on the standard triaxial test (Bishop and

Henkel, 1962). However, if the material is assumed to be isotropic, as

is usually done, then the number of independent variables in the yield

function reduces from six to three; i.e., the three stress invariants or

three principal stresses replace the six independent components of a.

In other words, material directions are not important, only the

intensity of the stress is. Therefore, by ignoring anisotropy, all that

the theoretician needs is a device, like the cubical triaxial device,

which can vary a,, 02, and 03 independently.

Another implication of the isotropy assumption is that stress data

can be plotted in a three dimensional stress space with the principal

stresses as axes. This stress space is known as the Haigh-Westergaard

stress space (Hill, 1950). Working in this stress space has the

pleasant consequence of an intuitive geometric interpretation for a

special set of three independent stress invariants. In order to see

them, the rectangular coordinate reference system (oa, Oz, a3) must be

transformed to an equivalent cylindrical coordinate system (r, 6, z) as

described in the following.

Figure 2.9 depicts a yield surface in Haigh-Westergaard (or

principal) stress space. The hydrostatic axis is defined by the line

01 = 02 = 03,

which is identified with the axis of revolution (z). For cohesionless

soils (no tensile strength), the origin of stress space is also the

origin of this axis. A plane perpendicular to the hydrostatic axis

called a deviatoric or octahedral plane and is given by

01 + 02 + a3 = constant.

When this constant is equal to zero, the octahedral plane passes through

the origin of stress space and is then known as a i plane.

If we perform a constant pressure test (paths TC or TE of Figure

2.4), the stress point follows a curve on a fixed deviatoric plane for

the entire loading. Such stress paths provide a useful method for

probing the shape and/or size of the yield surface's r-plane projection

for different levels of mean stress. Polar coordinates (r, e) are used

to locate stress points on a given deviatoric plane.

By elementary vector operations, the polar coordinates r, 8, and z

can be correlated to each of the stress invariants /JJ, 6 and II, which

were previously defined in equations,, and

Yield Surface


von Mises

Hydrostatic Axis
-(-, = 7Z 0"3 )

Deviatoric Plane
( (I + O2+ 0-": Constant)

Hydrostatic Point
( = 02 = C3 )

Figure 2.9 Yield surface representation in Haigh-Westergaard
stress space

respectively. A measure of the shear stress intensity is given by the


r = /(2J2) (

from the hydrostatic point on the octahedral plane to the stress point.

The polar angle shown in Figure 2.9 is the same as the Lode angle

6. It provides a quantitative measure of the relative magnitude of the

intermediate principal stress (a,). For example,

a2 = a3 (compression tests) -+ = +300

aO = oa (extension tests) 6 = -300


oa + 03 = 2 a2 (torsion tests) e = 00.

Lastly, the average pressure, an important consideration for

frictional materials, is proportional to the perpendicular distance "d"

from the origin of stress space to the deviatoric plane;

d = Ii//3, (

where Ii is the first invariant of a.

For isotropic materials, the yield function (equation may

therefore be recast in an easily visualized form (Figure 2.9)

F(I1, /J,, 6) = 0. (

Some of the more popular failure/yield criteria for isotropic soils and

metals are reviewed in the following.

The much used Mohr-Coulomb failure criterion (Coulomb, 1773) for

soils is usually encountered in practice as

(01 a,) = sin 0 = k, (
(oa + 03)

where is a constant termed the angle of internal friction. The symbol

"k" is used as a generic parameter in this section to represent the size

of yield surfaces. This criterion asserts that plastic flow occurs when

the shear stress to normal stress ratio on a plane reaches a critical

maximum. If the equations which express the principal stresses in terms

of the stress invariants (equation are substituted into

equation, the Mohr-Coulomb criterion can be generalized to

(Shield, 1955)

F = II sin 4 + /J2 { sin e sin cos 6 } = 0. (
3 73
A trace of this locus on the r plane is shown in Figure 2.9. The

surface plots as an irregular hexagonal pyramid with its apex at the

origin of stress space for non-cohesive soils.

Also depicted in this figure are the well-known Tresca and Mises

yield surfaces used in metal plasticity. Mises (1928) postulated a

yield representation of the form

F = /J, k = 0, (

and physically, this criteria can be interpreted to mean that plastic

flow commences when the load-deformation process produces a critical

strain energy of distortion (i.e., strain energy neglecting the effects

of hydrostatic pressure and volume change).

Tresca (1864), on the other hand, hypothesized that a metal will

flow plastically when the maximum shear stress on any plane through the

point reaches a critical value. In the Mohr's circle stress

representation, the radius of the largest circle [(or 0,)/2] is the

maximum shear stress. Replacing the principal stresses with the stress

invariants gives the following alternative form for the Tresca


F = -1 /JJ [ sin (e + 4 i) sin (9 + 2 w) ] k = 0,
7 3 3
which, upon expansion of the trigonometric terms, simplifies to

F = /J2 cos 6 k = 0. (

A noticeable difference between the Mises or Tresca criterion and

the Mohr-Coulomb criterion is the absence of the variable It in the

former. This reminds us that yielding of metals is usually not

considered to be dependent on hydrostatic pressure, as the experiments

of Bridgman (1945) have demonstrated.

Drucker and Prager (1952) modified the Mises criterion to account

for pressure-sensitivity and proposed the form

F = /J2 k = 0. (

To match the Drucker-Prager and Mohr-Coulomb yield points in compression

space (o2 = 03), one must use

k = 2 sin 0 (
/3 (3 sin 0)
but, to obtain coincidence in extension space (oa = 02),

k = 2 sin 0 (
/3 (3 + sin 0)
must be specified. Although the development of the Drucker-Prager yield

function was motivated mainly by mathematical convenience, it has been

widely applied to soil and rock mechanics. However, there is

considerable evidence to indicate that the Mohr-Coulomb law provides a

better fit to experimental results (see, for example, Bishop, 1966).

Scrutiny of sketches of the previously defined yield surfaces in

principal stress space (see Figure 2.9) reveals that they are all "open"

along the hydrostatic stress axis. Therefore, for an isotropic

compression path, no plastic strains will be predicted. This

contradicts the typical behavior observed along such paths, Figure 2.2.

Recognizing this deficiency, Drucker et al. (1957) capped the Drucker-

Prager cone with a sphere to allow for plastic yielding for generally

outward but non-failure loading paths. The equation for the spherical

cap (of radius k) centered on the origin of stress space can be derived

by rearranging equation,

F( ) = "ijij k2 = I 2 I1 k2 0. (

As a result of the development of more sophisticated testing

devices, sensing equipment, and data capture units, more reliable and

reproducible stress-strain data is becoming available. This has quite

naturally led to the development of many new mathematical

representations of yielding in soils. Most notably, Lade and Duncan

(1975), using a comprehensive series of test data obtained from a true

triaxial device (Lade, 1973), have suggested that failure is most

accurately modeled by the function

F = (I3/I,) (Ii/Pa) k = 0, (

where 13 is the third stress invariant defined in equation, pa

is the atmospheric pressure in consistent units, and m is a constant to

model deviation from purely frictional response. A spherical cap was

subsequently added by Lade (1977) to "close" this "open-ended" function

along the hydrostatic axis.

Another recent proposal, based on a sliding model, was put forward

by Matsuoka and Nakai (1974). They defined the spatial mobilized plane

as the plane on which soil particles are most mobilized on the average

in three dimensional stress space. Only for special cases when any two

of the three principal stresses are equal does this criterion coincide

with the Mohr-Coulomb criterion. Based on the postulate that the

shear/normal stress ratio on the spatial mobilized plane governs

failure, Matsuoka and Nakai have derived the following failure


F = I[ 1, IZ 9 13] k = 0. (
9 I,

The mobilized plane concept is essentially a three-dimensional extension

of the Mohr-Coulomb criterion that takes into account the relative

weight of the intermediate principal stress.

Even more recently, Desai (1980) has shown that the Mises, Drucker-

Prager, Lade, and Matsuoka surfaces are all special cases of a general

third-order tensor invariant polynomial he proposed. Using statistical

analyses, he found that the failure criterion

F = [I2 + (Ii 131/3)] k = 0 (

gave the best fit to experimental data sets on Ottawa sand and an

artificial soil. Research in this field is presently very active, and

as more high quality data becomes available, it is anticipated that even

more proposals for failure/yield functions will emerge in the near


2.7.3 Incremental Plastic Stress-Strain Relation, and Prager's Theory

A material at yield signals the onset of plastic strain, and this

section describes the computation of the resulting plastic strain

increment. By definition, plasticity theory excludes any influence of

the rate of application of the stress increment on the predicted plastic

strain increment, and as will be shown later, this leads to restrictions

on the possible forms of the stress-strain relation.

In analogy to the flow lines and equipotential lines used in

seepage analysis, the existence of a plastic potential, G, in stress

space can be postulated such that (Mises, 1928)

de. = A G, A > 0 (
ii -a

where A is a scalar factor which controls the magnitude of the generated

plastic strain increment, and G is a surface in stress space (like the

yield surface) that dictates the direction of the plastic strain

increment. More specifically, the plastic strain increment is

perpendicular to the level surface G( ij) = 0 at the stress point.

To get a better grasp of equation, the soils engineer may

think of the function G as a fixed equipotential line in a flow net

problem. The partial derivatives 3G/aoij specify the coordinate

components of a vector pointing in the direction perpendicular to the

equipotential. This direction is, in fact, the direction of flow (along

a flow line) of a particle of water instantaneously at that spatial

point. Supplanting now the spatial coordinates (x,y,z) of the seepage

problem with stress axes (ax, ay oz), while keeping the potential and

flow lines in place, illustrates the mathematical connection between the

movement of a particle of water and the plastic deformation of a soil

element. The plastic geometrical change of a soil element is in a

direction perpendicular to the equipotential surface G(oij) = 0. At

different points in the flow problem, the particles of water move at

speeds governed by Darcy's law; therefore, it is possible to construct a

scalar point function which gives the speed at each location. In an

equivalent manner, the scalar multiplier A in equation

determines the speed (or equivalently, the magnitude of the incremental

deformation) of the soil particle at different locations in stress

space. For example, the closer the stress point is to the failure line,

a larger magnitude of A (with a corresponding larger magnitude of dE )

is expected. Therefore, in the crudest sense, the two elements of

plasticity theory which immediately confront us are: a) the

specification of the direction of the plastic strain increment through a

choice of the function G(oij), and b) the computation of the magnitude

of dE There are, of course, other important questions to be answered,

such as "What does the subsequent yield surface look like?", and these

will be treated in later sections and chapters.

Mises (1928) made the assumption that the yield surface and the

plastic potential coincide and proposed the stress-strain relation

de?. = A 3F (

This suggests a strong connection between the flow law and the yield

criterion. When this assumption is made, the flow rule (equation is said to be associated and equation is called the

normality rule. However, if we do not insist upon associating the

plastic potential with the yield function (as suggested by Melan, 1938),

the flow rule is termed non-associated. The implications of the

normality rule, it turns out, are far reaching, and as a first step to

an incisive understanding of them, Prager's (1949) treatment of the

incremental plastic stress-strain relation will be summarized.

The first assumption is designed to preclude the effects of rate of

loading, and it requires the constitutive equation

dEp = d-p (Y, do, 3n)

to be homogenous of degree one in the stress increment do. Recall that

homogeneity of order n ensures that

de = dp (t, A do, q) = An dp (ot, do, n), (

where A is a positive constant.

A simple example will help clarify this seemingly complex

mathematical statement. Suppose an axial stress increment of 1 psi

produced an axial plastic strain increment of .01 %; this means that if

A is equal to 2, and n = 1, the stress increment of 2 psi (A x 1 psi)

will predict a plastic strain increment of .02% (A x .01%). Ideally

then, the solution should be independent of the stress increment,

provided the stiffness change is negligible over the range of stress

spanned by the stress increment.

The simplest option, which ensures homogeneity of order one, is the

linear form

de.j = D do (
ij ijk1 kl'
where D is a fourth order plastic compliance tensor, the components of
t t
which may depend on the stress history o the strain history the

fabric parameters, etc., but not on the stress increment do. This is

referred to as the linearity assumption.

The second assumption, the condition of continuity, is intended to

eliminate the possibility of jump discontinuities in the stress-strain

curve as the stress state either penetrates the elastic domain (i.e.,

the yield hypersurface) from within or is unloaded from a plastic state

back into the elastic regime. To guarantee a smooth transition from

elastic to elastic-plastic response and vice-versa, a limiting stress

increment vector, do tangential to the exterior of the yield surface

must produce no plastic strain (note: the superscript "t" used here is