UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Nonlinear temperature and consolidation analysis of gassy soils Vaziri-Zanjani, Hans Hamid 1986

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-UBC_1986_A1 V39.pdf [ 21.58MB ]
JSON: 831-1.0062911.json
JSON-LD: 831-1.0062911-ld.json
RDF/XML (Pretty): 831-1.0062911-rdf.xml
RDF/JSON: 831-1.0062911-rdf.json
Turtle: 831-1.0062911-turtle.txt
N-Triples: 831-1.0062911-rdf-ntriples.txt
Original Record: 831-1.0062911-source.json
Full Text

Full Text

NONLINEAR TEMPERATURE AND CONSOLIDATION ANALYSIS OF GASSY SOILS by HANS HAMID VAZIRI-ZANJANI A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA June, 1986 © HANS HAMID VAZIRI-ZANJANI, 1986 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the requirements f o r an advanced degree a t the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head o f my department or by h i s or her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . Department of C i v i l Engineering The U n i v e r s i t y of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 October 1 5 , 1986 ABSTRACT A study is undertaken to formulate and solve the equations governing the time dependent response of gassy soils. The formulations have been implemented into a finite element program capable of analyzing stress state and flow phenomenon under a variety of boundary conditions. The validity of such a program has been verified by comparing its results with closed form solutions. This program has then been applied to simulate the processes involved in depleting o i l sand reservoirs in order to give some insight into the mechanism causing fluid flow and sand production. The theory developed couples the effects of both stress and flow. It takes account of the varying permeability and compressibility of the pore fluid, and the nonlinear stress-strain behaviour of the soil. A hyperbolic model is employed to represent the soil's stress-strain characteristics and a distinction in behaviour is made between shear and tensile failures. Various schemes are proposed for transferring loads which violate the yield criterion and in ensuring that a reasonable behaviour is modelled during failure. In order to compute the change in pore pressure in gassy soils, both under undrained and transient states, the concept of a homogenized compressible phase is introduced which is used to treat a multiphase soil system as a two phase material. Such a hypothesis is found to be highly akin to the procedure normally followed in finite element analysis since i t replaces the compressibility of fluid and solid phases by one phase which occupies the entire soil volume. Assuming that the gases are present only in the form of bubbles within the liquid phase, the compressibility of the fluid phase is obtained by giving due consideration to the mixture of - i i -liquid and gas phases using Boyle's law and Henry's law and taking account of the surface tension effects. Under undrained conditions the pore pressure is computed by invoking volumetric compatibility between the soil skeleton and the compressible phase. Under transient conditions, the pore pressure is calculated by using Biot's theory of consolidation and modifying i t to account for a soil with an incremental stress-strain law and a compressible fluid phase. Formulations are derived to compute the change in pore pressure and effective stress as a result of changes in temperature and a methodology is proposed for implementing these effects into f i n i t e element analysis. Various numerical techniques are incorporated for increasing the accuracy, efficiency, and stability of the finite element procedure. The computer program based on these formulations is verified by comparing the computational results with known solutions for several problems. Application of the finite element program to analyze the problem of unloading a cavity in o i l sand reservoir has revealed that the principal factor leading to fluid production is the compressibility attained by the pore fluid as a result of the gas evolution. It is also demonstrated that large movements develop around uncased wellbores when the internal fluid pressure is reduced below the in-situ pore pressure. The factor which governs the overall stability of the reservoir is the flow rate which is a function of the pore pressure gradient, soil strength properties, perme-ability, and the volume of evolved gases. - i i i -II TABLE OF CONTENTS Page Abstract i i List of Tables v i i List of Figures v i i i Symbols x i Acknowledgements xvi CHAPTER 1 - GENERAL INTRODUCTION AND SCOPE 1 1.1 O i l Sand Characteristics 3 1.2 Scope 6 CHAPTER 2 - EFFECTIVE STRESS ANALYSIS OF GASSY SOILS 9 2.1 Introduction 9 2.2 Literature Review 10 2.3 Fundamental Laws Describing the Behaviour of Liquid/Gas Mixture 15 2.3.1 Gas Laws 15 2.3.2 Rate of Gas Diffusion 17 2.3.3 Surface Tension and Capillary Effects 17 2.3.4 Compressibility of a Miscible Gas/Liquid Mixture ... 20 2.4 Computation of Pore Pressure Under Undrained Conditions ... 23 2.4.1 Analytical Rationale 26 2.5 Summary 29 CHAPTER 3 - REPSONSE OF UNSATURATED SOILS SUBJECTED TO THERMAL EFFECTS 32 3.1 Introduction 32 3.2 Effect of Temperature on O i l Sand 34 3.3 Prediction of Pore Pressure and Volume Changes Due to Temperature Variations 39 3.3.1 Literature Review 39 3.3.2 Theoretical Analysis 44 3.3.2(a) Free Case (Aa=0) 48 3.3.2(b) Restrained Case (A6=0) 50 3.4 Discussion on the Application 51 3.5 Summary and Conclusions 53 CHAPTER 4 - STRESS-STRAIN-STRENGTH MODEL 54 4.1 Introduction 54 4.2 Nonlinear Elastic Stress-Strain Formulation 56 4.3 Limitations of the Incremental Elastic Models 59 4.4 Shear-Volume Coupling 60 4.5 Stress Transfer Concept 62 4.6 Methodology for Representing the Behaviour at Failure 64 4.7 Summary and Conclusions 66 - iv -TABLE OF CONTENTS (Continued) Page CHAPTER 5 - NONLINEAR FINITE ELEMENT CONSOLIDATION ANALYSIS OF UNSATURATED SOILS 68 5.1 Scope 68 5.2 Consolidation Phenomenon 68 5.2.1 Introduction 68 5.2.2 Physical Description 69 5.2.3 Theoretical Considerations 70 5.2.4 Literature Review 73 5.3 The Proposed Approach for Coupled Consolidation Analysis of Unsaturated Soils 77 5.3.1 Introduction 77 5.3.2 Representation of Permeability 78 5.3.3 Analytical Formulation 80 5.4 Finite Element Equations for a Soil with an Elastic Skeleton 84 5.5 Extension of Finite Element Equations for a General Incremental Law 90 5.6 Discussion 95 5.7 Summary and Conclusions 99 CHAPTER 6 - DESCRIPTION OF SOME ASPECTS OF THE FINITE ELEMENT PROCEDURE EMPLOYED 100 6.1 Introduction 101 6.2 Selection of the Elements 105 6.3 The Adopted Procedure for Generation of Stiffness Matrix Triple Product 108 6.4 The Solution of the Simultaneous Equations 109 6.5 Description of the Element Order Optimizer Incorporated ... 110 6.6 The Proposed Technique for the Nonlinear Analysis 110 6.6.1 Introduction I l l 6.6.2 The Basic Euler Procedure 114 6.6.3 Modified Euler Method 115 6.6.4 The Improved Modified Euler Method 116 6.7 Summary and Conclusions 119 CHAPTER 7 - VERIFICATION AND APPLICATION OF THE FINITE ELEMENT PROGRAM 119 7.1 The Aspects Checked by Cheung 119 7.2 One Dimensional Consolidation Analysis of Unsaturated Soils 121 7.2.1 Introduction 121 7.2.2 Finite Element Representation 123 7.2.3 Analyses and Results 123 7.2.4 Discussion 130 7.3 Two Dimensional Consolidation Analysis 130 - v -TABLE OF CONTENTS (Continued) Page 7.3.1 Introduction 130 7.3.2 Finite Element Idealization 131 7.3.3 Analyses and Results 132 7.4 Verification of the Temperature Analysis 135 7.5 Application to Selected Geotechnical Problems 138 7.5.1 Model Description 140 7.5.2 Factors Affecting the Fluid Production 144 7.5.3 Factors Affecting the Reservoir Stability 147 Introduction 147 Analyses and Results 148 Discussions 158 7.6 Summary and Conclusions 159 CHAPTER 8 - CONCLUDING REMARKS 162 REFERENCES 168 APPENDICES A -: Some General Descriptions and Basic Definitions 180 B - Equations Defining the Compressibility of Gases and Liquids 183 C - Finite Element Implementation of the Undrained Analysis .... 187 D - Derivation of the Stiffness Matrix for Soil Elements in terms of effective stresses 190 E - Finite Element Procedure for Thermal Analysis 193 F - Method for Inducing Coaxiality Between Strains and Stresses 196 G - Methods for Yield Surface Correction 198 H - The Semi-Empirical Permeability Formulation 206 I - The Front Solution Method 212 J - Derivation of the One-Dimensional Consolidation Solution for Unsaturated Soils 222 K - Validation of the Finite Element Procedure for Thermal Effects 229 L - Program User Manual 'ENHANS' 234 M - Mesh Generator User Manual 'GILGA' 338 N - Listing of the Finite Element Code 359 - vi -LIST OF TABLES Page Table 1.1 L i s t of Active and Deactivated Nodes 218 1.2 Process of Marking the Last Appearance of the Nodes 218 1.3 Order of the Nodes Occupying the Vacated Positions 218 - v i i -LIST OF FIGURES Page Figure 1.1 Composition of In-Situ Oil Sand 4 3.1 Thermal Expansion of Oil Sand Samples. (After Scott and Kosar, 1982) 36 3.2(a) Thermal Expansion of Sand Grains and Soil Struture in Drained Oil Sand Tests (After Scott and Kosar, 1982) 38 3.2(b) Thermal Volume Change of Soil Structure (After Scott and Kosar, 1982) 38 3.3(a) Thermal Volume Change of Clay Particles and Soil Structure (After Scott and Kosar, 1985) 40 3.3(b) Thermal Volume Change of Clay Soil Structure (After Scott and Kosar, 1985) 40 3.4 Dynamic Viscosity and Hydraulic Conductivity of Oil Sands (After Scott and Kosar, 1984) 41 3.5 Conceptual Model of Soil Constituents 45 4.1 Permissible and Non-Permissible Stress States 62 5.1 Elastic and Hydraulic Boundary Conditions for a Body of Volume V 80 5.2 Location of Displacement and Pore Pressure Nodes 97 6.1 Simple Incremental Approximation 113 6.2 Modified Euler Scheme 114 6.3 Improved Modified Euler Scheme 116 7.1 Stresses and Displacements Around Circular Opening in an Elastic Material (After Cheung, 1985) 120 7.2 Stresses and Displacements Around Circular Opening in an Elastic-Plastic Material (After Cheung, 1985) 122 7.3(a) Comparisons of Predicted and Observed Pore Pressure (After Cheung, 1985) 124 7.3(b) Comparisons of Predicted and Observed Strains (After Cheung, 1985) 125 7.4(a) Stresses Around a Wellbore in an Elastic-Plastic Material (After Cheung, 1985) 126 - v i i i -LIST OF FIGURES (Continued) Page Figure 7.4(b) Stresses Around a Wellbore in an Elastic-Plastic Material (After Cheung, 1985) 127 7.5 Finite Element Mesh 128 7.6(a) Degree of Consolidation for a 1-D Soil Element 129 7.6(b) Dissipation of Excess Pore Pressure at the Middle of a 1-D Soil Element 130 7.7 Finite Element Mesh of a Footing on a Thin Layer 131 7.8(a) Settlement of a Circular Footing on a Layer of Finite Thickness 133 7.8(b) Degree of Settlement of a Circular Footing on a Finite Layer 133 7.9(a) Settlement of a Strip Footing on a Layer of Finite Thickness 134 7.9(b) Degree of Settlement of a Strip Footing on a Finite Layer 134 7.10 Finite Element Mesh Subjected to Non-Uniform Temperature Field 138 7.11 Stress and Displacement Variation with Distance 139 7.12 Geometry of the Problem 141 7.13 Finite Element Mesh of a Wall Through a Layer of Oil Sand .. 143 7.14 Variation of Flowrate With Time 146 7.15 Unloading Sequence 150 7.16 Stress Distributions Following a Well Pressure Reduction to 2000 kPa. (c' = 35 kPa, H = 0.3) 152 7.17 Stress Distributions Following 47 Days Drainage at a Well Pressure of 2000 kPa. (c* = 35 kPa, H = 0.3) 153 7.18 Stress Distributions Following a Further Reduction of the Well Pressure to 500 kPa. (c' = 35 kPa, H = 0.3) 154 7.19 Stress Distributions Following 130 Days Drainage at a Well Pressure of 500 kPa. (c' = 35 kPa, H = 0.3) 155 7.20 Total Fluid Production 156 - ix -LIST OF FIGURES (Continued) Page Figure 7.21 Total Volume of Displaced Sand 157 E.l Equivalent Representation of Stresses and Strains Due to Thermal Effects 193 G.l Mohr-Coulomb Failure Surface 198 G.2 Correction of Stresses Normal to the Yield Surface 199 G.3 Incorrect and Corrected Mohr Circles for the Plane Strain Case 202 G.4 Incorrect and Corrected Mohr Circles 202 G. 5 Uncorrect and Corrected Mohr Circles 203 H. l Variation of Permeability With Void Ratio (After Lambe and Whitman, 1969) 210 I. 1 Finite Element Mesh 216 1.2 Arrangement of the Terms in the Available Core Space 221 J.l Isochrones in a Bed of 1-D Consolidating Soil 226 K.l(a) Axisymmetric (Triaxial Test) 230 K.l(b) One-Dimensional 230 K.l(c) Fully fixed 230 - x -LIST OF SYMBOLS English Letter B Bulk modulus B Compressible phase bulk modulus cp B^  Fluid bulk modulus B Gas bulk modulus g B^  Liquid bulk modulus B Solids bulk modulus s B , Solid skeleton bulk modulus sk B Tangent bulk modulus t ( c Cohesion intercept c Coefficient of consolidation v C Compressible phase compressibility C^  Fluid compressibility C Gas compressibility g C Liquid compressibility Cg Solids compressibility C , Soil skeleton compressibility S iC D' Drained constraint modulus e Void ratio e. Volume of dissolved gas "g e ^ Volume of free gas e Total gas volume e„ Volume of liquid £ e Volume of solids s - xi -e , Total soil volume sk E Young's modulus E^ Initial Young's modulus E^ Tangent Young's modulus G Shear modulus G^  Tangent shear modulus H Henry's solubility constant k Coefficient of permeability k Bulk modulus number B kg Young's modulus number m Bulk modulus exponent m Coefficient of volumetric compression v n Young's modulus exponent n porosity P Absolute pore pressure P^  Atmospheric pressure Q Quantity of flow q Critical flow rate c r c Capillary radius r Capillary radius at saturation c s R Ratio of constraint modulus to fluid bulk modulus R^  Failure ratio S Degree of saturation S Total boundary surface S, Portion of the boundary surface subjected to applied d displacements S^  Lowest degree of saturation at which liquid begins to flow S^  Portion of the boundary surface subjected to no flow - x i i -Portion of the boundary surface subjected to applied po pressure Portion of the boundary surface subjected to applied tractions Time Absolute temperature Surface tension Time factor Excess pore pressure Capillary pore pressure Compressible phase pore pressure Pore fluid pressure Pore gas pressure Liquid/gas saturation pressure Average degree of consolidation Superficial fluid velocity Volume of the body Volume of liquid Volume of voids Weight of gas Weight of dissolved gas tters Integration rule parameter Coefficient of the volumetric thermal expansion of the gas Coefficient of the volumetric thermal expansions of the liquid Coefficient of the volumetric thermal expansions of the solids - x i i i -a Coefficient of structural volume change st y Unit weight Y shear strain Fluid unit weight e Normal strain e volumetric strain v (e ) Volumetric strain of the fluid phase v cp (e ) , Volumetric strain of the soil skeleton v sk \i Viscosity v Poisson's ratio o Total normal stress a' Effective normal stress a Mean normal stress m T Shear stress <|> Angle of internal friction Matrix Notation B Matrix of shape function derivatives B Matrix of displacement shape function derivatives u Bp Matrix of pore pressure shape function derivatives D,D. ., , Matrix of elastic coefficients in the generalized Hooke's law i j k l F,F^ Vector of body forces k^j Coefficients of permeability in the generalized Darcy's law K Stiffness matrix Incremental tangent stiffness m Array containing unity and zero - xiv -n i Normal unit vector N Shape function N u Displacement shape function N P Pore pressure shape function P Vector of excess pore pressure q Vector of nodal pore pressures r Array of residual errors T i Applied tractions u i Components of the displacement vector V i Components of the superficial velocity vector w Array of weighing functions 6 Vector of nodal displacements Kronecker delta e ,e i j Array of strain components a ,a i j Array of total stress components a', a i j Array of effective stress components - xv -ACKNOWLEDGEMENTS The work described in this thesis was made possible by an award from the National Sciences and Engineering Council Scholarship (NSERC), University of British Columbia Fellowship, and by a grant from Alberta Oil Sands Technology and Research Authority (AOSTRA). I am deeply indebted to my supervisor Professor P.M. Byrne who had the vision to suggest that research with respect to developing a numerical tool for solving reservoir engineering problems using a geomechanical approach would not only be academically interesting and challenging but would also have significant practical values. Throughout my research program he gave me a judicious blend of guidance, free rein, and encouragement which provided the ideal setting for completing this study. I am thankful to Professor D.L. Anderson for his interest and pain-staking task of checking the technical aspects and overall structure of this thesis. In my view, his sharp remarks have considerably improved the credibility of this report. At the outset of this research, several people were consulted to ensure that my efforts toward developing a finite element code would be steered in the right direction. Foremost, I would like to express my gratitude to Dr. A.M. Britto of the University of Cambridge for gener-ously furnishing me with a copy of their finite element program. This program provided the framework for accommodating the developments described herein. Also to be mentioned are Professor J.M. Duncan of Virginia Polytechnic Institute, Professor J. Booker and Dr. J. Carter of the University of Sydney, and in particular Dr. B. Simpson of Ove Arup and Partners (London). - xvi -Mrs. Kelly Lamb is responsible for the typing of this dissertation. Her professional attitude, precision and remarkable efficiency coupled with a very pleasant personality provided a l l the necessary ingredients to turn an ordeal into a happy event. I truly appreciate her efforts. This is an opportunity to express my deepest gratitude to my kind and loving parents for the care, encouragement and support that they have given me throughout my lif e . My brother Reza, who commenced his postgraduate course simultaneously with me, made my time at university far more reward-ing and enjoyable. Since mid-1985 when I left the university to begin my employment, Reza devoted a considerable amount of time and energy in expediting the finalization of this dissertation. His concern and enthusiasm will always be remembered. In dedicating this thesis to my family, Homa, Hose, Alex and Reza, i t is appropriate to include my fianc£, Rachel, who wisely agreed to keep away whilst my study was in progress. In my entire l i f e I have never lost sight of God's role in giving me the previlege to continue with my educational goals. The following motto closely reflects my sense of responsibility in attempting to return the fruits of his benevolence: "What we are is God's gift to us. What we become is our gift to God." - xvii -1. CHAPTER 1  GENERAL INTRODUCTION AND SCOPE Soils are particulate materials consisting of assemblies of mineral grains, with liquids and gases f i l l i n g the pore spaces. The mechanical properties of soils derive directly from the interactions of these phases with each other and with applied potentials (e.g., stress, hydraulic head, temperature). Soil mechanics is the study of the fluid conductivity, volume change, and strength behaviour of an assemblage of particles and is normally divided into three broad groups: dry soil, wet soil with no flow or constant flow, and wet soil with transient flow. The study of soil mechanics has very largely concentrated on the behaviour of fully saturated soils, that is, soils with only one pore fluid which may be either air or water. However, there are geographical and practical reasons for often using soil in a state in which the voids are fil l e d partly with gas and partly with liquid. Many natural soils contain large amounts of pressure dissolved gases in their pore fluids. For instance deep sea sediments may contain dissolved and free methane, the groundwater in geothermal areas may be charged with carbon dioxide, and petroleum reservoirs characteristically contain pore fluids in which complex gaseous mixtures are present in dissolved and free states. There are also some engineering projects which involve soils that cannot be considered to be in a totally saturated state. For example, in the construction of embankment dams or formation for roads or any other construction involving the placement and compaction of a f i l l material, i t would be reasonable to assume that the saturation state is not 100% throughout the soil. Another important, albeit not frequently encountered, 2. class of engineering problems which necessitate the consideration of soil in an unsaturated state are those of o i l sands subjected to unloading or reduction of their in-situ stresses. Such a process leads to evolution of large quantities of dissolved gases which will significantly lower the saturation level in the soil. In fact, whenever the confining pressure of any soil is reduced to below the liquid/gas saturation level of the fluid occupying its void spaces, the exsolution of some dissolved gases should be expected. The problems associated with depleting o i l sand reservoirs form an important area of investigation in this thesis. The significance of the problem is related to the vast volumes of known reserves which have been left relatively unexploited. Estimated in-place volume of heavy hydrocarbons in o i l sand deposits throughout the world approach 3000 billion barrels which is nearly equivalent to the total discovered conventional medium and light gravity reserves in-place in the world (Agar et al. (1983)). More than 90% of known hydrocarbon reserves occur in o i l sand deposits located in Alberta, and Eastern Venezuela (Demaison (1977)). Approximately 4% of Alberta o i l sand reserves are buried at depths less than 50 m which are economically recoverable by surface mining techniques. The remaining 96% is only exploitable by in-situ extraction procedures, generally requiring some form of heating as the extremely high viscosity of the crude bitumen in o i l sands makes conventional recovery by pumping impractical. One of the objectives in this study Is to develop a numerical procedure for simulating the processes involved in o i l recovery. Analysis can then be performed to determine the mechanisms governing the stability of the reservoir and the factors which result in the fluid production. 3. In the next section a brief introduction to the behaviour of o i l sands is given. Some basic definitions of the terms frequently used in this thesis such as soil skeleton, fluid phase, saturation, and liquid/gas saturation are presented in Appendix A. 1.1 Oil Sand Characteristics Oil sand may be considered as a three phase system comprising of a dense interlocked skeleton of predominantly quartz sand grains whose void spaces are occupied by bitumen, water, and gases. An illustration of o i l sand structure was produced by Dusseault (1977) and is shown here in Fig. 1.1. It can be observed that the bitumen occupies a large portion of the pore space, being separated from the solid grains by a thin film of water. Since o i l and water form continuous phases, gases (mostly methane and carbon dioxide) can only exist In the form of discrete bubbles (free gas) with some quantity dissolved in both the bitumen and water phases. In-situ o i l sand is a very dense, uncemented fine grained sand exhi-biting high shear strengths and dilatancy, and low compressibility charac-teristics compared to normal dense sand of similar mineralogy. The effec-tive permeability of o i l sand is very low because of the extremely high viscosity of the bitumen. This causes the soil to respond in an undrained manner to construction loading. Another unusual characteristic of o i l sand is its response during load removal. An unloading process that causes the level of confining stresses to decrease below the gas/liquid saturation pressure will result in gas exsolution. Due to the low permeability of o i l sands to gases, the evolution of gases are likely to occur under undrained conditions which will consequently induce a change in the volume of soil matrix. 4. Silty layer, largely oil free Micaceous, partings Matrix containing a high proportion of oil Intimate grain-to-grain contact Occluded gas bubble Surrounding water sheath, containing the majority of fines in the sand layers (after Dusseault, 1977) Fig. 1.1 Composition of In-Situ O i l Sand 5. The production and expansion of gas under undrained conditions causes the pore fluid to become more flexible, thus increasing its tendency to maintain higher pore pressures. For problems involving unloading of o i l sand, the relatively stiffer soil matrix will predominantly experience most of the stress change and responds with a corresponding reduction in effec-tive stresses until a tendency for developing negative stresses arises. In general since soils cannot sustain negative stress they will develop a substantial Increase in extensibility which results in stress changes being transferred to the now stiffer fluid phase. The physical consequences of such a process are a significant increase in volume and a marked reduction in shear strength. When the decrease in total stress occurs slowly over a period of time, the o i l sand may respond in a drained manner. In this case, the pore pres-sure change and resulting pore volume change will not disturb the soil fabric. It is worth pointing out that a finite time is required for the gas to completely evolve from the solution. Consequently o i l sand exhibits a short term undrained (i.e. before any gas exsolution) and long term undrained (i.e. after f u l l production of gases) behaviour. The period of time required for distinguishing these phases is very short relative to the time which takes for the excess pore pressures to dissipate. In dealing with long term problems, therefore, the assumption is made that gas evolu-tion occurs instantaneously thus ignoring the rate of gas buildup from the analysis. In general, however, such an assumption is considered to be inapplicable to cases concerned with the instantaneous response of gassy soils. Sobkowicz (1982), and Sobkowicz and Morgenstern (1984) have produced a comprehensive treatment of gas exsolution phenomenon, at both theoretical and experimental levels. 6. 1.2 Scope This thesis is concerned with the study of instantaneous, long term, and time dependent response of soils subjected to various boundary condi-tions such as pressure, displacement, temperature, and flow. Formulations have been developed to describe the mechanical behaviour of soil under these conditions which are subsequently implemented into a finite element program. Within the context of this work, soil is regarded as a multiphase system (i.e., composed of solids, liquid, and gas phases) with non-linear stress-strain and flow characteristics. Central to the whole approach adopted in this study is that the response of soil is governed by the principle of effective stress. Thus, one of the key objectives is to develop a suitable effective stress rela-tionship for partly saturated soils. This requires consideration of some of the laws governing the behaviour of gases and their interactions with solids and liquids. Discussions of these issues form the first segment of Chapter 2. In the remaining parts of this chapter a new scheme is presented for computing pore pressures under undrained conditions which is highly applicable to unsaturated soils and easily adaptable to the procedure normally followed in the course of finite element analysis. The basic requirement of this technique is a knowledge of the fluid compressi-bility whose derivation is also outlined in this chapter. The fundamental premise behind the proposed approach is that the gas remains in an occluded form within the liquid phase, thus the gas and liquid phases travel together as a single entity, herein referred to as the fluid phase. One of the important issues addressed in this thesis is the computa-tion of soil response upon exposure to temperature variations. This is of particular interest in problems dealing with the extraction of bitumen from o i l sands by the process of steam injection. Chapter 3 is exclusively dedicated to discussing these points and i t commences with a general intro-duction to the effects of temperature on o i l sands followed by development 7. of expressions for changes in pore pressure and effective stress as a consequence of thermal effects under different boundary conditions. The f i n i t e element technique for incorporating these formulations is subsequently outlined. Chapter 4 is concerned with the stress-strain model used. Since the present study is principally aimed at the behaviour of granular soils, a hyperbolic model is employed to describe its load-deformation response. For the particular application of the model to the problem of unloading a cavity in o i l sands, additional considerations must be given in ascertain-ing that the correct behaviour at failure is maintained. A crucial consideration is in ensuring that at failure, particularly in tension, the elements do not violate the yield criterion and redistribute the stresses in order to maintain a state on the failure envelope. Schemes for perform-ing such a task are also included in Chapter 4. The major thrust of this thesis is presented in Chapter 5 which deals with the time dependent or transient response of soils. The objective here is to develop suitable formulations for finite element analysis of consoli-dation in gassy soils. The general expressions derived take account of the nonlinear fluid and soil characteristics. The formulations developed couple the effects of both stress and flow. The numerical procedure adopted In this research is specifically designed to achieve the most efficient method of performing accurate finite element analysis. The aspects pertaining to the development of such tech-nique are covered in Chapter 6. It is shown in this chapter that the main source of inaccuracies with the solution of finite element analyses stem from the assumed variation of the unknowns within the elements and the method employed in solving the equations. Schemes are then presented which 8. eliminate some of these errors and greatly enhance the efficiency of the numerical algorithm, with particular emphasis on techniques for minimizing the computational effort associated with forming and solving the governing non-linear equations. The intentions in Chapter 7 are twofold: to verify the validity of the finite element computer program developed and to demonstrate its applica-tion in solving some difficult geotechnical problems. The former objective is satisfied by comparing the computed results with the available closed form solutions for some particular features of the program such as consoli-dation and temperature analyses. The latter is dealt with by applying the program to analyse the soil response around deep wellbores in o i l sands. By assuming reasonably realistic boundary conditions, several problems are analyzed in an effort to assess the stability, and to identify the factors which result in the flow of fluid from these materials. The concluding remarks are made in Chapter 8 along with some comments on the aspects which warrant further investigation. 9. CHAPTER 2 EFFECTIVE STRESS ANALYSIS OF GASSY SOILS 2.1 Introduction The most fundamental concept in soil mechanics is that the response of soils are governed by the effective stresses. This can be formulated as (Terzaghi, 1923): o» = o - u f (2.1) where a and a' are the total and effective stress components and u^ is the pore fluid pressure. In the case of fully saturated soils, the pore fluid pressure is that of the liquid which occupies the entire void space. In the case of unsaturated soils, however, the pore fluid pressure is a function of both the liquid pressure (u ) and the gas pressure (u ). In order to establish *• g a unified approach for determining pore pressure, regardless of its satura-tion state, the assumption can be made that the mixture of gas and liquid components can be replaced with a fluid phase having an equivalent compressibility. Under isothermal conditions, i t can then be stated that: A U f = B f ( A e v ) f (2.2) where is the bulk modulus of the f l u i d phase, which in the case of a saturated s o i l i s that of the l i q u i d , and (Ae^)^ is the change in volumetric strain experienced by the fluid. For unsaturated soils, there-fore, a knowledge of the fluid compressibility is required for computing pore pressures. The considerations which must be given in deriving a suit-10. able expression for the fluid modulus is described in this chapter. Such an expression is then used to compute pore pressure under undrained conditions by making use of the fact that the fluid volumetric strain can be directly obtained from a knowledge of the volumetric strain developed in the soil. A suitable scheme has been proposed for this purpose which is akin to the procedure followed in the finite element analysis. The method for computing pore pressure under transient conditions is covered in Chapter 5. Prior to proceeding with these issues, a review is presented of the attempts made in describing and formulating the factors which must be considered in the effective analysis of unsaturated soils. 2.2 Literature Review The significance of the principle of effective stress in soil mechanics has sometimes led to its consideration as a physical law and attempts have been made to examine the principle theoretially. These examinations of the principle of effective stress, and of the physical meaning of effective stress, encounter difficulties when forces at the interparticle contacts have to be expressed in the form of a stress (Atkinson and Bransby (1978)). In fact Terzaghi (1936), does not assign any physical meaning to the effective stress other than 'it has its seat exclusively in the solid phase of the soil' and, second, that he limits the principle to 'all measureable effects'. Indeed, Terzaghi arrived at the principle of effective stress from the results of laboratory experiments on saturated soils under moderate stress levels and i t is evident that he viewed the principle as a working hypothesis that was sufficiently accurate for practical engineering purposes. Detailed theoretical and experimental 11. examinations of the principle of effective stress have been published by Skempton (1960) and by Bishop et al. (1965) and no conclusive evidence has yet been found which invalidates Terzaghi's original postulate for saturated soils at normal levels of pressure. For an unsaturated soil, the pore fluid consists of gas and liquid. The effective stress equation then involves liquid pressure (u^)> Sas pressure (u ), and surface tension between gas and liquid. For a material with a two phase pore fluid where both gas and liquid are present having separate pore pressures u and u respectively, Bishop (1959) suggested g * that the pressure in the liquid phase would act over a reduced area so that the effective stress would become: t °' = ° " u g " X (u £ - u g) (2.3) Others have looked at more restrictive cases where, for example, the pressure in the pore gas was assumed to be constant at atmospheric (the gas voids being assumed continuous and open to the atmosphere) so that u =0; S and have used different symbols for different terms (e.g., Aitchison (1960), Jennings (I960)). There appears, however, to have been some general agreement on the form of this effective stress equation (Aitchison and Bishop (I960)) in the aftermath of the 1960 conference on pore pressure and suction in soils. But criticism has also been levelled at i t from two directions: on the one hand because of its supposed incompleteness, and on the other from suggestions of its irrelevance to the understanding of the behaviour of partly saturated soils. For instance, Lambe (1960a) suggests that the parameter x includes structural effects which may greatly limit the use of measured x values: x I s certainly not a given soil property 12. dependent only on degree of saturation. Lambe (1960b) considers the forces acting across a potential shear surface and implies an effective stress relationship in the following form: a' = a - u a - u„aA - R + A (2.4) g g I i where R and A are repulsive and attractive pressures acting over the plane and a and a denote the areas of a gas and liquid intersected per unit g i area of shear plane. The measurement of the parameters required in Lambe's expression severely limits its practical application. Donald (1963) introduces the surface tension into the x parameter defining: f(S ) x - S + - — (2.5) * r u - u„ v ' g * where is the degree of saturation of the pores, and f(S f) incorporates the surface tension effects. Sparks (1963) also claims that Equation (2.3) is incomplete due to the omission of surface tension effects. By considering a typical section through an element of soil that cuts through interfaces between gas and liquid on which surface tension forces act and by summing components of interparticles forces (including repulsive forces such as electrical and attractive forces such as van der Waal's between particles), pore pressure forces and surface tension effects he suggests: a' = a - X.u - A,u + AT 1 g 2 i 3 s (2.6) 13. where T is the surface tension at the gas liquid interface. If the pores s are fil l e d with liquid Aj = Xg = 0, and i f filled with gas * 2 = * 3 = 0. Sparks employed a model composed of uniform spheres arranged in various packings, and having various contact angles between the liquid and the soil particles. On the basis of these studies, he concluded that x ^ y assume values greater than one or less than zero. Sparks also made up an extensive set of charts to be used to determine Aj, A2 and Ag for different packing conditions and contact angles. For the special case of pores containing occluded gas bubbles (i.e., 0.8 < S < 1.0) Sparks recommends the following relationship: a' = a - (1 - A)u £ (2.7) where A is the effective soil particle contact area per unit area. For most soils i t is reasonable to assume that the area of contact between the solid particles is small compared to the area of contact between the solid particles and the liquid. Thus by setting A=0, Equation (2.7) reduces to the familiar form of: a» = o - nl (2.8) Therefore according to Sparks' studies, the same effective stress equation which applies to saturated soils can also be used for unsaturated soils where the gas in the voids is in the form of occluded bubbles. By considering the fluid phase to be composed of air and water, Fredlund (1973) argues that partly saturated soils should be regarded as four phase materials with the air-water interface (referred to as 14. 'contractile skin') being considered as an independent phase. Fredlund and Morgenstern (1977) and Fredlund (1979) cite references for evidence that the contractile skin has differing properties from those of the continuous materials; and has definite boundary surfaces. Appropriate stress state parameters for the four phase unsaturated soil are derived by supposing that the soil can be considered as a mixture of two phases that come to equilibrium under applied stresses (the soil particles and the contractile skin) with two phases that flow under applied pressures (the air and the water). By conceptually re-arranging a soil element into a portion containing the contractile skin and air, Fredlund and Morgenstern (1977) conclude that the general overall equilibrium equation contains three independent sets of surface tractions: i) (a - u w) i i ) (u - u ) w a i i i ) u w where u and u are the pressures in the air phase and water phase respect-a w , ively. It is further stated that u can be eliminated when the soil particles w are assumed incompressible. A slightly different analysis produces a different pair of combinations of a, u , u as the pair of independent a w stress parameters and, in fact, Fredlund (1979) concludes that the pair (o - u ) and (u - u ) is the most useful pair for application to earth a a w pressure problems. Fredlund (1979) has presented an appropriate conceptual basis for study of unsaturated soils, however, the number of parameters required in 15. describing the behaviour of partly saturated soils makes i t unreasonable to expect that a simple treatment will apply widely. The brief review of some of the past work presented in this section serves to indicate that i t is not possible to derive a unique relationship between stresses and pore pressures for unsaturated soils which can be broadly applied to a l l soils and under a l l conditions. The successful application of any formulation requires certain idealizations (e.g., as suggested by Sparks (1963)) or some degree of empiricism (e.g., as suggested by Bishop (1959)). 2.3 Fundamental Laws Describing the Behaviour of Liquid/Gas Mixture 2.3.1 Gas Laws If the gas and liquid are immiscible, the physical laws governing compression of gas will determine the fluid compressibility, whereas i f the gas is soluble to some extent in the liquid, the fluid compressibility will also be influenced by the solubility relationship. The basic laws govern-ing the volume and pressure relationships are Boyle's law and Henry's law. Boyle's law states that (Sisler et al. (1953)) under constant tempera-ture conditions, the density of the gas is proportional to the absolute pressure, i.e.: P e = w R T (2.9) g g where P is the absolute pressure T is the absolute temperature R is the universal gas constant e is the volume of gas 8 w is the weight of gas 16. For the undrained case, where the weight of gas does not vary, Eq. (2.9) can be written as: P e = k (2.10) g where k is a constant. Henry's law states that (Sisler et al. (1953)) the weight of gas dissolved in a fixed quantity of a liquid, at constant temperature, is directly proportional to the absolute pressure of the gas above the solution, i.e.: "dg P° = — (2.11) wj P1 dg where w^  is the weight of gas dissoyled at pressure P and superscripts 0 and 1 denote i n i t i a l and final conditions respectively. Henry's law implies that the volume of dissolved gas in a fixed quantity of liquid subjected to constant temperature and to a confining pressure P is constant when the volume is measured at P. Thus e, = He (2.12) dg £ ' where H is the Henry's constant which is temperature dependent and over a wide range of pressure is also pressure dependent, e^ is the volume of liquid. Since the volume of dissolved gas remains constant, i t is also possible to combine the free and dissolved gas and apply Boyle's law to the entire volume (Fredlund (1976)), i.e., (e, + e° ) P° « (e. + e* ) P 1 (2.13) dg fg' dg fg' v ' where e ^ denotes the free portion of the gas phase. 17. 2.3.2 Rate of Gas Diffusion The amount of gas going in and out of solution is time dependent and is either ignored or taken into account, depending upon the engineering problem being considered. The rate at which gas goes into solution is described by Fick's law of diffusion. The driving force is a concentration or density difference between the free gas and the gas entrapped in the liquid (Fredlund, 1976). The time 't' necessary to dissolve the gas for a unit change in pressure can be determined from the following expression (Beck, 1963): r 2 t = — (2.14) where D^  is the diffusion coefficient and r is the radius of the bubble. The value of D^  for most gases is in the order of 10~5 cm2/sec hence for small gas bubbles, time for diffusion is quite fast and for many practical situations may be assumed to occur instantaneously. Sobkowicz (1982) in studying the behaviour of gassy soils, has shown that a very small, but finite, time is required for complete exsolution of gases. For practical purposes, particularly in comparison with the time necessary for the consolidation process; i t is reasonable to assume that evolution of gases occurs instantaneously. Such an assumption has been employed in the present study. 2.3.3 Surface Tension and Capillary Effects Liquids have cohesion and adhesion, both of which are forms of molecu-lar attraction. Cohesion enables a liquid to resist tensile stress, while adhesion enables i t to adhere to another body. As stated by Daugherty and 18. Franzini (1965), the attraction between molecules forms an imaginary film capable of resisting tension at the interface between two immiscible liquids or at the interface between a liquid and a gas. The liquid property that creates this capability is known as surface tension. Capillarity is due to both cohesion and adhesion. When the former is of less effect than the latter, the liquid will wet a solid surface with which It is in contact and rise at the point of contact; i f cohesion predomin-ates, the liquid surface will be depressed at the point of contact. For an unsaturated soil, the pore fluid consists of gas and liquid. The effective stress equation involves liquid pressure (u ), gas pressure JC (u ), and surface tension (T ) between gas and liquid. As a result of 8 s surface tension, the gas pressure and the water pressure of a liquid-gas mixture are not the same. Kelvin's equation for a single capillary tube, relates the difference in the gas and liquid pressure to the surface tension and effective radius of the gas-liquid interphase (r). Schuurman (1966) has used such a relationship to develop an expression for fluid compressibility of unsaturated soils containing occluded gas bubbles (i.e., 0.8 < S < 1.0). As a result of surface tension, the gas and liquid pressures in a system differ by an amount equal to the capillary pressure (u c). Their relationship can be expressed as u„ = u - u (2.15) A g e In a mixture of gas and liquid, the bubbles are not usually the same size and shape. Denoting the average capillary radius in the mixture by r , then: 19. T u c r (2.16) Schuurman (1966) took surface tension into account by assuming that the air bubbles were a l l of equal size and spherical in shape. By this assumption, the value of r in Eq. (2.16) would be equal to the radius of c the air bubbles. Taking the i n i t i a l radius of air bubbles as r Q , the radius after a change in void ratio can be expressed as: The accuracy of this expression depends on how closely the idealiza-tion of uniform spherical bubbles and the assumption that the total number of gas bubbles do not change during the compression of the soil reflect the reality. In fact, the situation is more complicated than is shown by the simple relationship of Eq. (2.16). This is because capillary pressure does not necessarily arise only from the bubbles. It can also develop from the other liquid surfaces exposed to gas at the specimen boundaries. It is, therefore, more logical to use an empirically determined average capillary radius to simulate capillary effect as suggested by Mitchell et al. (1965). Their study has shown that reasonable agreement with test results can be obtained i f the following expression for the average capillary radius is used: (2.17) S -r c r cs (1 - ) (2.18) 20. where r is the capillary radius at saturation, and cs is the lowest degree of saturation at which the liquid begins to flow freely. may be determined from permeability tests for samples with different degrees of saturation. 2.3.4 Compressibility of a Miscible Gas/Liquid Mixture By applying Boyle's law to the total volume of gas in the system: fg dg a g fg dg a g where P is the atmospheric pressure and subscripts 'fg' and 'dg' denote a free and dissolved components of the gas. u° and u 1 are the i n i t i a l and 8 8 final gas gauge pressures respectively. Since u = u - u (2.15) bis I g c and 2T u = — - (2.16) bis c r c From Eq. (2.19) i t follows that: u l = — — ± £ — — 2 £ — i - - — - (2.20) (el +e, ) rc fg dg 21. According to Sparks (1963), for the case where the gases exist as occluded bubbles, i t can be reasonably assumed that fluid pressures and liquid pressures become identical. Following a similar procedure as outlined in Appendix B, the fluid compressibility can be expressed as: °f " " e dP or 1 de de c f = - 7 ( - d f + d P i ) <2'21> where P is the absolute pore fluid pressure. Since (P = u_ + P ) and (u = u ); then Eq. (2.21) can be written as: r 3. r Jc 1 d e f e d e£ cv = - - ( - j - ^ - + -j-*-) (2.22) f e du du„ Differentiation of Eq. (2.20) yields the following expression: duj (e° + e j )P + (e2 + e, )u° 2T dr _ J _ = £& dg a fS dg g + _ J . . _ £ _ (2.23) del ( e i +e, ) 2 r 2 de] fg fg dg' c fg By using the following relationships: e £ = e(l - S) and e, = SeH fg dg Eq. (2.23) can be written in the following form: dul ( e n - Sften + HS.en)(u° + P ) 2T dr S, —L = - l i °-l ° ° g 2 l + _ ± . _ £ . J . (2.24) d er g < el " S l e l + H S l e l > 2 r c 2 d S l e l 22. where e Q and are the i n i t i a l and final void ratios SQ and Sj are the i n i t i a l and final degrees of saturation As discussed in Section 2.3.3, i t can be assumed that: Therefore: S - S r c = r cs ( T^s7 ) dr r Hsf = T ^ S - <2'25> Substituting Eq. (2.25) into Eq. (2.24) results in: d u I ( e 0 " S 0 e 0 + H S 0 e 0 ) ( u g + V 2 T s rcs S l — — = — — B — + — - > — ^ — . — (2.26) de1 ( e - S e + H S e ) 2 r 2 (1-S ) e fg v 1 1 1 I K c f 1 Substituting Eq. (2.26) into Eq. (2.22) and using the definition that: du„ - u C = - S e C then the following expression for the fluid compressibility can be obtained: s i e i c 1 l i e c i + * (2.27) f e 0(e 0-S 0e 0«S 0e 0)(«J« ) 2Tg r Sj e Q & . e 0 (ej-Sjej-WSjej)2 r 2 (1-S ) ^ 23. For most practical applications this expression can be simplified to: m (1 - S + SH) ( 2 > 2 g ) (u° + Pfl) This is because generally the liquid phase has a very small compressi-b i l i t y ( i . e . , C * 0.0), surface tension effects are small (T « 0.0) and K, S for moderate levels of load application Sj " SQ and a e Q. The expression for fluid compressibility as shown by Eq. (2.28) is the same as that presented by Skempton and Bishop (1954) and Bishop and Henkel (1962). Equation (2.27) represents a rigorous expression for the compressi-bil i t y of the fluid whereby provisions are made for liquid compressibility, and gas compressibility as a function of dissolved gases, and surface tension effects. In the normal course of analysis, i t is expected that the stress changes are imposed in small enough increments such that e 1 and e Q, and similarly Sj and SQ are approximately equal. Eq. (2.27) can then be reduced to: C - + SC, (2.29) * (u° + P ) r * _S §_ _ 2 T -(1-S+HS) s 2 (1-Sf) r 2.4 Computation of Pore Pressures Under Undrained Conditions In soil mechanics considerable attention has been directed towards the prediction and measurement of the pore pressure set up by the changes in the state of stress under undrained conditions. The conventional approach 24. has been to compute changes in total stress using a total stress-strain law and then use some empirical equation to compute the pore pressure change. The most widely used equation is that proposed by Skempton (1954): Au = B[Aa3 + A(Aa - Aag)] (2.30) where A and B are pore pressure parameters derived from triaxial test results. The B parameter is a function of the relative stiffness of the fluid phase as compared to the skeleton phase and varies between zero and unity. The A parameter represents the shear volume coupling or dilatant characteristics of the skeleton. Christian (1968) adopted an effective stress approach and proceeded with the calculation of pore pressure subject to the constraint of no change in volume. In the technique he proposed, the conventional stiffness matrix is modified by the addition of a pore pressure term large enough to prevent volume change together with a new equation expressing zero volumetric strain. Naylor (1973) proposed a rational approach for determining pore pressures under undrained conditions. The method follows the principle of effective stress and is equally applicable to saturated or unsaturated soils irrespective of any assumed stress-strain behaviour. Particularly important, the proposed scheme provides a numerically stable procedure for the analysis of problems in which the effective stresses go to zero and a l l the load has to be carried by the pore fluid. The author has extended Naylor's formulation to allow for the gassy nature of the pore fluid (expressed in pore fluid compressibility). As is shown later, the proposed technique allows a unified approach to drained and undrained analysis (i.e. both types of analysis can be performed using the same data base) and is highly adaptable to the procedure followed in finite element computations. The principal idea behind this approach is that under undrained conditions there is no movement of pore fluid and solids relative to the skeleton, consequently the skeleton and its comprising components undergo the same deformations and their volumetric strains can be equated. Such an equation can then be used to compute the pore pressures (i.e. Eq. (2.2)), provided that the bulk moduli of a l l components which make up the soil system can be replaced with a phase having similar compressibility characteristics. Such a phase, herein referred to as the 'Homogenized Compressible Phase', can be loosely thought of as a pore fluid which not only occupies the voids but also extends to account for the contribution of space taken up by the solid grains. With this terminology, the three phase soil matrix can be reduced to the conventional two-phase system comprising of a saturated material with a homogenized compressible phase for which the principle of effective stresses can be applied in the following form: Aa' = Aa - Au (2.31) cp where the subscript cp denotes the hypothetical compressible phase. Now i f the compressible phase undergoes a volumetric strain (Ae ) due to the r «- o \ v cp application of stresses under undrained conditions, then: Au = B (Ae ) (2.32) cp cp v cp where B denotes the bulk modulus of the homogenized compressible phase. According to the foregoing discussion, i f C (=1/B ) is used to cp cp define the compressibility of this homogenized material, the following relationship can be stated: C = n Cc + (l-n)C (2.33) cp r s where n is the porosity and is the compressibility of the pore fluid which includes the effects of gas and liquid compressibilities and also the surface tension as expressed by Eq. (2.27). Equation (2.33) in terms of bulk moduli becomes: cp f s 2.4.1 Analytical Rationale Based on the effective stress principle, the pore fluid and soil skeleton contributions to stress can be separated and written as: {Ao} = {Ao'} + {m}Au (2.35) cp where {a} and {a'} are vectors representing the total and effective stress states and {m} is a vector indicating that only the direct stresses parti-cipate in the effective stress law. The components of these vectors are: T {a} = { o a a T x T } x y z xy yz zx {m}T = { 1 1 1 0 0 0 } The incremental effective stresses are related to the incremental strains by the relationship: {Aa1} = [D']{Ae} (2.36) where {e} is a vector representing the strain components of the soil skeleton i.e.: T {e} = {e e e Y y y } x y z xy yz zx and [D'] consists of parameters which represent the effective stress-strain behaviour of the material and can take any form, e.g. elastic or elasto-plastic. Similarly, total stress and strain changes are related by: {Aa} = [D]{Ae} (2.37) Now i f the element of soil undergoes volumetric strains {Ae } . as a v sk consequence of the imposed stresses {Aa}, i t can be stated that under undrained conditions these strains are equal to those developed in the homogenized compressible phase, thus: {Ae } = {Ae } . (2.38) v cp v sk From definition of volumetric strain: {Ae } . = Ae + Ae + Ae (2.39) v sk x y z Alternatively , the above can be expressed as: 2 8 . ( A O v sk {n}T{AE} (2.40) Hence v cp {m}T{Ae} (2.41) The conversion of Eq. (2.39) to the form presented by Eq. (2.40) i s the paramount feature of the present technique. This step allows the same l inking vector ( I . e . , {m}) that was used to Indicate the normal components of the pore pressure vector ( i . e . , {u}, as shown by Eq. (2.35)) to be used to serve the same function for the strain components. Substitution of Equations (2.36), (2.37) and (2.42) into the def in i -t ion of effective stress, i . e . Eq. (2.35), introduces {Ac} into each of the three terms. Since the components of { A e } are independent of each other, they can be cancelled thus leading to the following expression: The above expression represents the governing relationship between the to ta l and effective stress stiffness matrices as a function of the equiva-lent bulk f lu id modulus or the compressible phase, i f the so l id partic les are considered to be compressible. Substitution of Eq. (2.41) into Eq. (2.32) yields: (2.42) [D] - [D'] + {m}{m}x B cp (2.43) It Is apparent from Eq. (2.43) that the total and effective stress stiffness terms are explicitly related by the compressible phase modulus. From the analysis point of view such a relationship provides a convenient way of converting between total and effective stress parameters. This means that both undrained and drained analyses can be carried out from the same data base (in terms of effective stresses) by simply specifying appropriate values for B . For instance in performing undrained analysis cp a value of B^  equal to the bulk modulus of the fluid would be used, but by simply setting B^  equal to zero a drained analysis can be simulated. It is also worth emphasizing that in Eq. (2.43) no restriction is placed on the possible forms which the [D'j or B can take. The same cp formulation can be applied to any material whether elastic, elasto-plastic, dilative, saturated, or unsaturated, etc. This is because the pore pressures are determined explicitly from the computed deformations, regard-less of how these movements are generated. The inclusion of the developed technique in finite element analysis requires no extra Information over the data which is already provided to perform the numerical procedure. In fact no additional computational effort is demanded other than using the strains which are normally worked out in the course of finite element analysis. Eq. (2.43) is used in the finite element formulation so that the pore pressure is not an unknown in the solution of the system of equations but is obtained later from Eq. (2.42). Implementation of this technique into finite element analysis is outlined in Appendix C. 2.5 Summary Soil is Inherently a multiphase system consisting of a mineral phase plus a fluid phase. The nature of pore fluid influences the magnitude of shear resistance in the soil and as such, its correct representation 30. becomes imperative in any rigorous study of the soil behaviour. The presumption that soil is composed of an incompressible fluid phase and incompressible solid particles neglects the presence of some of the physical parameters which affect the mechanical characteristics of the material. Even soils that are initially fully saturated with water will not always remain so. For instance, Okumura (1977) reported on the disturbance of samples taken just below the sea floor but subjected to i n i t i a l pore pressures in excess of 100 m of overlying water. These marine sediments were of low permeability and saturated with water and dissolved air. Although the solubility of air in water is small, he observed that the large decreases in total stress on retrieval of the sample led to substantial decreases in effective stress. The earlier sections in this chapter discuss the physical and theo-retical concepts which are pertinent to the analysis of soil as a multi-phase material. A review of the previous work on this topic has revealed that a practical method of modelling the soil characteristics entails some degree of Idealization and empiricism. The facets which are considered for this purpose are the constitutive laws related to the behaviour of gases, compressibility characteristics of the various components forming the soil structure and their interactions, and surface tension effects. By making the assumption that the gases remain occluded, i t is hypothesized that a multiphase soil system can be phenomenonlogically treated as a saturated material with an equivalent fluid phase whose compressibility is composed of the relevant contributions of each phase (i.e., solids, liquids, and gases). One of the applications of the formulation derived for the fluid modu-lus is in computing pore pressures under undrained conditions for soils 31. w i t h any c o m p r e s s i b i l i t y c h a r a c t e r i s t i c s o r a t any s t a t e o f s a t u r a t i o n . T h i s c a n be c o n v e n i e n t l y c a r r i e d o u t by r e l a t i n g t h e p o r e p r e s s u r e s d e v e l -oped t o t h e v o l u m e t r i c s t r a i n s and t h e e q u i v a l e n t s t i f f n e s s o f t h e f l u i d . S i n c e t h e magnitude o f v o l u m e t r i c s t r a i n s a r e g o v e r n e d by t h e chosen s t r e s s - s t r a i n r e l a t i o n s h i p , t h e method r e s u l t s i n changes i n p o r e p r e s s u r e w h i c h a r e c o n s i s t e n t w i t h t h e g e n e r a l a s s u m p t i o n s employed t h r o u g h o u t t h e a n a l y s i s . One o f t h e most a t t r a c t i v e f e a t u r e s o f t h e p r o p o s e d t e c h n i q u e i s t h a t i t f u r n i s h e s t h e r e s u l t s w i t h o u t any e x t r a c o m p u t a t i o n a l e f f o r t o r a need f o r a d d i t i o n a l s o i l p a r a m e t e r s . 32. CHAPTER 3 RESPONSE OF UNSATURATED SOILS SUBJECTED TO THERMAL EFFECTS 3.1 Introduction Temperature variations have a significant influence on the behaviour of soils. In general temperature affects pore pressures and interparticle forces, and induces changes in volume. Temperature can also significantly alter some of the engineering properties of the soil, for instance i t may change the hydraulic conductivity or result in gas exsolution which will increase the compressibility. Some of these effects have been considered by Campanella and Mitchell (1968), Lambe (1961), Henkel and Sowa (1963), Scott (1963), Mitchell and Campanella (1963), and Mitchell et al. (1968). A knowledge of changes in volume, pore pressure, and stress as a result of variation in temperature is of practical interest in many geo-technical activities such as the engineering projects involving power plants, underground power cables, liquefied natural gas storage and pipe-lines, and also in conducting laboratory tests on soil samples. Amongst the more common field problems, the extraction of bitumen from o i l sands by the process of heat injection is probably the most important example of a case where an understanding of soil's response to thermal effects is greatly needed. Although o i l sands have been found in abundance, they have hitherto been relatively unexploited as the extremely high viscosity of the crude bitumen makes conventional recovery by pumping impractical. In-situ extraction methods generally involve heating the o i l sand with pressurized steam or by In-situ combustion. Such a process decreases the viscosity of the bitumen and thus enhances the rate of flow. But the temperature 33. Increase also tends to increase the pore fluid pressure which, under undrained conditions, result in a reduction In effective stresses and hence a loss in shear strength. This In turn will cause some shear deformations that can adversely affect well casings, shafts and tunnels. Moreover, volume changes in the affected strata will occur as a consequence of temperature increase. With the propagation of heat through the o i l sand, the thermal expansion will produce ground movements towards underground utili t i e s or the surface. Thermal expansion may further cause differential displacements which could result in the propagation of fractures, threaten-ing the stability of well casings and shafts (Scott and Kosar (1982), Agar et al. (1983)). The severity of the problems can be further compounded i f the change in temperature causes the dissolved gases to come out of solution, thus increasing the expansion of the sand matrix. The determination of the thermal expansion of the soil under investi-gation, the pore fluid pressure generation, and the fluid flow is therefore of paramount importance to the f u l l understanding of the processes taking place during in-situ heat simulation and bitumen withdrawal. In this chapter formulation for the magnitude of the excess pore pressure developed due to temperature change for the case of a three phase soil system is derived. For consistency with the general approach adopted in this thesis, the solid particles and the liquid component of the sand matrix are considered to have finite compressibilities. Presented also, are the theoretical expressions for the volume change due to temperature variation under both drained and undrained conditions. The finite element scheme used for computing changes in stress and pore pressure as a result of thermal effects is included in Appendix E. 34. Since i t is expected that the greatest application of these formula-tions will be for the problems involving heat induction into o i l sands, the fir s t section is devoted to describing the general response of o i l sand to thermal effects. 3.2. Effect of Temperature on Oil Sand An in-depth examination of the effect of temperature on oil sand is beyond the scope of the present study, however, i t is beneficial to elabor-ate on some of the engineering parameters which are more significantly influenced by temperature variations. Oil sand subjected to increases in temperature will experience volume increases, the magnitude of which are dependent on the amount of pore fluid drainage which is permitted. Assuming that the o i l sand under in-situ conditions is fully saturated, the application of a rapid increase in temperature will result in the development of a significant positive pore pressure due to the greater volumetric expansion of the pore fluid relative to the mineral solids. Provided that the boundary conditions are such that a hydraulic gradient is set up, then the fluid will drain from the soil at a rate governed by its hydraulic conductivity. Based on physical reason-ings, drained thermal expansion can be considered to represent a lower bound to thermal expansion behaviour while the upper bound occurs under undrained conditions. The change in volume is composed of the expansion of solid particles and soil skeleton together with the pore fluid. The latter plays the pre-dominant role under undrained conditions, being responsible for over 90 percent of the overall volume change even before any gas exsolution. Should gases evolve from solution, additional large volume changes will occur irrespective of the drainage conditions. 35. Fig. 3.1 taken from Scott and Kosar (1982) illustrates the appreciable difference in volumetric expansion between drained and undrained cases for a typical sample of oi l sand. Another component of volume change which develops in any but the densest arrangement of the soil particles, can be attributed to the reorientation of the sand grains during the change in temperature. This occurs because an increase in temperature causes a decrease in the fric-tional resistance and the shearing strength of individual interparticle contacts. Consequently, there is a partial collapse of the soil structure and a decrease in void ratio until a sufficient number of additional bonds are formed to enable the soil to carry the same effective stress at the higher temperature. This phenomenon was originally recognized by Campanella and Mitchell (1968) who further quantified its magnitude in terms of change in temperature and a parameter o known as the coefficient S u of structural volume change of s o i l structure due to change in temperature. Findings from extensive studies conducted at the University of Alberta (e.g. Scott and Kosar (1982,1985)) into the behaviour of o i l sands subjec-ted to temperature variations reveals that a is not a constant soil s t property but varies with the state of void ratio and temperature. The reason for such a behaviour is that once the o i l sand matrix undergoes some heat induced compaction i t will achieve a tighter packing particle arrange-ment, leaving very l i t t l e room for further grain slippage. Since the greatest structural reorientation occurs in soils having their loosest possible state, one can deduce that such an effect is not a dominant feature of o i l sand's characteristics which are generally very tightly TEHPERRTURE (DEGREES CELSIUS) Fig. 3.1 Thermal Expansion of Oil Sand Samples. (After Scott and Kosar, 1982) 37. packed and under high confining stresses. Figs. 3.2(a) and 3.2(b), repro-duced from Scott and Kosar (1982), clearly demonstrates the volumetric response of o i l sand and the corresponding variation of a due to S t increases in temperature. Fig. 3.2(a) represents the. change in volume under drained conditions for samples under different confining pressures and back pressures, and having different i n i t i a l densities. When the responses are compared with that of the solid particles alone, i t can be observed that the densest sample exhibits the smallest amount of structural rearrangement. Unfortunately, due to Insufficient back pressure in the case of the dense sample (I.e., 2000 kPa), gases commenced to exsolve at temperatures over 180°C which has somewhat complicated the general appear-ance of this diagram. The main observation from Fig. 3.2(a) is that beyond a temperature of approximately 120°C no further reorientation of soil particles occur as the expansions follow parallel lines with that of the solids. Fig. 3.2(b) isolates the contribution of volume change due to the grain slippage for the tests shown on Fig. 3.2(a) (i.e., the net difference between the observed change in volume and that of the solids alone). It can be visual-ized that the temperature induced changes in the soil structure Is greatest at approximately 50°C, changing sign beyond this point for the dense samples (i.e., becoming expansive) and reducing in rate for the sample with the lowest density. These tests provide a clear understanding of the physical meaning of a and show that this parameter is not a unique soil property as is the coefficient of material expansion for example. Indeed a is highly dependent on the soil type. For instance the tests carried out on the i l l i t e clay by Campanella and Mitchell (1968) illustrate that in contrast 38. to 0 2S 50 75 100 125 150 175 200 TEMPERATURE (DECREES CELSIUS) Fig. 3.2(a) Thermal Expansion of Sand Grains and Soil Struture in Drained Oil Sand Tests (After Scott and Kosar, 1982) E 1 1 1 1 1 1 ' 1 ' 1 ' 1 ' 1 ' v w - H u c d In Mini mil ctructire VV \ bulk demlty - 2.Q na/cu • \ \ \ cry demlty - ISO rto/ai. «^— v—aol l atrueturt for drained Icrx bulk denalty - ZC21 r V c u . • dYy density -1721 Mc/cu. • mil itructure for drained te»t - bulk deraity • 1968 no/aj . a dry denalty - IS2I lig/cu. • -mil m n m far &~tlrcd test ' -sulk d m l t y - 1390 rvai. • dry denalty - 1840 r V c u . • I . I . I i . i . i i 0 25 50 75 100 125 150 175 200 lE I l tHHIURE DECREES m T i l IS) Fig. 3.2(b) Thermal Volume Change of Soil Structure (After Scott and Kosar, 1982) . 39. to o i l sand, the particle rearrangement tends to expand the soil structure up to temperatures of about 50°C and then to cause contraction until a temperature of approximately 120°C, before its effect is completely diminished. Figs. 3.3(a) and 3.3(b) depict these observations. Another property of o i l sand which is significantly affected by temperature is its hydraulic conductivity. As mentioned in the introduc-tory part of this chapter, the bitumen in o i l sands is essentially immobile at in-situ temperatures, however, its viscosity is significantly lowered at elevated temperatures which results in an increase in soil permeability. This aspect of o i l sands has also been studied at the University of Alberta and Fig. 3.4 reproduced from Scott and Kosar (1984) shows the variation of bitumen viscosity and o i l sand permeability with temperature. Also shown on this figure is the variation of water viscosity and the permeability of clay studied by Campanella and Mitchell (1968). The comparison elucidates the unusual response of o i l sands to temperature which is not manifested by non-oil-comprising soils. 3.3 Prediction of Pore Pressure and Volume Changes Due to Temperature  Variations 3.3.1 Literature Review A theoretical analysis for predicting volume changes in saturated soil subjected to temperature variations and pore pressure changes due to undrained heating has been proposed by Campanella and Mitchell (1968) and examined further by Mitchell (1976). The same formulations have been used by Agar et al. (1983) and by Scott and Kosar (1982, 1985) as part of their studies on the response of o i l sands. The formulation proposed by Campanella and Mitchell is valid for pressures above the gas/liquid saturation pressure (i.e., the pore fluid pressure sufficient to prevent A O . 8>S 0) c CO u 0) 6 3 i — I O > 0.5 0.4-0.2 0.0 clay p a r t i c l e s S o i l structure clay p a r t i c l e s & s o i l structure 50 100 Temperature °C 150 F i g . 3.3(a) Thermal Volume Change of Clay Particles and So i l Structure (After Scott and Kosar, 1985) 0.2 . 0.1-Temperature °C F i g . 3.3(b) Thermal Volume Change of Clay Soi l Structure (After Scott and Kosar, 1985) TEMPERATURE (degrees C e l s i u s ) Fig. 3.4 Dynamic Viscosity and Hydraulic Conductivity of Oil Sands (After Scott and Kosar, 1984) 42. gas exsolution) and also for the assumption that the solid particles are incompressible. Although in their analysis provisions are made for solid particle compressibility, such a feature is not treated in a consistent manner, and their formulation is in fact only strictly true i f the solids compressibility is ignored. Later in this chapter, an approach which correctly allows for the solids flexibility as well as the inclusion of gases is presented. Harris and Sobkowicz (1977) developed a closed form solution for the change in pore pressure in terms of changes in total stress and tempera-ture, assuming a linear constitutive relationship for the soil skeleton. The model is intended to be applied to 1-D consolidation or 2-D plane strain problems. Their mathematical model, however, neither considers the compressibilities of the liquid and solid particles nor any change in volume due to structural reorientation of the particles as observed by Campanella and Mitchell (1968). Another limitation of this model is that i t may overpredict the change in volume under drained conditions since i t does not account for the quantity of fluid draining. This is an important consideration particularly in problems involving heat flow. In such cases, since the effect of temperature is to significantly increase the rate of flow of o i l , i t is then likely that the rate of change of temperature to be of the same order of magnitude as the drainage of fluid pressure from the o i l sand. Disregarding the change in volume due to drainage results in incorrect estimation of the expansion of oil sand and hence the change in the compressibility and shear strength. Byrne and Janzen (1984) extended Harris and Sobkowicz formulation to allow for the nonlinear stress-strain relationship of the sand skeleton as well as the shear induced volume change. The analysis procedure involves 43. an effective stress approach in which the stresses in the sand skeleton are computed using a finite element approach and the pore fluid pressures are computed from the gas laws together with volume compatibility between fluid and skeleton phases. At the beginning of each load or temperature incre-ment the i n i t i a l pore fluid pressure is checked against the gas/liquid saturation pressure. For situations when the former exceeds the latter, i t is tacitly assumed that no gases are present and the computation of excess pore pressures will follow the conventional technique based on Skempton's pore pressure parameters. Otherwise, excess pore pressures are computed using the expression developed from consideration of gas laws. In either case, the finite element analysis requires a knowledge of the change in total stress or the change in volumetric strain within the element. But since neither of these values are initially known, an iterative approach must be adopted which unnecessarily results in an escalating computational cost. An alternative scheme which avoids the aforementioned complications has already been described in Chapter 2. Charlwood et al. (1980) used a modified form of Byrne and Janzen's computer model 'INCOIL' for application to mine assisted in-situ methods of o i l extraction. Their model includes transient heat transfer effects. In this model, the water and bitumen are both considered to be incompressible. In fact, this is not a reasonable assumption for o i l sands because the skeleton compressibility is of the same order of magnitude as its liquid constituent. Furthermore, their proposed formulation neither considers the thermal expansion of pore fluids nor differentiates between drained and undrained expansions (i.e., the change in volume due to drainage). This relatively brief review of the previous work serves to indicate that although a satisfactory working framework has been established, 44. further developments can be made to remove some of the limitations. In the next section, a new formulation and method is proposed for computing the response of soils to thermal effects. 3.3.2 Theoretical Analysis In discussing the physical concepts associated with the response of soils to changes in temperature, i t was stated that excess pore pressures will be developed under undrained conditions because of the different expansibility characteristics of the soil matrix and its comprising com-ponents. The governing criterion for undrained conditions is that the sum of separate volume changes of the soil constituents due to any temperature and pressure change must be equal to the change in volume of the soil skeleton due to both temperature and pressure changes. By invoking such strain compatibility, theoretical expressions can be developed for relating the pore pressures to the external boundary changes. Herein, derivation of the pertinent formulations will be outlined. In order to give wide generality to the analytical treatment which follows and make i t particularly applicable to situations where high stress conditions or significant temperature changes are encountered, i t is considered that the solid particles have finite compressibilities and the voids in the soil skeleton are fill e d with both compressible liquid and gas. Additionally the assumption is made that the pores of the porous medium are interconnecting and the solid material forming the skeleton is elastic and isotropic, and the liquid in the voids is linearly compres-sible. The fluid compressibility is of course nonlinear (as indicated by Eq. 2.31) and the soil skeleton is modelled as a nonlinear elastic material 45. with dilative characteristics as will be described in some detail in Chapter 4. Consider an element of soil composed of compressible solids, liquid and gas phases as illustrated in Fig. 3.5. AT Proportion by Volume Ao Fig. 3.5. Conceptual Model of Soil Constituents By subjecting this element of soil to a change in temperature (AT) and a change in total stress (Ao) (which may arise i f the elements' displace-ments are restrained) a change in pore pressure (Au^) will be developed. Under these conditions, the overall change in volume is composed of those associated with a change in pore f l u i d (Au^), a change in the effective stress (Aa-Au), and the change due to thermal expansion. In order to account for the contribution of the first two factors, similar to Bishop (1973), the problem can be analyzed in two stages as suggested by the following diagram. 1 Aa AUj Auf Aa = A U f Auf AUf ^ ,+ i i Aa i i Auf 46. (Aa-Auf) (Aa-Auf) AUf=0 (Aa~Auf) Problem Application of an a l l round stress Aa giving rise to a change in pore fluid pressure (Auf). Stage (a) Application of a total stress and pore pres-sure change equal in magnitude to Au^. (Aa~Auf) Stage (b) Application of a total stress change of magni-tude (Ao-Auf) and zero change in pore fluid pressure. For stage (a) the following equations can be written: i) The decrease In volume of the pore fluid Ae = e Au^  (3.1) i i ) The decrease in volume of the solid grains Ae = C Au, s s f (3.2) i i i ) The decrease In volume of the soil skeleton due to the compression of the solids Ae , = (1+e) C Au£ (3.3) sk s f For stage (b) the relevant equations are: iv) A decrease in the soil skeleton 47. Ae , = (1+e) C , (Ao-Au ) (3.4) s k sk r where C , is the compressibility of the soil skeleton under drained sk conditions. v) A decrease in the volume of solid grains* Ae = C (1+e) (Ao-AuJ (3.5) s s r Next the changes In volume caused by the thermal expansion must be considered. These are: Ae = -e a AT (3.6) g g g Ae = -e a AT (3.7) Ae = -e a AT (3.8) s s s Ae = -e , a AT (3.9) sk sk s In addition, as i t has been physically observed, a change in tempera-ture may induce a change In interparticle forces or frictional resistance that necessitates some reorientations or movements of soil grains to permit Since the pores are assumed to be statically random in distribution, the area porosity on a plane intersecting the element at any level or orien-tation is equal to the volumetric porosity (n). The average normal stress in the solid material of the skeleton in any plane intersecting the element is thus: (Ao-Au) — . d - n ) 48. the soil structure to carry the same effective stress. Denoting this component of volume change by Ae , we can write: Ae st = e AT (3.10) where a can be positive or negative depending on the boundary conditions, s t state of temperature, thermal history, and the soil type. Now that the changes in volume due to each effect for every component of the soil matrix have been calculated, we can proceed by considering two cases: (a) Free case where the soil element is subjected to a change in temperature only (i.e. with no change in the external loads). (b) Restrained case where the element of soil is subjected to a change in temperature and a change in boundary pressure. The latter results from the constraint of thermal movements. 3.3.2(a) Free Case (Aa = 0) For this situation, as there is no change in the applied load, then the change in effective stress is: Under undrained conditions, volumetric compatibility can be invoked between the soil skeleton and its comprising constituents. This condition requires that: Aa1 = -Au (3.11) Ae sk = Ae + Ae. + Ae + Ae s I g st (3.12) 49. After making the relevant substitutions: e . e , e e„ e e . sk . sk . ._ . . g , £ , s sk. — Au - — Au - e g k a s AT = Au + — + j - - — ) s sk g £ s s - AT(e a + e„a„ + e a ) + e , ct AT (3.13) g g £ £ s s sk st Upon rearrangement of the terms, the following expression for the excess pore pressure can be obtained: AT(a e + a.e„ + a e - a e , - a e.) A u = g 8 L i L_S I-**:  s t sk n 14) flU e /B + e0/B„ + e /B + e ,/B - 2e ,/B g g I % s s sk sk sk s Eq . (3.14) provides the expression for predicting pore pressures under undrained conditions with no change in boundary stresses. The formulation is valid for unsaturated as well as saturated soils and i t takes account of solids and liquid compressibilities. It can be shown that for saturated soils (i.e., e =0), Eq. (3.14) simplifies to: n AT(a^-ag) - ct g t AT A u = 1/B v + n/B0 - (l+n)/B ( 3 * 1 5 ) sk. £ s where n Is the porosity; n = e /e , . £ sk Eq. (3.15) is different from the one proposed by Campanella and Mitchell (1968) since in their analysis they neglected the contribution of Eq. (3.3), i.e., the change in skeleton volume due to the change in pore fluid pressure (Au). 50. By assuming that the solids are incompressible relative to the other phases, Eq. (3.15) reduces to the following form: n AT(ct -a ) - a st AT Au = 1/B8k + n/B£ (3.16) Eq. (3.16) is Identical to the one suggested by Campanella and Mitchell which confirms the validity of the proposed approach for the case of saturated soils with incompressible solid particles. 3.3.2(b) Restrained Case (A6 = 0) In this case, the soil element under consideration is totally restrained from deforming. Such a constraint can arise i f fixity is enforced by prescribed boundary conditions. The main reason, however, for studying this particular situation is that i t forms the basis of the thermal elastic approach adopted in the finite element method described in Appendix E. Under these circumstances the change in temperature results in some change in pore pressure as well as changes in effective stress caused by the boundary stresses. The assumption of no volume change in this case implies that: Ae + Ae„ + Ae + Ae s I g st = 0 (3.17) and Ae sk = 0 (3.18) After making the relevant substitutions: 51. e e e e AuCr5- + r*- + - r 8 - ) - AT(e a + e n a „ + e a - e . a ) + ^  Ao' = 0 (3.19) VB B„ B v s s H SL gg sk st' B s £ g 0 0 s ! s k A u + ! ^ A a . - e s k a g AT =0 (3.20) s sk Eqs. (3.19) and (3.20) can be solved for the unknowns Au and Ao' to give: Bsk AT(e a + e.a„ + eot - e . c t - e . a -—) s s 11 gg skst sk s B Au = — ^ — (3.21) < e s / B s + V B * + e 8 / B g " - w f } B Aa' = B .a AT - Au (3.22) s i c S B s Eqs. (3.21) and (3.22) provide theoretical expressions for the change in pore pressure and effective stress when the soil is restrained from any change in volume. Moreover, under drained conditions when there is no change in the excess pore pressure, effective stresses can be explicitly determined from Eq. (3.22) as a function of change in temperature. The finite element procedure employed for incorporating the above formulations are described in Appendix E. 3.4 Discussion on the Application A rigorous application of the formulations developed for the change in pore pressure and effective stress requires data on a number of parameters. These are the bulk moduli, coefficients of thermal expansion, and the pro-portions in volume occupied by the skeleton, solids, liquid, and gas phases. In addition, the value of the structural reorientation coefficient ( a g t ) must also be known. Although these parameters may initially appear to complicate the use of the proposed formulae, for many engineering problems their application is quite straight-forward. This is so, because 52. for any problem Information on the porosity or void ratio is assumed to be available. The soil skeleton bulk modulus is also one of the key para-meters which is generally known. For almost a l l the conventional soil types, the solids compressibility can be reasonably assumed as negligible compared with the other soil components. Apart from o i l sands, in the case of other saturated soils which make up the majority of practical problems, the liquid compressibility can also be excluded from the calculations as its value is at least two orders of magnitude less than the soil skeleton. Indeed most liquids (e.g., water and bitumen) tend to have similar compres-sibil i t y characteristics and their values can usually be found in standard handbooks. Likewise, the coefficients of volumetric thermal expansion for many materials are readily available in the relevant publications. For the case of partly saturated soils, the bulk modulus of the gas, Bg, (assuming an ideal gas) is actually equal to the absolute pore pressure (see Appendix B). Furthermore, using the gas laws, i t is easy to show that the coefficient of thermal expansion, ctg, for an ideal gas is equal to the inverse of absolute temperature. Thus: Bo = (Pa + u) and 1 a8 = T Therefore, whether the soil is saturated or not, the practical appli-cation of the proposed formulations should not pose any difficulty. The only parameter which does, however, require specific measurement is the so-called coefficient of structural reorientation due to thermal effects 53. (a ^ ) . As has already been pointed out, this parameter is a function of soil type and boundary conditions and varies with the state of temperature. It is, therefore, not possible to relate its magnitude to other soil properties and its proper determination calls for appropriate laboratory tests. Fortunately, o i l sands being very densely packed under in-situ conditions, would be minimally affected by such a phenomenon, thus relaxing the requirement for its accurate measurement. No doubt with the accumu-lated experience which is rapidly growing in this area, i t will be possible to develop some general guidelines in estimating such a parameter for practical purposes. 3.5 Summary and Conclusions Heating of soils results in volume increases, and pore fluid pressure increases i f no drainage Is allowed. Theoretical expressions, in terms of elastic and thermal properties of the soil, have been derived for these changes by maintaining compatibility between the soil skeleton and its comprising constituents. The proposed formulations are applicable to both saturated and unsaturated soils under any condition of drainage. The theoretical analysis does not place any restriction on the stress-strain behaviour of the soil and properly accounts for the compressibility of the solid particles as well as its liquid and gas counterparts. Provisions have also been made for the component of volume change due to reorientation of soil particles as a result of thermal effects. The magnitude of this component, however, requires some form of experimental measurement as its behaviour cannot be theoretically linked to other soil properties. 54. CHAPTER 4  STRESS-STRAIN-STRENGTH MODEL 4.1 Introduction In any analytical study of the behaviour of a soil mass subject to loading or any other exciting agent, i t is necessary to select a suitable mathematical model to describe the elemental response of the soil to applied stress. Such models are termed stress-strain laws, though in fact they may involve relations between quantities other than stress and strain (e.g. stress rate and strain rate). A complete representation of soils is comprised of a stress-strain relationship and a failure criterion. No real soil can be accurately modelled by any stress-strain law, partly because of the complexity of its behaviour and partly because of its variability in the ground. The inherent nonlinear and inelastic behaviour of soils makes i t necessary to use complex theories (as compared to linear elasticity) for its modelling. The validity of any conclusions that may be drawn from analytical studies of soil behaviour are, of course, dependent upon the degree to which the mathematical model approximates real soil response. However, i t should be noted that in a l l models the real behaviour must be somewhat idealized in order to make the mathematical analysis tractable. It is not the intent of this chapter to review the existing models as their rigorous treatment has been excellently covered in several state-of-the-art papers, namely Ko and Sture (1981), Christian and Desai (1977), Naylor (1978), and Wood (1984). A thorough description of the most popular and so called 'proven' models with an assessment of their predictability is 55. covered in the North American Workshop on Generalized Stress-Strain and Plasticity for Soils (1980). Herein, attention is only focussed on describing the model incorporated in the general program developed for this study. Since the opportunity for comparing the adopted model with its rivals does not arise, a discussion section is included so that the model can be viewed from a more realistic perspective. Wherever appropriate, limitations or other difficulties inherent in the selected scheme, will be clearly pointed out. Prior to discussing the chosen model, i t is appropriate to l i s t the features which an incremental elastic model should possess in order to represent the observed properties of granular soils. These are as follows: 1. The bulk stiffness B increases as particles are pressed closer together, i.e. as mean stress increases or void ratio decreases. 2. The shear modulus G also increases with tighter packing, but more significantly, reduces with distortion. 3. A Mohr-Coulomb or similar type of failure criterion should be satis-fied. This implies that the tangential shear modulus tends to zero when yielding occurs. 4. On unloading, there is an abrupt increase in stiffness. 5. Following an i n i t i a l degree of contraction, dense sands tend to become dilatant whereas loose sands exhibit a contractive behaviour when subjected to shearing. The dilatancy depends on confining stress as well as density, becoming suppressed as this stress increases. 56. 6. The behaviour is anisotropic and stress path dependent and also influ-enced by the rotation of the principle stresses. Ideally one would strive to include as many of these facets of soil as possible in one general model, however, such a mathematically complex model would be difficult to use in engineering practice. This is because the laboratory techniques or alternative means of producing the necessary information to activate the f u l l capability of these models are not conven-tionally available. For this reason an incremental elastic hyperbolic model is adopted as i t appears to provide the optimum balance between a theoretically acceptable and a practically viable model in dealing with the majority of geotechnical problems. 4.2 Nonlinear Elastic Stress-Strain Formulation The work to be described throughout the remaining parts of this chapter is based on the premise that the soil behaviour is governed by the effective stresses. Therefore, although effective stress parameters are implied, the conventional prime symbol is omitted for the sake of tidier presentation. Assuming the material to be isotropic and piecewise elastic, the stress-strain relationship over each load increment can be described in terms of any pair of elastic parameters (E and v) or (B and G). Fundamentally though, the latter pair describe the soil behaviour in a more realistic manner since their variation can be independently controlled (e.g., at failure the shear modulus can be reduced to almost zero while the bulk modulus can be maintained at a high value) . In terms of these parameters the stress-strain relationship can be written as: 57. °a = ( B t - i V ekk 6 i j + 2 G t e i j <4-x> where the subscript t refers to the tangent values. The shear modulus is somewhat difficult to measure accurately, therefore, in practice one obtains G through its correspondence with Young's modulus E which is also a measure of the soil response to distortion. This relationship is given by: where G " 2CFvT ( 4 - 2 ) - - (1 - — ) 2 u 3 B ; For an isotropic material 1, Eq. (4.1) can be written in an incremental form: A o i j = ( B t - T G t > A ekk 6 i j + 2 G t A e i j ( 4 ' 3 ) Konder (1963) observed that for many soils the stress-strain curve under triaxial testing conditions resembled a hyperbola. Duncan and Chang (1970) extended Konder's work to propose the following expression for the tangent modulus E : R (1 - sin<(,)(a -a.) a. f 1 3 i-> , / ~ t 1 " 2c cos<{, + 20 sin* ^ ^ a ^ T ^ ( 4 , 4 ) E t a The five constants appearing in the above formulation are: 1 Similar relationships in terms of incremental stress vector and the total stress vector only prevail for the particular case when the stress Increment vector is parallel to the total stress vector (e.g., Naylor and Pande (1981)). This condition is satisfied i f there is no rotation of the principal stresses. 58. k^, = modulus number n = modulus exponent = failure ratio which is the ratio of the maximum deviator stress from the test to the ultimate value predicted from the hyperbola. c and o> are the Mohr-Coulomb strength parameters and P is the a atmospheric pressure. The most up-to-date method of determining these parameters from triaxial test can be found in Duncan et al. (1980). In order to allow for the observed nonlinear failure envelope tangent to the Mohr circles, the angle of friction a) Is considered to vary in accordance with: °3 o> = o>1 - Ao> log (_) (4.5) a where = the value of the friction angle at a confining stress of 1 atm. A<(> = the reduction in friction angle for a ten-fold increase in confining stress. The tangent bulk modulus Bfc at any stress level is given by: »t - V . < r > " ( 4 - 6 ) a where kg = bulk modulus number m = bulk modulus exponent a = mean normal effective stress which can also be replaced by the m minor principle stress 59. k can be determined from the measured values of volumetric strain B e and a as described by Byrne (1983). v m The complete stress-strain behaviour of the soil is therefore defined by E and B. These parameters in turn depend on the soil type and the level of stress within the soil and are specified in terms of seven soil constants. Tables of computed soil constants in terms of SPT blow counts are presented by Byrne and Cheung (1984) for sands and by Byrne et al. (1982) for clays. 4.3 Limitations of the Incremental Elastic models The most serious limitation of elastic models of soil behaviour is that they do not adequately account for dilatancy of real soils when they are sheared. Because of the basic nature of the elastic stress-strain relations, there is no coupling between shear stresses and volumetric strains. A corollary of this is that in undrained situations, pore pressures are poorly estimated. Another shortcoming of the selected model is that i t assumes the principal axes of strain increment always coincide with the principal axes of stress increment. Although this is a perfectly valid assumption in the elastic phase, at stress levels approaching failure; the strain increment axes coincide more closely with the principal axes of stress as is given by plasticity theory. This limitation is particularly important for conditions where an element of soil approaches failure through a reduction in the magnitude of the mean normal stress, with constant deviator stress. Under this condition the hyperbolic relationships predict expansive strains in a l l directions, whereas the strain in the direction of the major principal stress would actually be compressive. 60. Other l i m i t a t i o n s such as the i n a b i l i t y to incorporate anisotropy, e f f e c t s of the p r i n c i p l e s t r e s s r o t a t i o n , stress path dependency, and above a l l , the relevance of the model to only the conventional t r i a x i a l test s t r e s s states are shared by almost a l l the other commonly used s t r e s s -s t r a i n models. One exception to t h i s i s the p l a s t i c i t y model developed by Prevost (1978), but that model i s only r e a l l y applicable to cohesive s o i l s under undrained conditions. As already pointed out, many of these l i m i t a -t i o n s are inconsequential from a p r a c t i c a l point of view since the means of obtaining reasonable parameters f o r t h e i r representation are beyond our present c a p a b i l i t i e s , at l e a s t w i t h i n the l i m i t s of conventional techniques of sampling and t e s t i n g . The f i r s t three d e f i c i e n c i e s of the model can, however, be r e c t i f i e d i f one adopts an e l a s t o - p l a s t i c model such as the one proposed by Lade (1972). The penalty f o r the a d d i t i o n a l elegance i s the extra cost involved i n conducting more elaborate laboratory tests as w e l l as the computer a n a l y s i s , since the non-associated flow r u l e used i n the c h a r a c t e r i z a t i o n of the s o i l leads to a nonsymraetric s t i f f n e s s matrix. It i s , however, possible to preserve a l l the convenience associated with the hyperbolic models and yet improve on some of the shortcomings l i s t e d above. In the sections that follow, methods for incorporating shear-volume coupling, and the load transfer from elements which are subjected to reduction i n the mean normal stress while i n a state of f a i l u r e or simply a t t a i n i n g stress states i n v i o l a t i o n of the y i e l d c r i t e r i o n , are discussed. A method f o r inducing c o a x i a l i t y between s t r a i n s and stresses, which i s an important feature of p l a s t i c behaviour, i s outlined i n Appendix F. 4.4 Shear-Volume Coupling According to the generalized Hooke's law, which forms the basis of the hyperbolic s t r e s s - s t r a i n r e l a t i o n s h i p s , changes i n shear stress do not cause changes i n volume, i . e . , s o i l cannot be modelled as a d i l a t a n t material. This significant shortcoming can, however, be remedied by the inclusion of an additional term in the stress-strain relation as suggested by Byrne and Eldridge (1982). Eq. (4.3) then becomes: A 0 i J = ( B t - I V A Ekk 6 i j + 2 G t A e i j - B t A E v 6 i j (4'7> where Ae* is the additional volume change due to dilatancy Ae* can be related to the shear deformation by means of the dilation v J angle 6 as proposed by Hansen (1958): Ae* - sine = — - (4.8) Ay ' where A Y is an appropriate chosen invariant of the strain tensor, i.e.: A Y = Aex - Aeg The dilation angle In turn can be computed using Rowe's (1971) stress dilatancy theory: sino) - sina) m cv sine = •= 2 — — ( 4 . 9 ) 1 - sin* sin* v ' m cv where and <|> is the angle of internal f r i c t i o n at the critical void ratio, cv o) is the mobilized angle of friction. 62. The formulation incorporated for this study uses a modified form of the Rowe's stress dilatancy theory as suggested by Hughes et al. (1977): •m m") (4.10) tan(4 5 + t a n ( 4 5 + | ) = | _ tan(45 + -|^) The above formulation is in good agreement with laboratory test data when <b >d> , but It overestimates volume contractions for values of 4> m cv m smaller than <f>cv« Thus, as recommended by Byrne and Janzen (1984), no dilation is assumed to occur when d> <d> Tm Tcv The Inverse form of Eq. (4.7) can be represented as: (Ao. . - -r—— Ao. , 6 . .) + 4 Ae* 6. . 2G i j 1+v kk i j 3 v i j (4.11) The finite element procedure for incorporating dilatancy effects has been described by Byrne and Eldridge (1982). 4.5 Stress Transfer Concept With the incremental elastic stress-strain model proposed, or indeed with any incremental solution technique, the possibility of certain elements violating the failure criterion before assuming a failed state is Fig. 4.1 Permissible and Non-Permissible Stress States 63. highly likely. Such a situation is schematically presented in Fig. 4.1. In this figure, the application of load increment causes the stress path to move outside the region bounded by the failure surface as indicated by o.j and a^. The occurrence of this violation is due to the fact that the element is considered elastic prior to failure and thus capable of accept-ing any imposed load. However, i f i t happens that its state lies just below the failure surface or the magnitude of the applied incremental load is too large, then i t is possible for the new state to move well into the illegal zone which clearly introduces errors in the analysis. One way of circumventing this problem is to control the applied loads such that the stress state of each affected element falls on the failure surface or very close to i t . But this may involve too many load increment steps throughout a load sequence because i t is not easy to assess the optimum magnitude a priori. With the incremental elastic analysis employed here another situation may arise that can lead to a stress violation. This occurs when an element whose state lies on the failure surface is subjected to unloading of normal stress. For this problem, the shear modulus is set to zero (or very close to zero) which means that the shear stress in that element will not change for a further increment of loading. Since the increment of load is such that the normal stress decreases, then, unless the strength envelope is horizontal, the predicted stress state will violate the failure criterion. In order to overcome the problems associated with these offending stress states, corrections must be applied to bring the stresses back to the failure surface after each load increment or iteration in the program. Several options exist for applying such a correction depending upon which path is considered to be the most appropriate one for tracking the stresses back along i t . In Appendix G three alternative approaches for determining 64. the correcting stresses are suggested. The method for transferring such overstress to adjacent elements and i t s f i n i t e element implementation are also presented. 4 . 6 Methodology for Representing the Behaviour at Failure The nonlinear and stress dependent stress-strain properties of soils are approximated by varying the values of both Young's modulus and bulk, modulus in accordance with the calculated stresses. In addition, the dilational plastic volume changes are incorporated i n the analysis by a method similar to that used to include thermal volume changes i n thermal ela s t i c i t y (see Appendix E). In following the modified Euler technique adopted herein (see Chapter 6 ) , any incremental changes to the boundary conditions (e.g. application of load, temperature, displacement, or pore pressure) is analysed twice. The f i r s t time one uses moduli values based on the stresses at the beginning of the increment, and the second time moduli values based on the average stresses during the increment are used. If the dilation or stress transfer options are invoked, then additional iterations w i l l be required to satisfy their respective convergence c r i t e r i a . The stress transfer is normally carried out i f the ratio for the developed deviator stress with that at failure for any element (or alter-natively at any integration point) exceeds unity by say more than 5%. This tolerance can be varied depending on the accuracy desired for the problem under consideration. The exceedence of this ratio i n any element signifies the occurrence of failure. For such a state, the assumption i s made that the shear modulus (and hence the Young's modulus) is reduced to a very small value (e.g., 0.1% of the bulk modulus), and the Poisson's ratio i s 65. increased to approximately 0.5 (0.499 is the upper limit used). In analyzing problems subject to severe loading conditions such as unloading of a wellbore (see Chapter 7), where a significant number of elements develop instability, care must be exercised to ensure that the failed elements retain a state on the yield envelope. Unless such a measure is taken, some elements which have presumably failed may attain stress states well below the failure envelope. This can happen as a result of a large number of load increments and load shedding iterations in a problem. It is, therefore, necessary to adjust the stresses in these elements by performing a reversal of the load shedding process while keeping their shear modulus close to zero. Unless this step is taken, the stresses in some of the failed elements can become far less than what they are able to carry and as a consequence greater loads will be transferred to the neighbouring elements. This may eventually spoil the entire analysis and produce unrealistic stress and pore pressure regimes in the plastic zone. So for every element failing during the load shedding process or during a loading cycle which is along the direction that initiated failure, the shear modulus will be maintained at a low value with stresses on or very close to the yield envelope. If, however, the direction of loading changes or as a result of drainage some of the elements which had previously failed wish to move below the failure envelope, then they are permitted to do so and regain their stiffness properties. Thus no restric-tion is placed in keeping an element in a plastic state throughout the analysis. Development of any stress less than -c cot<J> implies failure in tension. Under these conditions, both bulk and shear moduli are reduced to very small values. The recommended scheme is to assign moduli values 66. which are relatively much smaller than the current fluid bulk modulus. For an element of soil which experiences a tensile failure, a l l the stresses will be set equal to -c cottfi and the residual stresses will be either transferred from or to the element in order to ensure that they retain such a value. In analyzing problems where the tensile zone actually breaks away from the soil mass (e.g., in the vicinity of open wellbores) the option is available to stop those elements from both regaining stiffness and attain-ing a stress state different from -c cot<f>. This is to simulate the condition that those elements are removed and thus unable to provide any further contribution to the analysis. Although the above mentioned approach ensures that the soil responds in a realistic manner in the direction of the loading which has caused the failure, i t does not properly represent the overall behaviour of the soil. This is because soils, being anisotropic in nature, have shearing resist-ance in other directions which the present isotropic model is unable to incorporate. Thus i t is essential to exercise some caution in examining the computational results pertaining to the failed zones in a problem and to apply good judgement in assessing the influence of any anomaly which might have led to incorrect modelling of the behaviour in those regions. 4.7 Summary and Conclusions The effective stress-strain constitutive law selected for the general program is a piecewise elastic hyperbolic formulation with the ultimate strength being defined by the Mohr-Coulomb failure criterion. Although this is not the most sophisticated model from a theoretical point of view, i t does have the advantage that its parameters have well understood physical significance and coefficients are easily evaluated from conven-67. tional laboratory tests. Furthermore, the widespread use of the model in many applications has culminated in a wealth of information on the range of values of the model constants for a variety of soil types under both drained and undrained conditions. In order to overcome some of the limitations of the model, additional features such as dilation and load transfer effects have also been incor-porated in the program. The option of re-distributing the over-stresses from the elements which are in a state of failure onto the adjacent elements is considered to constitute an extremely important part of any soil model. This is particularly true in situations where the response of soil in the failed zones is of interest. For instance in considering the behaviour of soil around excavated tunnels or shafts where the effective stresses are reduced to very low values, the use of the stress transfer option results in significantly different magnitudes and modes of movements and also different levels of pore pressures and stresses as compared to the predictions using the normal course of analysis. Some of the examples presented in Chapter 7 reflect the importance of such a treatment. 6 8 . CHAPTER 5 NON-LINEAR FINITE ELEMENT CONSOLIDATION ANALYSIS OF UNSATURATED SOILS 5.1 Scope This chapter is concerned with the time dependent or transient response of soils. The objective here is to develop suitable formulations for finite element analysis of consolidation in soils which may be in a partially saturated state. This task is pursued by treating the soil and fluid as having a nonlinear behaviour. The former is achieved by express-ing the stress-strain relationship with a hyperbolic elastic model while the compressibility characteristics of the latter is governed by the state of ambient pressure, degree of saturation, and void ratio which may vary throughout the consolidation process. The nonlinear flow of pore fluid is modelled by expressing the permeability as a function of void ratio, degree of saturation, and coefficient of viscosity