FINITE ELEMENT ANALYSIS OF CREEP PROBLEMS IN SOIL MECHANICS by JOHN J . EMERY B.A.Sc, U n i v e r s i t y o f B r i t i s h Columbia, 1966 A THESIS SUMBITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n the Department o f C i v i l E n gineering We accept t h i s t h e s i s as conforming to the r e q u i r e d standard: THE UNIVERSITY OF BRITISH COLUMBIA June, 1971 In present ing th i s thes is in pa r t i a l f u l f i lmen t o f the requirements for an advanced degree at the Un ive rs i t y of B r i t i s h Columbia, I agree that the L ib ra ry sha l l make i t f r ee l y ava i l ab le for reference and study. I fu r ther agree that permission for extensive copying o f th i s thes is fo r s cho la r l y purposes may be granted by the Head of my Department or by h is representa t i ves . It is understood that copying or pub l i c a t i on o f th i s thes is f o r f i nanc i a l gain sha l l not be allowed without my wr i t ten "permission. Department of Civ i l Engineering The Un ivers i t y of B r i t i s h Columbia Vancouver 8, Canada Date July 5, 1971 i i ABSTRACT Application of the f in i te element method of stress analysis to problems in soi l mechanics has enabled the engineer to gain much information on the deformations and stresses in earth structures under the assumption that the properties of soi l are independent of time. However, both laboratory studies and f i e ld observations indicate that time is a very important factor in the behaviour of cohesive so i l s . The f in i te element method is extended to deal with typical problems in soi l mechanics in which the time dependence of the mechanical pro-perties of soi l is considered. From creep studies reported in the l i terature , stress-strain-time relationships that have been developed for essential ly uniaxial constant stress creep tests are extended to the multiaxial changing stress condition that is applicable in s i tu . The f in i te element method is used to examine problems where the soi l is assumed to be linear v iscoelast ic . The correspondence rule of l inear v iscoelast ic i ty is used as part of this development. This is an idealized case and wi l l generally be limited to qualitative studies. Then, the incremental i n i t i a l strain f in i te element method is developed to deal with general soi l creep as described by empirical relationships. A cumulative creep law based on the strain-hardening rule is adopted in this .analysis. Creep rupture is introduced into the general incre-mental solution method by reducing the st iffness of elements of the earth, structure that have fa i led . The fai lure criterion used in the analysis is based on the total elapsed time from the beginning of creep and the creep strain rate. Some typical examples of problems in soi l mechanics are examined to i l lust rate the methods developed. A C K N O W L E D G E M E N T S In presenting this thesis, the author wishes to record his indebtedness to the Shell Oil Company of Canada for a Shell Centennial Fellowship that provided financial support for much of this work, to the National Research Council of Canada for a Bursary during the f i r s t year of graduate study, and to the University of Br i t ish Columbia for providing the excellent computing f a c i l i t i e s that made much of the work possible. j The author would also l ike to thank the following individuals: - His supervisor, W. D, Liam Finn, for his continued encouragement and guidance at every stage of the study. - R. G. Campanella and D. L. Anderson for their suggestions and advice. - His colleagues in the Soil Mechanics Graduate Laboratory, both past and present, for friendly crit icisms and discussions. - Ed Goronsay for his generous assistance in preparing the final drawings. - Dinah Lindsay for preparing the manuscripts and final thesis and for her interest in the work. And f ina l l y , and personally most important, for their under-standing and encouragement, this thesis is dedicated to my parents. V TABLE OF CONTENTS Pa^e 1 IMPORTANCE OF CREEP AND CREEP RUPTURE IN SOIL MECHANICS 1 1-A Introduction 1 1-B Creep and Creep Rupture 1 1-C Examples of Creep and Creep Rupture 2 1-D Factors Complicating the Analysis of Soil Creep 4 1- E Purpose and Scope 5 2 DEVELOPMENT OF CREEP AND CREEP RUPTURE RELATIONSHIPS 7 2- A Introduction 7 2-B Typical Creep Behaviour of Cohesive Soils Observed in Laboratory Tests 7 2-C Methods of Representing the Time-Dependent Deformation of Cohesive Soils 11 2-D Extension of the Uniaxial Relationships to the Multiaxial Case 22 2-E Creep Rupture Cr i ter ia 35 3 LINEAR VISCOELASTIC STRESS ANALYSIS USING FINITE ELEMENTS AND THE CORRESPONDENCE RULE 46 3-A Introduction 46 vt Page 3-B Correspondence Rule 47 3-C Material Characterization 48 3-D Linear Viscoelastic Example Problem 51 3-E Finite Element Method 54 3-F Methods of Representing {Q(t)}n 57 3-G Evaluation of the Required Transforms 61 3-H Evaluation of Time-Dependent Stresses 62 3- 1 Accuracy Check on the Finite Element Iterative Method 63 4 NONLINEAR CREEP ANALYSIS USING FINITE ELEMENTS AND THE INCREMENTAL INITIAL STRAIN METHOD 72 4- A Introduction 72 4-B Finite Element Solution Using the Incremental Init ial Strain Method 76 4- C Accuracy Check on the Finite Element Solution of Nonlinear Creep Problems 79 5 TYPICAL EXAMPLE PROBLEMS 86 5- A Introduction 86 5-B Problems Involving Soils Considered to be Linear Viscoelastic 86 5-C General Creep Problems 108 5-D Creep Rupture 117 6 CONCLUSIONS BIBLIOGRAPHY APPENDIX A' - BRIEF DEVELOPMENT OF THF, FINITE ELEMENT METHOD APPENDIX B - OPTIONS AND FEATURES OF THE PROGRAMS B-l Options and Features of the Axisymmetric and Plane Strain-Plane Stress Finite Element Programs B-2 Plane Strain-Plane Stress Program Input and Capabilities APPENDIX C - EQUATIONS REQUIRED FOR THE DETERMINATION OF THE COMPONENTS OF THE CREEP STRAIN INCREMENT LIST OF TABLES Table 5-1 Parameters for Rheological Models Used in Figure 5-1 5-2 Parameters for Rheological Models Used in Figures 5-2, 5-3, and 5-10 5-3 Parameters for Rheological Models Used in Figures 5-5, 5-6 and 5-7 5-4 Parameters for Rheological Models Used in Figure 5-8 (Shale) 5-5 Creep Rupture for Axisymmetric Problem LIST OF FIGURES Typical Creep Behavior of Soils Observed in Laboratory Tests Strain Rate Equation Parameters Proposed Rheological Models of Soils Common Rheological Models Cumulative Creep Rules (Steel at High Temperatures) Schematic Diagram I l lustrat ing the Relation Between Types of Tests When the Effect of Previous Strain Rate is Negligible (Steel at High Temperatures) Axial Strain Rate-Axial Strain Curves for Triaxial Tests on Normally Consoli-dated Haney Clay Axial Strain Rate-Time Results for Un-drained Triaxial Creep Tests on Normally Consolidated Haney Clay Line of Transient Minimum Strain Rates for Undrained Triaxial Creep Tests on Haney Clay Total Rupture Life for Laboratory Creep Tests on Cohesive Soils Time to Rupture-Strain Rate Curve for Haney Clay Time to Rupture-Strain Rate Relationship Example Problem—Thick-Walled Cylinder Piecewise Linear Approximations for u(t) Finite Element Meshes for the Thick-Walled Cylinder Comparison of Iterative Finite Element Solution with Closed-Form Solution. Thick-Walled Cylinder, Maxwell Model in Shear, Elast ic in Dilatation Comparison of Iterative Finite Element Solution with Closed-Form Solution. Thick-Walled Cylinder, Standard Linear Solid in Shear, E last ic in Dilatation Comparison of Iterative Finite Element Solution with Closed-Form Solution. Thick-Walled Cylinder, Burgers Body in Shear, Elast ic in Dilatation Approximation of a Smooth Stress-Time Curve by a Series of Steps Schematic Diagrams Showing How to Deter*-mine the Increment of Creep Strain for Changing Stress Conditions Plane Strain Quadrilateral Finite Element Idealization for the Thick-Walled Cylinder Comparison of Radial Stresses in the Thick-Walled Cylinder for Closed-Form and Finite Element Solutions of the Creep Problem Comparison of Radial Creep Deforma-tions for the Thick-Walled Cylinder for Various Finite Element Solutions Uniform Continuous Loading on a Laterally Constrained Soil Layer Circular Surcharge Loading on a Soil Layer Surface Displacements with Time for a Circular Surcharge Loading on a Soil Layer Finite Element Idealization of a Long Slope Horizontal Creep Displacement and Velocity for the Long Slope in Figure 5-4 Displacements with Time at the End of a Lonq Slope with a Gravity Loading, Clay " Horizontal Force with Time on a Rigid Retaining Wall from Node 1 to Node 5 of Figure 5-6 Displacements with Time at the End of a Long Slope with a Gravity Load-ing, Shale Finite Element Idealization of a 45° Slope Edge Displacements with Time for a 45" Slope with a Distributed Load Horizontal Creep Displacement and Velocity for a 'Glacier ' Represented by the Long Slope in Figure 5-4 Edge Displacements with Time for a 45° Excavated Slope. Creep Given by Singh and Mitchel l 's Empirical Relationship Finite Element Idealization for Exca-vation of a Slurry Trench Edge Displacements with Time after Excavating a Slurry Trench. Creep Given by Singh and Mitchel l 's Empirical Relationship Triangular Finite Element Division of a Plane Stress or Plane Strain Region Quadrilateral Finite Element Made up from Four Triangular Elements Finite Element Idealization for Axisym-metric Problems Simulation of Slope Excavation 1 CHAPTER 1 IMPORTANCE OF CREEP AND CREEP RUPTURE IN SOIL MECHANICS 1-A INTRODUCTION Limit equilibrium methods are generally used for the s tab i l i t y analysis of earth structures such as embankments and natural or excavated slopes. These l imit equilibrium methods do not take into account the actual deformation and stress behaviour throughout the earth structure. Further, and s t i l l more important in many problems, the time dependency of the deformations, stresses and shear strength is not considered at a l l . However, both laboratory tests and f i e ld observation indicate that time is a very important factor in the behaviour of s o i l . Application of the f in i te element method of stress analysis to problems in soi l mechanics (62) has enabled the engineer to gain much information on the deformations and stresses in earth structures under the assumption that the properties of soi l are independent of time. However, this method has not been extended to creep and creep rupture problems in so i l s . In this thesis the f in i te element method wil l be extended to deal with problems in soi l mechanics in which the time dependence of the mechanical properties of soi l is considered. The soi l properties used wi l l be adopted from a survey of lab-oratory and f i e ld observations of creep and creep rupture. 1-B CREEP AND CREEP RUPTURE Creep in soi l mechanics is generally considered to be the slow, more or less continuous deformation or flow of soi l which takes place under constant stresses that result from gravity or external loadings (30, 78). Most writers (86) have assumed that soi l creep occurs at constant volume. Such a restr icted 2 definit ion of soi l creep wil l not be adopted herein, and creep wi l l be con^ sidered to be the total time-dependent deformation which includes both the distortional and the volumetric components. Soil creep may eventually lead to rupture and consequently to embankment fai lures or landslides. Such creep rupture is usually defined as the fai lure of the soi l under a shear stress that is less than the shear strength measured in standard laboratory tests. While the importance of surface creep (26) is recognized, i t is the deep-seated continuous creep that can lead to excessive deformation or rupture that is considered here. 1-C EXAMPLES OF CREEP AND CREEP RUPTURE The importance of creep and creep rupture in earth structures is best i l lustrated by examples from the l i terature. Haefeli (30), in an excellent review of the creep process. in_so i l s , snow and ice , describes the fai lure of an arch bridge between two creeping embankments. In this case, the remedial design required an estimate of the pressure due to creep on the bridge so that the size of a horizontal brace between the embankments could be determined. After the bridge was reconstructed, the magnitude and changes with time of the creep pressure were determined by measuring changes in the length, of the brace. However, these measurements were complicated by the fact that creep within the concrete i t s e l f must also be considered. This combined soi1-structure creep problem is typical of many cases in which the foundation must be designed to resist iarge lateral creep pressures. Without methods for the analysis of creep, such designs must be based on high safety factors or the experience of previous fa i lures. Terzaghi (78), in a review of general landslide and creep problems, gives several case studies of large scale creep movements at the bottom of shallow 3 erosion valleys. The change in stresses which in i t iated the creep were pro-duced by the local removal of load during the formation of erosion valleys. Beneath the valley bottom the clay was 'squeezed1 up into valley bulges and the strata were intensely folded, broken and sheared. Such creep movements could have a profound effect on structures bui l t in valley bottoms. Similar bulging due to creep from unloading is sometimes observed in the excavation of open-pit mines and deep foundation s i tes . . Landslides in soi l ; and rock are usually the most obvious and destructive end result of creep leading to creep rupture. Examples of time-dependent deformation leading to slope failures in many different countries have been reported by Ter-Stepanian (76), Suklje (74), Saito (64), Henkel (33), Bjerrum (4) and others (69, 86). Perhaps one of the most disastrous landslides that was preceded by creep leading to rupture was the Vaiont Reservoir sl ide in . Italy. This s l ide , in which i t appears that the creep was in i t ia ted by runoff from heavy ra infa l l swelling clay layers in rather incompetent rock, involved a great mass of earth and rock - more than 240 mil l ion cubic meters - which s l i d into the reservoir, sent water over the dam, destroyed the dam, and k i l l ed 2,117 people in the valley below. As is usually the case for land-slides preceded by measurable creep, the rate of creep was observed to increase about one week before the sl ide (37). , The typical features of slope failures preceded by time-dependent deform-ations may be summarized as follows: they often involve relat ively f l a t clay slopes that form the sides of erosion valleys or cuts; the force that , in i t iates creep after a slope has been stable for many years is often due to groundwater changes or the removal of toe material by erosion or excavation; and, there is usually a noticeable increase in creep rate before the fai lure (86). The many examples found in the l i terature where creep and creep rupture 4 have lead to fai lure and damage show the importance of being able to determine the potential time-dependent behaviour of natural and excavated slopes. This knowledge wil l give some indication of the long-term s tab i l i t y of the slope to the designer and allow remedial measures, i f required, to be in i t ia ted . 1-D FACTORS COMPLICATING THE ANALYSIS OF SOIL CREEP Attempts to analyse the time-dependent behaviour of earth structures, are complicated by several very important factors. Most of these factors involve the problem of relating constitutive relationships determined from idealized laboratory tests and conditions to the actual in s i tu state. The examples presented have generally involved a wide range of cohesive soi ls in the plane strain condition where the creep was in i t ia ted by a reduction of lateral . pressure and hence a decrease in the minor principal stress. An example of a case in which the creep is in i t iated by loading leading to an increase in the major principal stress is the construction of dams or f i l l s on clay layers. The stresses wil l be changing continually with time throughout the earth structure undergoing creep. This is not the case in the typical constant stress t r iaxia l creep tests and therefore tr iaxia l creep tests wil l not be generally applicable. Therefore, plane strain creep tests in which various stress paths leading to creep rupture can be simulated are generally required (45). Further complications that must be considered involve the consolidation history and drainage conditions in s i tu (73). Most sedimentary clay strata are consolidated under one-dimensional or KQ strain conditions. This is an anisotropic condition, whereas laboratory test specimens are usually con-solidated isotropical ly . In addit ion, these clays may often be over-consolidated due to the removal of load by the melting of an ice cover or the 5 removal of surface material, etc. Thus, the applicable drainage conditions and consolidation history must be considered when applying the results of laboratory creep.and creep rupture tests to f i e ld studies. Since creep is basically a long-term phenomenon, the undrained condition is not applicable, except possibly for the shearing of clay layers at great depth or during the rupture stage when fai lure is fa i r l y rapid. To date, most laboratory creep studies have involved undrained tests on saturated undis-turbed or remolded clays. These studies are only applicable to the analysis of time-dependent deformation occurring at constant volume. This l imitation on the available data must be recognized in any creep analysis and is •, discussed further in Chapter 2. There are also less general compositional and environmental factors that are applicable to the analysis of soi l creep (34). Some of these composi-tional and environmental factors are: the type and amount of clay and non-clay materials involved; the pore-water chemistry; the degree of saturation and method of compaction for embankment so i l s ; the water content, structure and sensit iv i ty of undisturbed clays; and the temperature characteristics in s i tu . Many of these additional factors must be considered in order to extend the results of creep tests to meet the f i e ld conditions. 1-E PURPOSE AND SCOPE 't An understanding of the factors that l imit the general a p p l i c a b i l i t y o f creep and creep rupture studies wil l assist in the development of meaningful stress-strain-time relationships for creep analysis. Some of the more, recent creep and creep rupture results from the l i terature are examined in Chapter 2. From the studies reported, stress-strain-time relationships that have been developed for essential ly uniaxial constant stress creep tests wil l be , 6 extended to the multi axial changing stress condition that is applicable in s i tu . Results from creep rupture tests and f i e ld observations of rupture wi l l be used to develop some basic creep rupture c r i te r ia for the analysis. To handle the rather complicated geometries, boundary conditions and loadings involved in the analysis of earth structures, the f in i te element method wi l l be used in conjunction with the relationships developed in Chapter 2. In Chapter 3, the f in i te element method is used to examine prob-lems where the soi l is assumed to be linear viscoelast ic. This is an idealized assumption and its appl icabi l i ty wi l l generally be limited to qualitative studies of soi l behaviour. In Chapter 4, the i n i t i a l strain f in i te element method wi l l be developed to deal with general nonlinear creep. The non-linear analysis presented in Chapter 4 also allows for the introduction of creep rupture. The methods and computer programs developed in Chapter: 3 and Chapter 4 are checked by comparing the results to closed-form solutions from the l i terature. In Chapter 5, some typical embankment and slope problems are examined using the methods developed in the previous chapters. The analysis is : also extended to a few related geotechnical problems in rock and ice mechanics. ,The general results and conclusions are summarized in Chapter 6 and suggestions for further research are presented. 7 CHAPTER 2 DEVELOPMENT OF CREEP AND CREEP RUPTURE RELATIONSHIPS 2-A INTRODUCTION An examination of available laboratory and f i e ld creep studies is necessary in order to formulate creep and creep rupture relationships for use in the analysis of creep in earth structures by f in i te element methods. The primary information required involves: 1. Determining appropriate stress-strain-time relationships from tr iaxia l creep tests. 2. Extension of the uniaxial , constant stress, deformation relationships developed from t r iax ia l creep tests to the multiaxial changing stress state that is representative of the actual f i e ld conditions. 3. Developing rea l i s t i c creep rupture c r i te r ia to predict when the soi l wi l l fa i l due to creep. The time-dependent behaviour of other engineering materials wil l be d is -cussed where i t assists in the development of the relationships for s o i l . 2-B TYPICAL CREEP BEHAVIOUR OF COHESIVE SOILS OBSERVED IN LABORATORY TESTS The general creep behaviour of cohesive soi ls (70, 71) is quite similar to that.observed for other materials. Although the compositional and environ-mental factors outlined in the previous chapter wi l l determine the speci f ic creep characteristics of a given s o i l , a typical laboratory creep curve for axial deformation under a constant deviatoric stress wi l l be of the form shown in Figure 2-1 a. When the constant deviatoric stress a. is applied to the soi l X PRIMARY JL SECONDARY *m~ TERTIARY * FAILURE CONSTANT DEVIATORIC STRESS arA IME, t FIG. 2-la TYPICAL CREEP CURVE SHOWING TYPES OF POSSIBLE BEHAVIOR IN THE TERTIARY STAGE. STRESS a < cc i-3 FAILURE CONSTANT DEVIATORIC d3 CONSTANT DEVIATORIC STRESS cr._ CONSTANT DEVIATORIC STRESS cr dl TIME , t FIG. 2-lb EFFECT OF STRESS LEVEL ON TYPICAL CREEP CURVES. FIG. 2-1 TYPICAL CREEP BEHAVIOR OF SOILS OBSERVED IN LABORATORY TESTS (70). 9 at time t=o, an instantaneous strain e 0 is observed. This strain is undoubt-edly made up of e last ic and plast ic components, but i t is usually considered to be an instantaneous recoverable e last ic stra in. At any later time t 1 , the total strain e wil l be given by the sum of the recoverable strain e'Q and the non-recoverable creep strain e C : e = + E C (2-1) Relaxation tests are not generally performed on soi ls so that the relative magnitude of e! is not known. 3 o The creep curve for most materials is usually divided into three regions, although the boundaries are often quite arbitrary and a cohesive soi l may not exhibit a l l of the stages represented. In region I (Figure 2-1 a), the strain rate decreases continuously and the strain-time curve is concave downwards: creep in this region is called primary or transient creep. The strain-time curve for many materials has an approximately constant slope in region II: creep in this region is referred to as secondary creep. For cohesive soi ls and most other materials, i t appears that true steady-state creep does not exist (45, 71, 73, 43) and so the term 'steady-state' wi l l not be adopted in this thesis for describing region II. F inal ly , region III of the creep curve corresponds to a f inal transition stage when the strain rate decreases to an imperceptible value or increases to end in creep, rupture (concave upwards): creep in this region is called tert iary creep. A set of typical creep curves for soi l is shown in Figure 2-1b where the constant deviatoric stresses are such that a .„ > a . > a . . For the lowest d3 d2 di stress c/j , l i t t l e effect is produced and the strain rate may not be measurable after some time. For the highest stress , the strain-time curve has a point of inf lect ion or minimum strain rate where the strain rate changes from a decreasing to an increasing function. This behaviour indicates the existence 10 of an upper y ie ld strength for cohesive soi ls (53, 81, 73). For stresses below this upper y ie ld strength, the strain rate wi l l eventually become neg-l ig ib le after large times. Stresses above the upper y ie ld strength wil l lead to an increasing strain rate which results in large deformations and, even-tual ly , creep rupture. (73, 45, 72) The elapsed time from the beginning of creep to the instant that the minimum strain rate is reached, and the value of this minimum strain rate, is of importance in the development of creep rupture c r i t e r i a . Many investigators consider the apparent cessation of creep after a long time at lower stress levels to be due to strain hardening of the so i l s . Existence of a lower bound stress below which creep does not occur, as suggested by Murayama and Shibata (54), has not been confirmed by other researchers (11, 73). The change from decreasing to increasing strain rate after some strain has accumulated at the higher stress levels is usually attributed to strain softening (72). Between the lower and higher deviator stresses, there is a broad range of applied stresses for which the strain rate continues to decrease over a long period of time and creep rupture does not occur. This is generally the range of stresses that is of concern to the soi ls engineer and most available stress-strain-time relationships (70, 71, 72) have been developed to cover this range. With the exception of the recently reported work by Bishop and Lovenbury (3), the laboratory studies of so i l creep have been of relat ively short duration when compared to the times involved in the f i e l d . Bishop and Lovenbury found in a series of drained creep tests last ing up to 3j years that the creep behaviour could be represented with reasonable accuracy by simple logarithmic or power laws for a short i n i t i a l period of 20 to 200 days, depending on the stress level and type of c lay, and then the strain r^te either increases or levels off after the i n i t i a l period of steadily decreasing 11 axial stra in. If the strain rate increases, a new, almost l inear relationship between log strain and log time is observed. The behaviour observed in these long-term creep tests is qual itat ively the same as that observed in the more usual short-term creep tests. 2-C METHODS OF REPRESENTING THE TIME-DEPENDENT DEFORMATION OF COHESIVE SOILS There are three general approaches that have been adopted to study the creep of cohesive soi ls and other engineering materials: the fundamental approach; the empirical approach; and the phenomenological or rheological model approach (32, 34). A generalized relationship developed from any of these approaches should satisfy the following c r i te r i a suggested by Singh and Mitchell (71): "(1) i t must be applicable to a reasonable range of creep stresses; (2) i t must describe the behaviour of a range of soi l types; (3) i t must account for both l inear and curved relationships between strain and logarithm time; and (4) i t must contain parameters that are easily deter-mined". A review of the l i terature indicates that there has been l i t t l e analysis and synthesis of the available laboratory data and f ie ld information for application to the study of time-dependent deformations in s i tu . Thus, in addition to Singh and Mitchel l 's c r i t e r i a , the fundamental cr i ter ion of being usable in the analysis of f i e ld problems must be added. 2-C-l Fundamental Approach ; In .the fundamental or micromechanistic approach, an attempt is made to relate the creep behaviour of the material to processes on the atomic or mole-cular scale. For instance, the so-called dislocation theory of plast ic flow (24),is used to explain and predict many of the mechanical properties of metals. Since experimental evidence indicates that creep involves thermally 12 activated processes (24) the absolute reaction-rate theory or rate process theory, which was employed by Eyring in 1936 (20) to predict the v iscosity, p last ic i ty and diffusion of viscous f lu ids , has been used to study the creep of many materials such as metals, p last ics , bituminous materials and cohesive soi ls (53, 55, 11, 49, 50). Rate process theory has been used by Mitchel l , Campanella and Singh (50) to develop a strain rate equation for the creep of cohesive soi ls directly from considerations of micromechanistic behaviour. At present, the strain rate equations based on rate process theory are not adequately developed for use in the analysis of creep problems and the various parameters are d i f f i cu l t to determine with conventional laboratory equipment (34). Nevertheless, this micromechanistic approach is extremely useful in contributing to a better understanding of the creep deformation of soi ls and provides insight into,the bonding mechanisms that contribute to soi l resistance (52). 2-C-2 Empirical Approach In the empirical approach, various parameters such as strain and strain rate are measured experimentally as a function of time, stress and temperature under controlled conditions. The results may then be used to develop functional relationships between the parameters that adequately describe and predict the material behaviour or may be used in tabular form in computer programs. This approach has been ut i l ized extensively for the presentation of data from creep experiments in metal (43), rock (32, 36), and cohesive soi ls (70, 71, 72). Singh and Mitchell (71, 72) have developed a generalized stress-strain-time function for cohesive soi ls which is based on the study of creep curves for many cohesive soi ls over a range of sustained deviatoric stresses. They observed that the relationship between strain rate and time for various stress 13 levels was represented by parallel straight lines on a log-log plot as shown in Figure 2-2a for stresses below the upper y ie ld and until the minimum strain rate was reached for stresses above the upper y i e ld . A l inear relationship as shown in Figure 2-2b was also observed between log strain rate and deviator stress level for points corresponding to the same time after load application at stress levels below the upper y i e l d , but high enough to cause measurable deformations. The relationships indicated by Figures 2-2a and 2-2b can be expressed in the following equation: i In e(t,D) = In A + «D - m In f t or e(t,D) = Ae' f t . ) m (2:2) (2:3) where e(t,D) = strain rate at time t , a function of the deviator stress D, A = projected value of strain rate at time t 1 S and D=o in Figure 2-2b, <* = value of the slope of the mid-range l inear portion of the plot of log strain rate versus deviator stress, a l l points corres-ponding to the same time after load appl i -cation. (Figure 2-2b), m = slope of the straight line on the log strain rate versus log time plot (Figure 2-2a). This simple three-parameter relationship was developed by Singh and Mitchell from their experimental work and data in the l i terature ; and they state (71): " i t appears valid irrespective of whether the clays are undis-turbed or remolded, wet or dry, normally consolidated or overconsolidated, or LOG TIME , t FIG. 2 - 2 a D E T E R M I N I N G m. PRINCIPAL STRESS DIFFERENCE, D FIG. 2-2b D E T E R M I N I N G A AND a. FIG. 2 - 2 STRAIN RATE E Q U A T I O N P A R A M E T E R S (34). 15 tested drained or undrained." The parameters A, <* and m can be readily deter-mined for any given soi l from creep tests on two identical samples subjected to different deviatoric stress levels. The f i c t i t ious nature of the parameter A is shown in Figure 2-2b. However, Singh and Mitchell feel that the para-meter A is meaningful in that i t indicates the order of magnitude of the creep rate for the particular cohesive s o i l , and reflects the so i l ' s structure, composition and stress history. The parameter oc indicates the stress level effect on the creep rate and, from analogy with rate process theory, i t may be hypothesized that = reflects the number of bonds per unit area resist ing the creep deformation. (52, 71) The parameter m provides a measure of the creep potential of a s o i l : soi ls with m<i eventually fa i l in creep rupture; soi ls with m=i seem to exhibit the same strength before and after creep; and, soi ls with m>i exhibit cessation of creep with time. The lower the value of m, the greater is the tendency of the soi l to accumulate strain with time. It has been observed that m is not unique for a given cohesive soi l and depends on the consoli -dation history (72). Goldstein (25) also developed a generalized three-parameter deformation relationship for constant volume unconfined compression tests that is similar in form to Equation 2:3. Singh and Mitchel l 's three-parameter relationship (Equation 2:3) can be rewritten in the form: e = A e ^ (^)m (2:4) where t 1 has been taken as unity, D = the normalized stress leve l , defined as the ratio of the deviator stress to strength at the 16 beginning of creep, /D , or i t may be the stress 3 a r-> m a x » J ratio *Vp where q is the deviator stress and p, the mean normal stress, « = the dimensionless parameter defined as the value of the mid-range l inear portion of the log strain versus normalized stress level or stress ratios plot , a l l points corresponding to the same time after load application. Equations 2:3 and 2:4 were developed empirically from data for a large number of tr iaxia l creep tests in which the sample was usually f i r s t iso-tropical ly consolidated and then axial ly loaded. Thus, i t is s t i l l necessary to check the effect of K0 consolidation, plane strain test conditions, and total stress paths of loading, on the form of the equations. Mallawaratchie (45) recently considered these factors in a series of creep rupture tests at stresses near the upper y ie ld stress. While these tests are mostly out of the stress range that the empirical equation was developed for, the creep curves were found to be of the same form as those for standard tr iaxia l creep tests. Bishop and Lovenbury (3) checked total stress paths of loading in the t r iaxia l creep tests that were referred to previously. The form of the creep curves was found to be the same for the increased major principal stress tests and decreased minor principal stress tests. On the basis of the available data from creep tests on cohesive s o i l s , i t appears that Singh and Mitchel l 's empirical equation is adequate to describe the creep deformation of a wide range of soi ls for various stress levels , time intervals and test conditions. The parameters m, A and « must be determined for each soi l and test condition but, as indicated, this requires only two creep tests. For these reasons, Singh and Mitchel l 's empirical deformation 1 7 relationship will be adopted in this study, with the realization that i t must be extended to describe the multiaxial changing stress state that is represen-tative of the actual f i e ld conditions. 2-C-3 Rheological Model Approach In the rheological model approach, a model made up of l inear or non-l inear springs, dashpots and sl iders is developed that is supposed to behave mechanically l ike the actual material. The mathematical description of the creep behaviour of the rheological model is then used in the stress analysis. Since the constitutive equations are developed for the model of the material, the adequacy of the model selected l ies in how well i ts behaviour represents the actual material behaviour. This approach has been applied extensively in the f i e ld of rheology to study the flow of rubber, p las t i cs , viscous fluids and geologic materials (32). A few representative examples of the many rheological models that have been,developed to describe the creep of cohesive soi ls are shown in Figure 2-3. Models such as those shown may be considered as either actually representing some of the characteristics of the microstructure of the soi l under consider-ation or merely as providing a useful analogy for determining appropriate forms of the stress-strain relationships. Most of the rheological models shown have l inear spring elements with e last ic moduli denoted by E.., K.. or G.. The dashpot elements may be linear (Newtonian) with a coeff icient of viscosity n - , or nonlinear (non-Newtonian) with the viscous flow obeying a hyperbolic law based on rate process theory or an exponential law based on empirical studies. . While the mathematical description of the various rheological models; appears adequate to describe the creep behaviour of a specif ic soi l over a particular stress range and time interva l , the models do not offer a T E i 18 ^0 3d A, , B. FIG. 2-3a MURAYAMA AND SHIBATA'S NONLINEAR MODEL (53) | G 2 ^ 2 FIG. 2-3b SCHIFFMAN'S LINEAR BURGERS BODY MODEL (67). 1^ G 2 ^2 63 [ f v W < — T ^4 G4 FIG. 2-3c LARA-TOMAS' LINEAR MODEL ( 54). S <—vVA-1^ 77T7T7T — ^2 FIG. 2-3d SHERIF'S NONLINEAR MODEL (69). 1 T7 FIG. 2-3e SINGH'S NONLINEAR MODEL (70). FIG. 2-3 PROPOSED RHEOLOGICAL MODELS OF SOILS. 19 generalized representation of the creep behaviour of cohesive soi ls (34, 68). Also, i t is now generally agreed that the creep of most soi ls is nonlinear (17, 34, 67, 68) so that l inear rheological models usually represent an idealized soi l behaviour. This nonlinear creep behaviour has been studied in tests by Kondner and Krizek (38) and more recently by Drescher (17). Drescher found that superposition, and hence l inear v iscoelast ic i ty , did not hold for his short-term uniaxial creep tests and suggested using a nonlinear integral rep-resentation of the creep behaviour as proposed by Green and Ri v i i n (27). At present, i t is very d i f f i cu l t to determine the parameters for such a represen-tat ion, and extension of the results to long-term creep behaviour is s t i l l required. Four simple l inear one-dimensional mechanical models that are often used to describe viscoelastic behaviour and are included in some of the models,for cohesive soi ls are shown in Figure 2-4. The defining one-dimensional stress-strain equations of these models are as follows: Maxwell Body or el astico-viscous material (Figure 2-4a); a + r]-o = £ (2:5) n l K i Kelvin Solid or firmo-viscous material (Figure 2-4b); a = kje + (2:6) Standard Linear Solid or 3-parameter so l id (Figure 2-4c); (kj + k 2) a + n 2 a = klk2 e + k ^ e (2:7) Burgers Body or 4-parameter f lu id (Figure 2-4d). " k 2 k 2 n 2 n 2 _ a + ( 1 + _ _ + _ ) - = ^ + ^ ( 2 : 8 ) In Equations 2:5 to 2:8, o is the total applied stress, e is the total s t ra in , n. and k. are the dashpot viscosit ies and spring stiffnesses respectively as k l V\ vWV—Tf FIG. 2-4a MAXWELL BODY (ELASTICO-VISCOUS MATERIAL). M/\A ^—H— FIG. 2-4b KELVIN SOLID (FIRMO-VISCOUS MATERIAL). -AVW k. I—)J\f^ j k 2 FIG. 2-4c STANDARD LINEAR SOLID (3-PARAMETER SOLID) -WW—3J-k l V\ r~WWV^—I ^2 FIG. 2-4d BURGERS BODY ( 4 - PARAMETER FLUID) FIG. 2-4 COMMON RHEOLOGICAL MODELS. 21 shown for each model in Figure 2-4, and the superscript dots indicate deriv-atives with respect to time. In order to obtain stress-strain relationships that more closely approx-imate the observed results , springs and dashpots can be added to the basic models in various series and parallel arrays such as in the Lara-Tomas. model (Figure 2-3c). This procedure leads to stress-strain relationships containing higher derivatives of stress and strain with respect to time. Alternatively, the stress-strain relationships can be generalized directly by adding higher time derivatives without considering a corresponding mechanical model. Thus, the general one-dimensional l inear differential viscoelastic stress-strain relationship can be expressed in the form: Pa = Qe (2:9) where P and Q are l inear dif ferent ia l operators given by: P P = T a m ^ — (2:10) m=o 9t q .m Q = i o b ^ ( 2 : 1 1 ) m=o at The practical consequence of such a generalized relationship is that more parameters, a m and b , are available to be f i t ted to the available experimental data. Even though l inear rheological models do not offer a generalized descrip-tion of the behaviour of cohesive s o i l s , they often provide useful approxi-mations.' The deformation and stress behaviour of a l inear viscoelastic b,ody can be qual itat ively examined for a wide range of models and parameters. This information is particularly valuable when laboratory test results are not available and the material properties must be assumed, or developed from f i e ld measurements (77, 86). Also, there are some well-developed methods available m 22 for the stress analysis of l inear viscoelastic bodies. For these reasons, the assumption of l inear viscoelastic soi l behaviour wil l be adopted in this study in addition to the previously discussed empirical approach. This idealized l inear v iscoelast ic i ty approach to solving problems where the time-dependent soi l behaviour is considered, is analogous to the idealized l inear e las t ic i ty approach that is usually used in soi l mechanics problems where time-dependent behaviour is not considered. 2-D EXTENSION OF THE UNIAXIAL RELATIONSHIPS TO THE MULTIAXIAL CASE The extension of the one-dimensional stress-strain relationships for. a constant deviatoric stress to the general three-dimensional state of stress and strain must now be considered. Extension of the one-dimensional l inear viscoelastic stress-strain relationship given by Equation 2:9 is quite straight-forward so this case wil l be examined f i r s t . 2-D-l Constitutive Equations for Isotropic Linear Viscoelastic Behaviour : The usual procedure in developing the constitutive equations of l inear e las t ic i ty or v iscoelast ic i ty is to decompose the stress and strain tensors a- • and e.j into their mean components and e^, and their deviatoric components, S.j . and e.y True uncoupling of the shear and dilatational effects is not really the case for soi ls but this usual assumption wil l be adopted in the present l inear viscoelastic analysis. Isothermal behaviour wil l also be , assumed since the temperature variation below the surface of earth structures wi l l be minimal. The stress-strain relationships for isotropic l inear e last ic behaviour are then given by: S^. = ZGR.. in shear, ( i , j = 1 , 2 , 3 ) (2:12) CTM = 3KeM in di latation (2:13) where G and K are the shear and bulk moduli respectively. 23 By analogy, the stress-strain relationships for isotropic l inear visco-e last ic behaviour are given in l inear differential form by: P'S,. = Q'e. . in shear (2:14) P"aM = Q"eM in dilatation (2:15) where P', Q 1 , P", Q" are l inear differential time operators with constant.co-e f f i c ients : P' = I a n ^ , Q' = [ b ^ , (2:16) n=o n at n=o n 3t P" an q " .n P" = I c n ^ , Q" = H „ ^ . (2:17) n=o at n=o a t " It can be seen that the one-dimensional l inear dif ferent ia l v iscoelastic stress-strain relationship given by.Equation 2:9 is a special case of Equations 2:14 and 2:15. The assumptions of uncoupling, isothermal conditions, and isotropic be-haviour are idealizations that have been made in order to allow for the stress analysis of soi l mechanics problems. With the development of more refined laboratory descriptions of the time-dependent soil properties i t should be possible to handle the coupled anisotropic case (2). ' The three-dimensional shear behaviour of the models shown in Figure 2-4 is then given by the following equations: Maxwell Body in shear (Figure 2-4a) 1 S. . + - i - S. . = e , , (2:18) 2nj i j 2G X i j i j Kelvin Solid in shear (Figure 2-4b) , S i j = 2 G i e 1 J + 2 n i 4 i J ( 2 : 1 9 ) Standard Linear Solid in shear (Figure 2-4c) 24 Burgers Body in shear (Figure 2-4d) S.. + f ^ + ^ + ^ ] § = 2 n e.. + 2 ^ e . . (2:21) 1J \G2 Gx G2; i j GjGg i j 'l i j G 2 i j where and are the viscosit ies and shear moduli respectively for each model (5). Similar expressions can be written to describe the dilatational behaviour for the various models. It is usual, but not necessary, to assume that materials are e last ic in dilatational behaviour. The relation between .and is then given by Equation 2:13 for l inear e last ic behaviour. The const i -tutive equations developed in this section wil l be used in Chapter 3 when a method for the viscoelastic stress analysis of earth structures is presented. 2-D-2 Constitutive Equations based on Singh and Mitchel l 's Empirical Relationshi p The development of three-dimensional stress-strain-time relationships based on Singh and Mitchel l 's nonlinear one-dimensional empirical relationship (Equation 2:3) wi l l require the following basic steps: 1. Application of the constant deviatoric stress relation-^ ship given by Equation 2:3 to the more usual case of changing stress conditions involved in f i e ld problems with mixed boundary condition. This requires the dev-elopment of a cumulative creep rule that is appropriate for cohesive so i l s . 2. Relating the nonlinear creep behaviour under multiaxial stresses to the behaviour observed in t r iaxia l creep tests on cohesive so i l s . 3. Determining the components of creep strain for any multi-axial relationship developed. 25 1. Cumulative Rule Very l i t t l e work has been reported on the formulation of cumulative rules to describe the nonlinear creep of cohesive soi ls under changing stresses with time. However, a number of cumulative creep rules have been proposed to des-cribe the time-dependency of the strain path for the creep of metals, plastics and concrete. The two basic rules that have been developed are the strain-hardening and time-hardening creep rules which are i l lustrated by the hypothe-t ical creep curves for a decreasing uniaxial load case shown in Figure'2-5. For the strain-hardening rule, the current strain rate depends only on the current creep strain and the current stress; the effect of previous strain rates being negligible. In terms of the creep strain-time curves (for a fixed temperature) a change in stress is represented by a horizontal l ine on Figure 2-5a from one constant stress curve to the new stress curve, corresponding to the same creep stra in. On the other hand, for the time-hardening rule, the current strain rate depends only on the elapsed time and the current stress. A change in stress wil l show up as indicated in Figure 2-5b, by a vertical transfer from one constant stress curve to the new stress curve, corresponding to the same elapsed time. Modifications and combinations of these two basic rules have been developed to represent the cumulative creep of many materials (61). Snead (73) and others (45, 82) have examined the appl icabi l i ty of some of these cumulative creep rules for cohesive so i l s . Generally, this research has been concentrated on studying the strain-hardening cumulative creep rule and, for this reason, this section wil l be focused on the strain-hardening rule. The strain-hardening cumulative creep rule for a uniaxial loading and isother-mal conditions is usually expressed in the functional form (61); e c = f(e c ,a ) : (2:22) FIG. 2-5a STRAIN-HARDENING CUMULATIVE CREEP RULE. TIME , t HOURS FIG. 2-5b TIME - HARDENING CUMULATIVE CREEP RULE. FIG. 2-5 CUMULATIVE CREEP RULES (47). (STEEL AT HIGH TEMPERATURE) 27 where e c = creep strain rate, e c = current creep s t ra in , a = current stress. Equation 2:22 is a particular form of the mechanical-equation-of-state for isothermal conditions developed by workers in metallurgy (24, 43, 52). For the strain-hardening cumulative creep rule to apply for any material, the effects of prior strain rate must be negl igible. This is usually checked by comparing the data for two, or preferably more, tests with different rate histories such as constant stress creep tests , incremental stress tests , and strain controlled tests. Consider the constant strain rate tensile test and the constant tensile stress creep tests for steel at high temperature shown in Figure 2-6. If the effects of prior strain rate are negl igible, then when any two of the three variables S (stress) , e c (creep strain) or e c (creep strain rate) are known or constant for the test , the third variable is uniquely deter-mined. For instance, when the creep stress level has reached a value of S 1 5 the creep strain e ' is known for the constant strain rate tensile test at a strain rate of e .' For the constant tensile stress creep test at a stress of S j , when the strain rate has a value of e c, the creep strain wil l also be;e as shown in Figure 2-6 i f the effects of prior strain rate are negl igible. It appears that for plastics and metals, where the temperature remains constant and metallurgical changes do not occur, that data obtained under one set of conditions such as a constant stress creep test or constant strain rate test can be applied to a situation in which the stress is varying (43, 52).; Snead (73) conducted a series of t r iax ia l tests on Haney clay to determine whether there is a functional relationship between current deviator stress, current strain and current strain rate for this particular clay that is indep-endent of the strain rate history. Results for constant stress creep rupture 28 FIG. 2-6 SCHEMATIC DIAGRAM ILLUSTRATING THE RELATION BETWEEN TYPES OF TESTS WHEN THE EFFECT OF PREVIOUS STRAIN RATE IS NEGLIGIBLE (43). (STEEL AT HIGH TEMPERATURE) 29 tests, incremental increasing stress tests , and constant strain rate tests were plotted for each consolidation history (normally consolidated, over-consolidation ratio = 2, and overconsolidation ratio = 6) as shown in Figure 2-7 for normally consolidated Haney clay. From Figure 2-7, i t can be seen that the values of the deviator stress at any strain and strain rate are very close, regardless of the type of t r iax ia l test. Snead concluded that a r e l -ationship does exist between the current deviator stress, current strajn and current strain rate ( i .e . = f(e, e)) that is independent of the strain rate history for the following conditions: the strain is continually increasing; the temperature remains constant; adequate time is allowed for pore pressure equalization within the test sample; the consolidation history of thejsamples is the same; the tests are undrained, since drainage introduces additional parameters; and, the soi l is not heavily overconsolidated so that the shear strains wil l be uniform throughout the sample. No attempt was made to evaluate the actual parameters for the deviator stress-strain-strain rate relationship, but examples were given to show how the behaviour of samples in one type of test can be evaluated from the results of another type of test. Snead feels that this relationship wil l hold for other saturated cohesive so i l s . Mallawaratchie (45) has examined the val idity of Snead1s hypothesis of a deviator stress-strain-strain rate relationship for plane strain and KQ t r i -axial creep and incremental loading tests on Haney clay. Since the samples are under a deviatoric stress at the end of KQ consolidation, there are two deviator stresses that can be considered: the additional deviator stress applied after consolidation; and, the total deviator stress on the sample during creep. Mailawaratchie1s results for undrained a increasing plane strain creep and incremental loading tests show a reasonable agreement for the total deviator stresses at any strain and strain rate. However, in terms of the 30 FIG. 2-7 AXIAL STRAIN RATE-AXIAL STRAIN CURVES FOR TRIAXIAL TESTS ON NORMALLY CONSOLIDATED HANEY CLAY. (73) 31 additional deviator stresses, differences of 20% and higher were found between the stresses. Similar strain-time and strain rate-time results were found for undrained a 3 decreasing plane strain creep and incremental loading tests. It should be noted that i t is very d i f f i cu l t to determine the appropriate strain rate in incremental loading tests and this d i f f i cu l t y wil l be reflected in any comparison. Mallawaratchie found that the stress-strain-strain rate relationship was not quantitatively independent of test type. A different quantitative relat ion-ship was observed for increasing a l 5 KQ t r i ax i a l , conventional t r i ax i a l , and plane strain creep tests. As mentioned previously, there is a fundamental difference between the in i t i a l conditions for the plane strain and KQ t r iax ia l samples and for the conventional t r iax ia l samples. After consolidation, the plane strain and KQ t r iaxia l samples, which have been KQ consolidated, are under a deviator stress; while the conventional t r iax ia l samples, which have been isotropical ly conso-l idated, are not under an i n i t i a l deviator stress. Also, during, creep,, the conventional t r iaxia l and KQ t r iax ia l samples wil l deform latera l ly with a 2 = a , while the plane strain samples wil l deform latera l ly with no longitud-inal strain and a2 t o3. Therefore, during creep, the modes of application and variation of stresses and strains are different for each type of test. Thus, in order to predict f i e ld behaviour using empirical laws and a cumulative creep rule, i t is important that the type of laboratory test used reproduces the f i e ld conditions since there does not at the present time appear to be a fundamental relationship between the various types of tests and test conditions. Vyalov and Meschyan (82) have experimentally and theoretically checked the appl icabi l i ty of the various cumulative creep laws for cohesive s o i l . From their research, i t was concluded that the strain-hardening cumulative creep law 32 best described the creep of soi ls subjected to variable stresses. The time-hardening and hereditary cumulative creep laws also gave fa i r l y reasonable results , especially for small stress changes. From the experimental evidence available to date, i t would appear that the strain-hardening cumulative creep law gives a reasonable description of how the creep strain rate should be determined when the stresses are changing during the undrained creep of a saturated cohesive s o i l . However, for creep problems in soi l mechanics, the f i e ld drainage conditions wil l be somewhere between the undrained and drained cases, and will generally be indeterminate. In the pre-sent development, i t wi l l be assumed that there is no volume change due to the creep strains. This assumption is made to fac i l i ta te the use of the strain-hardening cumulative creep law in it's present total stress form. With further laboratory research on creep in terms of effective stresses and on the effect of drainage conditions on the stress-strain-strain rate relationship for various types of tests, i t may be possible to include time-dependent volume changes in a nonlinear creep analysis. 2. Multiaxial Stress-Strain Relationship The nonlinear creep relationships for cohesive soi ls that have been d is -cussed are based on tr iaxia l tests. In order to use these relationships in two-dimensional and three-dimensional soi l mechanics problems, i t is necessary to determine some generalized strain increment and generalized stress that can be used in the creep analysis. Unfortunately, no work has been reported on the formulation of generalized relationships for the nonlinear creep of so i l s . However, there is some data available on the multiaxial creep of metals at elevated temperatures indicating that the classical p last ic i ty laws can be used to relate the multiaxial creep behaviour to the uniaxial creep behaviour (43, 46). Although the physical processes involved in creep are probably different 33 than those involved in plast ic flow, i t is usual to assume that the same r e l -ationships hold for creep strains as for plast ic strains. The nonlinear creep strains wil l be considered to be non-recoverable. The current state of creep strain in a body is an accumulation of creep strain increments, each of which was governed by the stresses prevail ing at the time at which the particular creep strain increment occurred. This concept, that the creep strain increment, and not the total creep s t ra in , is related to the stresses, is usually referred to as an incremental theory and corresponds to the flow theory of p las t i c i ty . The usual generalized strain increment adopted in isotropic work-hardening p last ic i ty and in the study of the creep of metals is the equivalent strain increment given by: C A C C~ _ i . c (2:23) r where Ae g = equivalent creep strain increment, r Ae -• = increments of the creep strain tensor, r A y o c t = o c t a h e d r a l creep shear strain increment. The prefix A wi l l be used to denote increments and the superscript C wil l be used to denote creep values throughout. The usual generalized stress adopted in p last ic i ty and in the study of the creep of metals is the equivalent stress a g given by: a = / - S. .S.. e 2 i j " l j (2:24) 3_ = fl T ° C t where S. . = the stress deviator tensor, T . = octahedral shear stress, oct 34 Thus, the equivalent creep strain increment and equivalent stress are proportional to the octahedral creep shear strain increment and the octahedral shear stress, respectively. For most problems i t is more convenient to use the equivalent creep strain increment and equivalent stress because they be-come numerically equal to the uniaxial creep strain and stress in the l imit ing C uniaxial case (43). The equivalent creep strain increment Aeg and theiequiv-alent stress aQ wi l l be used in the nonlinear uniaxial creep relationships for cohesive soi ls so that the nonlinear creep behaviour under multiaxial stresses can be determined for soi l mechanics problems. In general, since creep strain increments are involved, solutions must be obtained by incremental methods. 3. Components of the Equivalent Creep Strain Increment C Once the equivalent creep strain increment Aee at any stage of creep has been determined by substituting the equivalent stress a g into the desired uni-axial creep relationship, i t is necessary to determine the components of the creep strain increment. Once again, the same relationships are assumed to apply for the components of the non-recoverable creep strain increment as for the components of the plast ic strain increment. It is assumed that the compon-ents of the creep strain increment are proportional to the instantaneous values of the deviatoric stresses for an isotropic material; i . e . , A e i j = V c ( 2 : 2 5 ) where Ac is the constant of proportionality for any increment. Equation 2:25 is analogous to the Prandtl-Reuss equation (46) of p l as t i -c i ty and states that the increments of creep strain depend on the current values of the deviatoric stresses and not on the stress increment required to reach this state. This equation also implies that the principal axes of the stress and of the creep strain increment tensors coincide. The assumption of no volume change due to the creep strains is also included in Equation 2:25 35 s i n c e : Equation 2:25 merely gives a relationship between the ratios of the com-ponents of the creep strain increment in the different directions and, in order to determine the actual magnitude of the increments, the expression for Ac must be determined. By expanding Equation 2:25 i t can be shown that: AC = f — (2:27) e • and, substituting this value of Ac into Equation 2:25 y ie lds : AeV. = -3 — S - s . . (2:28) 1J 2 ° e C With the equivalent creep strain increment Ae g determined at any stage of r creep, the components of the creep strain increment A e . . can be determined. 13 An incremental solution of the f in i te element formulation of nonlinear creep problems wil l be developed in Chapter 4. The solution method that is presented is based on the strain-hardening cumulative creep rule, the multi-axial stress-strain relationship in terms of the equivalent creep strain incre-ment and equivalent stress, and the expression for the components of the creep strain developed in this section. The solution method can also be extended to handle other cumulative creep laws and other generalized incremental creep strain-stress relationships. 2-E CREEP RUPTURE CRITERIA Hirst and Mitchell (34) have updated Mitchel l , Seed and Paduana's (51) extensive review of the l i terature discussing the influence of creep on the shear strength of so i l s . Without detai l ing these excellent reviews, the general conclusions that were reached can be summarized as follows: 36 1. "In general, saturated soft sensitive clays and over-consolidated clays lose strength during undrained creep. On the other hand, the effects of creep on the strength of part ia l ly saturated and compacted soi ls have been observed to be quite variable. Strength gain, no change in strength and strength loss have been observed during creep of these so i l s . 2. . . . strength loss is usually attributed to shear strains that remold the specimen, and strength gain is postulated as result ing, in part, from the stru-ctural alterations caused by shear strains. 3. No mathematical expression relating strength change during creep to composition or environment was ev i -dent from the data analysed." While mathematical expressions have not been developed for the strength change during creep, Saito and Uezawa (65), Saito (64), Snead (73), Singh and Mitchell (72) and Mallawaratchie (45) have al l developed relationships for the creep rupture l i f e of cohesive so i l s . Also, Sherif (69) and Snead (73) have developed expressions for the upper y ie ld strength of particular cohesive so i l s . i 2-E-l Creep Rupture Life of Cohesive Soils Some of the terms used to describe the creep rupture of cohesive soi ls can be explained with the aid of Snead's log strain rate - log elapsed time curves for undrained tr iaxia l creep tests on normally consolidated Haney clay shown in Figure 2-8 (73). The point of transient minimum strain rate has been marked on each creep curve and i t can be seen that a l l of the transient mini-mum strain rates nearly l i e on a straight l ine for this plot. All of the, 37 FIG. 2-8 AXIAL STRAIN RATE-TIME RESULTS FOR UNDRAINED TRIAXIAL CREEP TESTS ON NORMALLY CONSOLIDATED HANEY CLAY. (73) 38 strain rate-time curves start below the transient minimum strain rate l ine and, i f they cross this l i ne , creep rupture wil l eventually occur. The lowest curve, with a deviator stress of 42.8 ps i , is unlikely to reach the transient minimum strain rate l ine and this sample wi l l not f a i l in creep rupture. This indicates the existence of an upper y ie ld strength for Haney clay which is in agreement with Murayama and Shibata (55) and Vialov and Skibitsky's (81) hypo-thesis that an upper y ie ld strength exists for s o i l s , below which creep rupture wil l not occur. For deviator stresses below this upper y ie ld stress, the strain rates eventually become negligible after long times i f a l l the test conditions remain the same. Snead found that a transient minimum strain rate line could be drawn for each series of tests with the consolidation histories l i s ted on Figure; 2-9. Figure 2-9 also gives the l ine of transient minimum strain rates, e^, for a l l of his creep rupture tests. Snead proposed that for Haney clay, this l ine of transient minimum strain rates exists independent of consolidation history, stress level and drainage conditions, as long as the total stresses remain constant. The per cent reduction in strength between the 'normal' reference strength and upper y ie ld strength is nearly the same (14.6 to 17.7 per cent) for al l of the consolidation histor ies. The time from the beginning of creep is usually not known for f i e ld cases involving natural slopes, and therefore any elapsed time relationship cannot be used to predict the creep rupture l i f e for such cases. On the basis of their t r iax ia l compression tests and the results of other researchers, Saito and Uezawa (65) proposed a l inear relationship between the logarithm of , 'secondary' strain rate, e^, and the logarithm of the total time to creep rupture, t r- This relationship, along with some test results from Snead, Sherif (69) and Campanella's (9) research, is shown in Figure 2-10. Saito 39 LINE OF TRANSIENT MINIMUM STRAIN RATES log (0t =.037-1.08 log|0 e m LEGEND HANEY CLAY CONSOLIDATION HISTORY x - NORMALLY CONSOLIDATED 0 - 77 = 2 A - T) = 6 n ~ 17 = 25 • -77 = 2 5 DRAINED TESTS ~ (77 = OVERCONSOLIDATION RATIO) UNDRAINED TESTS I0W IO IO IO IO TIME OF TRANSIENT MINIMUM STRAIN RATES, t MINUTES FIG. 2-9 LINE OF TRANSIENT MINIMUM STRAIN RATES FOR UNDRAINED TRIAXIAL CREEP TESTS ON HANEY CLAY. (73) 40 FIG. 2-10 TOTAL RUPTURE LIFE FOR LABORATORY CREEP TESTS ON COHESIVE SOILS. 41 and Uezawa's relationship, which was obtained for different types of s o i l s , stress levels , consolidation histories and drainage conditions can be approx-imated by: t r = (2:29) e s where t = total rupture l i f e (in minutes), e s = 'secondary' strain rate (in per cent per minute). Snead (73) found that in common with most cohesive s o i l s , Haney clay does not exhibit a constant secondary strain rate. However, Snead suggested that Saito and Uezawa's secondary strain rate can be approximated by the transient minimum strain rate, e , for Haney clay. The resulting l inear relationship between log transient minimum strain rate, e , and the log total rupture l i f e , t , is also shown on Figure 2-10. Snead suggested that a similar relationship between t and e exists for other so i l s , r m Snead proposed a method of predicting the time to rupture during the ' ter t iary ' creep stage when the strain rate is increasing. Figure 2-11 shows the log time to rupture, t^ r , (elapsed time from the instant considered until fai lure) - log strain rate, e, curve for a creep test on normally consolidated Haney clay. On this plot , the test progresses from right to l e f t , as ind i -cated, so that there are two occasions when the strain rate may have a certain value. Thus, the minimum time l e f t until rupture wi l l occur can be determined for any strain rate during the increasing strain rate stage to the lef t of the transient minimum strain rate. Results from al l creep rupture tests on Haney clay when the strain rate is increasing are shown on Figure 2-12. From these results , Snead concludes that for Haney clay there is a l inear relationship between the logarithm of the strain rate and the logarithm of the remaining time to rupture when the strain rate is increasing, which is independent pf stress leve l , consolidation history and drainage conditions. The equation of 10* 5? 10 •UJ UJ < cc < cc I -<n _i < x < 10 \ N RUPTUI NORM HAN ALU EY f CONSOL CLAY .IDA- ED rARi r OF TEST 42 I00 10' IO5 FIG. IO1 IO2 IO3 TIME TO RUPTURE , MINUTES 2-11 TIME TO RUPTURE-STRAIN RATE CURVE FOR HANEY CLAY. (73) IO2 10 3 10 4 TIME TO RUPTURE , MINUTES FIG. 2-12 TIME TO RUPTURE-STRAIN RATE RELATIONSHIP. (73) 43 this relationship i s : log t t r = .23 - log 1 Q £ or t (2:30) - L l tr " : where t t = time to rupture from time t (in minutes) e = accelerating strain rate at time t (in per cent per minute) ' Singh and Mitchell (72) have also observed that for a given cohesive s o i l , the value of et which causes creep rupture is fa i r l y independent of deviatoric stress level so that (^t)^a i-iu appears to be a material property. At any given time, an element of soi l under a sustained deviatoric stress wil l have an instantaneous value of et which may increase (m<i in Singh and Mitchell 's empirical equation), remain constant (m=i), or decrease (m>i) as time elapses. When the value of et reaches (et) f a . j i for m<i, creep rupture may be ant i -cipated. It can be seen that this method of predicting the time of slope fai lure is very similar to Snead's method. Mallawaratchie (45) has examined the effect of type of creep test and consolidation history on the creep rupture l i f e of Haney clay. The lines of minimum strain rates for the increasing and a g decreasing plane strain creep rupture tests were found to be very close and can be taken as one straight l ine . The resulting l inear relationship between logarithm of total rupture l i f e and logarithm of minimum strain rate for al increasing and a 3 de-creasing plane strain creep tests taken together is shown on Figure 2-10 along with the other total relationships. Results for four KQ t r iaxia l creep tests performed by Mailawaratchie^are also shown on Figure 2-10. From these tests i t can be seen that the use of data from conventional t r iaxia l creep tests may be extremely unconservative at 44 lower times since the predicted time to fai lure may be 3 to 6 times greater than that obtained from the data for the more representative plane strain creep tests. The available data for K0 t r iaxia l creep tests indicates a l inear relationship between log transient minimum strain rate and log total rupture l i f e that is parallel to and l ies below the l inear relationship for conventional t r iaxia l creep tests. However, the comparisons shown in Figure 2-10 may not be valid in that the strain rate used in the tr iaxia l creep tests is the axial strain rate, whereas the strain rate used in the plane strain creep tests is the vertical strain rate. Mallawaratchie's computations show that the octahedral shear strain at the minimum strain rate varies from 3.4 to 7.65% for samples in conventional t r iaxia l creep tests compared to 0.9% for plane strain creep tests and 0.45% for KQ t r iax ia l creep tests. The corresponding axial (vert i -cal) strains are: 2.4 to 5.4% for conventional t r iax ia l creep tests; 0.6% for plane strain creep tests; and, 0.3% for KQ t r iaxia l creep tests. Vaid (79) has obtained similar values for the axial (vertical) strain at fai lure based on a maximum deviator stress fai lure condition for conventional t r i a x i a l , plane s t ra in , and KQ t r iaxia l tests run at constant strain rates similar to the values of the minimum strain rate in the creep tests. Therefore, the strains at the minimum strain rate appear to depend on the particular type of creep test. If a minimum strain rate is to be used in a total rupture l i f e relationship such as Equation 2:30, i t is extremely important that the test conditions ref lect the actual in s i tu conditions. It should be noted that the ratio of the total rupture l i f e to the time until the minimum transient strain rate is reached is about 2.5 to 3 for plane strain creep tests, 3 to 4 for conventional t r iaxia l creep tests , and 4 to 5 for Kn t r iaxia l creep tests (45). Thus, most of the total rupture l i f e is 45 spent during the tertiary creep state at progressively increasing strain rates. In the nonlinear creep analysis, a total rupture l i f e relationship wi l l be used in conjunction with Singh and Mitchel l 's (^)f a-j-| u r e condition. For the multiaxial stress cases, the principal shear strain rate wil l be used in the various rupture l i f e relationships. The evaluation of upper y ie ld strengths for cohesive soi ls as suggested by Sherif (69) and Snead (73) may allow the determination of working stress levels at which creep rupture wil l not occur. If i t is known that the,stresses in an earth structure are below the upper y ie ld strength, then creep rupture is not anticipated and i t is not necessary to consider a time to potential rupture. It is hoped that further research in this area wil l lead to £he determination of upper y ie ld stresses for a variety of cohesive soi ls for in s i tu conditions. The solution methods for creep analysis developed in Chapter 4 are f lexible enough to handle such a fai lure c r i te r i a when i t is developed further. 46 CHAPTER 3 LINEAR VISCOELASTIC STRESS ANALYSIS USING FINITE ELEMENTS AND THE CORRESPONDENCE RULE 3-A INTRODUCTION The constitutive equations for a cohesive soi l assumed to be l inear visco-elast ic developed in Section 2-D-l are a system of ordinary l inear differential equations which are to be solved in conjunction with the kinematic equations, equilibrium equations, and boundary conditions for a particular problem. The solution of these equations is often simplif ied by using the Laplace trans-formation, which enables many e las t i c i t y solutions to be used in l inear visco-elast ic stress analysis. Consider a constant axial compressive stress a to be applied at time t=o to an undisturbed axial test specimen that is e last ic in material behaviour. The axial strain e is then given by: •; e " E (3:1) where E is Young's modulus. If the test specimen behaves as a Kelvin Solid (Equation 2:6), the axial strain e(t) is given by: ( K l E ( t ) + n xe(t) = a(t) (3:2) where Kl is the spring stiffness and r\1 is the dashpot viscosity for the model. Taking the Laplace transform of Equation 3:2 y ie lds : . 7 ^ - s ( K l ° V ) = ( 3 : 3 ) where s is the Laplace transform parameter and bars are used to indicate transformed expressions (60). The solution of Equation 3:3 for e(t) is found from a table of inverse Laplace transforms (44). -Kit e(t) = f- (l - e n l ) (3:4) 47 It can be seen that Equation 3:3 corresponds to Equation 3:1 in which the stress has been replaced by its Laplace transform and E by the transformed parameters — ^ , where P and Q are l inear dif ferent ia l operators. This P(s) elast ic-viscoelast ic analogy allows many uniaxial v iscoelast ic i ty problems to be solved by the following correspondence rule: i f , in the solution for a problem in e last ic i ty in which the stress is applied to the body at time t=o when i t is undisturbed, the stress is replaced by i ts Laplace transform, and E by — ^ , the Laplace transform of the solution to the corresponding visco-P(s) e last ic problem is obtained. In order to use this correspondence between e last ic and viscoelastic. solutions, the solution to the e las t i c i t y problem is required. However, for many problems in soi l mechanics the closed-form e las t i c i t y solution is not available. For such problems, the f in i te element method is part icularly useful. An outline of the history of the f in i te element method in applied mechanics and soi l mechanics wil l not be attempted here as there are many excellent reviews available in the l i terature (88, 62). 3-B CORRESPONDENCE RULE The elast ic-viscoelast ic correspondence, which has been outlined, is the most commonly used method for obtaining the displacements and stress d i s t r i -butions in bodies that are assumed to be l inear viscoelast ic. Usually referred to as the correspondence rule, i t appears to have been f i r s t formu-lated by Alfrey (1) using the separation of variables technique, and by Read (63) using Fourier integral transforms. It was subsequently shown by Lee (41) that the correspondence rule could be derived by the use of the Laplace trans-form. This was deduced for isotropic materials, but later extended by Biot (2) for anisotropic materials. Biot also indicated that this method of 48 solution is applicable to variational methods of approximate e last ic analysis. Schapery (66) has also shown the basis from variational principles and thermo-dynamics, and developed methods that can be used with anisotropic and inhomo-geneous materials. Some very useful approximate inversion techniques also developed from this work. The general development of the correspondence rule wil l not be given here since i t is covered in detail by many authors such as Corneliussen and Lee (16) and Flligge (21). The use of the correspondence rule in conjunction with f in i te element solutions is suggested by Zienkiewicz and Cheung (87) and a method of solution outlined by Webber (83). 3-C MATERIAL CHARACTERIZATION The isothermal stress-strain relations for a l inear viscoelast ic material developed in Section 2-D-l have been given in the l inear dif ferent ia l operator form which can be associated with a rheological model represented by a combin-ation of springs and dashpots. Writing these stress-strain relations out: [ a p , ^ r + + a„] S^ .^ . t ) = [ b q l ^ + ... + b 0 ] e ^ x ^ t ) * (3:5) and, . p " q" [V'~^ + + c ° ] ^ V * * = [ V V + + d o ] ( 3 : 6 ) dt d t where the material constants ap,, b^,, Cp„, d^,,, are determined from.experi-mental measurements, for example, by f i t t i ng the stress response to specified bulk and shear inputs, or by extending uniaxial test data in the case of a cohesive s o i l . The Laplace transforms of Equations 3:5 and 3:6 are: [b , s q ' + + b n ] 1 J k [a p,sP + . . . + aQ] 1 J k 49 *"(x*'s) = i X l ^ l — T 7 T e - - ( x k ' s ) ( 3 : 8 ) [ C p i i S K + ••• + c 0 J which again reveals the correspondence of the viscoelastic stress-strain r e l -ations in the transform plane to the e last ic stress-strain relations. The purpose of material characterization is to determine the nature of G(s) and K(s) or E"(s) and v"(s) in a convenient form that can then be used in the stress analysis. The algebraic relations in s given by Equations 3:7 and 3:8 become identical with their e last ic counterparts i f : £ M , 2G (3:9) | M - SK : (3:10) Using the relationships: G = » K = - 7 - ^ — r (3:11) 2 ( l + v ) 3 ( l - 2 v ) v ' r 9KG 3K-2G , 0 n o x then, 3K+G * v = IK^G ( 3 : 1 2 ) so that F(s) and \T(s), in terms of quotients of operator polynomials, become: r < s > • 2p#Tw ( 3 : 1 3 ) »<S> - PF^F (3:14) The correspondence rule can now be stated as follows: If the solution of an e las t i c i t y problem is known, the Laplace transform of the time-dependent functions in the corresponding v iscoelast ic i ty problem may be found by replacing the e last ic constants K and G according to Equations 3:9 and 3:10, or by replacing the e last ic constants E and v according to Equations 3:13 and 3:14, and the actual time dependent loads by their Laplace 50 transforms. Thus, the transformation has reduced derivatives and integrals with respect to time to algebraic expressions of the transform parameter s. The expressions that result are analogous to the constitutive, equilibrium and kinematic equations, and boundary and i n i t i a l conditions that govern the be-haviour of an e last ic body of the same geometry as the viscoelastic body. If the solution to this associated e last ic problem can be determined, the solution to the time-dependent viscoelastic problem can be obtained by operating upon the stresses and displacements of the associated e last ic problem with inverse Laplace transforms. Often, i t is assumed that viscoelastic materials are essential ly incom-pressible, which for small strains is equivalent to assuming an in f in i te bulk modulus or a Poisson's ratio of one-half. Viscoelastic shear behaviour is s t i l l permitted, therefore: K(s) = » v(s) = \ E(s) = 3G(s) (3:15) Because the precise measurement of compressibility is d i f f i c u l t , a second approximation is possible. A f in i te value of the bulk modulus is permitted, but any time-dependency is neglected. The actual time-dependent bulk modulus is replaced by an average (elastic) value. Again, viscoelastic shear behaviour is permitted: 3K - 2G(s) _ 9K G(s) K(s) = K v(s) = — =7-^ E(s) = „ e ^ . (3:16) e 6K e + 2G(s) 3Kg + G(s) With detailed materials information, the assumption that both the bulk and shear moduli are viscoelastic can be used in the analysis. Equations 3:13 and 3:14 are then applicable. The viscoelastic material is s t i l l assumed to be l inear , isotropic and homogeneous. It should be noted that the final step of inverting the transformed 51 solution may often prove to be extremely d i f f i c u l t and may require approximate techniques. Tables (44) are available that allow for the inversion of the transformed expressions that result from assuming rheological models with up to 7 or 8 elements. A model with this number of elements should be more than adequate to describe the behaviour of most l inear , isotropic viscoelastic materials. 3-D LINEAR VISCOELASTIC EXAMPLE PROBLEM Unfortunately, there are no typical 'simple' soi l mechanics problems for which closed-form analytical solutions are available to i l lust rate the use of the correspondence rule. However, the l inear viscoelastic plane strain ana-ly t ica l solution for the radial displacement and stresses of the thick-walled cylinder shown in Figure 3-la can be found quite readily using the correspon-dence rule (21). This example can also be visualized as representing the possible behaviour of a thick concrete tunnel l iner in a very soft clay. From e las t i c i t y , the radial displacement of the cylinder is given by: u(r) Pb2 3r l a' (a 2-b 2) |^K+2G 2G r (3:17) and the stresses by: Pb2 a 2 -b 2 Pb 2 a 2 -b 2 1 + (3:18) The stresses are independent of the moduli, and hence the same for both the e last ic and l inear viscoelastic materials. Applying the correspondence rule to Equation 3:17 for the radial displacement, and replacing K and G by their equivalent viscoelastic moduli: FIG. 3-la CROSS SECTION OF THICK-WALLED CYLINDER PROBLEM. t i l K MAXWELL SHEAR MODEL ELASTIC DILATATIONAL MODEL cr + p,cr =q,6 S j j ( t ) + ^ Sij(t) = 27 ? l e j j ( t ) cr.. = 3K£.. II II FIG. 3-1 b MODELS OF VISCOELASTIC RESPONSE FIG. 3-1 EXAMPLE PROBLEM-THICK-WALLED CYLINDER (21) 53 u(r,s) = Pb2 (a 2-b 2) , P'(s) a 2 2Q"(s)P(s) + P"(s)Q'(s) Q'(s) r 3rP ' (s)P"(s) (3:19) Considering the linear viscoelastic material to follow the Maxwell law in shear and to be e last ic in dilatation as shown by Figure 3:1b: n l c J i j + G7 Sij = ^ i e i j 3Ke, in shear, in d i la tat ion, (3:20) (3:21) (3:22) n n where n is the viscosity. Then, in terms of the Laplace parameter: P"(s) = i , Q"(s) = 3K, P'(s) = l + p ^ , Q'(s) = q xs where P, = f - and ql = 2 n 2 . l The internal pressure is assumed to be a step function loading: (3:23) where H(t) is the Heaviside step function. The Laplace transform of Equation 3:23 is given by: P(t) = PQH(t) P(s) = Substituting Equations 3:22 and 3:24 into 3:19. (3:24) u(r,s) P o b 2 i (a 2-b 2) s 3r(i+p 1s) 1+PjS g 2 6K(i+p1s) + q 2s q xs r (3:25) and taking the inverse of Equation 3:25, the time-dependent radial displacement is given by: u(r,t) = Inv (u(r,s)) = p o b 2 (a 2-b 2) f r 2K l -^ p ^ -6Kt e6Kp 1+q 1 (t+p x) (3:26) Values of u(r,t) are then readily tabulated for this viscoelastic model by using a computer program. This same problem wil l be examined later using the f in i te element iterative method for an accuracy check. 54 3-E FINITE ELEMENT METHOD The solution to the e last ic problem for f in i te elements has been derived in Appendix A-l using the displacement method and is given by: [K] N {Q>N = {P>N (3:27) where the element stiffness matrix [K]^ = / [A~ 1]J [BlJ [Dl N [B] N [A~ 1 ] NdV, in which [ A _ 1 ] N is the transpose matrix, [B] N is the strain-displacement matrix, [D]^ is the constitutive matrix, {Q}^ is the vector of unknown nodal displace-ments and {P}^ is the nodal loading vector, for the element. Combining al l of the elements for the particular structure y ie lds : [K] {Q} = {P} . (3:28) where [K] is the structural st iffness matrix assembled from al l N elements, {Q} is the vector, of unknown nodal displacements and {P} is the nodal loading vector for the structure. The unknown nodal displacements for the e last ic problem are then given by: {Q} = [K - 1 ] {P} (3:29) where the st i f fness matrix has been inverted. The structural st iffness matrix [K] contains Young's modulus and Poisson's ratio as terms within each coefficient of the matrix. Application of the correspondence rule to Equation 3:29 requires the replacement of these e last ic moduli by their viscoelastic equivalents as outlined in the previous section. Unfortunately, the solution of the resulting set of transformed equations.and their subsequent Laplace transform inversion is extremely d i f f i cu l t for other than very simple cases. However, the terms of [K]^ can be sp l i t as in the following manner suggested by Webber (83): [B]J [D]N [B] N = F ^ ] + F 2[X 2] : (3:30 where [X 2 ] , [X 2] do not contain moduli, and ? 1 , F 2 are functions of the moduli. 55 Then, [ K ] N = / [ A _ 1 ] J ( F^XJ + F 2[X Z]) [A - 1 ] dV (3:31) . . . . . . V = F J C K J ^ + F 2 [K 2 ] N (3:32) where [K,] = / [ A - 1 ] J [ X 1 l [ A" : l3 N dV and [ K 2 ] N = / [A " 1 ] J [ X 2 ] [A ' 1 ] N dV v V For the assembled structure, considered homogenous and is t rop ic , (F jK j ] + F2[K2]){Q} = {P} (3:33) where [K^, [K ] are matrices of constants for the structure and F l 5 F2 are functions of E and v. It'should be noted that the coefficients of the matrices of constants and [K 2 J N for any element can be a ratio of the coefficients for any other element. Thus, elements can be considered to have different stiffnesses without changing the functions of the moduli, Fx and F 2. In this manner, nonhomogeneous elements can be assembled for the structure to give Equation 3:33. Equation 3:33 can be rewritten as: [K x] {Q} = ^ [K 2] {Q} (3:34) A procedure using matrix iteration that has been outlined by Varga (80) wil l be adopted to solve Equation 3:34. Consider {Q}} as a f i r s t approximation for the unknown nodal displacement vector. Then, a new value of the nodal d is -placement vector that is closer to the correct vector can be determined by substituting {Q}1 in the right side of Equation 3:34: i . L" K I ]W 2 = 7 7 - F7 [ K 2 ^ < 3 : 3 5 ) The new value of the nodal displacement vector {Q}2 is then compared with the assumed vector {Q} l S and, i f the agreement is not close enough, a new vector is determined to substitute in the right side of Equation 3:34 as follows: {Q}3 .= (Q}2 + A({Q}2 - {Q}^ (.3:36) 56 where A is an over-relaxation factor. These iterations on {Q} are carried out until the mean square error is less than some acceptable value. This matrix iterative procedure is usually given in the following form: (3:37) where the subscript on {Q} refers to the cycle of i terat ion. Applying the correspondence rule to Equation 3:37 and substituting the equivalent viscoelastic moduli for E and v in F 1 and F 2 in accordance with the correspondence rule gives': i F 2(s) [K1]{Q}r { P } " f ^ y C K 2 ] { Q } n n + i . F 2(s) T (*\ L 2-where bars are used to indicate transformed expressions, nodal load vector to be applied at t=o, {P} = {P}Q H(t) where H(t) is the Heaviside step function. Then, l {P} =. {P} o s (3:38) Considering the (3:39) (3:40) and Equation 3:47 becomes:. {P}0 F (s) s Fi(s) 2 n Fi(s) L " 2 Taking the Laplace inverse of 3:41 gives the following expression.for the viscoelastic displacements, l (3:41) [K1]{Q(t)> n+i = Inv sF^s ) {P>, [K 2] Inv F 2(s) _ FTTiT^n (3:42) where the time-dependent representation of {Q (t)> n is required to determine'{Q}. It is not desirable to invert [K : ] since this requires a large amount of computer storage and time. However, Equation 3:42 can be rewritten as: [K1](Q(t)} n+i = {W} (3:43) where {W} Inv ( s F ^ s ) {P} 0 - [K 2] Inv F 2(s) _ J- {Q} F x(s) " 57 (3:44) and since [K x] i s a banded symmetric matrix, the square root (Choleski) method of solution (60) can be used where {W} i s treated as a new load case for each i t e r a t i o n . Thus, only one determination of the lower transpose of [K x] i s required for the problem. The matrix i t e r a t i o n procedure i s repeated for each time step required in the solution of the problem u n t i l the f i n a l desired total time is reached. This leads to the problem of the representation of the time-dependent displacements {Q(t)> n required to determine {Q>n for ite r a t i o n s at each solution step. 3-F METHODS OF REPRESENTING {Q(t)> n Consider a displacement u(t) which is one of the components of {Q(t)}. The f i r s t piecewise l i n e a r approximation of u(t) employed was a sum of step functions (44) as shown in Figure 3-2a. However, this approximation of u(t) was found to be i n f e r i o r to the piecewise l i n e a r approximation of u(t) as a sum of ramp functions. This approximation is shown in Figure 3-2b, where the i n i t i a l time increment Atx is to be very small so that the i n i t i a l displace-ment increment AU x is representative of the instantaneous response at time t=o, For this piecewise l i n e a r approximation of u(t) as a sum of ramp functions, u(t) = I I j= i AU. A U J H(t-t j._ i) - W{t-t-) AU + AU. H(t-t.) jA ( t-tj^) Hft-t .^) - J (t-tj) H(t-tj) + MjHlt-t j ) "J J (3:45) Taking the Laplace transform of Equation 3:45 gives: FIG. 3-2 PIECEWISE LINEAR APPROXIMATIONS FOR u(t). 59 oo AU. J=l J -1 _ x S i -t,-S e J - e J 4 + >^ t,- - 1 + At.e J -J s I AU- - t . S -t-S _ J _ L ( e J i _ e J ) At. c2 0=1 J S (3:46) Substituting this expression for the components u(s) of {Q} into the right term of Equation 3:44, fir Inv F 2(s) Since, ^ ( s ) Inv {Q} = Inv (e J 1 - e J ) F,(s) l s 2 p i ( s ) (3:47) -t,s f(s) = f ( t - t j ) H(t-tj) (3:48) Equation 3:47 then becomes: T 2 ( s ) Inv F ^ s ) {Q} oo {AU} . = I 1 3=1 A t j F 3 ( t - t . x) H(t-t. x) - F 3 ( t - t . ) H(t-t.) where F (t) = Inv F,(s) ^ ( s ) (3:49) (3:50) and {AU} is the' vector of the nodal displacement increments. Equation 3:44 can then be expanded for time step M to: CK 1](Q(t M)} n+i = F j t M ) { P } Q - [K 2] M-l {AU}. y J .L At. J = l J. W W " F 3 < W this term for only M>i {AU}^ ~AtT F 3 ( A t M ) - F 3(o) (3:51) a l l 60 M-i {AU}. For time step M, I A f J is known and (Au>M is required such that the difference { Q ( t M ) > -At. J=i J 'M-l I {AU}. +' {AU} M U = l J is within acceptable l imits . Since changes for each time step, i t would usually be necessary to store M-i sets of {AU}. However, this large storage requirement can be much reduced. ! F 2(s) F, and F0 are differential operators so that .— , . w i l l have the general 1 1 s 2 F 1 (s ) form: F 2 ( S > C / + . . . . + C X S + C Q S ^ s ) S 2 ( d n S n + + d x S + dQ) (3:52) where c and d are constants. When m^(n-i), which is the case for the rheo-logical models considered in Section 2-D-l, Equation 3:52 can be readily inverted by partial-fraction methods using the Heaviside expansion theorems (60), giving an expression of the general form: Inv ( F 2(s) -t -t -t = a n t + b n + b.e / T l + b,e / x 2 + - . . . + b . e / x i (3:53) s 2 F 1 (s ) where a 0 , b., are constants for the model selected. Considering the f i r s t 4 terms of Equation 3:53 for convenience: - t M-i {AU}. y — i •_ At-J-1 J [ " W W " F 3 ( V V M-i {AU}. • L At. J=l J *M a 0 A t j + b i e (e t J" 1 / t 1 ^ 7 , e l ) + b2e (e V i A , 2 _ t . , j / n ao { Q ( W } + J e t . , M-i {All} . J * 1 ; T . y 1 fe - 1) 61 (3:54) Therefore, i t is necessary to only store a Q{Q(t M_ 1)} and I sets of M-i {Au}. t j - i / t j / . b. I , J (e T l - e T l ) , where I depends on the model selected and j=i J M-i {Q(t M .)} = I {AU} . . The expressions for F (t) and F (t) are now required. 3-G EVALUATION OF THE REQUIRED TRANSFORMS The inverse transforms required in the solution method are given by Equations 3:50 and 3:52. Fj and F 2 f o r the plane s t r a i n and axisymmetric cases are given by: , E eK + 2G Ev and ( l + v ) ( l-2v) 3K - 2G 6 K + 2G ( i + v ) ( i-2v) = V = (3:55) (3:56) Then, replacing K and G by t h e i r v i s c o e l a s t i c equivalents (Equations 3:9 and 3:10), l _ 3P~" (s )P ' ( s ) s F ^ s ) s(2Q"(s)P'(s) + P"(s)Q'(s)) F 2(s) _ P'(s)Q"(s) - P"(s)Q'(s) (3:57) (3:58) s 2 F ! ( s ) s 2(2P'(s)Q"(s) + P"(s)Q'(s)) The expressions for F 3 ( t ) and F 4 ( t ) w i l l then depend on which reheological model is substituted into Equations 3:57 and 3:58. Considering the same ex-ample of the thick-walled cylinder that was solved a n a l y t i c a l l y e a r l i e r and substituting the Laplace transformed terms for a material that follows the Maxwell model in shear and is e l a s t i c in d i l a t a t i o n , as given by Equation;3:22, into Equations 3:57 and 3:58 y i e l d s : 62 F,(t) = Inv ( 3 K P l - q i ) (sKp^qj) 3K 3KPi-qi + s 2 ( 6K + s ) I 6Kp 1+q 1 and -6Kt t - ^ ( i - e 6 K P i + q i ) F,(t) = Inv 3P, ( 6 K P l + q i ) — + s s( 6K + s) v6Kp.+q. ' I 1 1 J (3:59) l 2K -6Kt 11 e 6 K P i + q i 6Kp i+q i (3:60) •t/ It can be seen that Equation 3:59 i s of the general form a t + b Q+ b^e i and the expressions derived in Equation 3:54 are v a l i d . The f i n i t e element i t e r a t i v e method given by Equation 3:51 yields the time-dependent displace-ments, and using the strain-displacement relationships given in Appendix A - l , the time-dependent strains are determined. A method of evaluating the stresses from these time-dependent strains is required. 3-H EVALUATION OF TIME-DEPENDENT STRESSES The element strains are a l i n e a r function of the nodal displacements that have been determined. Thus, the time-dependent shear stresses are given by: M-i {Ae..}„ S i i ( t ) = I -It-1 1 3 1=1 A t £ F5 ( tM " W " ^ M " V {Ae,,} i j J M At, M F 5 ( A t M ) - F 5 ( 0 ) where Inv Q'(s) s 2p'(s) (3:61) for the pa r t i c u l a r reheological model selected. (3:62) 63 For the models selected, the d i l a t a t i o n a l response is assumed to be e l a s t i c , and therefore independent of time. Then, at any time, the stresses are given by: ,,(t) = S,,(t) + \ 5,,c (3:63) ' U l " " 3 i j kk Considering the material that follows the Maxwell model in shear and is e l a s t i c in d i l a t a t i o n again, Fc(t) = Inv s P l ( ^ + s ) Pi - t / (3:64) = q x ( i - e ^ l ) (3:65) From the expressions for the time-dependent displacements, strains and stresses that have been developed, i t is possible to check the accuracy of the f i n i t e element i t e r a t i v e method. 3-1 ACCURACY CHECK ON THE FINITE ELEMENT ITERATIVE METHOD An axisymmetric f i n i t e element program and a plane strain-plane stress f i n i t e element program were written to solve problems in which the s o i l ' s be-haviour can be considered to be time-dependent. The assumption of i s o t r o p i c , l i n e a r v i s c o e l a s t i c material behaviour and the solution method that has been developed in this chapter is one of the several solution options that are permitted in these general programs. This l i n e a r v i s c o e l a s t i c solution portion of the plane strain-plane stress f i n i t e element program i s outlined in the f i r s t part of the s i m p l i f i e d description of the program inputs and c a p a b i l i t i e s given in Appendix B-2. The axisymmetric f i n i t e element program is quite s i m i l a r , but does not have a l l the solution options of the plane strain-plane stress f i n i t e element program, which has wider application in s o i l mechanics problems. The various options incorporated in the programs, such as i n i t i a l 64 load conditions, are given in Appendix B - l . The fundamental difference between the l i n e a r v i s c o e l a s t i c solution portion of the f i n i t e element programs developed here and the usual f i n i t e element programs available in the l i t e r a t u r e is the generation of each element s t i f f n e s s matrix and the overall structure s t i f f n e s s matrix as two separate matrices as outlined previously. Also, the i t e r a t i v e method of solution i s used to solve the basic s t i f f n e s s relationship (Equation 3:43) at each desired time step to determine the time-dependent displacements, strains and stresses. Each problem i s solved for the assumption of l i n e a r e l a s t i c material behaviour in a standard f i n i t e element option of the program before s t a r t i n g the l i n e a r v i s c o e l a s t i c solution option of the program. The e l a s t i c moduli are selected to give the same displacements as the time t=o displacements for the visco-e l a s t i c problem. This e l a s t i c solution for the nodal displacements i s then used as the f i r s t assumed value of the nodal displacements in the i t e r a t i v e solution for the f i r s t time step of the v i s c o e l a s t i c problem. Thus, the gener-ation of the various s t i f f n e s s matrices is checked and a large number of i t e r -ations for the f i r s t time step is avoided. At each successive time step of the v i s c o e l a s t i c solution, the nodal displacements determined for the previous time step are used as the i n i t i a l values for the i t e r a t i o n cycle. The l i n e a r visco-e l a s t i c solution option has been limited to compressible triangular elements. However, with some modifications, quadrilateral and incompressible elements (85, 12) can be included in the l i n e a r v i s c o e l a s t i c solution option. The f i r s t step in the checking of the l i n e a r v i s c o e l a s t i c solution option was to determine the value of the over-relaxation factor A in Equation 3:36 that gives the solution at each time step in the least number of i t e r a t i o n s . The mean square error in any i t e r a t i o n cycle is defined by the following equation: 65 M I (u. 1(n+i) MSE = i=i (3:66) where MSE stands for the mean square err o r , u. are the individual nodal displacement, and the subscript n refers to the i t e r a t i o n number. Iterations on the nodal displacement vector {Q} are performed until the mean square error i s less than some s p e c i f i e d value. The f i n i t e element i d e a l i z a t i o n of the thick-walled cylinder shown in Figure 3-3a was used in the l i n e a r v i s c o e l a s t i c solution option of the axi-symmetric f i n i t e element program to determine the value of A that reduced the mean square error the quickest. This i d e a l i z a t i o n was used for two basic reasons: the number of unknown nodal displacements, 26, is quite small and hence the computation time per step i s small; and, this same i d e a l i z a t i o n with only the dimensions changed has been used by others (28, 29, 75) to check creep programs. At the same time, a check was made with the closed form solu-tion to determine the size of the time steps for the piecewise l i n e a r approx-imation of {Q} that could be used to reach the f i n a l desired total time. > It was found through a series of t r i a l s that an over-relaxation factor A of 1.60 to 1.65 optimized the i t e r a t i o n procedure f o r the material assumed to follow a Maxwell model in shear and to be e l a s t i c in d i l a t a t i o n . An accuracy check on this p a r t i c u l a r test case with the closed-form solution indicated that time steps could be taken such that the new displacement increment was 2 to 4 times greater than the previous displacement increment with reasonable accuracy. However, the use of smaller i n i t i a l time steps gave somewhat better agreement for the same s p e c i f i e d mean square error. With the over-relaxation factor determined, accuracy checks were conducted AXIS OF CYLINDER 26 NODES 24 ELEMENTS 4f f j f y - f i f FIG. 3-3a AXISYMMETRIC FINITE ELEMENT IDEALIZATION < 240 " ^ < 360 " > FIG. 3 - 3 b P L A N E STRAIN FINITE E L E M E N T IDEALIZATION. FIG. 3-3 FINITE ELEMENT MESHES FOR THE THICK-WALLED CYLINDER. 67 for various material assumptions. This was done by comparing the axisymmetric f i n i t e element results and the plane s t r a i n f i n i t e element results with the closed form results for the same material assumptions. The plane s t r a i n f i n i t e element i d e a l i z a t i o n of the cylinder problem i s shown in Figure 3-3b, and this representation, with the large increase in unknown nodal displace-ments, required a much greater solution time than the axisymmetric i d e a l i z a t i o n required. Three test cases are given in Figures 3-4 to 3-6 with the material be-haviour assumption and parameters l i s t e d . These parameters were selected., to give the same time t=o response for each case and might be representative of the values for a s t i f f clay (32). The l e f t portion of each figure i s the i n i t i a l 200 minutes of creep on an expanded time scale. It can be seen from Figures 3-4 to 3-6 that there is excellent agreement in the solutions for a s p e c i f i e d mean square error, MSE, of 1 X 10" 6 in the i t e r a t i o n s , with the accuracy decreasing s l i g h t l y with increasing time. Approximately 2 to 10 i t e r a t i o n cycles were required for each time step of the f i n i t e element s o l -utions depending on the change in displacements. When the s p e c i f i e d mean square error i s decreased to 1 X 10" 1 0, the results a l l f a l l together on the same l i n e as the closed-form solution. However, approximately 5 to 50 i t e r -ation cycles are now required to achieve the s p e c i f i e d mean square error. The cost for the plane-strain solution is about $0.25 per time step. As a n t i c i -pated, the materials that follow the Maxwell model and Burgers Body in shear behaviour (Figures 3-4 and 3-6) show typical unlimited viscous deformation. The material that follows the Standard Linear S o l i d in shear behaviour shows the limited deformation typical of v i s c o e l a s t i c s o l i d behaviour. It can be concluded that the solution methods developed in this chapter r e s u l t in the desired accuracy on the basis of the comparisons with closed-form solutions. >-cr < o Z O ca cc UJ z I - o> Z UJ UJ X S O UJ z 3-_) a. < a < a : 0.75 -0.5 0.25 100 200 TIME , t MINUTES ( FIRST 200 MINUTES) - 20 3 >-DC < a z z> o CO MODEL PARAMETERS v = 0.25 E = I.25x IO4 PSI K=8333 PSI G, = 5000 PSI 77 = l.25x I0 6 PSI-SEC X O LEGEND CLOSED —FORM SOLUTION AXISYMMETRIC IDEALIZATION PLANE STRAIN IDEALIZATION 2 000 4000 6 000 TIME , t MINUTES 8000 10 000 FIG. 3-4 COMPARISON OF ITERATIVE FINITE ELEMENT SOLUTION WITH CLOSED-FORM SOLUTION. THICK-WALLED CYLINDER, MAXWELL MODEL IN SHEAR, ELASTIC IN DILATATION. CO 1.0 >-cc < o z o CQ cc UJ 5 w *I UJ £ c_> < _l CL CO < Q < CC 0.75 0.5 0.25 100 200 TIME , t MINUTES (FIRST 200 MINUTES) 1.0 >-cc <t o z z> o CQ CC UJ 0.75 < i -z UJ Ul a < CO UJ X o 0.5 ^ 0.25 -CL co o _ l < < MODEL PARAMETERS i/ =0.25 E = l.25x IO 4 PSI K = 8333 PSI G, = 5000 PSI G 0 = 5 x l 0 4 PSI ^2 = 5 x l 06 PS I-SEC X o LEGEND CLOSED-FORM SOLUTION AXISYMMETRIC IDEALIZATION PLANE STRAIN IDEALIZATION 2000 4000 6000 T IME, t MINUTES 8 000 10 000 FIG. 3-5 COMPARISON OF ITERATIVE FINITE E L E M E N T SOLUTION WITH CLOSED-FORM SOLUTION. THICK-WALLED CYLIN DER, STANDARD LINEAR SOLID IN S H E A R , ELASTIC IN DILATATION. ^ MODEL PARAMETERS. v =0.25 E = I.25xl04 PSI K = 8333 PSI TIME, t MINUTES TIME , t (MINUTES) ( FIRST 200 MINUTES) FIG. 3-6 COMPARISON OF ITERATIVE FINITE ELEMENT SOLUTION WITH CLOSED-FORM SOLUTION. THICK-WALLED CYLINDER, BURGERS BODY IN SHEAR, ELASTIC IN DILATATION. 71 Examples of the further use of the methods developed in this chapter to study some typical s o i l mechanics problems in which the s o i l ' s behaviour is assumed to be l i n e a r v i s c o e l a s t i c are given in Chapter 5. As a precaution, the e f f e c t of other material behaviour assumptions and parameters on the solu-tion accuracy i s always checked with the simple closed-form example used here. This allows the determination of the size of the time steps and the value of the mean square error that give a reasonable solution accuracy for each problem. 72 CHAPTER 4 NONLINEAR CREEP ANALYSIS USING FINITE ELEMENTS AND THE INCREMENTAL INITIAL STRAIN METHOD 4-A INTRODUCTION In general, creep problems are solved by incremental procedures in which the f i n a l solution is determined by solving a number of l i n e a r i z e d problems dis t r i b u t e d between the i n i t i a l and f i n a l desired time. For cohesive s o i l s , r the components of the creep st r a i n increment AE.. for each time interval are assumed to be given by Equation 2:28 which i s repeated here for convenience: where a = the equivalent stress (— T .) e ^ V oct The assumption of no volume change due to creep strains i s included in Equation 4:1. During creep, the stresses, a.., at any point in an earth structure w i l l 3 be changing with time and can be represented by the typical smooth stress curve r shown in Figure 4-1. The creep s t r a i n increments Ae.. depend on the time and 13 stress history and w i l l be assumed to follow the strain-hardening cumulative creep law that is shown schematically in Figure 4-2. The analysis i s simpli-f i e d by considering incremental creep strains in a small time interval At and by considering the typical stress curve in Figure 4-1 to be made up of many l C the equivalent creep s t r a i n increment (~&YQCt)> which is determined by substituting aQ into the desired uniaxial creep relationship. the stress deviator tensor. FIG. 4-1 APPROXIMATION OF A SMOOTH STRESS-TIME CURVE BY A SERIES OF STEPS. TIME CASE !•• DECREASING STRESS CREEPCURVEIORSTRES? oTsTclwr INTERVAL CREEP CURVE_FORJTRESS OF FIRST INTERVAL Us STRAIN AT END OF SECOND TIME INTERVAL STRAIN AT END OF FIRST TIME INTERVAL CASE ||: INCREASING STRESS FIG. 4-2 SCHEMATIC DIAGRAMS SHOWING HOW TO DETERMINE THE INCREMENT OF CREEP STRAIN FOR CHANGING STRESS CONDITIONS (49). 74 incremental steps. Each step consists of a constant stress period At followed by an instantaneous increment of stress Ac. .. The error introduced through keeping the stresses constant in the time interval At decreases as t he time interval i s decreased. At the beginning of each time i n t e r v a l , the stresses, e l a s t i c strains and creep strains w i l l be known by virtue of the calculations for the preceding time i n t e r v a l . The e l a s t i c solution of the p a r t i c u l a r problem i s usually the s t a r t i n g point. Thus, the incremental procedure can be summed up as five basic steps: 1. Obtain the value of the new equivalent stress cQ from the stresses a.. . determined at the end of the preceding time i n t e r v a l . C 2. Calculate the new equivalent creep st r a i n increment Ae g by substituting the equivalent stress determined in Step 1 into the desired uniaxial nonlinear creep relationship. Since the equivalent stress w i l l usually change for each time i n t e r v a l , the strain-hardening cumulative creep law shown in Figure 4-2 is used to determine the f i c t i t i o u s time t f at the beginning of the i n t e r v a l . Then t ^ , At and a g can be used in the creep relationship to determine Ae^. e 3. Using Equation 4:1 with the stresses determined at the end of the preceding time i n t e r v a l , the equivalent stress aQ determined in Step 1 and the equivalent creep s t r a i n C increment Ae g determined in Step 2; calculate the creep r s t r a i n increments Ae.. for the i n t e r v a l . r 4. Addition of the creep s t r a i n increments Ae.. to the creep 75 strains at the beginning of the time interval gives the creep strains at the end of the i n t e r v a l . 5. Using the creep s t r a i n increments determined in Step 3, together with the constitutive equations, kinematic equations, equilibrium equations, and boundary conditions for the p a r t i c u l a r problem; calculate the increment of stress ha., at the end of the time i n t e r v a l . This procedure is repeated for each time i n t e r v a l . An obvious advantage of this incremental solution procedure i s that a nonlinear e l a s t i c material can be considered by a piecewise l i n e a r i z a t i o n of the s t r e s s - s t r a i n curves at each time i n t e r v a l . If desired, the solution accuracy can be improved by iterat i o n s using the stresses determined at the end of the time interval or the average stress during the time interval to s t a r t through the steps again. Generally, the time intervals can be selected small enough to y i e l d the desired accuracy without i t e r a t i v e methods. It would appear that an analytical formulation of the incremental solution of creep problems described above was f i r s t presented by Mendelson, Hirschberg and Manson (47) in 1959. Their work was an extension of the method of successive approximations used to compute the p l a s t i c flow in plates, cylinders, and spheres (48, 31) and was limited to examining the uniaxial case of a f l a t plate and the bi a x i a l case of rotating discs. Several examples of the ap p l i -cation of the incremental solution method are discussed in Lubahn and Felgar's book (43). Lin (42) has also used this method to analyse the bending of a plate displaying nonlinear s t r a i n hardening creep. By tr e a t i n g the creep str a i n increments as i n i t i a l loads in an e l a s t i c problem for each time i n t e r -v a l , the numerical solution for the stress increment in Step 5 was much si m p l i f i e d over the method given by Mendelson, Hirschberg and Manson. 76 Except f o r a few special cases with simple boundary conditions, the deter-mination of the stress increment for each time interval i s extremely d i f f i c u l t in any analytical formulation of the creep problem. However, the f i n i t e element formulation of creep problems, which allows more complicated problems to be considered, has been investigated by many researchers (59, 28, 29, 75) in the l a s t ten years. Most of these researchers used d i r e c t i t e r a t i o n proce-dures in which the f i n a l solution is found d i r e c t l y by i t e r a t i v e techniques. The f i n i t e element formulation of the nonlinear creep problem and incremental i n i t i a l s t r a i n procedure used by Greenbaum (28) and Greenbaum and Rubinstein (29) has the advantage of being more general and providing a description of the intermediate creep stages. Greenbaum and Rubinstein developed t h e i r f i n i t e element method to study the nonlinear creep of axisymmetric structures made of concrete. The writer extended t h e i r basic method to solve plane . s t r a i n and plane stress problems in addition to axisymmetric problems,-parti-c u l a r l y problems involving earth structures. Since the programs were developed, Sutherland (75) has indicated the extension of Greenbaum and Rubinstein's work to handle plane st r a i n and plane stress problems. The programs discussed in this chapter include a large number of solution options that are applicable to s o i l mechanics problems (such as creep rupture) that have not been considered by workers in applied mechanics. 4-B FINITE ELEMENT SOLUTION USING THE INCREMENTAL INITIAL STRAIN METHOD , ; The increment of total s t r a i n Ae. . for any time interval is the sum of E r the increment of e l a s t i c s t r a i n A e . . plus the increment of creep s t r a i n A e . - . A e i j = A e i j + A E i j : ( 4 : 2> where the superscript E denotes the increment of e l a s t i c s t r a i n . General 77 e l a s t i c behaviour i s assumed so that the relationship between the increment of stress La., and the increment of e l a s t i c s t r a i n Ae^ is given by: Aa.. = D. .. .(Ae. . - Ae?.) (4:3) U ijk£ v i j U , where P^-^ is the material tensor. The solution to the nonlinear e l a s t i c problem for the f i n i t e element for-mulation and a piecewise l i n e a r i z a t i o n of the material tensor has been derived in Appendix A-2 using the displacement method. When the expression for the increment of stress ACT.. . given by Equation 4:3 is substituted into this solu-tion in Appendix A-3, the f i n i t e element solution method for one time interval of the nonlinear creep problem i s given by: [K]N{AQ}N = {AP}N + {AL}N ; (4:4) where the element s t i f f n e s s matrix [ K ] N = / [A ~ 1 ] J [ B ] | J [ D ] N [ B ] N [ A " 1 ] N dV, VN in which [ A - 1 ] ^ is the transpose matrix, [B]^ is the strain-displacement matrix, [D]^ is the cons t i t u t i v e matrix, {AQ}^ i s the vector of unknown nodal displace-ment increments, {AP}^ i s the vector of nodal load increments, and {AL}^ i s the vector of creep s t r a i n nodal 'load' increments, for the element. The vector of creep s t r a i n nodal load increments is found for each , element from: {AL) N = / [A-' ] i J [B ] J [D] N {A e C} NdV (4:5) V N r where the vector of creep s t r a i n components {Ae }^ i s found from Equation 4:1 and the integration i s over the volume of the element. For nonlinear e l a s t i c materials, the co e f f i c i e n t s of the constitutive matrix [D]^ w i l l depend on the state of stress and s t r a i n in the element for each time i n t e r v a l . The c expressions for the components of {Ae }^ for plane s t r a i n , plane stress and 78 axisymmetric problems are given in Appendix C. {Ae }^ w i l l be calculated at the beginning of each time interval on the basis of the stresses determined at r the end of the previous interval and {Ae }^ i s then assumed to be constant for the time i n t e r v a l . Combining a l l of the elements for the p a r t i c u l a r structure y i e l d s : [K]{AQ} = {AP} •+'{AL} (4:6) where [K] i s the structural s t i f f n e s s matrix assembled from a l l N elements, {AQ} is the vector of unknown nodal displacement increments, {AP} i s the vector of nodal load increments and {AL} is the vector of creep s t r a i n nodal load increments, for the structure. For any time i n t e r v a l , the vectors {AP} and {AL} are known and Equation 4:6 can be solved for the unknown nodal di s -placement increments {AQ} by the square root (Choleski) method of solution (60) which takes advantage of the fact that [K] i s a banded symmetric matrix. From these values of the nodal displacement increments determined for the time i n t e r v a l , the increment of total s t r a i n {Ae}^ can be calculated for each ele-ment using the strain-displacement relationship. The increment of e l a s t i c E C str a i n {Ae }^ is then found from Equation 4:2 since {Ae }^ is known and assumed constant for the i n t e r v a l , and the increment of stress {Aa}^ for each element i s found from Equation 4:3. This process is then continued for each time-interval until the f i n a l desired time i s reached. For l i n e a r e l a s t i c materials, the structural s t i f f n e s s matrix [K] does not change with time and the square root method of solution is extremely e f f i c i e n t since the same lower transpose is used for each time i n t e r v a l . However, for nonlinear e l a s t i c materials, the structural s t i f f n e s s matrix [K] changes with the s t r e s s - s t r a i n history of each element and a new lower transpose is required for each time i n t e r v a l . As a consequence, the solution time is increased by a factor of about 4 for the nonlinear e l a s t i c material. 79 4-C ACCURACY CHECK ON THE FINITE ELEMENT SOLUTION OF NONLINEAR CREEP PROBLEMS The solution method that has been outlined f o r problems in which the s o i l exhibits nonlinear creep is one of the options permitted in the axisymmetric f i n i t e element program and the plane strain-plane stress f i n i t e element pro-gram. This general creep solution portion of the plane strain-plane stress f i n i t e element program is described in the second part of the b r i e f description of the program inputs and.capabilities given in Appendix B-2. The options and features of the axisymmetric f i n i t e element program and plane strain-plane stress f i n i t e element program are given in Appendix B - l . The fundamental difference between the creep solution portion of the f i n i t e element programs developed here and the usual f i n i t e element programs available in the l i t e r -ature is that incremental quantities are dealt with in a series of l i n e a r i z e d solutions for the various time i n t e r v a l s . Some of the important advantages of this incremental approach w i l l be i l l u s t r a t e d by the example problems in Chapter 5. The thick-walled cylinder problem that was previously used to check the li n e a r v i s c o e l a s t i c solution option of the f i n i t e element programs w i l l be used to check the nonlinear creep option. However, the radial dimensions have been modified as shown in the plane s t r a i n quadrilateral f i n i t e element i d e a l i z a t i o n for the thick-walled cylinder given in Figure 4-3. The plane s t r a i n closed form solution of this problem based on the assumption that a l l the creep strains are non-recoverable was developed by Odquist and Hult (56) and graphs showing the solution for p a r t i c u l a r creep relationships are given by Greenbaum (28) and Greenbaum and Rubinstein (29). The creep relationship selected to check the programs i s given by: Eg = 6.4 X 10" 1 8a£' , ft (4:7) 7 2 ELEh,eNrs , 6 3 NODES / V = 0 . 4 9 ° ' 6 "RAD/US Q 2 5 " R 4 0 / U S ™ ' < * - W 4 L L E 0 ^ j J E « r 81 Equation 4:7 i s a p a r t i c u l a r form of the power law f o r the creep s t r a i n , e£ = A a " t m (4:8) e e where A, n and m are constants for the p a r t i c u l a r material and stress and temperature range. A power law of this form has been used by many researchers to describe the creep of various materials such as: concrete at higher stresses and temperatures (5, 22); some types of rock, usually at depths where the pressure and temperature is quite high (36); and g l a c i e r ice (15,58). It is necessary in the incremental procedure for solving nonlinear creep problems to insure that the time intervals selected are small enough so that the solution process i s stable. From the work of other researchers (28, 29, 42, 75) and the writer's experience, i t is known that the incremental proce-Q dure becomes unstable i f the equivalent creep s t r a i n increment Ae g for any element during the time interval exceeds the equivalent e l a s t i c s t r a i n for that element. To obtain the desired accuracy in the solution, i t would appear that the time interval should be selected so that the maximum ratio of the equivalent creep s t r a i n increment to the equivalent e l a s t i c s t r a i n for any element is less than about ~, i.e. Max f. C e E £e V (4:9) Also, a further l i m i t a t i o n is made on any time interval (j+i) that: A t j + 1 < 1 .2 At. (4:10) which has been suggested by Greenbaum and adopted in the programs. For the f i r s t time i n t e r v a l , the e l a s t i c or nonlinear e l a s t i c solution i s known for time t=o and is used as the s t a r t i n g point to determine the f i r s t time i n t e r -val A t r The time interval for further steps i s also found automatically by f i r s t determining A t j + 1 from Equation 4:10 and then checking the ra t i o of the strains given by Equation 4:9 for each element. I f necessary, At. i s reduced 82 unti l the maximum rat i o of the strains is less than The creep relationship given by Equation 4:7 i s assumed to be appropriate for describing the creep of concrete at higher temperatures and pressures (5, 22). Also, the strain-hardening cumulative creep law is generally assumed to be applicable for concrete (28) under such conditions. The concrete i s con-sidered to be l i n e a r e l a s t i c with the moduli indicated on Figure 4-3 and a l l the creep st r a i n is considered non-recoverable. From Equation 4:8 and the strain-hardening cumulative creep law,,the equivalent creep s t r a i n increment for a time interval At i s : Ae^ = Ac£ [ ( t f + A t ) m - t j ] (4:11) where the f i c t i t i o u s time t f i s determined at the beginning of the interval as shown in Figure 4:2 and i s given by: O1 e e (4:12) For the p a r t i c u l a r values of A, m and n given in Equation 4:7, Equation 4:11 becomes: A E ^ = 6.4 X 10" a e At (4:13) and i t can be seen that when the creep relationship i s l i n e a r in time, that the equivalent creep s t r a i n increment w i l l be the same for both the s t r a i n and time-hardening cumulative creep laws. The radial stresses in the thick-walled cylinder for the creep r e l a t i o n -ship given by Equation 4:13 are shown in Figure 4-4. After about 2 hours of creep, a steady state stress d i s t r i b u t i o n has been reached where the cylinder creeps at constant stress. It can be seen that the e l a s t i c solutions and the steady state creep solutions determined by using the axisymmetric program and the plane s t r a i n program (for both triangular and quadrilateral elements) agree - 500 LEGEND CLOSED —FORM SOLUTION (28) 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 RADIUS, r INCHES FIG. 4-4 COMPARISON OF RADIAL STRESSES IN THE THICK-WALLED CYLINDER FOR CLOSED-FORM AND FINITE ELEMENT SOLUTIONS OF THE CREEP PROBLEM. oo c o 84 very closely with the closed-form solutions for the radial stresses shown by Greenbaum. (The agreement with Greenbaum's own axisymmetric f i n i t e element solution, while not shown, is also excellent.) The comparison of the radial creep deformations for the various f i n i t e element solutions shown in Figure 4-5 indicates an excellent agreement with Greenbaum's axisymmetric f i n i t e ele-ment solution. It should be noted that Greenbaum used a higher value of Poisson's r a t i o (0.499) than in the present study (0.49) but a s i g n i f i c a n t difference between the various solutions does not show up. It can be concluded that the solution methods discussed in these sections give very close agreement with closed-form solutions and other f i n i t e element solutions. Some examples of the analysis of creep problems in s o i l mechanics using the methods developed in this chapter and Singh and Mitchell's empirical creep relationship are given in Chapter 5. Creep rupture w i l l also be i n t r o -duced into the analysis by reducing the s t i f f n e s s of any element i n which creep rupture has occurred. 2.0 0 20 40 60 80 100 120 TIME , t MINUTES FIG. 4-5 COMPARISON OF RADIAL CREEP DEFORMATIONS FOR THE THICK-WALLED CYLINDER FOR VARIOUS FINITE ELEMENT SOLUTIONS. oo CHAPTER 5 TYPICAL EXAMPLE PROBLEMS 5-A INTRODUCTION Some typical problems in s o i l mechanics that involve creep and creep rupture are examined in the following sections to demonstrate the various methods of analysis. The cases examined are divided into three groups: problems involving s o i l s considered to be l i n e a r v i s c o e l a s t i c ; problems involving s o i l s described by Singh and Mitchell's empirical creep r e l a -tionship; and problems involving creep rupture. In addition, a case involving rock and a case involving g l a c i e r ice are described. 5-B PROBLEMS INVOLVING SOILS CONSIDERED TO BE LINEAR VISCOELASTIC The problems described in this section are primarily intended to indicate some qua l i t a t i v e results f o r cohesive s o i l s considered to be li n e a r v i s c o e l a s t i c in shear behaviour and e l a s t i c in d i l a t a t i o n a l beha-viour. The various rheological models used in each problem have the same i n i t i a l e l a s t i c response. 5-B-1 Uniform Continuous Loading on a Laterall y Constrained Soil Layer The plane s t r a i n f i n i t e element i d e a l i z a t i o n for a layer of s o i l that is l a t e r a l l y constrained i s shown in Figure 5-la. A uniform continuous surface loading of 1,000 psf i s applied to the s o i l and the time dependent stresses and deformations of the s o i l layer are examined for the various rheological models representing the soil'sshear behaviour given in Table 5-This problem is very s i m i l a r to the v i s c o e l a s t i c cylinder in a r i g i d die that Flugge (21) has examined since no l a t e r a l deformation i s allowed at I T 50' 87 66 NODES 100 ELEMENTS SCALE l"=33.3' FIG. 5-la PLANE STRAIN FINITE ELEMENT IDEALIZATION. UJ u i c_> < _ l C L co -0.1 -0.2 u i UJ u. -0.3ft o r-cc U l > - 0.4 -0.5 1500 CO CO Ul cc I-co u. _ l CO < C L o N CC o X 1000 500 MATERIAL PROPERTIES GIVEN IN TABLE 5-| •X-=5F 10 20 30 TIME , MINUTES 40 X - - X X o — o — o LEGEND BURGERS BODY IN SHEAR. STANDARD LINEAR SOLID IN SHEAR. MAXWELL BODY IN SHEAR. SOIL IS ELASTIC IN DILATATION. =<•> -0-V E R T I C A L STRESS LEVEL =<?> 50 -<5> 10 20 30 T IME , MINUTES 40 50 FIG. 5- lb VERTICAL DISPLACEMENT OF SURFACE AND STRESSES IN THE SOIL LAYER. FIG. 5-1 UNIFORM CONTINUOUS LOADING ON A LATERALLY CONSTRAINED SOIL LAYER. TABLE 5-1 PARAMETERS FOR RHEOLOGICAL MODELS USED IN FIGURE 5-1 BURGERS BODY IN SHEAR 7.2 x l O ^ l b / f t 2 3.0 x l O ' M b-min / f t 2 7.2 x 1 0 5 1 b / f t 2 1.2 x 1 0 5 l b-min / f t 2 STANDARD LINEAR SOLID IN SHEAR . Gn 7.2 x l O ' n b / f t 2 G 2 7.2 x 1 0 5 l b / f t 2 ri2 1 -2 x 1 0 5 l b-min / f t 2 MAXWELL BODY IN SHEAR Gi 7.2 x l O ^ l b / f t 2 m 3.0 x l O ^ l b-min / f t 2 ELASTIC IN DILATATION n i G 2 112 V K E 0.25 1.2 x 1 0 5 l b / f t 2 1.8 x 1 0 5 l b / f t 2 89 any boundary. The parameters for each rheological model were developed from Cakmak and Yau's (8) Burgers Body representation of the short term response of pavements under a i r c r a f t wheel loadings. The d i l a t a t i o n a l behaviour of the s o i l i s considered to be e l a s t i c . From Table 5-1, i t can be seen that each rheological model has the same i n i t i a l response. Upon loading, the v e r t i c a l s t r e s s , a , i s equal to the applied v e r t i -cal pressure, p, (1,000 psf) and the horizontal s t r e s s , a x , i s less than a for the s o i l with the three d i f f e r e n t shear behaviour models as shown in Figure 5-lb. With increasing time, a increases asymptotically to o (and x y hence p) for the Burgers Body and Maxwell Body models since they exhibit f l u i d flow. When a =CT for these two models, the shear stress i s zero x y throughout the s o i l layer and, as expected, no further v e r t i c a l deformation or change in stresses occurs with increasing time a f t e r about 10 minutes for this hydrostatic state of stress. The s o i l with Standard Linear Soli d beha-viour shows v i s c o e l a s t i c s o l i d response so that a x and the v e r t i c a l deforma-tions stop increasing af t e r about 3 minutes and never reach the values asso-ciated with the v i s c o e l a s t i c f l u i d behaviour of the Burgers Body and Maxwell Body. Since the s o i l i s l a t e r a l l y constrained, the l i m i t a t i o n on the f l u i d response was anticipated and this p a r t i c u l a r example served as a further check on the f i n i t e element solution method. Also, the change in volume of the s o i l due to v e r t i c a l deformations was compared with the total volumetric s t r a i n and they were found to be equal as expected. 5-B-2 C i r c u l a r Surcharge Loading on a Soil Layer The axisymmetric f i n i t e element i d e a l i z a t i o n of a c i r c u l a r surcharge loading over part of a s o i l layer i s shown in Figure 5-2a and the parameters for each rheological model representing the s o i l ' s shear behaviour for this 41 90 k - 15' -I RADIUS A 70" $ 114 N O D E S 181 E L E M E N T S FIG. 5-2a AXISYMMETRIC FINITE ELEMENT IDEALIZATION. 0 . 2 5 U l UJ 0 .20 u . H EN 0 . 15 Ul o LA 0 .10 CL in o _ i 0 . 0 5 < a h-cc 0 UJ > M A T E R I A L P R O P E R T I E S G I V E N IN T A B L E 5 - 2 L E G E N D B U R G E R S BODY IN S H E A R x- - x- , _ x STANDARD L I N E A R SOLID IN S H E A R . o _ o M A X W E L L B O D Y IN S H E A R l , l : l 2 0 0 0 4 0 0 0 T I M E , FIG. 5-2b VERTICAL DISPLACEMENT AT CENTERLINE OF LOADED AREA 6 0 0 0 M I N U T E S 8 0 0 0 1 0 0 0 0 FIG. 5-2 CIRCULAR SURCHARGE LOADING ON A SOIL LAYER. 91 problem are given in Table 5-2. These parameters were adopted from the l i t e r a t u r e by taking the long term response of the models from f i e l d studies reported by Yen (86) and the short term response from Lara-Thomas' laboratory data (40). Once again, the rheological models each have the same i n i t i a l response and the s o i l i s considered to be e l a s t i c in d i l a t a t i o n . In Figure 5-2b, the v e r t i c a l displacement at the centerline of the load-ed area is shown. For the Burgers s o i l and the Maxwell s o i l , the v e r t i c a l displacement increases continually with time. However, a decrease in rate was noted at higher times due to the stresses in the s o i l below the loaded area slowly moving closer to a hydrostatic condition. On the other hand, for the Linear Soli d s o i l , the v e r t i c a l displacement and stresses stopped chang-ing a f t e r about 200 minutes. The surface displacements with time across a v e r t i c a l section of this axisymmetric problem are shown in Figure 5-3. Since the Burgers Body con-tains the additional coupled spring and dashpot of the Linear S o l i d , the" displacements for the Burgers s o i l are somewhat greater than the displacements for the Maxwell s o i l . This difference i s equal to the v i s c o e l a s t i c s o l i d re-sponse of the Linear S o l i d s o i l . The displacements below the c i r c u l a r loaded area and the surface heave around the loaded area show up c l e a r l y in Figure 5-3. As shown in Figure 5-2b and 5-3, the Linear S o l i d s o i l shows very l i t t l e time dependent deformation from the i n i t i a l e l a s t i c configuration. 5-B-3 Creep of a Long Slope ; The f i n i t e element i d e a l i z a t i o n of a long s o i l slope with the low6r end fixed and the basal plane fixed to bed rock is shown in Figure 5-4. Ter-Stepanian has given measured values for the horizontal creep displacements of such a long slope determined by digging a test p i t at Sochi on the Caucasus TABLE 5-2 PARAMETERS FOR RHEOLOGICAL MODELS USED IN FIGURES 5-2, 5-3, and 5-10. BODY IN SHEAR 1.44 x 10 51b/ft 2 9.16 x 10 8lb-min/ft 2 1.90 x 1 0 6 l b / f t 2 5.35 x 10 8lb-min/ft 2 STANDARD LINEAR SOLID IN SHEAR Gx 1.44 x 10 51b/ft 2 G 2 1.90 x 1 0 6 l b / f t 2 n 2 5.35 x 10 8lb-min/ft 2 MAXWELL BODY IN SHEAR G2 1.44 x 1 0 5 l b / f t 2 m 9.16 x 10 8lb-min/ft 2 ELASTIC IN DILATATION v 0.43 K 9.81 x 10 51b/ft 2 E 4.115x 1 0 5 l b / f t 2 BURGERS Gi ni G 2 2000 PSF — X X — o o X-o • + 0.1 LEGEND ELASTIC SOLUTION BURGERS BODY IN SHEAR STANDARD LINEAR SOLID IN SHEAR MAXWELL BODY IN SHEAR (3250 MIN.) 0 10 20 30 40 50 60 70 DISTANCE FROM CENTERLINE OF LOADED AREA , FEET FIG. 5-3 SURFACE DISPLACEMENTS WITH TIME FOR A CIRCULAR SURCHARGE LOADING ON A SOIL LAYER. FIG. 5-4 FINITE ELEMENT IDEALIZATION OF A LONG SLOPE. to 4 ^ 95 Coast (77). It i s not known how far up the slope that the creep well was placed or the end conditions for the p a r t i c u l a r slope. However, the lower end probably ended in a valle y f l o o r or gentle slope, and i t i s f e l t that the representation shown i s reasonable. The cohesive s o i l was founded on bed rock and the water table was located at this basal l e v e l . A Maxwell Body was used to represent the s o i l ' s time dependent shear behaviour and the s o i l was assumed to be e l a s t i c in d i l a t a t i o n . The v i s c o s i t y of the dashpot in the model was taken from Ter-Stepanian's analysis of his f i e l d measure-ments made over a period of about 1.55 years. Young's modulus for the; s o i l was taken as 300 times the undrained shear strength for this s o i l which would seem reasonable from the work by Bjerrum (39). The various parameters for the s o i l are given in Table 5-3. It was found in the plane s t r a i n f i n i t e element analysis of this prob-lem that the stresses and displacement p r o f i l e s were quite uniform for sec-tions through the middle portion of the slope from about node 61 to 89 in Figure 5-4. The f i n i t e element solution for the horizontal creep p r o f i l e after 1.55 years was found for several d i f f e r e n t cases. These r e s u l t s , along with Ter-Stepanian's measured horizontal creep p r o f i l e are shown in Figure 5-5a. From Figure 5-5a, i t can be seen that the f i n i t e element solution for a slope angle a of 23°, which is the slope angle for Ter-Stepanian's f i e l d case, and a Poisson's r a t i o of 0.49, gives good agreement with the measured creep p r o f i l e after 1.55 years. The e f f e c t of changing the slope angle and Poisson's r a t i o shows up c l e a r l y in the creep displacement p r o f i l e s and velo-c i t y p r o f i l e s of Figure 5-5. The creep displacements are quite sensitive to changes in these parameters. A further example was considered in which the lowest layer of elements was made ha l f as s t i f f as the upper layers of ele-TABLE 5-3 PARAMETERS FOR RHEOLOGICAL MODELS USED IN FIGURES 5-5, 5-6 and 5-7. MAXWELL BODY IN SHEAR Gi 8.25 x l O ^ l b / f t 2 m 2.24 x 10 9lb-min/ft 2 STANDARD LINEAR SOLID IN SHEAR Gi 8.25 x lO'Mb/ft 2 G 2 8.25 x .lOMb/ft 2 n 2 2.24 x 10 9lb-min/ft 2 ELASTIC IN DILATATION v 0.49 K 4.1 x 1 0 6 l b / f t 2 E 2.46 x 1 0 5 l b / f t 2 Unit Weight 127 l b / f t 3 97 LEGEND T E R - S T E P A N I A N ' S M E A S U R E D V A L U E S ( 7 7 ) . a = 2 3 ° • F I N I T E E L E M E N T A N A L Y S I S , v = 0 . 4 9 , a = 2 3 ° O o O F I N I T E E L E M E N T A N A L Y S I S , v = 0. 4 9 , a = 2 6 . 5 ° X — X — X F I N I T E E L E M E N T A N A L Y S I S , v = 0 . 4 0 , a = 2 3 ° • • - - • F I N I T E E L E M E N T A N A L Y S I S , S O F T L O W E R L A Y E R EL0WER LAYER = 0 5 EUPPER LAYERS ' V = 0 4 9 i a = 2 3 M A X W E L L B O D Y IN S H E A R , E L A S T I C IN D I L A T A T I O N . M A T E R I A L P R O P E R T I E S G I V E N IN T A B L E 5 - 3 H = 11.5 F E E T C R E E P D I S P L A C E M E N T , F E E T F I G . 5 - 5 a H O R I Z O N T A L C R E E P A F T E R 1.55 Y E A R S ( 4 . 9 x l 0 7 S E C O N D S ) . 11.5 J— UJ UJ u . „ 8 . 6 3 Ul tn < CO 5.75 O CC u . 1- 2 . 8 7 X e> Ul X 0 FIG. 5 0.5 1.0 1.5 C R E E P V E L O C I T Y , F E E T / M I N U T E x l O " 6 F I G . 5 - 5 b ' S T E A D Y - S T A T E * H O R I Z O N T A L C R E E P V E L O C I T Y . 5 H O R I Z O N T A L CREEP D I S P L A C E M E N T AND V E L O C I T Y FOR T H E LONG S L O P E IN F I G U R E 5-4. 98 merits. The e f f e c t of this change shows up c l e a r l y in the increased creep displacement and v e l o c i t y in the basal zone of the slope. In each case, over half the creep displacement occurred in the lower t h i r d of the s o i l and there was l i t t l e creep in the top quarter of the s o i l except for measured seasonal surface creep in the f i e l d creep p r o f i l e . This is the general type of creep behaviour that has been observed i n other f i e l d measurements of long slopes (86). For creep problems involving long slopes, there is a basal zone of active material that contributes the majority of the creep deformation. The material above this zone moves e s s e n t i a l l y as a r i g i d body. This shows up p a r t i c u l a r l y well for the slope with the s o f t e r basal layer. Thus, measurements of r e l a t i v e surface creep may prove quite misleading and the actual creep displacement p r o f i l e should be measured from bed rock i f creep rates are required f o r a long slope. 5-B-4 Creep Displacements at a Vertical Face The free end displacements associated with a 20 foot deep excavation in a long 2 to 1 s o i l slope are shown in Figure 5-6 along with the dimensions and boundary conditions of the problem. The f i n i t e element i d e a l i z a t i o n of the problem i s s i m i l a r to that shown i n Figure 5-4 for the long slope. The slope is considered to be creeping under the influence of gravity a f t e r the toe i s excavated. Such a removal of material by r i v e r erosion or man-made excavations often i n i t i a t e s the creep of slopes that have been stable for some time. Two shear behaviour models are considered for the s o i l , the Maxwell Body and the Standard Linear S o l i d , with parameters given in Table 5-3. The s o i l i s assumed to be e l a s t i c in d i l a t a t i o n . For the Maxwell s o i l , the free dashpot in the rheological model allows unlimited viscous flow. On the other hand, a s o i l represented by the Linear S o l i d , or 3-parameter visco-99 LEGEND O O o MAXWELL BODY IN SHEAR. X - - X - - X STANDARD LINEAR SOLID IN SHEAR. (ELAPSED TIME AFTER 'LOADING' IS GIVEN IN MINUTES ON EACH NODE I 0 0.6 1.2 HORIZONTAL DISPLACEMENT FEET FIG. 5-6 DISPLACEMENTS WITH TIME AT THE END OF A LONG SLOPE WITH A GRAVITY LOADING , CLAY. 100 e l a s t i c s o l i d does not exhibit unlimited, non-recoverable viscous flow since there are no free dashpots in this rheological model. This difference be-tween v i s c o e l a s t i c s o l i d behaviour and v i s c o e l a s t i c f l u i d behaviour shows up c l e a r l y for the creep displacements shown in Figure 5-6. Approximately 200,000 minutes af t e r the gravity loading i s applied, the face of the v e r t i c a l excavation for the Linear S o l i d s o i l has stopped creep-ing, while the Maxwell s o i l slope continues to displace in a viscous manner for times greater than 200,000 minutes. Most of the creep deformation occurs in the lower t h i r d of the excavation which i s what would be anticipated from the previous problem. After 200,000 minutes (approximately 139 days), the total displacements for the Linear S o l i d s o i l are about twice the i n i t i a l e l a s t i c displacements and no further movement takes place. After 200,000 minutes, the total displacements for the Maxwell s o i l are about eight times the e l a s t i c displacements and increasing at a r e l a t i v e l y constant rate. This amounts to a total movement of about one foot outward and one-half foot down-ward at the top of the v e r t i c a l face. Such large creep movements would cer-t a i n l y indicate potential f a i l u r e of the excavated slope and require that remedial action be taken for a s o i l that has such v i s c o e l a s t i c f l u i d behaviour. In contrast, the creep movements for the s o i l considered to be a v i s c o e l a s t i c s o l i d are quite small and are probably acceptable. 5-B-5 Heavy Retaining Wall at the End of a Creeping Slope If a r i g i d , or very heavy retaining wall is b u i l t to l i m i t the creep deformations of the excavated face i n the previous problem, the total h ori-zontal force exerted by the s o i l on the wall builds up as shown in Figure 5-7. As expected, the increase in horizontal force from the i n i t i a l value is much lower for the Linear Soli d s o i l than for the Maxwell s o i l , and the force 101 (O Q Z O OL X 1-Q O O O u_ cc UJ CL UJ o cc o 100 000 90 000 80 000 70000 60 000 LEGEND O •o« •X-O MAXWELL BODY IN SHEAR X STANDARD LINEAR SOLID IN SHEAR. SOIL IS ELASTIC IN DILATATION. MATERIAL PROPERTIES GIVEN IN TABLE 5-3 < 2 50000& 40 000 J 3 -O O •X-O -x 0 40 000 80 000 120 000 TIME , MINUTES 160 000 200 000 FIG. 5-7 HORIZONTAL FORCE WITH TIME ON A RIGID RETAINING WALL FROM NODE I TO NODE 5 OF FIG. 5-6. 1Q2 stops increasing at an elapsed time of about 200,000 minutes af t e r a f a i r l y rapid r i s e in the f i r s t 20,000 minutes. The increase in pressure i s less than 15% for this v i s c o e l a s t i c s o l i d behaviour. On the other hand, the force on the retaining wall for the Maxwell s o i l increases continually with time, but the rate of increase is slowing down since the wall acts as a constraint. For this v i s c o e l a s t i c f l u i d behaviour, the increase i n pressure i s about 60% in 200,000 minutes. This v i s c o e l a s t i c f l u i d behaviour i s typical of the type of behaviour t that would be expected as a mass of s o i l behind a heavy retaining wall re-laxes and the pressure in the s o i l builds up towards the hydrostatic state. Of course, as the horizontal force builds up, the wall w i l l tend to y i e l d and the pressure build up w i l l not be as great as f o r the r i g i d wall consi-dered. It i s possible in the programs to allow the wall to creep a l s o , and for this case, as for a y i e l d i n g w a l l , the. pressure build up would be much less for both models considered. By making a judgment on the type of s o i l behaviour anticipated in the f i e l d where time dependent movements might cause problems, i . e . v i s c o e l a s t i c s o l i d or v i s c o e l a s t i c f l u i d behaviour, estimates of the potential forces on the wall due to s o i l creep can be made. From these q u a l i t a t i v e r e s u l t s , some provision can be made in the design to handle any anticipated increase in horizontal forces i f they are s i g n i f i c a n t . 5-B-6 Creep Displacements at a Vertical Face in Rock To i l l u s t r a t e the creep of another material, the free end displacements for the same ve r t i c a l excavation in a long 2 to 1 slope are shown in Figure 5-8 for shale. The shale slope was considered for the two cases of Burgers Body, or four-parameter f l u i d , shear behaviour and Standard Linear S o l i d shear behaviour. The parameters for these models were taken from Hardy's ; 103 LEGEND • • • B U R G E R S BODY IN S H E A R X - - X - - X S T A N D A R D L I N E A R S O L I D IN S H E A R . ( E L A P S E D T I M E A F T E R ' L O A D I N G ' IS G I V E N IN M I N U T E S O N E A C H 0 10 2 0 H O R I Z O N T A L D I S P L A C E M E N T x I O " 5 F E E T FIG. 5-8 DISPLACEMENTS WITH TIME AT THE END OF A LONG SLOPE WITH A GRAVITY LOADING, S H A L E . 104 survey of the creep of geologic materials (32) and are given in Table 5-4. Once again, the zone of active basal material shows up c l e a r l y and the creep behaviour for the shale slope is very s i m i l a r to that of the clay slopes discussed previously. Of course, the amount of creep movement is much re-duced for the s t i f f e r shale. Rock i s often considered to be a l i n e a r visco-e l a s t i c material for design purposes where time dependent effects are impor-tant (36), and the various methods of analysis developed should be useful in many rock mechanics problems. 5-B-7 Creep of a 45° Slope With a Distributed Load The f i n i t e element i d e a l i z a t i o n of a 45° slope i s shown in Figure 5-9. The creep displacements with time for this slope with a uniform distributed load of 2,000 psf on the crest from node 107 to 145 are shown in Figure 5-10. Burgers Body shear behaviour and Standard Linear S o l i d shear behaviour were considered with the parameters for each rheological model given in Table 5-2. Gravity stresses were not considered in this problem. The Burgers s o i l shows viscous flow while the Linear S o l i d s o i l has stopped creeping after about 2,000 minutes. After 2,000 minutes of creep, the displacements under the load for the Linear S o l i d s o i l are s l i g h t l y greater than the e l a s t i c values, and for most of the slope, the difference between the e l a s t i c response and the creep response of this v i s c o e l a s t i c s o l i d i s n e g l i g i b l e . Even at 2,000 minutes, there i s already a large difference between the e l a s t i c response and the creep response for the Burgers s o i l . The creep d i s -placement pattern shows up c l e a r l y for the exaggerated displacement scale used in Figure 5-10. After 12,000 minutes have elapsed, there i s a very n o t i -ceable bulging of the slope face and some heave at the base of the slope. The total creep displacements under the distributed load are now about four times TABLE 5-4 PARAMETERS FOR RHEOLOGICAL MODELS USED IN FIGURE 5-8 (SHALE) • BURGERS BODY IN SHEAR 2.08 x 1 0 8 l b / f t 2 3.50 x 1 0 n l b - i T i i n / f t 2 7.45 x 1 0 8 l b / f t 2 2.08 x 10 9lb-min/ft 2 ni G 2 r i 2 STANDARD LINEAR SOLID IN SHEAR G : 2.08 x 1 0 8 l b / f t 2 G 2 7.45 x 1 0 8 l b / f t 2 Ti2 2.08 x 10 9lb-min/ft 2 ELASTIC IN DILATATION v 0.3 K 4.51 x 1 0 8 l b / f t 2 E 5.40 x 1 0 8 l b / f t 2 Unit Weight 150 l b / f t 3 185 NODES 316 ELEMENTS SCALE I" =20' FIG. 5-9 FINITE ELEMENT IDEALIZATION OF A 45° SLOPE. o F R E E V E R T I C A L L Y F I X E D H O R I Z O N T A L L Y , 108 the e l a s t i c values. For this s o i l , with v i s c o e l a s t i c f l u i d behaviour, the displacements w i l l become very large with elapsed time. 5-C GENERAL CREEP PROBLEMS Before examining the creep of s o i l s described by Singh and Mitchell's empirical relationship, an example of the creep of g l a c i e r ice w i l l be gi ven. 5-C-l Creep of Glacier Ice The horizontal creep displacements and creep v e l o c i t y for a 'glacier' represented by a long 1 to 2 slope 20 feet high are given in Figure 5-11. The g l a c i e r i s assumed to have the same boundary conditions and a si m i l a r f i n i t e element representation as the long slope in Figure 5-4. A uniaxial creep relationship for glac i e r ice that i s nonlinear in stress was taken from Col beck and Evans' f i e l d measurements of the creep of g l a c i e r ice (15) and this relationship, and the properties of the i c e , are shown on Figure 5-11. Paterson (58), in his book on the physics of g l a c i e r s , has indicated that such a creep relationship i s appropriate to describe the creep of many glaciers and that f o r the multiaxial state of s t r e s s , the equivalent stress and equivalent creep s t r a i n increment can be used in such a uniaxial creep relationship. It should be noted that the creep relationship i s l i n e a r in time and steady-state creep movements were found a f t e r about two months of creep under the gravity stresses. The steady-state creep veloc i t y and creep displacements for one month ^ are shown in Figure 5-11. These values were taken in the middle part of the gl a c i e r where the stresses and creep p r o f i l e s are uniform for any section through the glacier. The lower zone of movement shows up very c l e a r l y with 109 LEGEND U N I A X I A L C R E E P R E L A T I O N S H I P IS £ c = 0 . 4 4 c r ' - 9 ( 6 ° IN Y E A R S - 1 , a IN B A R S ) E = I . 3 4 x l 0 8 P S F v = 0 . 4 0 y = 5 5 L B / F T 3 • C R E E P D I S P L A C E M E N T X X X C R E E P V E L O C I T Y C R E E P V E L O C I T Y , F E E T / M I N U T E x I O " 7 2 0 LU LU li_ LU CO < m o cc X 15 10 uJ o x 10 15 } v \ 1 1 ll V m C R E E P D I S P L A C E M E N T IN I M O N T H , F E E T x I 0 " 2 FIG. 5-11 HORIZONTAL CREEP DISPLACEMENT AND VELOCITY FOR A 'GLACIER' REPRESENTED BY THE LONG SLOPE IN FIGURE 5-4. n o over half the movement taking place tn the lower quarter of the glacier and the upper half moving essential ly as a r ig id body. This behaviour is very similar to that observed in the clay and rock slopes considered to be l inear v iscoelast ic . Field observations reported by Paterson (58) indicate that this is the type of creep prof i le observed in glaciers where s l id ing on a basal plane is not too great. 5-C-2 Axisymmetric Problem An i n i t i a l problem was studied to check the use of the incremental solution method for a cohesive soi l described by Singh and Mitchel l 's empi-r ical creep relationship. The axisymmetric problem shown in Figure 5-2a was considered for the case where lateral deformations are permitted and the fu l l width of the surface is uniformly loaded. This f in i te element idealization then corresponds to a t r iaxia l specimen with fr ict ionless end plates. The state of stress and strain is uniform throughout the soi l with the vertical stress, a . equal to the applied pressure, p, and the other stresses, a , a., z r 6 and x equal to zero, rz n The cohesive soi l considered is London Clay with the parameters for Singh and Mitchel l 's empirical creep relationship taken from Bishop and Loven-bury's long-term creep tests (3). For this data, the uniaxial creep relation-ship is given by: : C = 2.66 x 10" 2 e ( 2 ' 8 4 x 10 D) _„ - 0.97 1 (5:1) C where £ is in %/day, D in psf and t in days. Bishop and Lovenbury conducted drained creep tests so that this data would not normally be used in the analy-sis since the cumulative creep law was developed from undrained test data. However, in this particular problem the stresses are not changing with time m and the cumulative creep law i s not required in the determination of the creep s t r a i n s . The e l a s t i c deformations and stresses were taken as the i n i t i a l conditions for the creep analysis and were considered to be applica-ble at time t=l day in the creep relationship. A pressure of 880 psf was applied since th i s corresponds to the stress level (19% of the reference strength) for one of the creep curves given by Bishop and Lovenbury (3). For t h i s problem, the equivalent stress aQ i s equal to the applied pressure of 880 psf. The stresses did not vary with time, and. af t e r 94 days of creep, the change in axial s t r a i n due to creep was • 0.135%. From the laboratory creep curve for this stress l e v e l , the change in axial s t r a i n measured for this time period i s about 0.14%. This close corres-pondence between the computed and experimental values of the change in s t r a i n provides a further check on the methods developed. 5-C-3 Excavation of a 45° Slope The s o i l considered in the following examples of the general creep analy-s i s of excavation problems i s undisturbed Redwood City Clay with Singh and Mitchell's uniaxial creep relationship shown on Figure 5-12. Parameters for the creep relationship were taken from the laboratory creep curves for un-drained tests on this p a r t i c u l a r s o i l given by Singh and Mitchell (72). The e l a s t i c deformations and stresses af t e r excavation were taken as the i n i t i a l conditions for the creep analysis of the two excavation problems examined. These i n i t i a l conditions after excavation were assumed to be applicable at time t=l minute in the creep relationship since information on the i n i t i a l condition of s t r a i n in the laboratory tests was not reported. In both prob-lems examined, creep during the period of excavation was not considered. A discussion of creep rupture for these problems i s given in the next section. WATER TABLE AT THE SURFACE » > a ZONE OF EXCAVATED DISPLACEMENT S C A L E LEGEND UNDISTURBED REDWOOD CITY CLAY TIME A F T E R EXCAVATION GIVEN BY • 0 MINUTES X X X 3 0 0 0 MINUTES o 0 o 2 0 0 0 0 0 MINUTES UNIAXIAL CREEP RELATIONSHIP IS: (I) 6 C = 5 . 8 3 x l 0 - 6 e ( 7 - , 4 x , 0 " 3 D ) ^ a 7 6 5 ( £ IN %/MINUTE, D IN PSF, t IN MINUTES) PARAMETERS FROM SINGH AND MITCHELL SOIL PROPERTIES : E = 288 000 PSF v - 0.45 /s = 125 L B / F T 3 K Q= 1.0 2 2 0 NODES 388 E L E M E N T S BOUNDARY CONDITIONS T H E SAME AS FOR FIGURE 5-9 FIG. 5-12 EDGE DISPLACEMENTS WITH TIME FOR A 45° EXCAVATED SLOPE. CREEP GIVEN BY SINGH AND MITCHELL'S EMPIRICAL RELATIONSHIP. 113 The creep displacements due to the excavation of a 45° slope are shown in Figure 5-12. The f i n i t e element i d e a l i z a t i o n of the part of the 45° slope not excavated is given in Figure 5-9. For this problem, the water table was considered to be at the ground surface and the c o e f f i c i e n t of earth pressure at rest, K q, was taken as 1.0. Since the s o i l was considered to be e l a s t i c during excavation, the excavation was accomplished in a single step. The boundary conditions for the excavated slope are the same as those shown for the 45° slope in Figure 5-9 and the various material properties used are indicated on Figure 5-12. The displacements in the f i r s t 3,000 minutes (2 days.) of creep were nearly as great as those in the following 200,000 minutes (139 days) of creep. After 200,000 minutes the increase in creep displacements was ne g l i g i b l e . This decrease in the rate of creep s t r a i n accumulation with time i s expected from the form of Singh and Mitchell's creep relationship. The exaggerated displacement scale shows the bulging at the base of ths slope very c l e a r l y . At the lower portion of the slope, the stresses are quite high, and since the creep strains vary nonlinearly with st r e s s , the creep displacements are ex-pected to be greatest in this zone of the slope. The creep displacements are more than twice the e l a s t i c displacements in this zone of the slope. Usually, during excavation, the slope would be cut to the f i n a l desired p r o f i l e and the e l a s t i c displacements are not of concern. However, the creep movements would result in up to another 0.25 feet of displacement on the lower slope face and toe of the slope and these movements could cause damage to structures placed close to the otherwise stable slope. 5-C-4 Excavation of a Slurry Trench The f i n i t e element i d e a l i z a t i o n for the excavation of a s l u r r y trench i s 114 shown in Figure 5-13 where the material above A-B-C is to be removed during construction. Creep displacements associated with the excavation of this sl u r r y trench in a s o i l with the same nonlinear creep relationship as the previous problem are shown in Figure 5-14. The dimensions of the trench, the s o i l properties and the boundary conditions for the excavation are also indicated on Figure 5-14. The water table was taken at ground level and K Q was taken as 0.5 for this excavation problem. Excavation was again accom-plished in a single step.. Two cases were considered: a) no s l u r r y in the trench; and, b) s l u r r y in the trench. a) No Slurry in the Trench Without the s l u r r y in the excavated trench, the e l a s t i c deformations are f a i r l y large. However, except for considerable bulging just above the bottom of the trench, the general increase in displacements due to creep i s quite small. The zone of bulging would be the zone of f a i r l y high stresses since the bottom of the trench i s not o f f e r i n g much support at this point. Considerable heaving occurred at the bottom of the trench from A to B, but for c l a r i t y , this is not indicated on Figure 5-14. After 18,000 minutes (12.5 days) the creep movements were ne g l i g i b l e with increasing time. The formation of a zone of bulging indicates the importance of getting the sl u r r y into place during excavation. This is c l e a r l y evident i n the next example when the trench is kept f i l l e d with a s l u r r y . b) Slurry in the Trench The s l u r r y was considered to have a unit weight of 70 l b / f t 3 , and an additional head of 2 feet by diking above the water table at the ground sur-face as i s often done in practice. The unit weight of the s l u r r y was taken from values reported for a number of s l u r r y trenches. With the s l u r r y in 5" 175 NODES 305 ELEMENTS FIG. 5-13 FINITE ELEMENT IDEALIZATION FOR EXCAVATION OF A SLURRY TRENCH. LEGEND U N D I S T U R B E D R E D W O O D C I T Y C L A Y C R E E P R E L A T I O N S H I P I S G I V E N O N F I G . 5 - 1 2 x - 0 — x -• _ 0 — • - -1 8 0 0 0 • • ->— * — ( 5 0 ' N O S L U R R Y I N T R E N C H x — x — x S L U R R Y I N T R E N C H ( E L A P S E D T I M E A F T E R E X C A V A T I O N I S G I V E N I N M I N U T E S O N E A C H C U R V E ) S O I L : E = 2 8 8 0 0 0 P S F , v = 0 . 4 5 / S = 1 2 5 L B / F T ? K 0 = 0 . 5 * S L U R R Y - y = 7 0 L B / F T 3 . Z O N E O F W A T E R T A B L E A T S U R F A C E 1 7 5 N O D E S 3 0 5 E L E M E N T S B I E X C A V A T E D M A T E R I A L 1 K - 1 2 0 * 1 0 ' B O U N D A R Y C O N D I T I O N S T H E S A M E A S F O R T H E 4 5 ° S L O P E . * 2 F T . S T A T I C H E A D O F S L U R R Y I N A D D I T I O N T O S L U R R Y I N T H E T R E N C H . B X x I 0 I x • 4-* 1 \ \ » • 1 • I I I I • • i • 0 i 1 1 8 0 0 0 • i • i • \ y <- 5 ' M 1 0 0 . 1 D I S P L A C E M E N T S C A L E F E E T > 0 . 1 FIG. 5-14 EDGE DISPLACEMENTS WITH TIME AFTER EXCAVATING A SLURRY TRENCH. CREEP GIVEN BY SINGH AND MITCHELL'S EMPIRICAL RELATIONSHIP. zl c n 117 place, the e l a s t i c deformations were reduced to about 25% of those without the s l u r r y in place. In addition, the creep displacements were n e g l i g i b l e at a l l times and there was no zone of bulging. This behaviour again shows how sensitive the creep relationship i s to stress l e v e l . The importance of keeping the s l u r r y in the trench during excavation i s shown c l e a r l y for this typical excavation problem. 5-D CREEP RUPTURE A total rupture l i f e relationship was used in conjunction with Singh and Mitchell's ( ^ t ) f a i ] u r e condition to study the potential creep rupture of the problems described in the previous section. This relationship states that creep rupture w i l l occur i f the value of the current s t r a i n rate multi-pli e d by the total elapsed time since i n i t i a t i o n of creep reaches a c r i t i c a l f a i l u r e value. For the multiaxial stress state, the principal shear s t r a i n rate was used in the uniaxial creep rupture relationship. Since the time from the i n i t i a t i o n of creep is known for the various problems examined, i t is possible to use such a total rupture l i f e relationship. However, for natural slopes, the time from the beginning of creep is not known and a total rupture l i f e relationship cannot be used to predict the creep rupture l i f e in such cases. For natural slopes, Snead's method (73) of predicting the time to rupture during the t e r t i a r y creep stage appears most promising. This method has not been considered in the present study. For the problems examined in this section, the cohesive s o i l was again considered to be undisturbed Redwood City Clay with the uniaxial creep, r e l a -tionship given on Figure 5-14. As Singh and Mitchell have indicated (72), a cohesive s o i l with m<l w i l l f a i l in creep rupture i f the value of (et) 118 reaches a f a i l u r e value, ( e t ) f a ^ u r g . For Redwood City Clay with m = 0.765, ( e t ) f a i l u r e i s given by Singh and Mitchell as (2.8 + 0.5)%. The lower value of 2.3% was adopted to examine the potential creep rupture of the various problems. 5-D-l Axisymmetric Problem The same i d e a l i z a t i o n of a t r i a x i a l specimen discussed in Section 5-C-2 was used to check the creep rupture relationship. The applied pressure, p, was increased until creep.rupture of the s o i l occurred. Since the state of stress and s t r a i n i s uniform throughout the s o i l , rupture should occur in every element at the same elapsed time i f (it) reaches the rupture value of 2.3%. The applied pressure (which is equal to the equivalent s t r e s s , a „ for this problem) and creep rupture behaviour are shown in Table 5-5. No rupture occurred until the applied pressure was 90% of the reference strength, and then as the pressure was increased, the time to rupture decreased rapidly from 207 minutes to less than 0.01 minutes. Failure occurred in a l l of the elements at the same time as anticipated. For a laboratory test at 90% of the reference strength, creep rupture occurred about 1,000 minutes a f t e r loading and the increase in str a i n during creep was about 9% (72). A labora-tory specimen at 85% of the reference strength did not show creep rupture. This comparison indicates that the solution methods developed are adequate to predict the potential creep rupture behaviour for problems in which the time from loading i s known. It i s important to note that as the applied pressure approaches that required to cause rupture, the permissible time interval f o r each solution step becomes very small since the s t r a i n rate i s large. This decrease in the time intervals can also be used as an indication of potential creep rupture. TABLE 5-5 CREEP RUPTURE FOR AXISYMMETRIC PROBLEM Applied Equivalent Percent of Pressure Stress, a Reference Creep psf psf e Strength* Behaviour 1300 1300 70 No Rupture 1470 1470 80 No Rupture 1660 1660 90 Rupture in 207 minutes (7% creep strain) 1750 1750 95 Rupture in 17 minutes (5% creep strain) 2210 2210 120 Immedi ate Rupture * Reference Strength for Undisturbed Redwood City Clay i s 1840 psf (72). (Shear strength measured in a standard laboratory test) 120 The s t i f f n e s s of each element was reduced a f t e r rupture to one-half i t s i n i t i a l value and creep was allowed to continue even though the s o i l mass had ruptured. Since the allowable creep s t r a i n increment (and hence time interva l ) i s computed automatically, this reduction in s t i f f n e s s did not cause any problems with solution s t a b i l i t y . 5-D-2 Excavation Problems The excavations for the 45° slope and the s l u r r y trench, without and with the s l u r r y in place, were checked for creep rupture during the creep analysis described in Sections 5-C-3 and 5-C-4. It was found that creep rupture did not occur for any of these excavation problems. To examine the solution method further, the value of the c o e f f i c i e n t of earth pressure at rest, K Q, was increased for the excavation of the s l u r r y trench without the sl u r r y in place. For K q = 1.0 and K q = 1.5, there was no creep rupture of the trench although the time intervals became f a i r l y small for K Q = 1.5. When K Q was taken as two, rupture occurred immediately in the zone of high stresses j u s t above B in Figure 5-14. Since the time intervals were extreme-ly small for this case, the propagation of the rupture zone was not traced. The various problems discussed in this chapter i l l u s t r a t e the use of the methods developed to consider problems in which time dependent s o i l be-haviour i s considered. Examples of the creep of rock and ice give some i n d i -cation of the extension of the methods of analysis to related geotechnical f i e l d s . While certain idealizations have been made throughout the problems, the solution methods developed are f l e x i b l e enough to handle more complicated problems as laboratory data and f i e l d measurements become available. 121 CHAPTER 6 CONCLUSIONS Application of the f i n i t e element method of stress analysis to non-time-dependent problems in s o i l mechanics has enabled the engineer to gain much information on the deformations and stresses in earth structures. However, a survey of laboratory studies and f i e l d observations indicates that time is a very important factor in. the behaviour of cohesive s o i l s and should be i n -cluded in the analysis of many problems. For this reason, the f i n i t e element method was extended to deal with problems in which the time dependence of the properties of s o i l is considered. The incremental i n i t i a l s t r a i n solution method developed for s o i l s described by a general creep relationship, and the i t e r a t i v e solution method developed for s o i l s considered to be l i n e a r visco-e l a s t i c , were checked by comparing computed results with analytical r e s u l t s , where available. In each case, excellent agreement between the results was found. At present, Singh and Mitchell's three-parameter empirical creep r e l a -tionship appears to o f f e r the most general description of the creep of cohe-sive s o i l s over a wide range of sustained deviator stress levels and test conditions. It i s very important that the parameters for such a general creep relationship be determined in laboratory tests that closely simulate the actu-al in s i t u s o i l conditions and loading. In p a r t i c u l a r , more plane st r a i n creep tests are required in which a range of stress levels are considered and various stress paths leading to rupture are studied. The extension of uniaxial constant stress creep relationships to the multiaxial changing stress state that i s representative of the f i e l d condi-122 tions was achieved by making three assumptions: 1) an equivalent stress and equivalent creep s t r a i n increment can be used in the uniaxial r e l a t i o n s h i p ; 2) Snead's hypothesis of a unique relationship between current deviator stress, current s t r a i n , and current s t r a i n rate holds; and, 3) the components of creep s t r a i n can be recovered from an equation analogous to the Prandtl-Reuss equa-tion of p l a s t i c i t y . Further laboratory work i s required to check the general a p p l i c a b i l i t y of these assumptions for various cohesive s o i l s . However, the incremental i n i t i a l s t r a i n solution method developed is extremely f l e x i b l e and can e a s i l y be modified to handle other and more complicated relationships i f they are developed. The excavation of a s l u r r y trench showed how sensitive the general creep relationship i s to stress l e v e l . With the s l u r r y in place, the creep deforma-tions were ne g l i g i b l e . However, without a s l u r r y in the trench, a zone of bulging formed with elapsed time in the area of high stresses near the bottom of the trench. This example i s typical of the type of problems considered with the incremental i n i t i a l s t r a i n solution method. For q u a l i t a t i v e studies, i t i s often convenient to consider a cohesive s o i l to be l i n e a r v i s c o e l a s t i c and to examine the s o i l ' s behaviour for a range of rheological models. Comparisons between v i s c o e l a s t i c s o l i d behaviour and v i s c o e l a s t i c f l u i d behaviour gave some indication of design problems that may occur in the f i e l d . This was shown for the case of a r i g i d retaining wall at the end of a slope. With elapsed time, there was a large increase in the horizontal force on the wall for a s o i l with v i s c o e l a s t i c f l u i d behaviour. The increase in force was n e g l i g i b l e for v i s c o e l a s t i c s o l i d behaviour. Also, f i e l d measurements can often be used to develop parameters for the rheological models. For instance, close agreement was found between the computed creep 123 displacement p r o f i l e and measured creep displacement p r o f i l e for a long slope. Creep rupture was considered in the various excavation problems examined by using a tot a l rupture l i f e r elationship. It was found that none of the ex-cavations would undergo creep rupture. By increasing the applied forces dur-ing excavation of a slurry trench, creep rupture was induced in the zone of high stresses where bulging due to creep was evident at lower stress l e v e l s . The four main areas in which the present research should be extended appear to be: 1) Consideration of creep during the actual excavation and con-struction stages of earth structures. Since loads can be applied at any time step in the incremental solution method for general creep, i t should be r e l a t i v e l y easy to make the necessary modifications in the programs. 2) Develop rheological models to describe the d i l a t a t i o n a l behaviour of s o i l s for use in the i t e r a t i v e solution method for s o i l s assumed to be l i n e a r v i s c o e l a s t i c . Also, additional shear behaviour models could be developed. 3) Introduce time-dependent volume changes into the incremental solution method for general creep. Since the creep problem is solved in a series of time steps, i t should be possible to con-sider volume changes at each step of the solution. Further laboratory research on creep in terms of e f f e c t i v e stresses, and on the e f f e c t of drainage conditions on the s t r e s s - s t r a i n - s t r a i n rate relationship, w i l l be required for this analysis. 4) Determination of an upper y i e l d stress relationship for a wide range of cohesive s o i l s for i n s i t u conditions. Such a f a i l u r e c r i t e r i a can be introduced into the incremental solution method to check on the potential creep rupture of earth structures. F i n a l l y , there are many problems in the related f i e l d s of rock and ice mechanics in which the various solution methods can be applied. It i s hoped that some of these potential applications w i l l be further developed. 125 BIBLIOGRAPHY 1. A l f r e y , T., "Non-homogeneous Stresses in V i s c o e l a s t i c Media", Quart. Appl. Math., Vol. 2, 1944, pp. 113-119. 2. Biot, M.A., "Dynamics of V i s c o e l a s t i c Anisotropic Media", Proceedings of the Fourth Midwestern Conference on Solid Mechanics, September 1955, pp. 94-108. 3. Bishop, A.W. and Lovenbury, H.T., "Creep Characteristics of Two Undisturbed Clays", Proceedings of the Seventh International Conference on Soil Mechanics and Foundation Engineering, Mexico Cit y , Vol. 1 , 1969, pp. 29-37. 4. Bjerrum, L., "The Third Terzaghi Lecture: Progressive Failure in Slopes of Overconsolidated P l a s t i c Clay and Clay Shales", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 93, Part 1, No. SM5, Proc. Paper 5456, September 1967, pp. 3-49. 5. Boley, B.A. and Weiner, J.H., Theory of Thermal Stresses, John Wiley and Sons, New York, 1960. 6. Brown, C.B. and King, I.P., "Automatic Embankment Analysis: Equilibrium and I n s t a b i l i t y Conditions", Geotechnique , Vol. 16, No. 3, 1966, pp. 209-219. 7. Bucher, Fe l i x , "Consolidation Theories and Fi n i t e Element Solutions", M.Eng. Report, University of B r i t i s h Columbia, December 1970. 8. Cakmak, A.S. and Yau, W.W.F., "Time-Dependent Concentrated Surface Load Moving with Diminishing Velocity on a Visc o e l a s t i c Half Space", Trans-actions of the Society of Rheology , Vol. 10, No. 1, 1966, pp. 85-95. 9. Campanella, R.G., "Effect of Temperature and Stress on the Time-Deformation Behaviour of Saturated Clay", Ph.D. Thesis, University of C a l i f o r n i a , Berkeley, 1965. 10. Chang, T.Y., Ko, H.Y., Scott, R.F. and Westmann, R.A. , "An Integrated Approach to the Stress Analysis of Granular Materials", Report on research conducted f o r the National Science Foundation, C a l i f o r n i a Institute of Technology, Pasadena, C a l i f o r n i a , 1967. 11. Christensen, R.W. and Wu, T.H., "Analysis of Clay Deformation as a Rate Process", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 90, No. SM6, Proc. Paper 4147, November 1964, pp. 125-157. 12. Christian, John T., "Undrained Stress Distribution by Numerical Methods", Journal of the Soil Mechanic's and Foundations Division, ASCE, Vol. 94, No. SM6, Proc. Paper 6243, November 1968, pp. 1333-1345. 13. Christian, John T. and Boehmer, Jan Willem, "Plane Strain Consolidation by F i n i t e Elements", Journal of the Soil Mechanics and Foundations Div-ision, ASCE, Vol. 96, No. SM4, Proc. Paper 7437, July 1970, pp. 1435-1457. 126 14. Clough, R.W. and Woodward, R.J., "Analysis of Embankment Stresses and Deformations", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 93, No. SM4, Proc. Paper 5329, July 1967, pp. 529-549. 15. Colbeck, S.C. and Evans, Roger J., "Experimental Studies Related to the Mechanics of Glacier Flow", The Trend in Engineering, Quarterly Journal of the University of Washington College of Engineering, Vol. 21, No. 2, April 1969, pp. 8-14. 16. Corneliussen, A.H. and Lee, E.H., "Stress Distribution Analysis for Linear V i s c o e l a s t i c Materials", Creep in Structures, Hoff, N.J. ( E d i t o r ) , Academic Press Inc., New York, 1962, pp. 1-20. 17. Drescher, A., "Nonlinear Creep of Cohesive S o i l " , Archiwum Mechaniki Stosowanej, Vol. 5, No. 19, 1967, pp. 745-763. 18. Duncan, J.M. and Dunlop, P., "Slopes in S t i f f - F i s s u r e d Clays and Shales", Journal of Soil Mechanics and Foundations Division, ASCE, Vol. 95, No. SM2, Proc. Paper 6449, March 1969, pp. 467-492. 19. Dunlop, P. and Duncan, J.M., "Development of Failure around Excavated Slopes", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 96, No. SM2, Proc. Paper 7162, March 1970, pp. 471-493. 20. Eyring, H., "Viscosity, P l a s t i c i t y and Diffusion as Examples of Absolute Reaction Rates", Journal of Chemical Physics, Vol. 4, April 1936, pp. 283-291. 21. Flugge, Wilhelm, Viscoelasticity , Bl a i s d e l l Publishing Company, Waltham, Massachusetts, 1967. 22. Freudenthal, A.M. and R o l l , F., "Creep and Creep Recovery of Concrete under High Compressive Stress", Proceedings, American Concrete Institute, Vol. 54, 1958, pp. 1111-1142. 23. Fung, Y.C., Foundations of Solid Mechanics , Prentice-Hall, Inc., Engle-wood C l i f f s , New Jersey, 1965. 24. Garafalo, F., Fundamentals of Creep and Creep Rupture in Metals, Macmillan Series in Materials Science, New York, 1965. 25. Goldstein, M., "Discussion", Proceedings of the Brussels Conference on Earth Pressure Problems, Vol. I l l , 1958, pp. 91-92. 26. Gray, Donald H., "Effects of Forest Clear Cutting on the S t a b i l i t y of Natural Slopes", Progress Report, ORA Project 01939, Department of C i v i l Engineering, University of Minnesota, September 1969. 27. Green, A.E. and R i v l i n , R.S., "The Mechanics of Nonlinear Materials with Memory", Archiv. Rat. Mech. Anal., Vol. 1 , 1957, pp. 1-21. 28. Greenbaum, G.A., "Creep Analysis of Axisymmetric Bodies", Ph.D. Thesis, Engineering Mechanics, University of C a l i f o r n i a , Los Angeles, 1966. 127 29. Greenbaum, G.A. and Rubinstein, M.F., "Creep Analysis of Axisymmetric Bodies using F i n i t e Elements", Nuclear Engineering and Design, Vol. 7, 1968, pp. 379-397. 30. H a e f e l i , R., "Creep Problems in S o i l s , Snow and Ice", Proceedings of the Third International Conference on Soil Mechanics and Foundation Engin-eering, Vol. 3, 1953, pp. 238-251. 31. Hanson, Roger S., "Practical Solution of P l a s t i c i t y Problems by the Incremental Theory of P l a s t i c i t y " , Ph.D. Thesis, Iowa State College, 1958. 32. Hardy, H.R., "Inelastic Behaviour of Geologic Materials, Part I", Canada Mines Branch Divisional Report FMP 66/51 -P, (Limited D i s t r i b u t i o n ) , 1966. 33. Henkel, D.J., "Investigation of Two Long-Term Failures in London Clay Slopes at Wood Green'and Northolt", Proceedings of the Fourth International Conference on Soil Mechanics and Foundation Engineering, Vol. 2, 1957, pp. 315-320. 34. H i r s t , T.J. and M i t c h e l l , J.K., "Compositional and Environmental Influ-ences on the Stress-Strain-Time Behaviour of S o i l s " , Report No. TE-68-4, Department of C i v i l Engineering, Institute of Transportation and T r a f f i c Engineering, University of C a l i f o r n i a , Berkeley, April 1968. 35. Ishihara, K., "Relations between Process of Cutting and Uniqueness of Solutions", Soils and Foundations , Vol. X, No. 3, September 1970, pp. 50-65. 36. Jaeger, J.C. and Cook, N.G.W., Fundamentals of Rock Mechanics, Methuen and Co. Ltd., London, 1969. 37. Kiersch, G.A., "Vaiont Reservoir Disaster", Geo Times, Vol. 9, No. 9, 1965, pp. 9-12. 38. Kondner, Robert L. and Krizek, Raymond J., "Creep Compliance Response of a Cohesive S o i l " , Journal of the Franklin Institute, Vol. 279, No. 5, May 1965, pp. 366-373. 39. Ladd, Charles C., "Stress-Strain Modulus of Clay in Undrained Shear", Design of Foundations for Control of Settlement, ASCE, 1964, pp. 127-161. 40. Lara-Tomas, Manrique, "Time-Dependent Deformation of Clay Soils under Shear Stress", Ph.D. Dissertation, Engineering Experiment Station, The Ohio State University, November 1961. 41. Lee, E.H., "Stress Analysis in V i s c o e l a s t i c Bodies", Quarterly Journal of Applied Mechanics , Vol. 13, 1955, pp. 183-190. 42. Lin, T.H., "Bending of a Plate with Nonlinear Strain Hardening Creep", Creep in Structures, Hoff, N.J. ( E d i t o r ) , Academic Press Inc., New York, 1962, pp. 215-228. 43. Lubahn, J.D. and Felgar, R.P., Plasticity and Creep of Metals , John Wiley and Sons Inc., New York, 1961. 128 44. McCollum, P.A. and Brown, B.F., Laplace Transform Tables and Theorems, Holt, Rinehart and Winston, New York, 1965. 45. Mallawaratchie, D.P., "Plane Strain Creep Rupture of a Saturated Undis-turbed Clay", M.A.Sc. Thesis, University of B r i t i s h Columbia, June 1970. 46. Mendelson, Alexander, Plasticity: Theory and Application, The Macmillan Company, New York, 1968. 47. Mendelson, A., Hirschberg, M.H. and Manson, S.S., "A General Approach to the Practical Solution of Creep Problems", Journal of Basic Engineering, ASME, December 1959, pp. 585-598. 48. Mendelson, A. and Manson, S.S., "Practical Solution of P l a s t i c Deformation Problems in the E l a s t i c - P l a s t i c Range", NACA Technical Bote 4088, 1958. 49. M i t c h e l l , J.K., "Shearing Resistance of So i l s as a Rate Process", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 90, No. SMI, Proc. Paper 3773, January 1964, pp. 29-61. 50. M i t c h e l l , J.K., Campanella, R.G. and Singh, A., "Soil Creep as a Rate Process", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 94, No. SMI, Proc. Paper 5751, January 1968, pp. 231-254. 51. M i t c h e l l , J.K., Seed, H.B. and Paduana, J., "The Creep Deformation and Strength Characteristics of So i l s under the Action of Sustained Stress", Report No. TE-65-8, Department of C i v i l Engineering, University of C a l i f o r n i a , Berkeley, 1965. 52. M i t c h e l l , J.K., Singh, A. and Campanella, R.G., "Bonding, Effective Stresses and Strength of S o i l s " , Journal of the Soil Mechanics and Found-ations Division, ASCE, Vol. 95, No. SM5, Proc. Paper 6786, 1969, pp. 1219-1246. 53. Murayama, S. and Shibata, T., "On the Rheological Characteristic of Clay, Part I", Bulletin No. 26, Disaster Prevention Research In s t i t u t e , Kyoto University, Japan, 1958. 54. Murayama, S. and Shibata, T., "Rheological Properties of Clays", Proceed-ings of the Fifth International Conference on Soil Mechanics and Found-ation Engineering, Paris, Vol. 1, 1961, pp. 269-273. 55. Murayama, S. and Shibata, T., "Flow and Stress Relaxation of Clays", Pro-ceedings of the Rheology and Soil Mechanics Symposium, Grenoble, 1964, pp. 99-129. 56. Odquist, F.K.G. and Hult, J . , Kriechfestigkeit Metallischer Werkstoffe, Springer Verlag, B e r l i n , 1962. 57. Pao, Y.H. and Marin, J . , "Prediction of Creep Curves from Stress Strain Data", Proceedings, ASTM, Vol. 52, 1952, pp. 951-957. 58. Paterson, W.S.B., The Physics of Glaciers, Pergamon Press, London, 1969.. 129 59. Percy, J.H., Loden, W.A. and Navaratna, D.R., "A Study of Matrix Analysis Methods f o r I n e l a s t i c Structures", A i r Force Systems Command, Technical Document, Report No. RTD-TDR-63-4032, October 1963. 60. Pipes, Louis A., Applied Mathematics for Engineers and Physicists , Second Editio n , McGraw-Hill Book Company, Inc., New York, 1958. 61. Poritsky, H., "Effect of Creep on Stresses in C y l i n d r i c a l S h e l l s " , Creep in Structures, Hoff, N.J. ( E d i t o r ) , Academic Press Inc., New York, 1962, pp. 229-244. 62. Radhakrishnan, N. and Reese, Lymon C., "A Review of Applications of the Fi n i t e Element Method of Analysis to Problems in Soil and Rock Mechanics", Soils and Foundations, Volume X, No. 3, September 1970, pp. 95-112. 63. Read, W.T., "Stress Analysis for Compressible V i s c o e l a s t i c Materials", Journal of Applied Phsyics , Vol. 21, 1950, pp. 671-674.. 64. Saito, M., "Forecasting the Time of Occurrence of a Slope F a i l u r e " , Pro-ceedings of the Sixth International Conference on Soil Mechanics and' Foundation Engineering, Vol. II, 1965, pp. 537-541. 65. Saito, M. and Uezawa, H., "Failure of Soil due to Creep", Proceedings of the Fifth International Conference on Soil Mechanics and Foundation Engin-eering, Vol. 1, 1961, pp. 315-318. 66. Schapery, R.A., "Irreversible Thermodynamics and Variational Principles with Applications to V i s c o e l a s t i c i t y " , Ph.D. Thesis, C a l i f o r n i a Institute of Technology, 1962. 67. Schiffman, R.L., "The Use of V i s c o e l a s t i c Stress-Strain Laws in Soil Testing", ASTM, Special Technical Publication No. 254, Papers on S o i l s , 1959 Meetings, 1959, pp. 131-155. 68. Scott, R.F. and Ko, Hon-Yim, "Stress-Deformation and Strength Character-i s t i c s " , State of the Art Volume, Seventh International Conference on Soil Mechanics and Foundation Engineering, Mexico, 1969, pp. 1-47. 69. Sherif, M.A., "Flow and Fracture Properties of Seattle Clays", Research Series No. 1, University of Washington Soil Engineering, January 1965. 70. Singh, Awtar, "Creep Phenomena in S o i l s " , Ph.D. Thesis, University of C a l i f o r n i a , Berkeley, September 1966. 71. Singh, A. and M i t c h e l l , J.K., "A General Stress-Strain-Time Function f o r S o i l s " , Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 94, No. SMI, Proc. Paper 5728, January 1968, pp. 21-46. 72. Singh, A. and M i t c h e l l , J.K., "Creep Potential and Creep Rupture of S o i l s " , Proceedings of the Seventh. International Conference on Soil Mechanics and Foundation Engineering, Vol. 1, 1969, pp. 379-384. 73. Snead, D., "Creep Studies on an Undisturbed Sensitive Clay", Ph.D. Thesis, University of B r i t i s h Columbia, Vancouver, 1970. \ 130 74. Suklje, L., Rheological Aspects of Soil Mechanics , Wiley-Interscience, London, 1969. 75. Sutherland, W.H., "AXICRP-Finite Element Computer Code for Creep Analysis of Plane Stress, Plane Strain and Axisymmetric Bodies", Nuclear Engineer-ing and Design, Vol. 11, 1970, pp. 269-285. 76. Ter-Stepanian, G., "On the Long Term Slope S t a b i l i t y of Slopes", Norwegian Geotechnical Institute, Publication No. 52, 1963, pp. 1-13. 77. Ter-Stepanian, G., "In Situ Determination of the Rheological Character-i s t i c s of Soils on Slope", Proceedings of the Sixth International Confer-ence on Soil Mechanics and Foundation Engineering, 1965, pp. 575-577. 78. Terzaghi, K., "Mechanism of Landslides", Application of Geology to Engin-eering Practice, The Geological Society of America (Berkey Volume), 1950, pp. 83-123. 79. Vaid, Y.P., "A Plane Strain Apparatus for S o i l s " , M.A.Sc. Thesis, Univer-s i t y of B r i t i s h Columbia, 1968. 80. Varga, R.S., Matrix Iterative Analysis, Prentice-Hall, Englewood C l i f f s , New Jersey, Chapter 3, 1962. 81. Vialov, S. and Skibitsky, A., "Rheological Processes in Frozen Soils and Dense Clay", Proceedings of the Fourth International Conference on Soil Mechanics and Foundation Engineering, Vol. I, 1957, pp'. 120-124. 82. Vyalov, S.S. and Meschyan, S.R., "Creep and Long-Term Strength of Soils Subjected to Variable Load", Proceedings of the Seventh International Conference on Soil Mechanics and Foundation Engineering, Vol. 1, 1969, pp. 423-431. 83. Webber, J.P.H., "Stress Analysis in V i s c o e l a s t i c Bodies using Finite Ele-ments and a Correspondence Rule with E l a s t i c i t y " , Journal of Strain Analysis, Vol. 4, No. 3, 1969, pp. 236-243. 84. Wilson, E.L., " F i n i t e Element Analysis of Two-Dimensional Structures", Structural Engineering Laboratory, Report 63-2, University of C a l i f o r n i a , Berkeley, C a l i f o r n i a , June 1963. 85. Wilson, E.L., "Structural Analysis of Axisymmetric S o l i d s " , J.A.I.A.A. , 1965, pp. 2269-2274. 86. Yen, Bing C , " S t a b i l i t y of Slopes Undergoing Creep Deformation", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 95, No. SM4, Proc. Paper 6675, July 1969, pp. 1075-1096. 87. Zienkiewicz, O.C. and Cheung, Y.K., The Finite Element Method in Structural and Continuum Mechanics , McGraw-Hill, London, 1967. 88. Zudans, Zenons, "Survey of Advanced Structural Design Analysis Techniques", Nuclear Engineering and Design, Vol. 10, pp. 400-440, 1969. 131 APPENDIX A BRIEF DEVELOPMENT OF THE FINITE ELEMENT METHOD This appendix gives the analytical background and development of the f i n i t e element method applied to plane str e s s , plane s t r a i n , and axisymmetric problems. The displacement method development used i s guided by the Principle of Minimum Potential Energy. The various assumptions made in the development are outlined as they occur. A-l GENERAL ASSUMPTIONS AND ANALYTICAL BACKGROUND a) Geometry , • In a purely analytical formulation of a stress analysis problem i t is usually required that the boundary shape be represented by a coordinate surface. However, with the use of the f i n i t e element method, special assump-tions about the geometry are not necessary and the boundary shape is not s i g n i f i c a n t . For plane stress and plane s t r a i n problems, the continuum i s subdivided into triangular or quadrilateral elements of unit thickness as shown in Figures A-l and A-2. The quadrilateral elements are each made up of four triangular elements sharing a common node. Axisymmetric bodies are subdivided into annular elements that are triangular in cross section as indicated in Figure A-3. The boundary for each body i s then simply approximated by the straight lines j o i n i n g the element nodes at the surface. b) Loading ( The loadings of slopes and embankments can r e s u l t from externally applied forces at a boundary such as foundation loads and hydrostatic pressures; body forces such as gravity forces and acceleration forces; and, internal stresses such as in s i t u stresses and thermal stresses due to temperature change. , Iso-132 FIG. A-l TRIANGULAR FINITE ELEMENT DIVISION OF A PLANE STRESS OR PLANE STRAIN REGION. + y 1,2,3,4 ( i , j , k , l ) Q U A D R I L A T E R A L NODE N U M B E R S . 5 ( m ) Q U A D R I L A T E R A L D U M M Y N O D E NO. © i ( D ^ © T R I A N G U L A R E L E M E N T S M A K I N G UP T H E Q U A D R I L A T E R A L E L E M E N T > FIG. A-2 QUADRILATERAL FINITE ELEMENT MADE UP FROM FOUR TRIANGULAR ELEMENTS. (85) 133 t i Trz (X r 2> crE(60) crz(82) -> + r I + u + «Y F I G . A - 3 a S T R E S S E S A N D S T R A I N S F O R A X I S Y M M E T R I C B O D Y A N A L Y S I S . 4Z / F I G . A - 3 b A C T U A L T R I A N G U L A R R I N G F I N I T E E L E M E N T . + Z ^ +• r F I G . A-3c T R I A N G U L A R F I N I T E E L E M E N T A P P R O X I M A T I O N . FIG. A - 3 F I N I T E E L E M E N T IDEALIZATION FOR AXISYMMETRIC P R O B L E M S ( 2 8 ) . 134 thermal conditions are assumed throughout so that thermal stresses are not considered in the development. Also, only s t a t i c problems are considered. c) Boundary Conditions The boundary conditions can be given by prescribed surface displacements and/or prescribed surface tractions. At points on the boundary where the com-(B) ponents of the externally prescribed displacement vector are u\ , i t is required that the components of the displacement vector, u^, evaluated at the boundary, s a t i s f y : u.(x.) = u j B ) ( X j ) i , j = 1,2,3 (A:l) At points on the bounding surface where tractions are prescribed, i t - i s required that the internal stresses a . . , evaluated at the surface, balance the components of the externally applied tractions T.. where n, are the components of the unit normal vector to the surface. d) F i e l d Equations For the i n f i n i t e s i m a l displacement case, i . e . |u. . t h e components of the s t r a i n tensor e, . are given in terms of the displacements by: ' i j ' 2 < U i , j + U J , i > < A : 3 ' For s t a t i c problems, the equilibrium equations are given by: a 1 d > . + F , - 0 (A:4) where F^ are the components of the known body forces. e) Constitutive Equations General e l a s t i c behaviour i s assumed so that the relationship between stresses and strains is given by: = Dijk£ ek£ < A : 5> 135 where D^-^ i s the material tensor, and the superscript E denotes e l a s t i c behaviour. For material nonlinearity, the solution i s based on a piecewise l i n e a r i z a t i o n of equation A:5 by applying the loading in increments. Thus, a sequence of l i n e a r e l a s t i c solutions is required. A-2 FINITE ELEMENT METHOD a) P r i n c i p l e of Minimum Potential Energy Rather than working with the boundary conditions, f i e l d equations, and constitutive relationship, which give a complete statement of the problem, i t i s usually f a r easier to use variational techniques and work with an integral formulation of the boundary value problem. The displacement method used here i s based on the minimization of the potential energy of the body. For s t a t i c equilibrium, the P r i n c i p l e of Minimum Potential Energy (23) states that: of a l l the compatible displacement f i e l d s that s a t i s f y the s p e c i f i e d boundary displacements, the displacement f i e l d that also s a t i s f i e s the condition of equilibrium ( i . e . the equilibrium equations and traction boundary conditions) minimizes the potential energy. This is given by: f /WdV - /T.u.dS - /F.u.dV V S V = 0 (A:6) where the terms in brackets represent the potential energy, "V, of the body, W i s the s t r a i n energy density, V is the volume of the body and Sj i s that part of the surface, S, where tractions are prescribed. This p r i n c i p l e is v a l i d for nonlinear, as well as l i n e a r constitutive relationships, as long as the body remains e l a s t i c . b) Formulation for the Linear E l a s t i c Case The l i n e a r e l a s t i c case w i l l be examined before developing the incre-mental approach. For the l i n e a r e l a s t i c case, the s t r a i n energy density i s 136 given by: w = i D i j w e i E j 4 ( A : 7 ) and the symmetry condition D...^ = Vj^l = Dk£ji = Dk£ij n o 1 d s- ^A:8^ The entire body being considered i s subdivided into volume elements ( f i n i t e elements) which are connected to each other at the nodal points. Considering the Nth triangular f i n i t e element of Figure A-l with volume V^, assume the displacement f i e l d is a function of the spatial coordinates: Ui = I * i c (x) «. (A:9) c=i where are the assumed functions of p o s i t i o n , n is the number of nodal dis-placement degrees of freedom (6 for a triangular element), and = are the un-known generalized coordinates. The spatial functions <i>.jc(x) are to be chosen so that compatibility between elements i s assured. Equation A:9 can be given in the usual matrix form: {u} = [$]{«} (A:10) The nodal displacements {Q} are equal in number to the generalized co-ordinates and are given by: {Q} = [A]{«> ( A : l l ) in which [A] is formed by successively introducing each nodal point coordinate into equation A:10. Thus, the matrix c o e f f i c i e n t s in equation A:11 are not a function of position. Solving for the generalized coordinates: ' {«} = [A _ 1]{Q} (A:12) From equations A:3 and A:10, a column matrix of the s t r a i n components {e} i s formed in terms of the generalized coordinates. ie} = [ B ] { c c } (A:13) where the c o e f f i c i e n t s of [B] are functions of position involving <j>. and t h e i r derivatives. Using equation A:12, the s t r a i n components are also given by: 137 {£} = {Q}T [ A _ 1 ] T [ B ] T (A:14) I f [D] i s the matrix o f the material property functions D ^ ^ , then, sub-s t i t u t i n g equations A:10, A: 12 and A:14 i n t o A:6, and using the symmetry property, A:7 g i v e s : / \ {Q>J [ A " 1 ] 1 [ B ] M T [ D ] M [B ] N [ A " 1 ] , {Q}M dV - / {Q}J [A" 1ij UjJ {T}, dS J N L - J N L - J N L " J N - / '(Q>I [ A - 1 ] ! rf{F}N dV N = 0 (A:15) where the summation i s c a r r i e d out f o r a l l the elements and [<}>s] represents the s p a t i a l functions s p e c i a l i z e d to surface p o s i t i o n s . Taking v a r i a t i o n s o f equation A:15 with respect to the nodal displacements y i e l d s the d e s i r e d s t i f f n e s s r e l a t i o n s h i p : / [ A " 1 ] 1 [ B ] J [ D ] N [ B ] N [ A " 1 ] N {Q}N dV - / [ A ' 1 ^ {T> N dS N or / [A-1]J[c(>]J{F}NdV [ K ] N W N = { P>N = 0 (A:16) (A:17) where the element s t i f f n e s s matrix [K] ^ i s given by: [ K ] N = / [ A - 1 ] ^ [ B ] J [ D ] N [ B ] N [A" 1],, dV (A:18) and the element nodal load, vector {P} N i s the two r i g h t hand terms o f equation A: 16 (P} N = / C A " 1 ] ! U l J {T} w dS + / [ A " 1 ] ! Ml {F}., dV JN s N " J H N N JN (A:19) 138 A-3 FORMULATION FOR THE INCREMENTAL APPROACH To handle the nonlinear e l a s t i c problem by the incremental approach, i t is necessary to expand the Pr i n c i p l e of Minimum Potential Energy to permit the determination of an incremental state superimposed upon a body with an exis t i n g equilibrium stress state (10). For this case, the potential energy, "V, of the body from Equation A:6 becomes: " V = V I } + AV = / (W ( I ) + AW) dV - / ( T . ( I ) + AT,) ( u - ( I ) + A U , ) dS V S T - / ( F . W + A F . ) ( u n . ( I ) + A U . ) dV (A:20) V where the superscript (I) and pref i x A indicate the i n i t i a l equilibrium values of the variables and the incremental change in these variables, respectively. Expanding this equation and truncating the Taylor series expression for the change in s t r a i n energy af t e r two terms gives: " V = V 1 ) + A V = "V ( I ) + / 3W 9e 1J .E (I) ' i j A E I J + 32W 8 e i j 9 E k £ . E E A e i j A e k £ .E (I) ' i j dV - / CT, Sx (I) A U . + A T . u . ^ + A T . A U . ) dS - / ( F . ^ A U - + A F , u . ( l ) + A F . A U , ) dV V (A:21) The incremental state must be such that the s t r a i n increment i s small enough to truncate the Taylor series expansion af t e r the second term. From the pr i n c i p l e of v i r t u a l displacements: / T. ( I )Au.dS S T / F . ^ A ^ . d V = 0 (A:22) 139 and, noting that; 8W 8e (I) .E (I) i j and 32W 9 £ i j 8 e k £ Ae .E (I) ' i j kl Aa (A:23) Equation A:21 reduces to: 1/ = "V^) + / 1 Aa. .AeE,dV - / AT . ( u . ( l ) + AU,) dS - / AF. ( u , ( l ) + AU.) dV V ST V (A:24) Substituting the l i n e a r incremental s t r e s s - s t r a i n law for A a - . , A a i j = D i j k - e A e k £ (A:25) into Equation A:24 and taking variations with respect to a l l admissable incre-mental displacement f i e l d s y i e l d s the desired variational p r i n c i p l e : 6"V = 6 {/ \ Dij^Ae^.Ae^dV - / AT.A^dS - / AF-Au.dV} = 0 (A:26) V T V Comparing this equation to Equation A:6, i t can be seen that the previous form-ulation for the l i n e a r e l a s t i c case may be used where the variables e^,, a.,, T.j, F.., u i are now replaced by increments of the variables Ae!^, Aa^, AT^, A F . , AU^ A-4 INTRODUCTION OF CREEP STRAINS The increment of total s t r a i n Ae.. can be considered to be made up of an increment of e l a s t i c s t r a i n , Ae^,, and an increment of non-recoverable creep s t r a i n , Ae,,. (In the previous two sections, e l a s t i c behaviour only was con-sidered and Ae,. = Ae...) The increment of creep s t r a i n i s assumed to be an i n i t i a l s t r a i n for any time interval At and to be constant during the interval 140 . E . . C o r A £ i j = Ae. . - Ae. • (A:28) S u b s t i t u t i n g Equation A:28 i n t o Equation A:26 y i e l d s : 6 { / \ D i j k £ ( A e i j A e k £ " A e i j A4 " A e i j A e k £ + A e i j A e k V d V " / A T i A U i d S / AF.AU..dV} = 0 (A:29) Then, s u b s t i t u t i n g the incremental form o f Equations A:10, A:12 and A:14 i n t o Equation A:29, and using the symmetry property A:7, g i v e s : I N \ <AQ>JJCA-1]JI:BDJJ[DDNCB:NCA-1DN{AQ1N - { A Q ^ C A - ^ J C B ^ E D l ^ A e 0 ] + I [ D ] N [ A e C ] [ A e C ] dV - / • U Q } J C A _ 1 ] J [ * S ] J { A T } N d S >T N / ( A Q ^ C A - ^ J C ^ ^ A F ^ d V = 0 (A:30) where the summation i s c a r r i e d out f o r a l l the elements and [<(» ] represents the s p a t i a l f u n c t i o n s s p e c i a l i z e d to surface p o s i t i o n s . Taking v a r i a t i o n s o f Equation A:30 with respect to the increment o f nodal displacements and assuming that the creep s t r a i n s are constant during the time increment y i e l d s the des i r e d s t i f f n e s s r e l a t i o n s h i p : I N [ A - ^ j E B ^ r D ^ t B ^ C A - ^ ^ A Q ^ - [ A ^ l J C B l j E D j ^ A e 0 ] / [A _ 1]J[<f> s]J{AT} NdS - / [ A - ^ j M j c A F ^ d V dV N N = 0 (A: 3.1) 141 or I [K] N{AQ} N = {AP}., + {AL} (A:32) N _ J where the element s t i f f n e s s matrix [ K ] M i s given by, (A:33) {AP) N = / [A-^JC+^JUD^S + / [A^lJwJuF^dV (A:34) and the creep s t r a i n nodal 'load' vector is given by, {AU N = / [A'^JCBlJCD^CAe^dV (A:35) The value of the c o e f f i c i e n t s of the material matrix [ D ] N f o r each element depends on the stresses and e l a s t i c strains for the increment. This section completes the development of the f i n i t e element methods required in the solution of s o i l mechanics problems in which the time-dependent behaviour of the s o i l i s considered. 142 APPENDIX B B-l OPTIONS AND FEATURES OF THE AXISYMMETRIC AND PLANE STRAIN-PLANE STRESS FINITE ELEMENT PROGRAMS For most s o i l mechanics problems, the i n i t i a l stress condition w i l l generally be quite complicated and depend on the history of the p a r t i c u l a r s o i l and what changes are being made by excavating the s o i l or building on the s o i l . Also, the load condition can be quite complex and even time-dependent. The f i n i t e element programs have been written with enough f l e x i b i l i t y to deal with many of these complicating factors, a) General Aspects of the Programs The l i n e a r v i s c o e l a s t i c solution option and the general creep solution option which are the main features of the programs have been described in Chapters 4 and 5. In general, most of the additional options and features have only been developed for the plane strain-plane stress f i n i t e element program which has wider application to s o i l mechanics problems than the axisymmetric program. The axisymmetric program has only been coded for triangular ring elements. However, the plane strain-plane stress program can handle the f o l l -owing element types or combinations of element types: triangular elements, quadrilateral elements made up from 4 triangular elements as given by Wilson (85); and, incompressible quadrilateral elements as given by Christian (16). The advantage of the quadrilateral elements i s that there i s some reduction in the number of unknown nodal displacements and usually a saving of computer time. For the incompressible (undrained) quadrilateral elements, the s t r a i n s are re-lated to the e f f e c t i v e stresses by e f f e c t i v e stress moduli which are related to the s t r e s s - s t r a i n relationship for total stresses. Christian's development is not repeated here since i t is readily available (16). The incompressible 143 quadrilateral elements have been checked against simple closed-form undrained problems and excellent agreement was found. Also, Bucher (7) has checked this portion of the plane strain-plane stress program against other f i n i t e element programs and found excellent agreement. Since incremental methods can be used throughout, i t should be possible to allow a volume change at each time in t e r -val and solve the time-dependent consolidation problem in conjunction with the creep problem (13, 7). However, at present this option has not been included in the programs. For both the axisymmetric and plane strain-plane stress programs, the nodal stresses are determined by Wilson's averaging method (84). This method gives more consistent results for the nodal stresses, p a r t i c u l a r l y for surface nodes, than ju s t averaging the stresses in the elements around the node. For the axisymmetric f i n i t e element program, the circumferential s t r a i n e A w i l l vary within each element as i t i s a function of the coordinates within e J the element. To simplify the various cal c u l a t i o n s , an average value of e Q i s used in a l l of the calculations. It was found that the mean of the circumfer-ential s t r a i n determined for each corner of the triangular element gave ex-c e l l e n t results when compared with closed-form solutions for l i n e a r visco-e l a s t i c and nonlinear creep problems. Other, more complicated methods of averaging have been suggested in the l i t e r a t u r e (28, 75). The other strains are constant throughout each element f o r the displacement method derivation used in Appendix A. It is possible at each solution step to modify the e l a s t i c moduli so that nonlinear e l a s t i c problems can be solved by using tangent moduli for each; s t r e s s - s t r a i n increment. Also, a l l loadings can be applied in steps to take advantage of the incremental solution method. These features are common to both programs. It is also possible to modify the geometry of the structure to 144 account for the displacements at each solution step in the plane strain-plane stress f i n i t e element program. b) I n i t i a l Stress Conditions for the Plane Strain-Plane Stress Program For most problems in s o i l mechanics, the i n i t i a l stress condition depends on the consolidation and subsequent loading history of the p a r t i c u l a r s o i l deposit. The i n i t i a l stresses for a deposit with a horizontal surface may be s p e c i f i e d for s o i l s consolidated under the one-dimensional or KQ consolidation condition. In terms of e f f e c t i v e stresses, these i n i t i a l values are given by: a'y = YZ - u : (B:l) and a'h = (B:2) where a'^ and are the v e r t i c a l and horizontal e f f e c t i v e stresses, respect-i v e l y , y i s the s o i l s unit weight, z i s the depth below the surface, u i s the pore-water pressure and KQ is the c o e f f i c i e n t of earth pressure at rest. To introduce these i n i t i a l stresses into any problem, i t i s only necessary to specify KQ and values of at various depths. The values of c/ and are then automatically computed at the mean height of each f i n i t e element. These values of and are then the i n i t i a l stresses (T^ = o) for the p a r t i c u l a r problem. The i n i t i a l stresses can also be computed in terms of total stresses i f desired. It i s also possible to consider the i n i t i a l stresses as those given by the self-weight of the s o i l or gravity stresses. For this option, i of the element weight i s assigned to each node and an e l a s t i c f i n i t e element solution is- used to determine the i n i t i a l stresses. c) Automatic Excavation and Construction of Earth Structures for the Plane Strain-Plane Stress Program The excavation of earth structures such as slopes or s l u r r y trenches can be simulated so that the i n i t i a l conditions for the s t a r t of the creep process 145 can be determined. This procedure has been developed from the work of Brown and King and other researchers (6, 18, 19). It has been shown (6, 35) that multi-step excavations can be simulated by a single step as long as the s o i l i s l i n e a r l y e l a s t i c . However, i f the s o i l is not l i n e a r e l a s t i c or exhibits time-dependent behaviour, the f i n a l state of stress can not be obtained unless the complete excavation history i s followed. The excavation process is best explained with the aid of Figure B - l . I n i t i a l l y , the stresses throughout the ground to be excavated are given by a . On the slope surface indicated i n Figure B-la these i n i t i a l stresses w i l l have a normal component on the horizontal portion and both a normal and shear com-ponent on the in c l i n e d portion. After excavation, the slope surface i s to be stress free. This is simulated a n a l y t i c a l l y by applying changes in stress to the excavated surface that are equal in magnitude and opposite in direction to the i n i t i a l stresses on the slope surface. The changes in stress applied at the surface induce changes in stress Aa throughout the excavated slope which are calculated by a f i n i t e element solution for each excavation step. The f i n a l stress condition throughout the slope at the end of any excavation step is then (a + Aa) as shown in Figure B-lc. This process is repeated i f more than one excavation step i s required. The excavation procedure described above is carried out automatically and the nodal displacement numbering i s changed to keep the band width at a minimum as elements are removed. Other researchers (18, 19) have reduced the s t i f f n e s s of the excavated elements to simulate the excavation process. However, this means that there are many additional unknown displacements to determine for these elements. For problems in which solutions at many time intervals are required, i t is far more e f f i c i e n t to completely remove these elements. The options available in the automatic excavation procedure are: S L O P E FIG. B - l a I N I T I A L S T R E S S E S IN T H E S L O P E T O B E E X C A V A T E D . A AAA A A A A A A 77777-^77777^77777^77777^7777^777, A a F I G . B - l b C H A N G E IN S T R E S S E S R E Q U I R E D T O S I M U L A T E E X C A V A T I O N \////=/////^E/77//^F-/77, /////=7/777^r77/77^//>77=/7/7±\ cr = aQ + A C T F I G . B - l c F I N A L S T R E S S E S IN T H E S L O P E . FIG. B-l SIMULATION OF SLOPE EXCAVATION (18). 147 1. I n i t i a l stress condition - none, gravity stresses, in s i t u KQ stresses. 2. External loads can be applied at any excavation stage and on the excavated structure. 3. The material properties and geometry of the problem can be changed at any stage. Nonlinear e l a s t i c materials can be represented by tangent moduli for each stress-s t r a i n increment. It is possible to consider the f a i l -ure of any element during an excavation step by reducing i t s s t i f f n e s s . 4. Incompressible quadrilateral elements can be used in an e f f e c t i v e stress analysis. With this range of available options, many s o i l mechanics problems can be examined for a variety of i n i t i a l conditions, loadings and material behaviour. In the present program, v i s c o e l a s t i c or nonlinear creep i s not considered to s t a r t until the excavation i s completed. However, the program i s f l e x i b l e enough to consider creep during the actual excavation period i f modifications are made. Automatic construction of embankments (6, 14) i s also possible with the plane strain-plane stress program. The procedure is quite s i m i l a r to that developed for excavating slopes. During construction, the f i n a l state of stress and displacement due to the weight of added s o i l w i l l depend on how the s o i l i s added, even for li n e a r e l a s t i c materials (6). The embankment is b u i l t up to i t s f i n a l desired p r o f i l e by adding layers of s o i l and determining the change in stress for the already constructed portion of the embankment due to the weight of the new s o i l layer. The s t i f f n e s s of the new s o i l layer i s con-sidered as part of the total s t i f f n e s s of the embankment up to that point in 148 the construction. The various options available for the analysis of the exca-vation of slopes are also available f o r the analysis of the construction of embankments. d) Creep and Creep Rupture The automatic excavation or construction of earth structures provides the i n i t i a l conditions for the s t a r t of many creep problems. However, natural slopes and many other typical s o i l mechanics problems can also be examined with the programs. With the many options available, a wide range of i n i t i a l stress and loading conditions or changes in these conditions can be considered. Creep rupture i s introduced into the creep analysis by reducing the s t i f f n e s s of any element in which creep rupture has occurred. With the f l e x i b i l i t y pro-vided in the programs, i t should be possible to incorporate more detailed creep and creep rupture relationships in the programs as they are made a v a i l -able in the l i t e r a t u r e . 149 B-2 PLANE STRAIN-PLANE STRESS PROGRAM INPUT AND CAPABILITIES C * PLANE STRAIN-PLANE STRESS CONSTANT STRAIN FINITE ELEMENT C * PROGRAM FOR EXCAVATING SLOPES, OR CONSTRUCTING EMBANKMENTS. C NONLINEAR (MATERIAL) FINITE ELEMENTS CAN BE USED. C TRIANGULAR AND OUADRI LATERAL ELEMENTS CAN BE USED. C INCOMPRESSIBLE OUADRI LATERAL ELEMENTS CAN BE USED. C LOAD CASES CAN BE APPLIED TO THE STRUCTURE AT AMY STAGE. C EXC/CON'STRESS STATE CAN BE USED AS INPUT FOR INCREMENTAL C CREEP SOLUTION OR LINEAR VISCOELASTIC SOLUTION. C (NOTE. ONLY TRIANGULAR ELEMENTS CAN HE USED FOR LINEAR C VISCOELASTIC SOLUTION.) C * NOTE. ALL NODES (INCLUDING DUMMY NODES) SHOULD BE C NUMBERED FOR MINIMUM BAND WIDTH. C ALL UNITS MUST BE COMPATIBLE. C PROGRAM USES FILES (TAPES) 2, 3, AND k. C ( F I L E S 2 AND k FOR SCRATCH, F I L E 3 FOR OUTPUT TO BE SAVED.) C * LIMITATIONS ON PRESENT PROGRAM CAPACITY. C DOUBLE PRECISION. C MAXIMUM NO. OF NODES IS 250. C MAXIMUM NO. OF ELEMENTS IS kOO. C HALF BAND WIDTH X NO. OF, UNKNOWNS LESS THAN 20000. C * DESCRIPTION AND FORMAT FOR INPUT DATA. C T I T L E (20Afc) C A DESCRIPTION OF THE PROBLEM. C NPP,NN,NE,NUE,NSS,NBB,NEN,NCI , NG I ,NLC,KV,NGC (1216) C NNTT,NC PR,K PC,J EU,NS TAT,NCRE,NS R,NS PO C NPP PROBLEM NO. C NN NO. OF NODES. (INCLUDING OUADRI LATERAL NODES.) C NE NO. OF ELEMENTS. C NUE =1 ELEMENTS HAVE THE SAME INITIAL PROPERTIES. C =2 UP TO 20 SETS OF INITIAL PROPERTIES. C =3 EACH ELEMENT HAS INITIAL PROPERTIES GIVEN. C NSS =1 PLANE STRAIN SOLUTION. C =2 PLANE STRESS SOLUTION. C NBB =1 FOR EXCAVATION PROPLEM. C =2 FOR CONSTRUCTION PROBLEM. C MEN =1 LINEAR ELASTIC MATERIAL. C =2 NONLINEAR ELASTIC MATERIAL. C NCI NO. OF CONSTRUCTION OR EXCAVATION INCREMENTS. C (INCLUDING NSTAT=1.) C NGI =0 NO GRAVITY LOADING CONSIDERED. C =1 GRAVITY LOADING CONSIDERED. (CONSTRUCTION) C =2 INS ITU STRESS UNLOADING CONSIDERED. C (EXCAVATION ONLY, NSTAT=1 FOR THIS CASE.) C =3 GRAVITY STRESS UNLOADING CONSIDERED. 150 C (EXCAVATION ONLY, NSTAT=0 FOR THIS CASE.) C NLC MO. OF LOAD CASES ON THE COMPLETE STRUCTURE. C KV =0 NO CREEP C =1 CREEP SOLUTION REQUIRED. C NGC =1 GEOMETRY NOT CHANGED. C =2 GEOMETRY CHANGED AFTER EACH INCREMENT. C ... NNTT =0 NO INFORMATION STORED ON TAPE FOR FUTURE USE. C =1 DISPL., STRAINS AMD STRESSES STORED Of! TAPE. C ( F I L E (TAPE) 3 IS USED FOR STORAGE.) C NCPR =0 ALL INFORMATION PRINTED OUT. C =1 INFO. FOR SPECIFIED NODES AND ELEMENTS PRINTED. C (APPLIES TO CREEP SOLUTIONS ONLY) C KPC =0 CREEP RUPTURE NOT CONSIDERED. C =1 CREEP RUPTURE CONSIDERED BY REDUCING STIFFNESS. C =2 CREEP RUPTURE CONSIDERED BY CHANGING MODULI. C = 3 ALL MODULI CHANGED FOR NONLINEAR AND RUPTURE. C JEU =0 NO ITERATION TO CORRECT MATERIAL PROPERTIES. C = MAXIMUM NO. OF ITERATIONS ALLOWED. C NSTAT =0 STATIC INSITU STRESSES NOT CONSIDERED. C =1 STATIC INS ITU STRESSES CONSIDERED. C NCRE =0 TOTAL CHANGE IM STRESS NOT REQUIRED FOR CREEP. C =1 TOTAL CHANGE IN STRESS REQUIRED FOR CREEP. C NSR =0 STATIC RUPTURE NOT CONSIDERED. C =1 STATIC RUPTURE CONSIDERED BY REDUCING STIFFNESS. C =2 STATIC RUPTURE CONSIDERED BY CHANGING MODULI. C =3 ALL MODULI CHANGED FOR NONLINEAR AND RUPTURE. C NSPO NO. OF TIME STEPS BETWEEN PRINT-OUT FOR CREEP C SOLUTIONS. C INCM =0 NO INCOMPRESSIBLE QUADS USED. C =1 INCOMPRESSIBLE QUADS USED. C (ND( I , J ) , J = 1 , 3 ) , I =1,NN) ( 1 8 U ) C ARE THE ALLOWABLE DISPLACEMENTS FOR EACH NODE. C I IS THE NODE NO. C J IS THE DIRECTION AND TYPE OF NODE. C (1 IS X DIR. CONTROL, 2 IS Y DIR. CONTROL, C 3 IS TYPE OF NODE.) C ND(I,1)=0 THEN NODE I IS FIXED IN X DIR. C ND(I,1)=1 THEN NODE I IS FREE IN X DIR. C ND(I,2)=0 THEN NODE I IS FIXED IN Y DIR. C ND(I,2)=1 THEN NODE I IS FREE IN Y DIR. C ND(I,3)=0 THEN MODE IS A REGULAR NODE. C ND(I,3)=1 THEN MODE IS A INCOMPRESSIBLE QUAD. C (DUMMY) NODE, (AND ND(I,1)=1, ND(I,2)=G). C ND(I,3)=2 THEN NODE IS A QUADRILATERAL (DUMMY) C NODE, (AMD ND(I,1)=0, ND(I,2)=0). C (SM(KK),KK=1,NN2) «NN2=2*NN' (12F6.0) C NODAL CO-ODS ARE TEMPORARILY STORED IN SM. C X AND Y CO-OD FOR EACH NODE IN ORDER. C (LEAVE CO-ODS FOR QUADRILATERAL (DUMMY) NODES AS C ZERO, THEY ARE CALCULATED IN THE PROGRAM.) C IF NUE=1, E l ( 1 ) , U I ( 1 ) , U N (1F12.0,1F6.3,1F6.0) 151 C E l ( 1 ) YOUNGS MODULUS. C U I ( 1 ) POISSONS RATIO. C UN UNIT WEIGHT. C IF NUE=1, N I I / N J J / N K K / N L L / N N Q / M E T (1216) C FOR EACH ELEMENT (TRIANGLES AND QUADS.) IN ORDER. C ... NLL AND NNQ = 0 FOR TRIANGULAR ELEMENTS. C NNQ IS QUADRILATERAL (DUMMY) NODE MO. C MET IS A DUMMY VARIABLE, LEAVE 0. C IF NUE=2, NEC (1216) C NEC NO OF SETS OF INITIAL PROPERTIES. C E l ( I ) , U I ( I ) , U N I ( I ) l=l,NEC (1F12.0,1F6.3,1F6.0) C E I. ( I ) YOUNGS MODULUS. C U l ( I ) POISSONS RATIO. C U N K I ) UNIT WEIGHT. C IF NUE=2, N I I / N J J / N K K / N L L / N N Q / M E T (1216) C FOR EACH ELEMENT (TRIANGLES AND QUADS.) IN ORDER. C MET IS THE INITIAL PROPERTIES INDICATOR. C IF NUE=3, N I I / N J J / N K K / N L L / N M O / M E T / E ( I ) , U ( I ) , U N C ( 6 1 6 / l F 1 2 . 0 / l F 6 . 3 / l F 6 . 0 ) C FOR EACH ELEMENT (TRIANGLES AND QUADS.) IN ORDER. C MET IS A DUMMY VARIABLE, LEAVE 0. C IF NSR=1 C RSSR (1F12.0) C REDUCTION IN ELEMENT STIFFNESS AFTER STATIC RUPTURE. C IF NCI GREATER THAN 1, (MD(I I ) , I I=1,MN) (1216) C MD(II) STAGE OF CONSTRUCTION 0,1,2...(NCI-1) C OR EXCAVATION (NCI-1)...2,1,0 C THAT THE NODE BELONGS TO. C (INCLUDE QUADRILATERAL (DUMMY) NODES.) C IF f!STAT=l (INS ITU STRESSES REQUIRED) C KK,NE2,CKO ( 2 I 6 / 1 F 6 . 0 ) C KK NO. OF LAYERS FOR INSITU STRESS CALCULATION. C NE2 LAST ELEMENT THAT INS ITU STRESSES REQUIRED FOR. C CKO COEFFICIENT OF EARTH PRESSURE AT REST. C ( U l ( J J ) / J J = 1 / K K ) (12F6.0) C Y CO-ODS AT BOTTOM OF LAYERS. C ( E l ( J J ) , J J = 1 , K K ) (12F6.0) C INITIAL INSITU VERTICAL STRESSES AT BOTTOM OF LAYERS. C * SUBROUTINE CNODE REQUIRES THE FOLLOWING DATA C IF NCI GREATER THAN 1. (INCLUDING NSTAT=1 CASE.) C N E1,N E 2,N L C E (NCI SETS) (1216) C NE1 FIRST ELEMENT OF STAGE. C NE2 LAST ELEMENT OF STAGE. C NLCE =0 NO EXTERNAL APPLIED LOAD CASE ON CON/FXC STEP. C =1 EXTERNAL APPLIED LOAD CASE ON CON/EXC STEP. C (INCLUDED IN THE TOTAL LOAD AT THAT STEP AND AS C PART OF THE TOTAL LOAD WHEN CON/EXC COMPLETE.) C IF NLCE=1 C NNL,NLI (1216) C NNL NO. OF NODES LOADED. 152 C NLI NO. OF INCREMENTS FOR LOAD. C J N U / F F ( 1 ) / F F ( 2 ) (MNL SETS) (116,2F12.3) C JMU NODE NO. C F F ( 1 ) X LOAD ON NODE. (POSITIVE TO RIGHT) C F F ( 2 ) Y LOAD ON NODE. (POSITIVE UP) C * FOR EACH EXTERNAL LOAD CASE. C ' NNL,NLI,NLS (1216) C NNL NO. OF NODES LOADED. C NLI NO. OF INCREMENTS FOR LOAD. C NLS =1 NOT ADDED TO TOTAL CON/EXC LOADING. C =2 ADDED TO TOTAL CON/EXC LOADING. C * SUBROUTINE NLOAD REOUIRES THE FOLLOW!NG DATA. C J N U / F F ( 1 ) / F F ( 2 ) (NNL SETS) ( 1 I 6 / 2 F 1 2 . 3 ) C JNU NODE NO. C F F ( 1 ) X LOAD ON NODE. (POSITIVE TO RIGHT) C F F ( 2 ) Y LOAD ON NODE. (POSITIVE UP) C * IF KV=1 (DATA FOR CREEP SOLUTION REQUIRED). C T I T L E (20AU) C A DESCRIPTION OF THE CREEP PROBLEM. C KVT,KVA,KVS,TDM ( 3 I 6 / 1 F 1 2 . 0 ) C KVT =1 INCREMENTAL CREEP SOLUTION. C =2 VISCOELASTIC SOLUTION USING TRANSFORM C METHOD AND NEW LOWER TRANSPOSE FOR EACH STEP. C =3 VISCOELASTIC SOLUTION USING TRANSFORM C METHOD AND MATRIX ITERATION WITH CHOLESKI. C =l» SPECIAL PROBLEM. C (MOTE. FOR LINEAR VISCOELASTIC, ONLY TRIANGULAR C ELEMENTS CAN BE USED.) C (NOTE. FOR LINEAR VISCOELASTIC, STRUCTURE MAY BE C NON-HOMOGENEOUS, BUT THE SAME LINEAR VISCOELASTIC C MODEL MUST BE APPLIED TO ALL ELEMENTS.) C KVA =1 GRAVITY LOADING STRESSES/LOADS ARE INITIAL C STRESSES/LOADS FOR CREEP. C =2 FINAL LOAD CASE STRESSES/LOADS ARE INITIAL C STRESSES/LOADS FOR CREEP. C =3 GRAVITY PLUS FINAL LOAD CASE STRESSES/LOADS ARE C INITIAL STRESSES/LOADS FOR CREEP. C =k UNLOADING STRESSES/LOADS ARE INITIAL C STRESSES/LOADS FOR CREEP. C =5 UNLOADING PLUS FINAL LOAD CASE STRESSES/LOADS C ARE INITIAL STRESSES/LOADS FOR CREEP. C (NOTE. FOR KVA = l , 3 , l | , 5 IF NCRE = 1, CHANGE IN STRESS C /LOADS DUE TO GRAVITY LOADING OR UNLOAD ING USED.) C (NOTE. FOR KVA=2,3,5 THIS MUST BE THE ACTUAL FINAL C (TOTAL) LOAD CASE SOLVED.) C KVS NUMBER OF SOLUTION STEPS REQUIRED. C TDM MAXIMUM TIME TO BE EXAMINED (INCREMENTAL CREEP). C * IF NCPR=1 (ONLY SPECIFIED DATA TO BE PRINTED). C (NNPR(I),1=1,NN) (1216) C ARE THE NODE PRINT OUT CONTROL. FOR EACH NODE IN ORDER, C =0 PRINT RESULT FOR NODE. 153 C =1 NO RESULTS PRINTED FOR NODE. C (N E PR(I),I=1,M E) (1216) C - ARE THE ELEMENT PRINT OUT CONTROL. C =0 PRINT RESULTS FOR ELEMENT. C =1 NO RESULTS PRINTED FOR ELEMENT. C * IF KPC=1 C RCCR (1F12.0) C REDUCTION IN ELEMENT STIFFNESS AFTER CREEP RUPTURE. C * IF KVT=1 NO FURTHER DATA IS REQUIRED. C * IF KVT=2,3 C ( U I ( I ) , I = 1 , K V S ) (6F12.0) C TOTAL TIME FROM THE START OF CREEP FOR EACH STEP. C NSVC,NOHOM (1216) C NSVC =0 LINEAR VISCOELASTIC STRESSES NOT REQUIRED. C =1 LINEAR VISCOELASTIC STRESSES REQUIRED. C NOHOM =0 HOMOGENEOUS STRUCTURE FOR LINEAR VISCOELASTIC. C =1 NON-HOMOGENEOUS STRUCTURE FOR LINEAR VI SC. C IF NOHOM=1 C ( F S K I ) / l = l / N E 2 ) (12F6.0) C RATIO OF THE INITIAL E VALUES OF EACH ELEMENT. C * IF KVT=3 C NJV,ERRV (1I6,1F12.8) C NJV MAXIMUM NUMBER OF ITERATIONS ALLOWED. C ERRV MAXIMUM ALLOWABLE ERROR. (SUM ERR**2/SUM DISP**2) C * CHANGES IN FORMAT FOR LARGE STRUCTURES. C FORMAT 0 (1X,5I6,6F6.2,1E12 . 1»,1F6.3,1F6.0,1F6.3) C ** SUBROUTINES AND SUB-PROGRAMS REQUIRED FOR INCREMENTAL CREEP C SOLUTION ARE NOT DESCRIBED HERE. C ** SUBROUTINES REQUIRED FOR LINEAR VISCOELASTIC SOLUTION C ARE NOT DESCRIBED HERE. 154 APPENDIX C C EQUATIONS REQUIRED FOR THE DETERMINATION OF THE COMPONENTS OF THE CREEP STRAIN INCREMENT This appendix gives the equations used to determine the components of the creep s t r a i n increment from Equation 4:1. r Ae^ " A s i i = J — S i i < C : 1> where A e . . = the components of the creep s t r a i n increment a = the equivalent stress (— T , . ) ^ 2 OCu C 1 c A e g = the equivalent creep s t r a i n increment (— A y o c t ) S.. = the stress deviator tensor. C-l EQUATIONS FOR THE PLANE STRAIN CASE 1 a e = i [ ( a x - ayy + ( a y - o J * + (a, - o,)* + G ^ ] 2 (C:2) ^ = ^ ^ x ' 2 + < A Ey' 2 + A e x A £ y + C A e e C f. C ) xy 2 1 2 F (C:3) A N D A E X = W ( 2 a x " a y " az) <C:4' e 17 Ae^ A e y = i ^ ( 2 a y - ° z - ^ < C : 5> Ae ^ ' A Y x y = 3 ^ f T x y ( C : 6 ) C-2 EQUATIONS FOR THE PLANE STRESS CASE 1 a = ( a 2 + a 2 - a a + 3 T 2 ) 2 (C:7) e x y x y xy' v ' 155 Ae g i s given by Equation C:3 and (C:8) Ae£ e J (C:9) 'xy A e e (C:10) C-3 EQUATIONS FOR THE AXISYMMETRIC CASE ae l = l~ [ ( a r - a Q ) 2 + (ff - a ^ 2 + (a, - a , ) 2 + 6 t ^ ] 5 v2 1 ( C : l l ) = / f C(Ae^) 2 + (Ae^) 2 + (Ae^) 2 + . i ( A y ^ ) 2 ] 5 (C:12) and = 2 a e { 2 a r ~ a e " a z ) (C:13) AeC = (lan - a - a ) 2 a e v 6 r z' (C:14) Z A e e = -— (la - a - cr) 2a e z r e ' (C:15) ' rz Ae£ = 3 T a e rz (C:16)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Finite element analysis of creep problems in soil mechanics
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Finite element analysis of creep problems in soil mechanics Emery, John Joseph 1971
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Finite element analysis of creep problems in soil mechanics |
Creator |
Emery, John Joseph |
Publisher | University of British Columbia |
Date Issued | 1971 |
Description | Application of the finite element method of stress analysis to problems in soil mechanics has enabled the engineer to gain much information on the deformations and stresses in earth structures under the assumption that the properties of soil are independent of time. However, both laboratory studies and field observations indicate that time is a very important factor in the behaviour of cohesive soils. The finite element method is extended to deal with typical problems in soil mechanics in which the time dependence of the mechanical properties of soil is considered. From creep studies reported in the literature, stress-strain-time relationships that have been developed for essentially uniaxial constant stress creep tests are extended to the multiaxial changing stress condition that is applicable in situ. The finite element method is used to examine problems where the soil is assumed to be linear viscoelastic. The correspondence rule of linear viscoelasticity is used as part of this development. This is an idealized case and will generally be limited to qualitative studies. Then, the incremental initial strain finite element method is developed to deal with general soil creep as described by empirical relationships. A cumulative creep law based on the strain-hardening rule is adopted in this .analysis. Creep rupture is introduced into the general incremental solution method by reducing the stiffness of elements of the earth, structure that have failed. The failure criterion used in the analysis is based on the total elapsed time from the beginning of creep and the creep strain rate. Some typical examples of problems in soil mechanics are examined to illustrate the methods developed. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-03-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050546 |
URI | http://hdl.handle.net/2429/33163 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1971_A1 E44.pdf [ 7.43MB ]
- Metadata
- JSON: 831-1.0050546.json
- JSON-LD: 831-1.0050546-ld.json
- RDF/XML (Pretty): 831-1.0050546-rdf.xml
- RDF/JSON: 831-1.0050546-rdf.json
- Turtle: 831-1.0050546-turtle.txt
- N-Triples: 831-1.0050546-rdf-ntriples.txt
- Original Record: 831-1.0050546-source.json
- Full Text
- 831-1.0050546-fulltext.txt
- Citation
- 831-1.0050546.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050546/manifest