UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
PLASTICITY THEORY FOR GRANULAR MEDIA By DEVO SEEREERAM A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1986 ACKNOWLEDGEMENTS 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 inhouse 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. AFOSR840108 (M.C. McVay, Principal Investigator), which made this study possible. TABLE OF C'rlTENT'i . PAGE ACKNOWLEDGEMENTS ......................................................i LIST OF TABLES ....................................................... vii LIST OF FIGURES. ................ .................................... viii KEY TO SYMBOLS .......... ..............................................xvi ABT A T .. .. ......... ... .... .................. ........... .. ...... xx CHAPTER 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 PRELIMINARY AND FUNDAMENTAL CONCEPTS 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 StressStrain 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 StressStrain 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 ElastoPlastic StressStrain Relation*.102 3 FR'FPOSED PLASTICITY THEORY FOR GRANULAR MEDIA 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 A STUDY OF THE PREVOST EFFECTIVE STRESS MODEL 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 CHAPTER PAGE APPENDICES PAGE A DERIVATION OF ANALYTICAL REPRESENTATION OF DILATION PORTION OF YIELD SURFACE ..............................................275 B COMPUTATION OF THE GRADIENT TENSOR TO THE YIELD SURFACE........280 C EQUATIONS FOR UPDATING THE SIZE OF THE YIELD SURFACE...........283 D PREDICTIONS OF HOLLOW CYLINDER TESTS USING THE PROPOSED MODEL286 E PREDICTIONS OF HETTLER'S DATA USING THE PROPOSED MODEL.........297 F COMPUTATION OF THE BOUNDING SURFACE SCALAR MAPPING PARAMETER B ....................................................301 G PREDICTIONS OF HOLLOW CYLINDER TESTS USING PREVOST'S MODEL.....303 LIST OF REFERENCES ..*................................ .............312 BIOGRAPHICAL SKETCH ................................................ 325 LIST OF TABLES TABLE PAGE 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 ReidBedford 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 ReidBedford Sand*..198 3.9 Model Constants Used to Simulate Pressuremeter Tests*.......199 4.1 Prevost Model Parameters for ReidBedford 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 LIST OF FIGURES FIGURE PAGE 2.1 Representation of plane stress state at a "point"............20 2.2 Typical stressstrain 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 stressstrain 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 Rateindependent idealizations of stressstrain response.....63 2.7 Two dimensional picture of MohrCoulomb failure criterion**..66 2.8 Commonly adopted techniques for locating the yield stress**..68 2.9 Yield surface representation in HaighWestergaard 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 elasticplastico..................... .107 viii 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 elasticplastic 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 qp stress space ................................................ 115 3.5 Axial compression stressstrain data for Karlsruhe sand over a range of porosities and at a constant confinement pressure of 50 kN/m2 ........................................ 116 3.6 Stressstrain response for a cyclic axial compression test on loose Fuji River sand....................................117 3.7 Medium amplitude axial compressionextension 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 stressstrain 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 MohrCoulomb 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 FIGURE PAGE 3.15 Trace of the yield surface on the triaxial qp plane*.......135 3.16 Stress state in "thin" hollow cylinder......................152 3.17 Saada's hollow cylinder stress paths in Mohr's stress space155 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 density175 3.30 Measured and predicted response for axial compression test (o3 = 50 kN/m2) on Karlsruhe sand at 92.3% relative density176 3.31 Measured and predicted response for axial compression test (03 = 50 kN/m2) on Karlsruhe sand at 99.0% relative density177 FIGURE PAGE 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 stressstrain response for type "A" loading path on loose Fuji River sand....................................183 3.38 Observed stressstrain 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 compressionextension 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 FIGURE PAGE FIGURE 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 ReidBedford 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 PAGE 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 pq (top) and Cpq subspaces (bottom) ............................................... .....242 4.3 Measured vs. fitted stressstrain 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 stressstrain 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 stressstrain 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 stressstrain 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 stressstrain 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 stressstrain response for DCR 15 stress path using proposed model............................286 D.2 Measured vs. predicted stressstrain response for DCR 32 stress path using proposed model ...........................287 xiii FIGURE PAGE FICURE D.3 Measured vs. predicted stressstrain response for DTR 58 stress path using proposed model............................288 D.4 Measured vs. predicted stressstrain response for DTR 75 stress path using proposed model ******.....................289 D.5 Measured vs. predicted stressstrain response for GCR 15 stress path using proposed model..........................**290 D.6 Measured vs. predicted stressstrain response for GCR 32 stress path using proposed model........................... 291 D.7 Measured vs. predicted stressstrain response for R 45 (or pure torsion) stress path using proposed model.........*292 D.8 Measured vs. predicted stressstrain response for GTR 58 stress path using proposed model ...........................293 D.9 Measured vs. predicted stressstrain response for GTR 75 stress path using proposed model ........................... 294 D.10 Measured vs. predicted stressstrain 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 stressstrain response for DCR 15 stress path using Prevost's model .....*................... 303 G.2 Measured vs. predicted stressstrain response for DCR 32 stress path using Prevost's model ................*** *...... 304 G.3 Measured vs. predicted stressstrain response for DTR 58 stress path using Prevost's model*..........................305 PAGE FIGURE PAGE G.4 Measured vs. predicted stressstrain response for DTR 75 stress path using Prevost's model...........................306 G.5 Measured vs. predicted stressstrain response for GCR 15 stress path using Prevost's model*..........................307 G.6 Measured vs. predicted stressstrain response for GCR 32 stress path using Prevost's model...........................308 G.7 Measured vs. predicted stressstrain response for R 45 stress path using Prevost's model ......................... 309 G.8 Measured vs. predicted stressstrain response for GTR 58 stress path using Prevost's model ...........................310 G.9 Measured vs. predicted stressstrain response for GTR 75 stress path using Prevost's model ..........................311 KEY TO SYMBOLS C dge, dop e p de de dee, dE dEkk e p dekk, de kk' kk ds do D r e eo E f(o) F(o) F (a) G g(e) II, 12, I, (Ii)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 respectively 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 loading 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 u K plastic modulus P K plastic modulus at conjugate point o p (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 /J2I, 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) xvii Q r R s S XN z Y Y1, Y2, Y3 r 6 60 6  e p , E Ekk e p Ekk' kk a 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 envelope deviatoric components of a slope of dilation portion of * yield surface at the origin of /J2I. stress space slope of radial line in /JJII 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 option reload modulus parameters for proposed cyclic hardening option Lame's elastic constant distance from current stress state to conjugate stress state 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 xviii 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 SMohrCoulomb friction angle or stress obiliquity ^c MohrCoulomb friction angle observed in a compression test (i.e., one in which a0 = 03) Oe MohrCoulomb 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) kk 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 PLASTICITY THEORY FOR GRANULAR MEDIA By Devo Seereeram May 1986 Chairman: Dr. Frank C. Townsend Major Department: Civil Engineering A special timeindependent or elasticplastic 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 stressdependent plastic moduli. Specific choices are presented for several sands of different origin and initial density. Based on wellknown 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 boundingsurface 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 multisurface representation was chosen because it is thought of as one of the most fully developed of the existing soil constitutive theories. CHAPTER 1 I[Tii FrDLJ:CTION 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 loaddeformation 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). 2 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 stressstrain 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 stressstrain 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 computerbased 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 stateoftheart numerical solution techniques. An allencompassing stressstrain 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 coursesin 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. Ninetytwo 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 multiphase 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 stressstrain relations in favor of the simpler elastic and quasilinear theories. 1.3 Approach Using concepts recently advanced by Drucker and Seereeram (1986), a new stressstrain 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 workhardening, 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) pressuresensitive isotropic/kinematic hardening theory is also undertaken. This model was chosen because it is thought of as the most complete analytical statement on elastoplastic anisotropic 'Irdr.,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 wellknown 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 stressstrain 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 stressstrain or strength constants, such as friction angle and angle of dilation (Rowe, 1962), considered fundamental by most geotechnical engineers. A number of hollow cylinder and solid cylinder test paths are used to demonstrate the predictive capacity of the simple "nonhardening" 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 loadunload 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 viceversa. 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. Althchi 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. CHAPTER 2 FRELIMINARf AND FULiNDAHIETAL CONCEPTS 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 engineeraccustomed to the many empirical correlations and conventional plane strain, limit equilibrium methods of analysisventuring into the field of generalized, elastoplastic stressstrain relations. The terms "generalized" and "elastoplastic" 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 computerbased 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 presentday 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 vectora frame indifferent quantitycomputed 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 2.2.1.1 is completely expressed as i' Q k i (2.2.1.2) 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 2.2.1.2 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 (2.2.1.3) k k m m Substituting the inverse of equation 2.2.1.2 (i.e., i =Q i') k kr ~r into equation 2.2.1.3 leads to t Q i' = t' i' k krr r r or (t' tk Qkr) i' = 0, r k kr ' from which we see S= t kr (2.2.1.4) 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 tensors. 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 2.2.1.4. These requirements for first order tensors (or vectors) can begeneralized and extended for higher order tensors. 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 (n) v = v [n] = v.n. (2.2.1.5) 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 (n) vector, t is associated with each direction n. Two important examples of this type of tensorthe stress tensor and the strain tensorare 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 (n) 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 (n) 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 2.2.1.5 is t = ; [n], (2.2.2.1) 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] (2.2.2.2) for arbitrary unit vectors n, and n2, and a[an] = a G[n] (2.2.2.3) for arbitrary unit vector n and real number a. Geometrically, equation 2.2.2.2 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 2.2.2.3 means that if the length of the vector n is multiplied by a factor a, then so is the length of the (n) 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 (2.2.2.4) 31 032 033 or alternatively, in the indicial notation, (n) t = o. n. (2.2.2.5) 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 2.2.1.4). For a pure rotation of the basis, the transformation formula is derived by employing a sequence of previous equations. Recall from equation 2.2.1.4 that t' = t Q r k kr' and by combining this equation with equation 2.2.2.5, we find that tr = jk nj Qkr (2.2.2.6) r 3k jkr* Furthermore, n in this equation can be transformed to n' resulting in t' = jk Qjs n' Qr (2.2.2.7) r jk s s kr The left hand side of equation 2.2.2.7 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. (2.2.2.8) pr p jk js s kr All the indices in equation 2.2.2.8 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 2.2.2.8 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 s (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' (2.2.2.9) sr jk is kr' or in tensor notation, o' = Q 9 9 (2.2.2.10) 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 A X 21 (Ty x) '22 (y) Representation of plane stress state at a "point" 0*22 (y ) 2, (xy) ({^ C a0 (0x) 012 ( Txy) Lj '3  0/ '21 ( yx ) Figure 2.1 L2 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 nonzero 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 stressstrain 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 twodimensional 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 ydirections 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 twodimensional 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 MohrCoulomb 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 framedependent 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 (n) ( = A n, (2.2.2.11) or in indicial notation, (n) ti = A n (2.2.2.12) 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 2.2.2.12 into equation 2.2.2.5 and rearranging leads to aji nj A ni = O. (2.2.2.13) 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 = (2.2.2.14) ij 0 if i j By writing out the terms in long form, one may easily verify that n. = 6.. n.. (2.2.2.15) Equation 2.2.2.15 can now be substituted into equation 2.2.2.13 to give .ji nj A 6ij nj = 0, 31 *J 13 3 (o.. A 6..) n. = 0. (2.2.2.16) For clarity, equation 2.2.2.16 is expanded out to (ol1 A) nl + o12 n2 + 013 n3 = 0 021 n" + (022 A) n. + 023 n3 = 0, (2.2.2.17) a31 n, + 032 nz + (o33 A) n3 = 0 which may be organized in the matrix form o 1A 012 013 n1 0 021 o22A 023 n2 = 0 (2.2.2.18) 031 032 o33A 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, (2.2.2.19) since n is a unit vector. Equation 2.2.2.16 has a nontrivial solution if and only if the determinant of the coefficient matrix in equation 2.2.2.18 is equal to zero (see, for example, Wylie and Barrett, 1982, p.188). That is, o11A o12 013 021 o22A 23 = 0 (2.2.2. 031 032 033A must be true for nontrivial 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, (2.2.2.2 12 = (011o22 + 022033 + 033011) + 0 + 31 + 02 = ( o ao I2 ) + 2, (2.2.2.2 13 ij 20) 1) ?2) 3) and 011 012 013 13 = 021 022 023 3 (2.2.2.24) 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 13must also be independent of the coordinate system. These are therefore invariant with respect to changes in the perspective of the observer and are the socalled 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 offdiagonal 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 follows: J1 = s1 + 322 + s33 = 1 22 1 + 33 I, Kk 3 3 3 and by recalling equation 2.2.2.22, it is clear that J, = 0. (2.2.2.25) From the last equation and equation 2.2.2.23, observe that the second invariant of the stress deviation (denoted by J2) is simply J2 = (s..s. ) + 2. (2.2.2.26) Denoting the third invariant of the stress deviation by J3, the cubic expression for the stress deviator s, in analogy to equation 2.2.2.21 for the stress tensor o, becomes A3 J2 A J3 = 0, (2.2.2.27) 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 2.2.2.27 is considerably easier than that of equation 2.2.2.21. It is therefore more convenient to solve for the principal values of s and then compute the principal values of a using the identities oi = si + p, 02 = s2 + p, and 03 = S3 + p. (2.2.2.28) The direct evaluation of the roots, A, of equation 2.2.2.27 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 choice, sin3e 3 sine + 1 sin 36 = 0. (2.2.2.29) Replacing A with r sine in equation 2.2.2.27 gives r3 sin3e J, r sine J3 = 0, which when divided through by r3 gives sin3e J_ sine J_ = 0. (2.2.2.30) r2 r3 A direct correlation of this equation with equation 2.2.2.29 shows that o r or r = 2 /J2, J3 = 1 sin 30, sin 3 = 4 sin 30 = 4 J3. (2.2.2.31) (2.2.2.32) Substitution of the negative root of equation 2.2.2.31 into equation 2.2.2.32 leads to sin 36 = [3/3 (J3//J23)], (2.2.2.33) 2 from which we find that e = 1 sin1 [3/3 (J3/JJz3)], (2.2.2.34) 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 2.2.2.33. 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), (2.2.2.35) s, = 2 /J2 sin(e), (2.2.2.36) 75 and, S3 = 2 /J2 sin(0 + 2 1). (2.2.2. 73 7 Finally, these relations can be combined with those of equation 2.2.2.28 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 I 6 = sin1 [oi + a3 2 02], 300 < 0 < 300. (2.2.2. 2 /(3 J2) 37) 38) 39) 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 (1949). 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, (2.2.3.1) 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)] (2.2.3.2) 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), (2.2.3.3) and x'(Q) = x (Q) + u (Q). (2.2.3.4) 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 deforminga rigid body motion. From equations 2.2.3.3 and 2.2.3.4, 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)], (2.2.3.5) and subtracting equation 2.2.3.2 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)]. (2.2.3.6) 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)]. (2.2.3.7) Substitution of this equation into equation 2.2.3.6 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)]. (2.2.3.8) Furthermore, we know approximately that [xr(Q) Xr(P)] = Lo nr, (2.2.3.9) where n are the components of the unit vector directed from P to Q; substitution of this relation into equation 2.2.3.8 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 (2.2.3.10) L1 Lo = [aur/3xslp ns u r/9xtP n + L2 2 n au r/ax nm]. (2.2.3.11) 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. (2.2.3.12) L2 Moreover, Lj L2 = L, Lo L, + Lo LL Lo Lo and with the that = L Lo L, Lo + 2 L9 Lo Lo = Li Lo [L1 Lo + 2] Lo Lo e ( e + 2 ), (2.2.3.13) assumption of small strain, e2 is negligible, which implies L L = 2 e. L20 (2.2.3.14) By equating the previous equation with equation 2.2.3.12, one finds that e = nr u r/ax I n (2.2.3.15) 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 ], (2.2.3.16) 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 (2.2.3.17) 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 offdiagonal 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' (2.2.3.18) 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 .). (2.2.3.19) 1J J1 2.3 StressStrain 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 loaddeformation response (or stressstrain 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 loaddeformation characterization of the solid media is usually expressed in the form of a constitutive law relating the forcetype 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 functionalit 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 stressstrain 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 nonlinear 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" stressstrain 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 rateindependent 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 stressstrain 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 weightperhaps 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 t 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 t 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 variablesexcept for a e and doas 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, rateindependent, incremental stressstrain 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 loaddeformation 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 socalled axioms or principles of constitutive theory. An axiom is a wellestablished 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 sequel. The Axiom of Causality states that the motion of the material points of a body is to be considered a selfevident, observable effect in the mechanical behavior of the body. Any remaining quantities (such as the stress) that enter the entropy production and the balance equationsi.e., the equations of conservation of mass, balance of momentum, and conservation of energyare 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 forminvariant 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 inequality. 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 loaddeformation response of loose and dense soil media subject to two conventional laboratory stress paths: hydrostatic compression, and conventional w I Z I u u Z u. F E3 0s 1 DENSE SAND OR OVERCONSOLIDATED CLAY 2 LOOSE SAND OR NORMALLY CONSOLIDATED CLAY 1 ^  2 o 2 PRINCIPAL STRAIN DIFFERENCE, r3 PRINCIPAL STRAIN DIFFERENCE, 1I 3 Figure 2.2 PLASTIC STRAIN [ AMOUNT OF COMPACTION ] VOLUMETRIC STRAIN,  = '6k V0 Typical stressstrain 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 stressstrain 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 qy Ox RTC Figure 2.3 Typical stress paths used to investigate the stress strain behavior of soil specimens in the triaxial environment Standard NAME OF TEST Designation DESCRIPTION Conventional Triaxial CTC Acrx = AOz = 0 A cy > 0 Compression 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>Aox (=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 X=0,z 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 e 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 simplifications. Three broad classes of continuum theories have evolved in the development and advancement of soil stressstrain models (Cowin, 1978): A C O irrev rev A elastic unloading O B B Ee SE rev E rev 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 stressstrain 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 NematNasser 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 (straininduced anisotropy), or the stresses which cause rearrangement of the microstructure (stressinduced 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 2direction and on the plane with axis 1 as its normal. A slight modification of the previous example permits a demonstration of a case of nfold axis symmetry or crossanisotropy. 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 specimen. Orthorhombic 0* 2 02 Materials = Materials properties properties (3) (2) Figure 2.5 Common fabric symmetry types (after Gerrard, 1977) Triclinic oi Monoclinic 02 Q2 n' fold ox symmetry 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 stressstrain 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 interparticle contact forces (Cambou, 1982); 4) the distribution of the interparticle 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 (NematNasser, 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 stressinduced 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 socalled inherent anisotropy (Casagrande and Carillo, 1944) of soil which differs from the stressinduced 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 stressstrain response of soil. In this section, the essential features of the three types of elasticitybased stressstrain 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 elasticplastic 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 unloadreload 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 oneway loading paths because even in its most complicated forms, elasticity theory may fail to predict critical aspects of stressstrain 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 singlevalued function of the strain tensor, ij fij ( kl) (2.6.1.1) 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' (2.6.1.2) 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. (2.6.1.3) 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' and 02 = a6, where a,, a2,. .., a6 are material constants (Desai and Siriwardane, 1984). 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 2.2.2.9) to its fourth order analogue produces C' =Q. Q. Q Q C (2.6.1.4) 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 (2.6.1.5) where r, p, and v are invariants. From the symmetry requirement, Cijkl Cijlk (2.6.1.6) or r 6 j 6kl + 6ik 6jl + 6i 6jk S6 ij 61k + 6il jk + v 6ik 6jl (2.6.1.7) and collecting terms, (P v) (6ik 6j 6 6j ) = 0, (2.6.1.8) which implies that v v. With this equality, equation 2.6.1.5 simplifies to Cjkl = r 6ij 61 + (6ik 6jl + 6il 6 ), (2.6.1.9) where r and i are Lame's elastic constants. The incremental form of the firstorder, isotropic, elastic stress strain relation is therefore doij r 6ij 6k1 + v (6ik 6j + 6l 6jk) ] dkl = 6j demm + 2 deij. (2.6.1.10) 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 (2.6.1.11) 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 3 (2.6.1.12) de.. = de + 1 dE, 6, di deij kk ij 3 into equation 2.6.1.10 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 2.6.1.11 in this expression shows that dsij/2 deij = G = j, (2.6.1.13: where G is the elastic shear modulus. Combining equations 2.6.1.12 and 2.6.1.13 gives a more familiar form of the isotropic, elastic stiffness tensor, namely Cjkl = (K 2 G) 6 j 6 + G (6 6k + 6 6k ) (2.6.1.14 ijkl ij 1 kl ik jl il jk 3 Many researchers have adapted this equation to simulate, on an incremental basis, the nonlinear response of soil; they have all essentially made K and G functions of the stress or strain level. Some of the betterknown 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 (2.6.2.1) 13 3Eij and conversely, E.. a 82 (2.6.2.2) 30 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, (2.6.2.3) 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 (2.6.2.4) 31, aj 3e ei 31, a = 1 6ij + 2 ij + 03 Eim Emj' (2.6.2.5) 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 2.6.2.3. 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 2.6.2.3. Green's method and Cauchy's method lead to the same form of the stressstrain 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 indepth discussion of these constraints in Eringen (1962). Also, detailed descriptionsincluding initialization proceduresfor 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 stressstrain relationship is usually written in the form do = f(o de), (2.6.3.1) 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) (2.6.3.2) 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 + 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 + 2 a, tr(a'2 de) a' + a,1 tr(a'2 dE) ',2, (2.6.3.3) where a' is the nondimensional stress o/2V (p being the Lame shear modulus of equation 2.6.1.10), 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 2.6.3.3. 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. (2.6.3.4) Comparing this equation with the first order Cauchy elastic model (equation 2.6.1.10) 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). 2 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. (2.6.3.5) 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 stressstrain relation is one way of compelling the stressstrain relation to be rateindependent; 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 wellestablished empirical stressstrain 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 pseudoelastic 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 elasticplastic constitutive equations in soil mechanics has grown in importance with the development of sophisticated computers and computerbased 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 stressstrain 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 stressstrain 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 rateinsensitive stressstrain idealizations. In particular, Figures (a) NONLINEARLY ELASTIC (C) NONELASTIC, OR PLASTIC (e) ELASTIC, PERFECTLY PLASTIC (b) LINEARLY ELASTIC 0" (d) RIGID, PERFECTLY PLASTIC a a (f) RIGID, (g) ELASTIC, WORK WORK HARDENING HARDENING Figure 2.6 Rateindependent idealizations of stressstrain response 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 response. 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 wellestablished forms of the yield function are previewed in the following. The MohrCoulomb 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 twodimensional yield curve of Figure 2.7 to a ninedimensional 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 ninedimensional equivalent of a stationary surface in stress space may be expressed as F(oij) = 0. (2.7.1.1) 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 DIRECTION OF 0, (a) YIELD LIMIT SURFACE = LINE ELASTIC OR RIGID REGION (b) o0 (C) 0,+ 03 2 Figure 2.7 Two dimensional picture of MohrCoulomb failure criterion 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 allimportant practical question still looms: How can we tell exactly where plastic deformation begins? Is the transition from elastic to elasticplastic response distinct? At least for soils, it is not that simple a task. The stressstrain 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 0.2%); 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 TaylorQuinney (1931) definition is used. This and some of the alternative definitions are illustrated in Figure 2.8. q / / / YIELD DEFINITIONS: / I MODERATE OFFSET / 2 TANGENT MODULUS = 3 TAYLORQUINNEY / 4 LARGE OFFSETa FAIL URE URE SHEAR STRAIN,z Figure 2.8 Commonly adopted techniques for locating the yield stress / / I / S/ I ( i i / I 4 L RGE OFFSEi, FAIL r Soil mechanicians will identify the TaylorQuinney 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 incrementfor sands, these surfaces have shapes much different from the typical failure or yield surfaces. 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 HaighWestergaard 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 HaighWestergaard (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 rplane 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 2.2.2.26, 2.2.2.39, and 2.2.2.22 Yield Surface Tresca von Mises Hydrostatic Axis (, = 7Z 0"3 ) Deviatoric Plane ( (I + O2+ 0": Constant) Hydrostatic Point ( = 02 = C3 ) Figure 2.9 Yield surface representation in HaighWestergaard stress space respectively. A measure of the shear stress intensity is given by the radius r = /(2J2) (2.7.2.1) 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 and 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, (2.7.2.2) where Ii is the first invariant of a. For isotropic materials, the yield function (equation 2.7.1.1) may therefore be recast in an easily visualized form (Figure 2.9) F(I1, /J,, 6) = 0. (2.7.2.3) Some of the more popular failure/yield criteria for isotropic soils and metals are reviewed in the following. The much used MohrCoulomb failure criterion (Coulomb, 1773) for soils is usually encountered in practice as (01 a,) = sin 0 = k, (2.7.2.4) (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 2.2.2.38) are substituted into equation 2.7.2.4, the MohrCoulomb criterion can be generalized to (Shield, 1955) F = II sin 4 + /J2 { sin e sin cos 6 } = 0. (2.7.2.5) 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 noncohesive soils. Also depicted in this figure are the wellknown Tresca and Mises yield surfaces used in metal plasticity. Mises (1928) postulated a yield representation of the form F = /J, k = 0, (2.7.2.6) and physically, this criteria can be interpreted to mean that plastic flow commences when the loaddeformation 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 criterion: 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. (2.7.2.6) A noticeable difference between the Mises or Tresca criterion and the MohrCoulomb 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 pressuresensitivity and proposed the form F = /J2 k = 0. (2.7.2.8) Ii To match the DruckerPrager and MohrCoulomb yield points in compression space (o2 = 03), one must use k = 2 sin 0 (2.7.2.9) /3 (3 sin 0) but, to obtain coincidence in extension space (oa = 02), k = 2 sin 0 (2.7.2.10) /3 (3 + sin 0) must be specified. Although the development of the DruckerPrager 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 MohrCoulomb 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 nonfailure loading paths. The equation for the spherical cap (of radius k) centered on the origin of stress space can be derived by rearranging equation 2.2.2.23, F( ) = "ijij k2 = I 2 I1 k2 0. (2.7.2.11) As a result of the development of more sophisticated testing devices, sensing equipment, and data capture units, more reliable and reproducible stressstrain 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, (2.7.2.12) where 13 is the third stress invariant defined in equation 2.2.2.24, 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 "openended" 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 MohrCoulomb 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 criterion: F = I[ 1, IZ 9 13] k = 0. (2.7.2.13) 9 I, The mobilized plane concept is essentially a threedimensional extension of the MohrCoulomb 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 thirdorder tensor invariant polynomial he proposed. Using statistical analyses, he found that the failure criterion F = [I2 + (Ii 131/3)] k = 0 (2.7.2.14) 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 future. 2.7.3 Incremental Plastic StressStrain 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 stressstrain 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 (2.7.3.1) 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 2.7.3.1, 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 2.7.3.1 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 stressstrain relation de?. = A 3F (2.7.3.2) Saij This suggests a strong connection between the flow law and the yield criterion. When this assumption is made, the flow rule (equation 2.7.3.1) is said to be associated and equation 2.7.3.2 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 nonassociated. 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 stressstrain relation will be summarized. The first assumption is designed to preclude the effects of rate of loading, and it requires the constitutive equation dEp = dp (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), (2.7.3.3) 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 (2.7.3.4) 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 stressstrain 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 elasticplastic response and viceversa, 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 