A STRESS-STRAIN MODEL FOR THE UNDRAINED RESPONSE OF OIL SAND by KA FAI HENRY/CHEUNG B.Sc, University of Manchester, 1983 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in FACULTY OF GRADUATE STUDIES Department of C i v i l Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA A p r i l 1985 © KA FAI HENRY CHEUNG » ?985 In presenting t h i s thesis i n p a r t i a l fulfilment of the requirements f o r an advanced degree at THE UNIVERSITY OF BRITISH COLUMBIA, I agree that the Library s h a l l make i t f r e e l y available f o r reference and study. I further agree that permission for extensive copying of this thesis f o r scholarly purposes may be granted by the Head of my Department or by his or her representative. I t i s understood that copying or publication of this thesis f o r f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of C i v i l Engineering THE UNIVERSITY OF BRITISH COLUMBIA 2075 Wesbrook Place Vancouver, B.C. CANADA, V6T 1W5 Date: A p r i l , 1985 ABSTRACT An e f f i c i e n t undrained model for the deformations analyses of o i l sand masses upon undrained loading i s presented i n t h i s thesis. An analysis which couples the s o i l skeleton and pore f l u i d s i s used. The s o i l skeleton i s modelled as a non-linear e l a s t i c - p l a s t i c i s o t r o p i c material. In undrained conditions, the constitutive r e l a t i o n s h i p s f o r the pore f l u i d s are formulated based on the i d e a l gas laws. The coupling between the s o i l skeleton and the pore f l u i d s i s based upon volume compatibility. The undrained model was v e r i f i e d with the experimental r e s u l t s and one dimensional expansion of s o i l sand cores. Comparisons between computed and measured responses are i n good agreement and suggest that this model may prove useful as a tool i n evaluating undrained response of o i l sand. The response of a wellbore i n o i l sand upon unloading was analysed using the developed model. Such analyses are important i n the r a t i o n a l design of o i l recovery systems i n o i l sand. iii TABLE OF CONTENTS Page Abstract i i Tabele of Contents i i i L i s t of Tables vii L i s t of Figures viii L i s t of Symbols x Acknowledgements CHAPTER 1 xiii Introduction 1.1 Introduction 1 1.2 Behaviour of O i l Sand 2 1.3 The Scope 5 1.4 The Organization of Thesis 6 CHAPTER 2 Review of Previous Work 2.1 Introduction 2.2 Mathematical Model for Undrained Behaviour 8 of O i l Sand CHAPTER 3 2.2.1 Harris and Sobkowicz 8 2.2.2 Dusseault 12 2.2.3 Byrne and Janzen 14 Stress-Strain Model 3.1 Introduction 17 3.2 Development of the Undrained Model 17 3.3 Constitutive Relations 3.3.1 Incremental Non-linear E l a s t i c S o i l s Model 22 3.3.2 Incremental Non-linear E l a s t i c Fluid Model 25 iv Page 3.3.2.1 General 25 3.3.2.2 Partly M i s c i b l e Gas/Liquid Mixture 27 a. Air/Water Mixture 28 b. Carbon Dioxide, Air/Water 30 Mixture c. CHAPTER 4 Gas/Bitumen and Water Mixture F i n i t e Element Formulations 4.1 Introduction 4.2 The Plane Strain Formulation 4.3 4.4 CHAPTER 5 36 4.2.1 The Constitutive Matrix [D] 36 4.2.2 The Strain Displacement Matrix [B] 40 4.2.3 The S t i f f n e s s Matrix [K] 44 The Axisymmetric Formulation 4.3.1 The Constitutive Matrix [D] 46 4.3.2 The Strain Displacement Matrix [B] 50 4.3.3 The S t i f f n e s s Matrix [K] 53 Load Shedding Formulation 55 Comparisons with E x i s t i n g Solutions 5.1 Introduction 5.2 Comparisons with Theoretical Solutions 5.3 33 60 5.2.1 E l a s t i c Closed Form Solutions 60 5.2.2 E l a s t i c - P l a s t i c Closed Form Solutions 61 Comparisons with Observed 5.3.1 Data One Dimensional Unloading of O i l Sand 68 Page 5.3.2 CHAPTER 6 T r i a x i a l Tests on Gassy Soils 72 Stresses Around a Wellbore or Shaft i n O i l Sand 6.1 Introduction 79 6.2 General Model Description 80 6.3 Theoretical Solutions for Stresses Around a Wellbore 6.3.1 6.4 80 Stresses 6.3.1.1 Stresses i n E l a s t i c Zone 83 6.3.1.2 Stresses i n P l a s t i c Zone 85 6.3.1.3 Radius of P l a s t i c Zone 87 6.3.2 Stability 89 6.3.3 Pore Pressure P r o f i l e 90 Comparisons of Predicted Response and Closed Form Solution 6.5 93 Analysis of Response of Borehole i n O i l Sand on Unloading 6.5.1 Undrained Response 6.5.2 Drained Response 6.5.3 Implications of Undrained and Drained Analysis 6.6 CHAPTER 7 93 100 104 Application of o i l recovery 104 Summary and Conclusions 108 Bibliography 110 Appendix A 114 Appendix B 116 Appendix C 119 VI Page Appendix D Appendix E Appendix F Appendix G Appendix H Appendix J . _, 136 vii LIST OF TABLES Table Page A comparison of computed and measured r e s u l t s (Test 11, Sobkowicz) 74 Vlll LIST OF FIGURES Figure _ __ 1.1 Schematic spring analogy for o i l sand 3.1 S t r e s s - s t r a i n curves f o r drained t r i a x i a l Page 3 tests on loose sand 24 3.2 Phase diagram f o r gassy s o i l 32 3.3 Phase diagram for o i l sand 34 4.1 Stresses associated with load shedding 56 5.1 thick wall cylinder 62 5.2 Stresses and displacements around c i r c u l a r opening in an e l a s t i c material 5.3 63 Stresses and displacements around c i r c u l a r opening in an e l a s t i c - p l a s t i c material 64 5.4a Model f o r one dimensional unloading of o i l sand 70 5.4b Response of o i l sand to one dimensional unloading 5.5 Comparisons of predicted and observed pore pressure 5.6 Comparisons of predicted and observed strains 77 6.1 Outline of the wellbore problem 81 6.2a F i n i t e element mesh for wellbore problem 82 6.3b Mohr Coulomb f a i l u r e envelope 84 6.3 Idealised flow to a wellbore 92 .... .. 71 76 LIST OF FIGURES - Continued Figure 6.4 Page Stresses around a wellbore i n an e l a s t i c - p l a s t i c material 6.5 Undrained 94 response of wellbore i n o i l sand on unloading 96 6.6 Pore pressure p r o f i l e around a borehole 101 6.7 Drained response of a wellbore i n o i l sand on unloading 6.8 Comparisons of undrained and drained 102 response of a wellbore i n o i l sand at a support pressure of 3500 kPa 6.9 Comparisons of undrained and drained 105 response of a wellbore i n o i l sand at a support of 3200 kPa pressure 106 X LIST OF SYMBOLS The following i s a l i s t of the commonly used symbols i n t h i s thesis. Multiple use of several symbols i s unavoidable because of the complexity of the formulations. The symbol w i l l be defined immediately i n the text where the use of symbols d i f f e r s from those l i s t e d below. SYMBOL MEANING 3 compressibility A change a stress e strain v Poisson's r a t i o (j) f r i c t i o n angle a temperature s o l u b i l i t y constants a r r a d i a l stress o"g tangential stress a z v e r t i c a l stress B bulk modulus e . void ratio E Young's modulus G shear modulus H Henry's constant k K permeability apparent bulk f l u i d modulus xi bulk f l u i d modulus n porosity P pressure q flow rate r radius (variable) R radius (constant) R r c w radius of p l a s t i c zone radius of wellbore S saturation T temperature u pore f l u i d pressure V volume xii Subscripts 1 final a air f pore f l u i d g gas i i n i t i a l ( i n t e r n a l i n Chapter 6) 0 o i l (outer i n Chapter 6) S s o i l skeleton v volumetric w water cb combined dg dissolved gas fg free gas Tg t o t a l gas Superscripts e elastic f final 1 initial p plastic xiii ACKNOWLEDGEMENTS I am greatly indebted to my supervisor, Professor P.M. Byrne, for his guidance, encouragement end enthusiastic interest throughout t h i s research. I would also l i k e to thank Professor Y.P.Vaid f o r reviewing the manuscript and making valuable suggestions. My colleagues, U. Atukorala, F. Salgado, J . She, H. V a z i r i and e s p e c i a l l y , C. Lum, shared a common active interest i n s o i l mechanics. I thank them a l l f o r t h e i r h e l p f u l discussions and constructive c r i t i c i s m s . Appreciation i s extended to Ms. S.N. Krunic for typing the manuscript and her patience during the preparation of t h i s t h e s i s . Support and assistance provided by the Natural Science and Engineering Research Council of Canada i s acknowledged with deep appreciation. A s p e c i a l thanks i s extended to my family and P r i s c i l l a , f o r their constant support and encouragement. 1 1.1 INTRODUCTION Many schemes for o i l recovery require open excavations, tunnels or wellbore i n o i l sand. As a r e s u l t , an accurate and e f f i c i e n t analysis of stresses and deformations around these openings i s becoming increasingly important. Mathematical models have been developed by Byrne et a l (1980), Dusseault (1979), and Harris and Sobkowicz (1977) to investigate these problems. An e f f i c i e n t undrained model f o r analysing the stresses and deformations around open excavations, tunnels and wellbores i n o i l sand i s presented i n t h i s t h e s i s . O i l sand i s comprised of a dense sand skeleton with i t s pore spaces f i l l e d with bitumen, water and free or dissolved gas. The presence of bitumen reduces the e f f e c t i v e permeability of o i l sand, hence undrained conditions occur on rapid unloading. Gas evolves from pore f l u i d s during unloading when the pore f l u i d pressure i s below gas saturation pressure. Because gas exsolution takes some time to occur, two undrained conditions a r i s e , (1) an immediate or short term condition i n which there i s no time f o r gas exsolution, and (2) a long term or equilibrium condition i n which complete gas exsolution has occured. Both of these conditions are considered i n this thesis. rate of gas exsolution i s not considered herein. The The drained analysis with pore pressures under steady-state conditions i s also addressed. The sand skeleton i s modelled as a non-linear e l a s t i c - p l a s t i c porous material. The f l u i d s s t r e s s - s t r a i n relationships f o r undrained condition are formulated on the basis of ideal gas law. For the undrained condition, the pore pressure changes ;.re computed from the constant of volume compatibility between the sand skeleton and the pore fluids. It i s assumed that the pore f l u i d pressures are known i n the drained analysis. 2 V a l i d a t i o n of the stress s t r a i n model i s made by comparing responses with t h e o r e t i c a l solutions and observed data. The observed data are the experimental results i n Sobkowicz's doctorate thesis. The stresses and deformations around a wellbore i n o i l sand upon unloading i s investigated using the new s t r e s s - s t r a i n model. 1.2 BEHAVIOUR OF OIL SAND O i l sand i s comprised of a dense, highly incompressible, uncemented, interlocked skeleton with pore spaces f i l l e d by water, bitumen, and dissolved or free gas. The interpenetrative structure leads to the low i n - s i t u void r a t i o and high shear strength. The response of o i l sand i s mainly governed by the rate of loading. Undrained conditions occur as a r e s u l t of: 1) low e f f e c t i v e permeability to pore f l u i d s 2) large amount of dissolved gas i n pore f l u i d s 3) rapid unloading The response of o i l sand may be physically modelled by a set of springs shown i n Figure 1.1. The o i l sand i s s p l i t into two load carrying components - s o i l skeleton and pore f l u i d s . The s o i l skeleton compressibility, 3 , characterizes the deformation of s o i l skeleton, g which r e s u l t s i n a change i n e f f e c t i v e s t r e s s . The pore f l u i d compressibility, 8^, characterizes the deformation of pore f l u i d s (because of free gas) which results i n a change i n pore pressure. 3 a + Aa X X X X 4 4 X 4 K=> O 6 § § f o § Water Soil Skeleton 5 \ § Stress = a'+ Aa § ) «o § Oil (bitumen) Pore Fluid Pressure = u + Au o § Gas Rigid Frame, No Lateral Yield a = a' + u Ao = Aa'+ A u A V = /(Aa, Au) Fig.1.1 - Schematic spring analogy for o i l sand ( After Dusseault , 1979 ) 4 Strain compatibility between the pore f l u i d s and s o i l skeleton controls the r e l a t i v e magnitudes of the changes i n pore pressures and e f f e c t i v e stress which together equal the change i n t o t a l s t r e s s . Stress changes (unloading) are shared between the sand skeleton and pore f l u i d s according to t h e i r c o m p r e s s i b i l i t i e s . When the pore pressure i s above the gas saturation pressure and the o i l sand i s 100% saturated, the stress changes w i l l be accommodated by the pore f l u i d s because t h e i r compressibilities are lower. Once the pore pressure drops below the gas saturation pressure, gas s t a r t s to evolve which increases the fluids' compressibilities. Therefore the sand skeleton becomes the less compressible phase and takes up the stresses rather than the pore fluids. When the e f f e c t i v e stress drops to zero, the s o i l skeleton i s very compressible r e l a t i v e to the pore f l u i d s ; hence any further decrease i n t o t a l stress i s e n t i r e l y accommodated by the pore pressure. The gas exsolution takes some time to occur, hence two undrained conditions (Sobkowicz, 1982) a r i s e . The expressions 'short term' and 'long term' w i l l be applied to these processes exclusively. 1) 'Short term' undrained i n which there i s no time for gas exsolution 2) 'Long term' undrained i n which equilibrium state has been reached, i . e . completion of gas exsolution. In the f i e l d , the unusal behaviour of o i l sand manifests i n a number of ways. They include: 5 1) Volumetric expansion of 5 to 15% occurred when core samples were l e f t i n an unconfined state, i . e . not retained by p l a s t i c core sleeves (Dusseault 1980; 2) Core samples spontaneously Byrne et a l 1980). s p l i t l o n g i t u d i n a l l y and perpendicularly to the core axis, effervescence was observed on several freshly recovered cores (Hardy and Hemstock 1963). 3) Retrogression of slopes i n o i l sand on rapid excavation. 4) O i l sand at the base of excavation subjects to softening and heaving, followed by settlement on reloading. When the decrease i n external stress occurs over a period of time, the evolved gas has time to drain o f f and e f f e c t i v e stress does not go to zero which results i n an undisturbed sand skeleton. i n - s i t u density and hence high shear strength i s retained. Its high Such situations can be seen on the exposures of o i l sand deposits along the Valley of Athabasca River where erosion (unloading) by the r i v e r has occurred over thousands of years. These o i l sand deposits are standing on steep stable slopes with slope angles i n excess of 60° and heights up to 60 metres, exhibiting high strength (Harris and Sobkowicz 1977). 1.3 THE SCOPE The purpose of this study i s to present a s t r e s s - s t r a i n model ( V a z i r i , 1985) f o r the deformation analyses of o i l sand masses, i . e . to explain the behaviour of o i l sand as described i n Section 1.2. Modelling the undrained response of o i l sand requires the pore f l u i d s pressure to be numerically evaluated. For the undrained condition, pore f l u i d pressures are computed from the constraint of volumetric compatibility between the sand skeleton and pore f l u i d phases. Compressibilities of sand skeleton and pore f l u i d s have to be evaluated i n the new s t r e s s - s t r a i n model. A non-linear s t r e s s - s t r a i n r e l a t i o n s h i p (Duncan et a l 1970) i s adopted for the sand skeleton. compressiblities The of pore f l u i d s are formulated on ths basis of ideal gas laws. The new s t r e s s - s t r a i n model i s incorporated into f i n i t e programmes (INCOIL, MHANS). element For the v a l i d a t i o n of the s t r e s s - s t r a i n r e l a t i o n s h i p , the computed results are compared with the t h e o r e t i c a l solutions and observed data. The validated programme was used to study the unloading response of stress and deformations around a deep wellbore i n o i l sand. 1.4 ORGANIZATION OF THE THESIS This thesis consists of seven chapters. A review of previous work on undrained models f o r o i l sand i s given i n Chapter 2. This chapter concentrates on the examinations of their c a p a b i l i t i e s and shortcomings. A s t r e s s - s t r a i n model which was developed by V a z i r i i s presented i n Chapter 3. Appropriate s o i l skeleton and pore f l u i d s constitutive relationships are also recommended i n t h i s chapter. Chapter 4 summarizes the f i n i t e element formulations used i n the development of the programme. V a l i d a t i o n of the developed model by comparing predicted response with existing t h e o r e t i c a l solutions, f i e l d and laboratory observations i s mentioned i n Chapter 5. The response of a wellbore o i l sand upon unloading i s investigated i n Chapter 6 using the validated f i n i t e element programme. A summary of work, and major conclusions are presented i n Chapter 7. 8 CHAPTER 2 2.1 : REVIEW OF PREVIOUS WORK INTRODUCTION Theoretical solutions f o r the undrained response of o i l sand were not a v a i l a b l e u n t i l 1977 because of i t s unusual behaviour. Due to the increase i n demand for construction i n o i l sand formation, such as open p i t mining, tunnels and deep shafts, a considerable amount of research has been done on this topic since 1977. These developed t h e o r e t i c a l r e l a t i o n s share the same basic approach of coupling the s o i l skeleton and pore f l u i d s together. Pore pressure changes are computed from the constraint of volumetric compatibility between the s o i l skeleton and pore f l u i d s . Harris and Sobkowicz's model i s capable of evaluating pore pressure change upon a stress or/and temperature change. This model was then extended by Byrne et a l and incorporated i n a f i n i t e element programme. Dusseault's approach i s r e s t r i c t e d to one dimensional problems. A c a r e f u l examination of these t h e o r e t i c a l models and t h e i r applications i s made and several shortcomings are also reviewed. 2.2 MATHEMATICAL MODELS FOR OIL SAND 2.2.1 Harris and Sobkowicz (1977) Harris and Sobkowicz presented a mathematical model to analyse the undrained response of o i l sand subjected to changes of stress and/or temperature. This i s the f i r s t a n a l y t i c a l model developed to explain the behaviour of o i l sand such as: 1) Movement and s t a b i l i t y of slopes and tunnels formed i n o i l sand. 2) Settlements or heave of structures placed i n o i l sand (e.g. hot o i l tank). 3) Heave at the_base of excavations i n o i l sand. The model can be explained by the same spring analogy as shown i n Figure 1.1. The response of o i l sand to changes i n stresses or/and to changes i n temperature may be computed by the following equations: a) One-Dimensional Analysis Au can be obtained from a quadratic equation L * Au + M * Au + N = 0 2.2 2 where L = 3 s Tl M = 3 (P. - Aoi) + n + Pi — (n H s i g T w w a 1 / N = P ± [n (1 - 1 v - A 3 ] 0 l s + n H ) oo' " P i AT ^ - ( n ^+ n^) 10 b) Two-Dimensional Analysis from s t r a i n compatibility again Aa 3 + A Au = (Ao ! - Ao ) - 3 n6 1+ * 2.3 f s and P Au z + Q Au + R = 0 2.4 where P = B s T1 Q = 8 ( P . - X) + n + P i (n H s i g T w w v R = P. [n l g L 1 (1 - — ) T + n H ) oo' — - X 0 ] - P, AT — s i T J X = A02 + Ap (Aai - Ao"2) and = volumetric B g n strain = s o i l skeleton compressibility = porosity of s o i l n Q = porosity of o i l n w = porosity of water n = porosity of gas Aa = change i n t o t a l stress ^ v (n a w w + n a ) 0 0 11 P = i n i t i a l pore f l u i d pressure Au = change i n pore f l u i d pressure Ao"i = change i n major p r i n c i p a l stress Ao"3 = change i n minor p r i n c i p a l stress T a = standard temperature (288 °K) = i n i t i a l pore f l u i d temperature °K Tl = f i n a l pore f l u i d temperature °K AT = change i n temperature = temperature s o l u b i l i t y constant for gas i n water = temperature s o l u b i l i t y constant for gas i n o i l H W = Pressure s o l u b i l i t y constants f o r gas i n water H q = Pressure s o l u b i l i t y constants f o r gas i n o i l Ap = s o i l skeleton dilatancy factor Harris and Sobkowicz examined the response of o i l sand during excavation and reloading of a square footing, and the behaviour of a tunnel excavation i n the same material. Their results w i l l allow an assessment of the a p p l i c a b i l i t y and shortcoming of the theory. 1) This solution..incorporated a linear constitutive relationship f o r the s o i l . As o i l sand behaves l i k e an e l a s t i c - p l a s t i c material, there w i l l be a p l a s t i c zone developed adjacent to the tunnel wall on unloading when the effective.stresses are such that f a i l u r e (Mohr Coulomb C r i t e r i a ) occurs. 2) The extent of the p l a s t i c zone i s only an approximation because r e d i s t r i b u t i o n of stresses has not been taken into account during the formation of p l a s t i c zone. 3) 2.2.2 An i t e r a t i v e procedure i s required to obtain Au DUSSEAULT Dusseault (1979) extended the one dimensional Skempton's B equation to analyse the behaviour of cohesionless amount of free or dissolved gas i n pore f l u i d s . materials with large He presented a more rigorous derivation for the compressibility of pore f l u i d s , and coupled the compressibility with that of s o i l skeleton ( e - a' relationship) i n Skempton's B equation. This model shares the same basic idea as the previous (Harris and Sobkowicz), which has two skeleton and pore f l u i d s . one load carrying components - s o i l The r e l a t i v e c o m p r e s s i b i l i t i e s of these components control the magnitude of changes in pore pressure and effective stress which together balance the change i n t o t a l stress i n an undrained state. The t r a d i t i o n a l Skempton's B equation i s 2.5 Aa 1 + n- 8^ = compressibility of pore f l u i d (water) 3 The = compressibility of s o i l extended one i s skeleton ^ = l / [ l + f(u,o)] 2.6 a+bAn(a-u)-e -e +H e -Hi e f(u,a) = . D [eg + e 3 + o o w w J J U a, b = void r a t i o - e f f e c t i v e stress relationship parameters e o', ew', e g = void r a t i o s : o i l», water,» e>gas H , H = Henry's constants : o i l , water o' w 3, 3 Q u, a w = compressibility : o i l , water = current values of pore pressure and t o t a l stress Dusseault applied this model to investigate the response of element of o i l sand upon unloading. and a deeper one (500 m). He examined a shallow case (15 His results w i l l allow an assessment of an m) the a p p l i c a b i l i t y and shortcoming of the theory: 1) In order for the solution to be numerically stable, reduction further i n t o t a l stress i s assumed to be e n t i r e l y taken by the pore pressure once e f f e c t i v e stress drops to zero. 2) This solution incorporated relationship for the 3) a non-linear constitutive soil. Accurate e - a' relationship or compressibility of o i l sand i s extremely d i f f i c u l t to obtain since they are very sensitive to sample disturbances. samples have been cored so f a r . No undisturbed o i l sand 4) 2.2.3. An i t e r a t i v e procedure i s required to get Au. Byrne and Janzen Byrne et a l (INCOIL, 1983) developed an incremental a n a l y t i c a l f i n i t e element method for predicting stresses and deformations i n excavations and around tunnels i n o i l sand using a nonlinear e l a s t i c sand skeleton with shear d i l a t i o n . An extension to Harris and Sobkowicz's model was used to evaluate the pore pressure change, Au. The f i n i t e element program, INCOIL, can handle both undrained and drained analyses. For the undrained condition, pore f l u i d pressures are computed from the constraint of volumetric compatibility between the sand skeleton and the pore f l u i d s . For the f u l l y drained condition, i t i s assumed that the pore f l u i d pressures are known. The general framework of the f i n i t e element model for o i l sand are as follows: [K] (fi) = (Af) - (K ) (Au) 2.7 w where u Ti | — * AT * (n a +n a ) - u fn * ( 1 - — ) T w w o o ' «oL g T ' v Au = u a L v a Ti L (n H +n H ) + n - Ae w w o o g v Derivation of equation 2.7 i s presented i n Appendix A Ae 1 v J 2.8 15 The e f f e c t i v e stress may (Ae) = [c] be evaluated from (A6) 2.9 (Aa') = [D»] (Ae) 2.10 where [c] = the matrix which depends on element geometry [D'J = the matrix property matrix ( a f f e c t i v e stress) [Aa'] = the change i n e f f e c t i v e stress vector [Ae] = the change i n s t r a i n vector [K] =» the element s t i f f n e s s matrix (A6) = the incremental element nodal deflections vector (Af) = the incremental element nodal forces vector (K ) = the pore pressure load vector n , n , n^ = porosity : water, o i l and gas phase a = temperature s o l u b i l i t y constant : water, o i l w w > H , W u Q a Q H Q = pressure s o l u b i l i t y constant : water, o i l = reference (atmospheric) pressure T a = reference temperature (288°K) T Q = i n i t i a l temperature (°K) Tl - f i n a l temperature (°K) AT = change i n temperature ( T i - T ) (°K) o = i n i t i a l absolute pore pressure U Q Ae = volumetric s t r a i n (compression p o s i t i v e ) They examined the response of c y l i n d r i c a l shaft i n o i l sand on reduction of support pressure and the response of an element of o i l sand to one dimensional unloading. The l a t t e r case i s to simulate core samples of o i l sand l e f t i n an unconfined state. Their results will allow an assessment of the a p p l i c a b i l i t y and shortcoming of the theory: 1) The predicted expansion of the core on unloading i s small compared with those measured i n the f i e l d which i s 5-15% when core samples of o i l sand are l e f t i n an unconfined state. 2) An i t e r a t i v e procedure i s required to obtain Au. CHAPTER 3 3.1 : STRESS-STRAIN MODEL INTRODUCTION It i s noted that the mathematical models described i n Chapter 2 have quite a few shortcomings. Therefore, a more sophisticated and e f f i c i e n t undrained model was developed (Naylor 1973, V a z i r i 1985) and i s incorporated i n a f i n i t e element programme. A t o t a l stress approach coupling the s o i l skeleton and pore f l u i d s i s used . Pore pressure changes are computed from the constraint of volume compatibility between the s o i l skeleton and pore f l u i d s under undrained conditions. Separate c o n s t i t u t i v e relationships f o r s o i l skeleton and pore f l u i d s are required i n the new undrained model. An incremental non-linear e l a s t i c and i s o t r o p i c s t r e s s - s t r a i n model as described by Duncan et a l i s adopted for the s o i l skeleton. Depending on the component of the pore f l u i d s , d i f f e r e n t formulations f o r the non-linear e l a s t i c and i s o t r o p i c s t r e s s - s t r a i n relationships of pore f l u i d are derived. 3.2 DEVELOPMENT OF UNDRAINED MODEL The e f f e c t i v e stress concepts (Terzaghi) i n conventional soil mechanics seem to be applicable to determine the shear strength of o i l sand (Hardy and Hemstock). However, the pore pressures of f l u i d s i n the o i l sand have to be numerically evaluated i n the e f f e c t i v e stress approach. Hence, a quantitative r e l a t i o n s h i p between the magnitudes of stress release and pore pressure i s required. When a saturated s o i l mass i s subjected to undrained loading, 18 stress change must be shared between the s o i l skeleton and pore f l u i d (Bishop and E l d i n 1950). A theoretical expression for the relationship between t o t a l stress change and r e s u l t i n g pore f l u i d pressure change was derived by Skempton (1954) which i s the Skempton's B equation. Using the f i n i t e element method, C h r i s t i a n (1968) introduced an e f f e c t i v e stress approach for s o i l s subjected to undrained loading, enabling r e s u l t i n g pore f l u i d pressure to be evaluated. Programmes incorporating this alternative are r e l a t i v e l y i n e f f i c i e n t . Naylor (1973) developed a more elegant approach which allows excess pore pressure to be computed e x p l i c i t l y i n terms of material skeleton s t i f f n e s s parameters and an independently s p e c i f i e d pore f l u i d stiffness. However, Naylor only considers s o i l s that are two-phase system - solids and water. Due to the presence of bitumen, free and dissolved gas i n o i l sand, the above mentioned approaches are inadequate to describe the undrained behaviour of o i l sand. Not u n t i l 1977, Harris and Sobkowicz developed an a n a l y t i c a l expression incorporating a linear constitutive r e l a t i o n s h i p f o r the o i l sand, to r e l a t e the change i n pore pressure to the change i n total stress. Dusseault (1979) extended Skempton's one dimensional B equation to model the equilibrium behaviour of o i l sand. Byrne et a l (1980) studied the behaviour of o i l sand by using f i n i t e element method. They extended Harris and Sobkowicz's model and incorporated i t into Christian's f i n i t e element formulation. As there are shortcomings of the approaches developed by Harris and Sobkowicz, Dusseault and Byrne et a l , a more sophisticated numerical model i s required. V a z i r i adopted Naylor's approach and extended i t to model the undrained behaviour of o i l sand using f i n i t e elements. The s t r e s s - s t r a i n model f o r o i l sand behaviour i s based on the following assumptions. 1) Volumetric change of s o i l skeleton i s governed by the e f f e c t i v e stress. 2) Liquids and gas i n the voids are at the same pressure, i . e . effect of surface tension between the pore f l u i d phases are neglected. 3) Free gas i n the pores behaves i n accordance with c l a s s i c gas laws with respect to pressure. Gas comes out from solution i n accordance with Henry's law. Henry solubility constants are constant. 4) The compressibility of s o i l grains i s n e g l i g i b l e , and has no contribution to the volume changes. 5) The gas Is i n the form of occluded bubbles inside the pore fluid. The t o t a l stress constitutive law may be written as: (Aa) = [D] (Ae) 3.1 where (Aa), (Ae) are the incremental t o t a l stress vectors and s t r a i n vectors respectively, and [D] i s the material property matrix ( t o t a l stress). Computation of element s t i f f n e s s matrix, assembly and solutions f o r displacements proceeds along the standard l i n e s . analysis yields a t o t a l stress f i e l d . The Since the s o i l skeleton and the pore f l u i d s deform together when conditions are undrained, strains - i n a macroscopic sense - are the same i n each phase. Thus i n addition to Equation 3.1 (Aa') = [D»] (Ae) (Au) = K = K (Ae a v n 3.2 + Ae 2 + A e ) 2 3 3 (Ae ) v' 3.3a 3.3b i n which prime means e f f e c t i v e , [D'] i s the material property matrix (affective stress), K i s the apparent pore f l u i d bulk modulus and (Au i s the pore pressure change vector. The actual pore f l u i d modulus, K^, i s related to the apparent one as f = — a n K K where n i s the s o i l 3.4 porosity. Derivations of equations 3.3 and 3.4 are presented i n Appendix B. Equation 3.3 may be expressed i n a form compatible with equation 3.1 and 3.2 as: (A° ) = [D ] (Ae) f where (Aa ) = [Au f 3.5 f AU AU 0 0 0] T Since the pore f l u i d cannot transmit shear, [D^] can be expressed i n terms of apparent bulk modulus, r [D ] = K f a F 3 0 3 °3 0 3.6 3 where I 3 i s as 3 x 3 matrix with a l l the elements equal to 1, whereas O3 are 3 x 3 n u l l matrices. The p r i n c i p l e of e f f e c t i v e stress may be used to r e l a t e the changes i n e f f e c t i v e stress and pore presusre caused by the applied loads to the corresponding change i n t o t a l stress (Ao) = (Aa') + (Aa ) 3.7 f Substituting from Equations 3.1, 3.2, and 3.3 into 3.7, yields: [D] = [D«] + [D ] 3.8 f The e l a s t i c material i s now considered to be two phase, with the s t i f f n e s s defined by e f f e c t i v e stress moduli, E, B and pore f l u i d apparent builk modulus K^. The material properties of [D'] and are read i n separately. They are combined i n the programme automatically using equations and 3.8. 3.6 Thereafter, computation of element s t i f f n e s s , assembly and solution of displacements and hence strains proceed along the standard lines. The e f f e c t i v e stress and pore pressure are obtained by Equations 3.2 and 3.5 respectively. This approach allows stresses, porepressure and deformation i n s o i l mass, with nearly incompressible to highly compressible pore f l u i d s , to be evaluated by f i n i t e elements f o r undrained conditions. Naylor studied the response of undrained t r i a x i a l test on clay using t h i s model. The end platens were r i g i d and rough. But the computed excess pore pressures near the centre of the sample were i n good agreement with the theoretical solutions (Cam-Clay). In the case of drained analysis, Equations 2.7 to 2.10 are used instead, assuming a l l the pore f l u i d pressures are known. 3.3 CONSTITUTIVE RELATIONS 3.3.1 Incremental Non-linear Elastic Soil Skeleton An incremental non-linear e l a s t i c - p l a s t i c and isotropic s t r e s s - s t r a i n s o i l model as described by Duncan et a l (1970) i s employed i n this thesis. In this approach, two independent e l a s t i c parameters are required to represent non-linear s t r e s s - s t r a i n and volume change behaviour. the Poisson's r a t i o , v. These are usually the Young's modulus, E, and The shear modulus, G, and the bulk modulus B, are the more fundamental parameters because they separate shear or d i s t o r t i o n and volume components of s t r a i n and would be the most desirable ones to use. However, the shear modulus i s d i f f i c u l t to obtain d i r e c t l y i n laboratory testings and f o r t h i s reason Duncan et a l (1980) used the Young's modulus and the bulk modulus as their parameters. two The Young's modulus i s very s i m i l a r i n character to the shear modulus as both are a measure of d i s t o r t i o n a l response. Therefore, E and B are used i n t h i s thesis. The stress-dependent E and B are usually obtained from laboratory tests. A t y p i c a l example for sand i s shown i n Figure 3.1. The d i s t o r t i o n a l response can be reasonably approximated by modified hyperbolas (Konder) and the volumetric response i n exponential form. They are expressed as follows: The tangent Young's modulus (l-sin4.) ( a i - 0 3 ) R E where t E ± = [l- =-j — 2oj sin<j> — L = Kg P a (p-) a 2 1 J E, i 3.9 3.10 n and <J. = <(>!- A* log (—) 3.11 The tangent bulk modulus B where E i t = K B a<F-> a P m = i n i t i a l Young's modulus = Young's modulus number n = Young's modulus exponent 3 ' 1 2 24 16 Fig.3.1 - S t r e s s - S t r a i n c u r v e s f o r d r a i n e d t r i a x i a l t e s t s on l o o s e sand ( A f t e r Byrne and E l d r i d g e , 1982 ) \ = bulk modulus number m = bulk modulus exponent P = atmospheric pressure a °{> R f = major and minor p r i n c i p a l e f f e c t i v e stress = failure ratio = f r i c t i o n angle at confining stress of 1 atm •l A<{>. = decrease i n f r i c t i o n angle f o r a tenfold increase i n confining stress The procedures for evaluating these parameters from laboratory tests are described i n d e t a i l by Duncan et a l (1980) and Byrne and Eldridge (1982). 3.3.2 Incremental Non-linear Fluid Modulus 3.3.2.1 General An incremental non-linear e l a s t i c and i s o t r o p i c s t r e s s - s t r a i n f l u i d model i s employed here. Since the f l u i d cannot transmit shear, the s t r e s s - s t r a i n relations are defined only by one e l a s t i c parameter, the bulk f l u i d modulus, K^, which i s a measure of volumetric response. Before the derivation of K^, Henry's law and Boyle's law must be mentioned because these laws govern the derivations. Pore f l u i d s may be immiscible, miscible, or a combination of the two. Examples of both miscible and immiscible f l u i d s w i l l be considered i n this thesis. That i s , water undersaturated with a i r and carbon dioxide, water and bitumen saturated with gas (methane, (X^)The f i r s t combination i s to describe the pore f l u i d s of the gassy s o i l s which Sobkowicz (1982) used i n his laboratory testings. combination i s to similate o i l sand pore f l u i d s . The second If both free gas and l i q u i d s are present i n the pore f l u i d s , and the gas i s soluble i n the pore l i q u i d i n a certain extent, the pore f l u i d compressibility w i l l be both pressure dependent and influenced by the s o l u b i l i t y relationship. Hence, Boyle's and Henry's laws are appropriate f o r describing these volume and pressure relationships. 1) Boyle's law (Laidler et a l ) : The volume of a free gas i s inversely proportional to the pressure applied to i t when the temperature i s kept constant. Mathematically, V a j 3.13 where V i s the volume, P i s the absolute pressure 2) Henry's law ( L a i d l e r et a l ) : The mass of gas m dissolved by a given volume of solvent at constant temperature, i s proportional to the pressure of the gas i n equilibrium with the solution. Mathematically, m (dissolved gas) = H * P 3.14 where H i s Henry's s o l u b i l i t y constant, P i s the absolute pressure. In other words, the volume of dissolved gas i s constant i n a f i x e d volume of solvent at constant temperature when the volume i s measured at P 27 (dissolved gas) = H * V (solvent) 3.15 Most gases obey Henry's law when the temperature i s not too low and the pressure i s moderate. If several gases from a mixture of gases dissolve i n a solution, Henry's law applied to each gas independently, regardless of the pressure of the other gases present i n the mixture. H i s both temperature and pressure dependent, p a r t i c u l a r l y for natural gases i n hydrocarbons (Burcik, 1956). Since the v a r i a t i o n of H on pressure i s not very s i g n i f i c a n t , i t i s assumed that H i s independent of pressure i n t h i s thesis. 3.3.2.2 Partly Mlsclble Gas/Liquid Mixture Definitions: 1) Pore f l u i d s compressibility, 3 3.16a where V i s the volume of pore f l u i d s , P i s the absolute pressure, assuming surface tension e f f e c t s are neglected, 3.16b 28 2) Compressibility of a gas g 1 d v 8 where V a) i s the volume of gas. Air/water Mixture Let the i n i t i a l volume of free gas and water i n a s o i l element be: V* and V fg w Thus the t o t a l volume of free and dissolved gas vi = Tg fg + H V w 3.18 For a change i n pressure, the volume of water i s assumed constant as the compressibility of water i s i n s i g n i f i c a n t compared with the pore gas. Applying Boyle's law to the .,1 t o t a l volume of gas (Fredlung 1973, Sobkowicz 1982) i n the element, f i V = V Tg Tg V V i * — P P 3.19 f Then the new volume of free gas after change i n pressure f f v = V fg Tg V V f - V dg V 29 P Tg ( V P. fg + V V dg dg ) * Pf " dg 3.20 V Change of free gas i n the s o i l element fg fg fg p fg - ( fg L V By + V dg' P ,!) dg' * P dg f +AP v fg ra» 3.21 definition g i fg AP (v* + fg i fg V, ) dg' 1 # v ±+ Also, by d e f i n i t i o n 1 V * f 1 P AP d f dP V 3.22 30 , dV. V V L dV dP dP P L f 1 - S + SH P, + AP J + AP + S e VwJ - w 3.23a As AP approaches zero, 3 f P + Sg w 3.23B and b) Carbon Dioxide, Air/water Mixture By d e f i n i t i o n , the compressibility of pore f l u i d s , 3 , f 3 f 1_ " V , = - — V ^ * f*f dP dV dV f—& + — dP dP r L J i n which there may have a i r and carbon dioxide as free gas, V* . g Consider the phase diagram i n Figure 3.2, the volume of solids i s assumed to be 1 u n i t , the volume of void i s e units according to d e f i n i t i o n of void ratios Therefore V* e - e fg w V " e 3.25a V* H e + H „ e dg _ a w co2 w V 3.25a = f f Substituting equations 3.25 into 3.24a, y i e l d s , , _ 1 r !, = — I f e L ( e - e ) + H e + H _ e w a w co2 w -i „ , r g e P + AP w w As AP approaches zero, L Q J „ , 3.24b 0 / 32 •+ •w Air & C02 v 9 Water Solids F i g . 3 . 2 - Phase diagram for gassy s o i l 33 . _ 1 r 3 = _ ^ f e a L ( e - e ) + H e + H _ e w a w co2 w , . -1 + 8 e J P w w - . . 3.24c J and K 1 - 2s £ c) Gas (Methane)/Bitumen and Water Mixture Again, by d e f i n i t i o n the compressibility of pore f l u i d s 1 8 = - — f V d f * dP P f i V V d v L f rr ^ ^ dv dP x dP J (V* + V* ) p , ° P^ + AP [ g v 8 - V 0 ww - V. B. ] b b 3.26a J Consider the phase diagrams i n Figure 3.3, the volume of solids i s assumed to be 1 u n i t , the volume of voids i s e units according to d e f i n i t i o n of void r a t i o s Therefore (e - e w - e ) b 3.27a F i g . 3 . 3 - Phase diagram for o i l sand where Hg/ > ^g/b W a r e Henry's s o l u b i l i t y constants of gas i n water and bitumen respectively. Substituting Equations 3.27 into 3.26a, and AP approaches zero, y i e l d s , B = I f e [ L (e - e - e, ) + H , e 2 * -^!_w P + H e, g/b_J> , + w w b b J 3.26b and CHAPTER 4 4.1 : FINITE ELEMENT FORMULATIONS INTRODUCTION Two types of formulations are presented herein which are suitable for modelling a variety of problems encountered i n practice. Depending on the nature of the problems, they w i l l f a l l into one of the following categories: 1) Plane s t r a i n - 2 Dimensional, e.g. tunnels, shaft, e t c . 2) Axisymmetric - 3 Dimensional, e.g. t r i a x i a l test, wellbore, e t c . The s o i l i s modelled by isoparametric q u a d r i l a t e r a l or triangular elements. Stress d i s t r i b u t i o n formulations are added to cope with problems where p l a s t i c zones are developed during loading or unloading. 4.2. THE PLANE SRAIN FORMULATIONS 4.2.1 The Constitutive Matrix [P] Plane s t r a i n problems are characterized by the following two properties: 1) no d e f l e c t i o n i n z d i r e c t i o n 2) f i r s t derivative of the x and y deflections with respect to z are zero. Therefore where u, v, w are the displacements i n x, y and z d i r e c t i o n s respectively with corresponding subscripts, 1, 2 and 3 r e s p e c t i v e l y . From the generalized Hooke's law for incremental e l a s t i c i t y Aei j = [Aon - v ( A 0 2 2 Ae 2 2 = ~ [AO"22 v ^ a i l + Ac )]/E 4.2a Ao" )]/E 4.2b 3 3 + 33 - v ( A o ^ + Aa 2)]/E Ae 3 3 = [Aa Ae 1 2 = Ao /G 4.2d Ae 2 3 = Aa /G 4.2e Ae 3 1 = Aa /G 4.2f 2 33 12 23 31 After substituting the conditions 4.2, i t follows 4.2c from equation 4.1b into equations that: Aoi 3 Ao2 3 = Ao = 0 4.3a = Aa 2 = 0 4.3b 3 1 3 where Aa^.. i s the incremental shear stress with the d i r e c t i o n indicated by the subscripts. With the above eliminations, the incremental stress and s t r a i n vectors become (Aa) = [ A a u Ao"22 Ao-12] 4.4a (Ae) = [ A e n Ae 4.4b Ae ] 2 2 12 The constitutive relations for a plane s t r a i n problem i n t o t a l stress analysis are written (Naylor, 1973) as: Ao r l-v n Ao" 2 1 2 + K a 0 l-v 0 0 d-2v) 2 (l+v)(l-2v) 2 Aa V 1 1 1 0 1 1 0 ll Ae 22 0 0 0 Ae 12 A e 4.5a or (Aa) = ([D'l + [D ]) f (Ae) 4.5b where E = tangent Young's modulus v = tangent Poisson's r a t i o K cL = apparent tangent bulk f l u i d modulus The constitutive relations can also be written i n another form which i s adopted i n INCOIL. equivalent * Aa 2 2 Ao 1 2 < B' + G' B' - G' 0 B B* - G' 0 1 - G 0 0 G' < + K 1 1 0 1 1 0 0 0 0 Ae • A £ n 22 Ae 1 2 or (Aa) - ([ ']) + [D ]) (Ae) n f Where •a G' i _ 3B 2(l+v) = E 2(l+v) B = tangent bulk modulus E = tangent Young's modulus v = tangent Poisson's r a t i o K = apparent tangent bulk f l u i d modulus Equation 4.6 can also be written as (Aa) = [D] (AS) 40 where [D] i s the constitutive matrix 4.2.2 The Strain Displacement Matrix [B] For isoparametric elements, the geometry (x,y) and displacement (u,v) are both expressed by te same shape functions and are approximated as: (*) = [N] («) 4.10 O 4.11 and - [M] («') where Nj 0 N Ni 0 (6) - [x! y! x 2 y 2 x 3 y 3 (6') = u 2 v 3 u 3 v 3 W to [uj V! 2 0 N N 2 3 0 0 N U k 3 0 X l + uit 0 NJ y j - v ] 4 i n which ( 6 ) i s the nodal coordinate vector and ( 6 ' ) i s the incremental nodal displacement vector and N x = (l-s) (l-t)/4 N 2 = ( l - s ) (l+t)/4 N 3 = (1+s) (l+t)/4 - (1+s) ( l - t ) / 4 x y = nodal coordinates i n x and y d i r e c t i o n s 1, i respectively u. v. = incremental nodal displacement i n x and y 2» 3 directions respectively. s,t are l o c a l coordinates The incremental s t r a i n vector can be expressed i n terms of displacement as follows: * 3u/3x Ae Ae 1 2 13 = 3v/8y 3u/3y + 3v/3x Substitution of u and v from equation 4.11 into equation 4.12 y i e l d 3N 3Ni Ae 11 3N 2 3x~ 3x~ 3N Ae 22 Ae 0 W 3Ni 3N 3y 3x W u l v l u 2 v 2 0 3N 2 2 3N\ 3 3No 3Nq 3N. 3N 3x 3y 3x 3y~ 3N1+: 4 37" I 4.10 3~ u v 0 3x~ 3y~~ 3Nj 13 3N^ 3 3 Equation 4.10 can also be written i n matrix notation as (Ae) - [B] (6) 4.13 where [B] i s the s t r a i n displacement matrix However, the shape functions for isoparametric elements are defined with respect to the l o c a l coordinates s and t and therefore cannot be d i f f e r e n t i a t e d d i r e c t l y with respect to the global x, y axes. In order to overcome this d i f f i c u l t y i t i s necessary to obtain the derivatives of the two sets of coordinates and t h i s can be achieved through the chain rule of p a r t i a l d i f f e r e n t i a t i o n . For p l a i n s t r a i n problems, the derivatives are related as 3_ 3s 3_ 3x 3_ 3t 3_ 3y 4.14 where 3x 3s 3y 3s M- i s c a l l e d the Jacobian matrix 3y 3t 3x 3t ) Hence the derivatives w*r»t» x and y can be expressed as derivatives w»r«t s and t as follows 3_ 3x 3_ 3y -1 3s 4.15 [j] 3t The s t r a i n displacement matrix i n Equation 4.13 are evaluated numerically, using Gaussian quadrature over q u a d r i l a t e r a l regions. quadrature rules are a l l of the form // f (s,t) ds dt - n E i-1 r. I K j-1 1 K. f ( s , t.) 2 1 J The where K^, K\ are weighting functions and s^, t^ are coordinate position within the element. A 2 x 2 Gauss quadrature i s used to evaluate the s t r a i n displacement. 4.2.3 The Stiffness Matrix [K] The s t i f f n e s s matrix for the force displacement relationship i s obtained by the p r i n c i p l e of v i r t u a l work. For a v i r t u a l nodal displacement vector ( 6 ) , the external work done, W(ext), by the g external force vector ( f ) caused by v i r t u a l displacements i s written as W(ext) - {6)1 ( f ) 4.16 The v i r t u a l s t r a i n vector caused by the v i r t u a l displacements vector i s written as (^) = [B] ( 6 ) 4.17 e Hence the internal work done, W(ext), caused by the v i r t u a l strain i s written as W(int) = J where ( A e ) (Ao) t dA 1 4.18a t = thickness of the element A = area of the element Substituting equation 4.17 into 4.18a y i e l d s W(int) = J [ « ] * [ B ] [D] [B] [6] t dA T A 4.18b Applying p r i n c i p l e of v i r t u a l work W(ext) = W(int) 4.19 (f) = / 4.20 Therefore [ B ] [D] [B] [6] t dA T A (f) = [K] (6) ' 4.21 T where [K] = / [B] [D] [B] t dA i s c a l l e d the s t i f f n e s s matrix f o r the element. The global s t i f f n e s s matrix, [K ], i s obtained by assembling a l l the element s t i f f n e s s matrices together. The procedure of assembling the element matrix i s based on the requirement of 'compatibility' at the element nodes. This means that at the nodes where elements are connected, the values of the unknown nodal degrees of freedom are the same for a l l the elements joining at that node. global force The vector, [ F ] , i s assembled by adding nodal loads of each of the elements sharing the node. procedure (e.g. Displacements are calculated using standard Gaussian elimination) to solve the simultaneous equations,[K ] [s] = [F],represented by the global s t i f f n e s s matrix and the force vector. Strains can be computed from equation 4.13 or 4.36 a f t e r knowing the elements nodal displacements. After solving f o r displacements and s t r a i n s , the e f f e c t i v e stress and pore pressure can be computed from (Aa') = [D'j (Ae) and (Au) = [D ] (Ae) f 4.3 THE AXISYMMETRIC FORMULATIONS 4.3.1 The Constitutive Matrix [D| Axisymmetric problems are characterized by the following properties: 1) Symmetry of both geometry and loading 2) Stress components are independent of*the angular (0) coordinates. Hence v = 0 4.22a A e 13 Ae 2 3 " A e 31 = Ae = 4.22b 0 = 0 3 2 4.22c where u, v, w are displacements i n the r , z and 0 directions with corresponding cubscripts, 1, 2 and 3 r e s p e c t i v e l y . From the generalized Hooke's law for incremental e l a s t i c i t y = [ l o - v (Ao 2 + Ao )]/E 33 4.23a 2 2 = [ACT + Aa )]/E n 4.23b 33 " [ "33 + Aa )]/E 4.23c Ae n Ae A e u 2 2 Ao 2 - v (A033 ~ v (Aa 22 n Aei3 = Ao /G 4.23d Ae 1 2 = Aa /G 4.23e Ae 2 3 = Aa /G 4.23f 13 12 23 After substituting the conditions from Equations 4.22b and 4.22c i n t o Equations 4.23, i t follows that Aai3 = A a i = 0 4.14a Aa 4.15a 3 where A o i j 2 3 = Aa = 0 3 2 i s the incremental shear stress with the d i r e c t i o n indicated by the subscripts. With the above eliminations, the incremental stress and s t r a i n vectors become (Ao) = [ A a n Aa 2 2 Aa 3 3 Aa ] 1 2 4.25 48 (Ae) = [ A e n Ae 2 2 Ae Ae 3 3 1 2 4.26 ] The constitutive relations for an axisymmetric problem i n t o t a l stress analysis are written (Naylor, 1973) as Aa Aa 2 2 Aa 3 3 Aa 1-v 11 (l+v)(l-2v) 1 2 V V 0 1-v V 0 V 1-v 0 0 0 1 1 1 0 Ae n 1 1 1 0 Ae 2 2 1 1 1 0 0 0 0 0 [D-] + + K l-2v 4.27a Ae 1 2 Also (Air) = [ D J (Ae) where E = tangent Young's modulus v = tangent Poisson's r a t i o K = apparent tangent bulk modulus 4.27b The constitutive r e l a t i o n s can also be written i n another equivalent form which i s adopted i n INCOIL. Ao ] B'+G' B'-G' B'-G' 0 A022 ] B'-G' B'+G' B'-G 0 ^33 ] B'-G' B'-G' B'+G' 0 Ao n 1 1 1 2 + K 1 1 1 1 0 1 1 1 1 0 1 1 1 1 0 0 0 0 0 0 A e ll I Ae 2 2 33 Ae j I 4.28a 1 2 or [Ao] - ([D«] + B' = [ D ] ) [Ae] f where 3 B 2(l+v) . _ J 2(l+v) B = tangent bulk modulus E = tangent Young's modulus v = tangent Poisson's r a t i o G , 4.28b K = a apparent tangent bulk f l u i d modulus Equation 4.28 may be written i n matrix notation as (Aa) = [D] (Ae) 4.29 where [D] i s the constitutive matrix. 4.3.2 The Strain Displacement Matrix [B] The isoparametric elements, the geometry ( r , z) and displacements (u, v) are both expressed by the same shape functions and are approximated as: c •> r - W = [ N ] 4.33 ( « ' ) 0 N 0 N , 0 N [N] 4.32 (6) x 0 2 N 3 0 N N q 0 z 3 r v 3 4 0 = (6) = (>! z l r (6') = [u! v x u 0 N , 2 2 z 2 r v 2 u 3 3 k Ul+ N Zlf b ] v„] i n which (fi) i s the nodal coordinate vector, ( 6 ' ) i s the incremental nodal displacement vector and the shape functions are the same as given i n Section 4.2.2. and r nodal coordinates i n x and y directions i» i Z respectively, incremental nodal displacements i n x and y u., v. J J directions respectively. The incremental s t r a i n vector can be expressed i n terms of displacements as follows: Ae s 11 i Ae 3v 3r 2 2 4.34 u r I Ae 33 Ae 3u 3r 1 2 ) 3v 3r 3u 3z Substitution of u and v from Equation 4.33 into Equation 4.30 y i e l d s 3N 3N X Ae 11 3r 3r~ 3Ni Ae 2 2 A e 33 Ae 1 2 0 3N 2 0 3z~ Nl — r 3Nj 3Nj 3z~ IF 0 TT~ 3N N — r 3N 0 N 0 2 3z~" 3N 3 3r~ 3N 3 3z~ 0 r 3 0 3N 2 0 2 0 3 3N 3 37" w 3 52 r 0 W vi j u | 2 3z~ 0 v u 3 v 3 0 r 31^ 4.35 2 UI+ 3r~ Equation 4.35 can also be written i n matrix notation as (Ae) = [B] (6) 4.36 where [B] i s the s t r a i n displacement matrix However, the shape functions for isoparametric elements are defined with respect to the l o c a l coordinates s and t therefore cannot be d i f f e r e n t i a t e d d i r e c t l y with respect to the global x, y axes. In order to overcome this d i f f i c u l t y i t i s necessary to obtain a r e l a t i o n s h i p between the derivatives of the two sets of and this can be achieved coordinates through the chain rule of p a r t i a l differentiation. ds 3_ 3r a_ at 3_ 9z _8_ 4.37 where [J] - 3r 3s 3z 3s 3r 3t 3z_ St Hence the derivatives wvt i s called the Jacobian matrix x and y can be expressed as derivatives w»r«t s and t as follows 3_ 3r 3_ 3z [j] 3_ 3s -1 4.38 3_ 3t A similar Gauss quadrature mentioned i n Section 4.2.2. i s employed to evaluate the above [B] matrix numerically. 4.3.3 The S t i f f n e s s M a t r i x [K] The s t i f f n e s s matrix for the force displacement relatinship i obtained by the p r i n c i p l e of v i r t u a l work. displacement vector (5) > e For a v i r t u a l nodal the external work done, W(ext), by the external force vector ( F ) caused by the v i r t u a l displacements i s written as W(ext) = [s) T e (f) 4.39 The v i r t u a l s t r a i n vector caused by the v i r t u a l displacements vector i written as (Ai) = [B] ( 6 ) 4.40 e Hence the internal work done, W(int), caused by the v i r t u a l s t r a i n i s written as W(int) = / where (Ae) T (Aa) dV 4.41a V = volume of the element Substituting Equation 4.40 into 4.41a y i e l d s W(int) = J Y [&f [ B ] [D] [B] [6] dV T e 4.41b Applying the p r i n c i p l e of v i r t u a l work Wext = Wint 4.42 Therefore (0 " / [ B ] [D] [B] T v (f) = [K] (6) [6] dV 4.43 4.44 where / [B] [ D ] [ B ] dV i s c a l l e d the s t i f f n e s s matrix for the element. The procedure of obtaining element stresses and s t r a i n s i s the same as described i n Section 6.2.3. 4.4 LOAD SHEDDING ( P L A I N S T R A I N ) (Byrne 1983) Problems arise when any element within the solution domain v i o l a t e s the f a i l u r e c r i t e r i a (Mohr Coulomb). That i s , for unloading of a shaft or tunnel, a p l a s t i c zone usually developes adjacent shaft. to the The extent of p l a s t i c zone and hence volume changes are only approximated since the stress r e d i s t r i b u t i o n has not been considered during the formation of the zone. The analysis predicts the stress path ABC unloading. But the stress state at C (Figure 4.1) c r i t e r i o n (Mohr Coulomb). instead of ABD on v i o l a t e s the f a i l u r e If load shedding technique i s used, the overstress can be distributed to the adjacent elements by applying an appropriate set of nodal forces described herein and brings the stress path BC back to the correct The overstress, Ax, BD. i n the element can be removed by subtracting the computed stresses by Aa^i, A a 2 2 and A a i 2 amount as shown i n Figure 4.21b. where Aa. . have the same notation as those i n Section 4.2 Aa (Ao) = Aa 12 4.3 Aa] 11 Aa 22 and = [T] Aa- Ax 4.45 13 56 Fig.4.1 - Stresses associated with load shedding where cos26 L_ 2 cos29 2 sin29 1 4 . cos29 2 2 sin26 (T) = cos26 2 i s called the 0 transformation 0 matrix A 0 3 = minor p r i n c i p a l stress = 0 Ao"i = major p r i n c i p a l stress = 0 1 - 0 3 tan A 13 T = 2 (45 + -|-) p r i n c i p a l shear stress = 0 The derivation of these stresses changes and the transformation matrix i s i n Appendix E, The removal of these overstresses can be achieved by applying a set of nodal forces which i s obtained by the principle of v i r t u a l work. vector. The incremental nodal force vector causes a v i r t u a l displacement Hence the external work done, W(ext), can be written as W(ext) = ( s ) ^ ( A f ) e 4.46 The incremental v i r t u a l s t r a i n vector caused by the v i r t u a l displacement vector i s Ui) e = [B] ( « ) 4.47 e 58 Therefore, the i n t e r n a l work, done, W(int), i s W(int) = /. (Ae) A. T (Ac) 6 t dA 4.48a 6 Substitution of Equation 4.43 to 4.48 y i e l d s W(int) - / Applying ( 6 ) [B] T ( A C ) t dA 4.48b the p r i n c i p l e of v i r t u a l work W(ext) = W(int) Hence , [Af ] g = t/ 4.49 [ B ] [ A a ] dA 4.50 T A e where [B] i s the s t r a i n displacement matrix i n Section 4.2.2. and [ A o ] i s the stress vector shown i n Equation 4.41. e The f a i l e d element w i l l have a stress change of Ao"n> A6 2 a n d 2 Ao"i2- However, the computed stresses may not l i e on the f a i l u r e envelope due to the a p p l i c a t i o n of the nodal forces. Therefore i t e r a t i o n s may be required to bring this to the assigned tolerance. The loading shedding technique presented herein gives the same results compared with INCOIL (a' = constant). m However, the number of i t e r a t i o n s required to bring the computed stresses back to the f a i l u r e envelope i s l e s s . With the incorporation of load shedding technique, the sand skeleton i s modelled as a non-linear e l a s t i c - p l a s t i c porous material. The s o i l skeleton i s coupled with the pore f l u i d s i n the undrained model. For undrained conditions, t h i s model allows stresses, pore pressure and deformations of o i l sand masses to be evaluated by f i n i t e elements. 60 CHAPTER 5 5.1 : COMPARISONS WITH EXISTING SOLUTIONS INTRODUCTION It i s important to check the v a l i d i t y of the developed model before any major application. Two types of comparisons are presented herein which are suitable for checking the s t r e s s - s t r a i n model of the s o i l skeleton and also the newly developed gas law model. For the v a l i d a t i o n of the a n a l y t i c a l procedure, the computed r e s u l t s are compared with e l a s t i c and e l a s t i c - p l a s t i c closed form solutions. The gas law model i s validated by comparing computed solutions with observed data, such as expansion of o i l sand cores and t r i a x i a l tests on gassy s o i l s . 5.2 COMPARISONS WITH THEORETICAL RESULTS 5.2.1 Elastic Closed Form Solutions The theory for e l a s t i c closed form solution was f i r s t developed by Timoshenko (1941). The plane s t r a i n solutions of stresses and displacements i n a thick wall cylinder are presented herein: Stresses a P. - b P i o 2 a r — —r~> Q (P - P ) a b i o 2 . — + —TTV ? b^ - a a ° 2 -T\—9 (b^ - a ) r z z (P - P ) a b 2 (b 2 b - a 2 2 2 - a ) r 2 2 5.1a z P. - b P 2 2 2 5 ' l b Displacement n N (1+v)N . _ /i (l-2v) (1+v) E i P P and a2p - 0 where ( - 1 ~ b P o')r ^TZTTI a' = o P r e s s u ( i~ o a b (b - a ) r P P z r e 2 v o n ) 2 5.1c 2 z t n e inner surface of cylinder = Pressure on outer surface of cylinder J E = Young's modulus v = Poisson's r a t i o a, b, r are defined i n Figure The 5.1 response of unloading a thick wall cylinder i s being investigated. The stresses and displacements predicted by the programme are i n remarkably good agreement with the closed form solution as shown i n Figure 5.2. 5.2.2 Hence the a n a l y t i c a l procedure i s validated Elastic-Plastic Closed Form Solution For a tunnel or shaft problem, the i n i t i a l state of stress w i l l be the same throughout the domain. As the support pressure of the tunnel drops, yielding w i l l occur i f the strength of the s o i l i s exceeded. Y i e l d i n g develops on the i n s i d e face f i r s t as a p l a s t i c annular zone and extends r a d i a l l y outward i f the support pressure i s reduced further. Hence theie e x i s t s a p l a s t i c and an e l a s t i c zone i n the domain concentrically. The solutions of stresses and displacements Fig.5.1 - Thick wall cylinder 63 o o O Radii (r/r ) 0 E = 3000 MPa v = 1/3 i n i t i a l stress : a = a, = 6000 kPa final stress : o = 2500 kPa inside radius : r„ = 1 m r r Fig.5.2 - Stresses and displacements around c i r c u l a r in an e l a s t i c material opening o o o 64 - Tl - J - J 2o •Q— 00 CO 0 0 0 j_> o > closed p o o o form program 0) Ed 1 1 1 _l 4 v- 1 6 Radii (r/r ) 1 i 8 10 0 = 3000 kPa E 1/3 40° i n i t i a l stress final stress inside radius V 0 1 L B e o a r r r 0 •= o = 6000 kPa = 1500 kPa = 1m e F i g . 5 . 3 - S t r e s s e s and d i s p l a c e m e n t s around c i r c u l a r elastic-plastic material opening i n an for the p l a s t i c and e l a s t i c zone w i l l be quite d i f f e r e n t . The stresses and displacements for the e l a s t i c zone are just an extension of Equation 5.1-by setting b to i n f i n i t y . They may be written as Stresses r = a P o + (P, -^V i - Po ) r" 5.2a 9 = a P - (P, - P ) —jr o i o' r 5.2b 2 a v / 1 2 a a v z Displacement E where r v i o P^ = Pressure on inner surface of cylinder P = Pressure on outer surface of cylinder E = Young's Modulus v = Poison's Ratio o J Different investigators, Gibson and Anderson (1951), Ladanyi (1963), Vesic (1972) and Hughes et a l developed closed form solutions of stresses and displacements for the p l a s t i c zone. Hughes et a l presented the more elegant solutions which are presented herein: A l l the sand i s assumed to f a i l with a constant ratio of p r i n c i p a l stresses, so that N = tan 2 (45 + |) 5.3 The equilibrium equation that must be s a t i s f i e d i s do a* - al -^+^ i-0 dr r 5.4 Substituting for a' from Equation 5.3, boundary conditions of a^ = integrating and using outer at r = R. a' in = (1-N) Jin (|) 5.5 R where a' r = r a d i a l stress at r within the p l a s t i c zone o = r a d i a l stress at the outer boundary, R, of the R r p l a s t i c zone = P (1 - sin<j>) p R - a [-2. ( l - s i n * ) ] I 5.6 l-sinfr 2 s i n * 5.7 Equation 5.5 governs the d i s t r i b u t i o n of the r a d i a l effective stress within the p l a s t i c zone. Continuity of stresses and displacements between the p l a s t i c and e l a s t i c zone must be maintained. into Equation 5.2c, contact i s Hence Equation 5.6 i s substituted the displacement, U at the e l a s t i c - p l a s t i c zone U R E" o R P S i n 5.8 * Hughes et a l also show that the displacements, u, within the p l a s t i c zone ^=(|) where j n + 1 n = tan 2 (^|) 5.9 (45 + u = d i l a t i o n angle The response of unloading a tunnel i s investigated. comparisons The i n Figure 5.3 show that the analysis of stresses and displacements i n e l a s t i c - p l a s t i c materials predicted by the programme i s i n good agreement with the closed form solutions. The minor discrepancies are due to the l i m i t of Poisson's r a t i o and the coarseness of the mesh. An upper l i m i t of 0.499 i s adopted i n the programme (MHANS) to maintain numerical s t a b i l i t y , whereas the actual value should be 0.5. The agreement i n displacements and the extent of the p l a s t i c zone also confirm that the load shedding technique i n the programme (MHANS) i s working properly. This i s a s a t i s f a c t o r y check o the programme i n drained analysis. Unfortunately, the load shedding technique i n INCOIL cannot b successfully tested. The element type i n this programme i s QM-6. Non-equilibrium of stresses a r i s e i n QM-6 elements at high Poisson's r a t i o (v > 0.4) i f the geometry of the elements i s non-rectangular. Since the elements i n a f i n i t e element mesh f o r modelling plane s t r a i n shaft problems are not rectangular, non-equilibrium of stresses arise before load shedding i s required. 5.3 COMPARISONS WITH OBSERVED DATA 5.3.1 One Dimenional Unloading of Oilsand This i s an opportunity f o r the new with some f i e l d data. gas law model to be checked Since the unloading i s 1-D, the v a l i d i t y of the schematic spring analogy model (Figure 1.1) can also be demonstrated. Unconfined oilsand core taken from d r i l l e d holes swells by 5 to 15% of the o r i g i n a l volume (Dusseault 1980, Byrne et a l 1980). Maximum expansion potential generally cannot be reached because expansion stops when there i s adequate intercommunication to permit flow of gas out of the sample. the s o i l f a b r i c . of gas voids This leads to disruption of This i s the equilibrium saturation point i n which gas becomes mobile and the maximum gas saturation value i s 15% (Amyx et a l 1960). However, the t o t a l amount of expansion i s impossible to predict and can only be measured for i n d i v i d u a l cores. The core l i n e r s are s p e c i f i c a l l y designed oversized to prevent the jamming of the core within the b a r r e l . Radial expansion i s assumed to be completed when the o i l sand core i s brought up from the d r i l l e d hole. The l i n e r s and steel containers are assumed to be r i g i d and f r i c t i o n l e s s so that only a x i a l expansion of the core i s allowed. Therefore, the oilsand core can be modelled as 1-D unloading after recovery. The i n i t i a l l y high stressed core specimen i s modelled as shown i n Figure 5.4a. V e r t i c a l stress i s reduced from 750 kPa to zero. stresses, pore pressure and displacements are shown i n Figure 5.4b. The It can be seen that the change i n t o t a l stress i s apportioned between the e f f e c t i v e stess and the pore pressure as suggested by the spring analogy model. Case A: Gas saturation pressure = 100 kPa. I n i t i a l l y pore pressure i s above gas saturation pressure and saturation remains 100%. Therefore loads come off from the pore f l u i d s while e f f e c t i v e stress remains f a i r l y constant because pore f l u i d i s the s t i f f e r phase as shown i n Figure 5.4b. As pore pressure drops below the gas saturation pressure, gas s t a r t s to evolve which causes the pore f l u i d to become f l e x i b l e . Stress change w i l l be taken up by the s o i l skeleton on further unloading u n t i l zero e f f e c t i v e s t r e s s . Further reduction on boundary load at zero e f f e c t i v e stress i s e n t i r e l y accommodated by the pore f l u i d s . Case Bt Gas Saturation pressure = 300 kPa. Because of gas exsolution the pore f l u i d s s t a r t as a f l e x i b l e phase so e f f e c t i v e stress drops to zero with no appreciable change of pore pressure on unloading. Reduction of boundary stress beyond this stage i s e n t i r e l y taken by the f l u i d phases as the s o i l skeleton e s s e n t i a l l y has no s t i f f n e s s at zero e f f e c t i v e stress as shown i n Figure 5.4b. It may be seen from figure 5.4b that there are no appreciable displacements when the e f f e c t i v e stress i s p o s i t i v e . The displacements e s s e n t i a l l y come from unloading at zero e f f e c t i v e stress. The t o t a l displacements upon t o t a l removal of v e r t i c a l stress l i e within the range of 5 to 15% of the o r i g i n a l length. 70 unload to 450 kPa 300 kPa •Im s • = 1565 n -= 1000 m - 65° ' 0.9 =• 0.243 = = A<f> = 0.5 0.25 o 22 f n o n == 0.045 w = 0.2 H o = 0.02 H w = 100% S R 5.4 - Model f o r one dimensional unloading of o i l sand - u = 100 kPa u = 300 kPa o bb w bb CO c '5o u _ \ j 150 g.5.4 300 L 450 i i 600 i External Load (kPa) i 750 900 - Response of o i l sand to one dimensional unloading 72 The r e s u l t s show that the new stress s t r a i n model has the excellent c a p a b i l i t y of predicting undrained response of o i l sand. 5.3.2 T r i a x i a l Tests on Gassy S o i l s (Sobkowicz 1982) Performance of gassy s o i l on laboratory t r i a x i a l tests are reported by Sobkowicz (1982). A review of his work shows that the gas law model i s generally correct. An appraisal of the predictive c a p a b i l i t i e s of the gas law model i s made by comparing predicted and laboratory observed response of the immediate pore pressure, the immediate (short term) B value, the equilibrium pore presusre, the equilibrium (long term) B value, and value of saturation and displacements. Since only test 11 i n Sobkowicz's thesis i s documented i n d e t a i l , a comparison between predicted and observed response for t h i s test i s made to evaluate the v a l i d i t y of the undrained model. Analysis i s performed using the same i n i t i a l condition as test 11: S = 99.75%, n = 32.28%, a = 1403.3 kPa, u = 652.3 kPa and 3 = 9E-6 k P a . s - 1 The compressibilkity of s o l i d i s comparable with that of water (3 = 5 x 10 -7 k P a ) so that the B (short term) value w i l l not equal to one even -1 for f u l l saturation. These components are converted into parameters to be read i n by the programme. The conversion and parameters are shown i n Appendix F. The unloading sequence i s shown i n Table 5.1. During any phase of the i s o t r o p i c unloading test, pore pressure responses are predicted from the knowledge of s o i l skeleton and f l u i d compressibilities which are a function or e f f e c t i v e stress, pore pressure, saturation and porosity. set For short term reponse, H i s to zero i n Equation 3.26b since there i s no time for gas exsolution. H . , ^ = 0.02 and H „ , ^ = 0.86 are used for the air/water co2/water equilibrium response when gas exsolution i s complete. The comparisons are summarized i n Table 5.1 and are presented graphically i n Figures 5.5 and 5.6. 1) They include: predicted undrained response by the present undrained model 2) measured undrained response (Test 11, Sobkowicz 1982) A careful examination of Table 5.1 and Figures 5.5 and 5.6 show that the predictive c a p a b i l i t y of the gas law model i s remarkably good, especially for the long term undrained response. The minor discrepancies are due to the loss of gas from the sample as the r e s u l t of gas d i f f u s i o n and leakage through the membrane. The observed immediate pore pressure are higher than the predicted values because of the time elapse (15 to 30 seconds) between reducing the t o t a l stress and taking the f i r s t reading. Thus, the predicted short term B i s always higher than the observed ones. It can be seen that the stress reduction i s apportioned between s o i l skeleton and pore f l u i d s , depending on their compressibilities. The sample i n Test 11 was i n i t i a l l y saturated with respect to a i r i n water and undersaturated with respect to carbon dioxide i n water. P On unloading, during the f i r s t few phases, as . < P < P . , _ , a small amount of gas exsolves for H co2/water air/water' 0 / (air/water) = 0.02. & The change of e f f e c t i v e stress and pore pressure TABLE 5.1a A Comparison of Computed and Measured Results (Test 11, Sobkowicz) Total Short Term B S a l w r a "t.» ori ( %,) Long Term B Phase Stress (kPa) Predicted A 1322.4 0.897 0.694 0.523 0.606 99.65 99.67 B 1220.5 9.843 0.69 0.477 0.433 99.50 99^54 C 1112.1 0.781 0.68 0.378 0.403 99.30 99.38 D 978.2 0.705 0.64 0.021 0.011 98.90 99.11 E 883.6 0.628 0.548 0.022 0.016 98.61 98.93 F 766.4 0.584 0.508 0.024 0.034 98.22 98.93 G 654.9 0.548 0.482 0.026 0.07 97.80 98.27 H 559.0 0.536 0.50 0.031 0.155 97.38 97.90 Jl 457.3 0.569 0.590 0.129 0.21 96.40 J2 1 1 1 1 95.15 J3 1 1 1 1 94.66 Measured Predicted Predicted: Results predicted by Programme Measured: Results measured i n Test 11 by Sobkowicz Measured Predicted Measured TABLE 5.1b A Comparison of Computed and Measured Results (Test 11, Sobkowicz) Measured Strains (%) Phase Horizontal Vertical Porosity (%) Volumeteric Predicted Measured A 0.112 E - l 0.112 E - l 0.336 E - l 32.30 32.30 B 0.275 E - l 0.275 E - l 0.825 E - l 32.33 32.33 C 0.490 E - l 0.490 E - l 0.147 32.38 32.36 D 0.929 E - l 0.919 E - l 0.2757 32.46 32.42 E 0.124 0.124 0.372 32.53 32.46 F 0.167 0.167 0.501 32.62 32.53 G 0.214 0.214 0.624 32.72 32.61 H 0.261 0.261 0.783 32.80 32.69 Jl 0.376 0.376 1.128 J2 J3 .5.5 - Comparisons of predicted and observed pore pressure 77 are roughly the same because the c o m p r e s s i b i l i t i e s of both s o i l skeleton and pore f l u i d s are comparable. s i m i l a r to those of unsaturated s o i l s . P This c h a r a c t e r i s t i c i s On further unloading, as P < „ , _ , a large amount of gas exsolves because the high s o l u b i l i t y co2/water' ° ° a J (H „ = 0.86) coz of carbon dioxide i n water. This causes a sudden increase i n f l e x i b i l i t y of the f l u i d phase and hence most of the load i s transferred to the s o i l skeleton. When the e f f e c t i v e stress i n the skeleton approaches zero, the f l u i d once again becomes the s t i f f e r phase, hence the B value rises to one. This i s the t y p i c a l behaviour of gassy s o i l on unloading. The predicted and measured displacements are i n remarkably good agreement. This indicates that the input parameters (Appendix f ) and the r a t i o of the parameters, and Cheung) are generally correct. = 0.6 K^, and n = 2m = 0.5 (Byrne CHAPTER 6 - STRESSES AROUND A WELLBORE OR SHAFT IN OIL SAND 6.1 INTRODUCTION The response of a wellbore i n o i l sand upon unloading i s considered because i t i s an important problem i n o i l recovery i n o i l sand. In general, knoweldge of the stress solutions around a borehole i s of great importance i n several s i t u a t i o n s : 1) borehole s t a b i l i t y 2) hydraulic f r a c t u r i n g 3) production or i n j e c t i o n A theoretical solution for stresses around a wellbore was developed by Risnes et a l (1982) and the equations are presented herein. Validation of the programme (MHANS) for drained analysis i s made by comparing computed response with the closed form solutions developed by Risnes et a l . A linear e l a s t i c - p l a s t i c constitutive r e l a t i o n s h i p i s used f o r the above v a l i d a t i o n . Upon v a l i d a t i o n of the programme, i t was used to study the behaviour of a wellbore i n o i l sand upon unloading. Undrained and drained analyses were performed to obtain the short term and long term response respectively. In the undrained analysis, the gas exsolution i s assumed to be very fast r e l a t i v e to the construction of wellbore. In the drained analysis, the pore pressure p r o f i l e i s estimated by using Dupuit's theory (Section 6.3.3). 6.2 GENERAL MODEL DESCRIPTION The wellbore under consideration i s supported by f l u i d pressure. As a model, a v e r t i c a l c y l i n d r i c a l hole through a horizontal layer of o i l sand i s considered. The geometry, f i n i t e element mesh, and i n i t i a l conditions of the problem are shown i n Figure 6.1 and 6.2. Loading and geometry are assumed to be symmetrical around the well axis. Only r a d i a l displacement a f t e r the i n i t i a l overburden loading are considered. These correspond to the assumption of axisymmetric and plane s t r a i n conditions. The sand formation homogeneous and i n i t i a l l y i s assumed permeable, i s o t r o p i c , f u l l y saturated. The material i s assumed e l a s t i c - p e r f e c t l y p l a s t i c and obeys Mohr Coulomb f a i l u r e c r i t e r i o n . Only stress solutions for > at the e l a s t i c - p l a s t i c boundary w i l l be investigated. 6.3 THEORETICAL SOLUTIONS FOR STRESSES AROUND A BOREHOLE A closed form solution for streses around a well, using a l i n e a r e l a s t i c - p l a s t i c s t r e s s - s t r a i n r e l a t i o n s h i p , can be obtained from Risnes et a l (1982). The derivations of the stress solutions follow that of Risnes et a l (1982) with two additional assumptions. 1) Insitu state of stress i s considered initially, 2) to be i s o t r o p i c i . e . a = o. = a . r o z The Mohr Coulomb f a i l u r e c r i t e r i o n i n a porous material is f = a{ - 2S tana - 0 3 tan a. = 0 2 6.1 Fig.6.1 - Outline of the problem 76 elements , 78 nodes Fig. 6.2 F i n i t e element mesh for wellbore problem where S i s the cohesion intercept (or apparent cohesion) a i s the f a i l u r e angle, i . e . + -|— <f>' i s the internal f r i c t i o n angle These symbols are also explained i n Figure 6.3.1 6.2b Stresses 6.3.1.1 Stresses In Elastic Zone The stresses around a hole i n an e l a s t i c thick wall cylinder, with porous material saturated with f l u i d , may be written as follows: R2, 'r = °ro < ro " r i > F ^ R T o + 0 o 0 [l i £n(R . -°> » "i ' *n(R /R.)o R^* ( P ^ o 2 i ' 2(l-v) R? °9 * °ro (^) ] + ( °ro " °ri> R ^ R T o _ P o [*n (^) ) i ; R i t 1 + 1-2* 2(l-v) i v - 1]) o i ' 6.3 Fig.6.2 b - Mohr Coulomb f a i l u r e envelope 6.4 The procedure for obtaining these stress solutions i s given i n Appendix G. 6.3*1.2 Stresses i n Plastic Zone As long as f = 0 (Equation 6.1), a p l a s t i c zone w i l l start to develop at the borehole wall, and then expanding i n size as support pressure i s decreased. Equation 6.1 w i l l apply within the p l a s t i c zone. If the stress state at the boundary between the e l a s t i c and p l a s t i c zone i s considered, the e l a s t i c stress solutions at this boundary are given by Equation 6.2 to 6.3, with R = r = R., a =a o i rc r and P = P.. With the assumption of no f l u i d flow ( i . e . P = P ) and c i c o' , r R Q » R^, the stress solutions from Equations 6.2, 6.3 and 6.4 may be written as: a rc = a = 2a a The zc = a 6.4 rc ro zo - a rc 6.5 6.6 e l a s t i c solutions (Figure 5.2) show c l e a r l y that the r a d i a l stress w i l l be the smallest at the boundary between e l a s t i c and p l a s t i c zones, and following the assumption (1) of i n i t i a l i s o t r o p i c stress, the state of stress at the e l a s t i c - p l a s t i c boundary w i l l be a < a < a . rc zc 9c n The stress solutions within the p l a s t i c zone may be derived by combining Coulomb f a i l u r e c r i t e r i o n (6.1) and the equation of equilibrium. do -* dr a -a + JL-± r =0 6 .7 The stress solutions for the p l a s t i c zone may be written as: For R. < r < R i c a r = P. + ir^r; Zn | - + ^- (2S tana - -^r-) i 2irhk R^ t 2Trhk [(I-)* - 1] 6.8 For R. < r < R i c °9 - p i + 2¥ht c * ir> F < 1+ n [(t + 1) (I-) - l] 11 R For R, < r < R. i + 2S t a n a - iSk> 6.9 87 °z " i 2^hk < * R7> I P + +X N + ( ° - zSfc t 2S tan ( t + L) <!/ 6.10 For R, < r < R b c a - (P + z «a- i £ n f-) + v R ' 2Trhk l i l - + Cl-y)(l-2v) 2irhk 1-v 6.11 <|)1 where t = tan a - 1, a = 45° + -j2 u = fluid viscocity The procedure to derive Equations 6.7 to 6.11 i s given i n Appendix H. 6.3.1.3 Radius of the P l a s t i c Zone Radius of the inner P l a s t i c Zone R, b At the boundary of inner and outer p l a s t i c zones, the tangential stress and v e r t i c a l stress given by equations 6.9 and 6.11 are equal. Setting Equations 6.9 equal to 6.11 and r = R^, y i e l d s . a x (^V + a i 2 = 0 6.12 88 where - £ (2S tana - a i a 2 * = ( 1 7 . v [(t + 1) - v (t + 2)] JiJL - (1^) 2irhk l - v ) ( 2 S t a n a ' 2$!k> ( 1 " zo (ff v P ) cr 2 v ) Radius of the E n t i r e P l a s t i c Zone There are two requirements that must be s a t i s f i e d at the boundary of e l a s t i c and p l a s t i c zones 1) Mohr Coulomb c r i t e r i o n must hold 2) Continuity of r a d i a l stress Inserting r a d i a l stress from Equation 6.8 and tangential stress from 6.3 into Mohr Coulomb f a i l u r e c r i t e r i o n 6.1, the r e s u l t i n g equation f o r the radius of p l a s t i c zone R^ i s b R x + b 6 C c 2 R £n ^ R + b 2 R R in Y~ • in jp- + b c i £ 7 n -2. + R in ^ R b + i 3 R b Q + 2 £ n b l + R 2 R _° + c £ b g n = |_ R + ,. 2 b R £ n _ | 0 6.13 89 where bx = (2S tan a - b • B - - ^r 1 2 Ci R o t " 1 2 2irhk b - - b^ R 6.3.2 9 = - C < + P + Q V + - 2 p i 1 a 2 5 tan uq t+2 i 2TThk t b 2 v = V 2uhk b5 8 o - ( P v 2 (1-v) b 2 1 3 * " 2(1=*) = R ^ J 3 R 2 Stability It i s noted that the r a d i a l stress component i n Equation 6.8 consists of two r-dependent terms, one logarithmic and one to the power 90 of t . The l a s t term w i l l become dominant when the exponent t has a value greater than about two. Setting ° C = 7-t (2S tan a - ^rtjr) Ri~ 2irhk 6.14 t If C i s positive, r a d i a l stress i n the p l a s t i c zone w i l l increase with r, and combined p l a s t i c - e l a s t i c solutions are possible. But when the flow rate q i s large enough to cause C to become negative, r a d i a l stress w i l l decrease with increasing distance r , and combined solutions are not possible. Hence, there exists a s t a b i l i t y c r i t e r i o n . C>0 6.15 with the l i m i t •X—T- 2irhk = 2S tan a 6.16 This study concentrates on o i l sand which has S = 0. If the wellbore i s supprted only by f l u i d pressure, equation 6.16 indicates that i n s t a b i l i t y arises when flow into the wellbore occurs. 6.3.3 Pore Pressure P r o f i l e When steady-state conditions around the wellbore have been reached, the pore pressure i n the s o i l elements may be estimated i f the piezometric surface i s known. Dupuit developed a theory which enables the quantity of steady-state seepage and the piezometric surface around a well to be evaluated. His theory i s based on three assumptions: 91 1) the hydraulic gradient i s equal to the slope of the free surface and i s constant with depth, 2) for small i n c l i n a t i o n s of the l i n e of seepage the streamlines may 3) be taken as h o r i z o n t a l . The permeability of the s o i l Is constant With the terminology i n Figure 6/2a. the flow when steady state conditions e x i s t i s given by: Q - w k (h J^J^ - h ) 2 2 2 L 6.17 w and the location of the free surface i s h - h! + . ,—rJin (—) £n(R/r ) r ' w w 2 2 2 h 2 = hi 1 2 / T i v 6.18 There Is a controversy about the location of the piezometric surface predicted by Dupuit's theory, especially i n the v i c i n i t y of the well. This i s because the surface of seepage i s omitted i n Dupuit's prediction. But the problem i s modelled as a disk of sand below the seepage surface (Figure 6.3b). disk. Hence equation 6.16 pressure profile. Only r a d i a l flow i s assumed i n the sand i s accurate enough to estimate the pore h, ^^— f Flow surface of seepage 7"~ F.F. M e s h Fig.6. 3 - Idealised flow to a wellbore 93 6.4. Comparisons of Predicted Response and Closed Form Solution The response of unloading a borehole, with l i n e a r e l a s t i c - p l a s t i c porous material, i s investigated. A drained analysis was performed. The i n i t i a l and f i n a l conditions of the problem are shown i n Figure 6.4(a). When the f l u i d support pressure i s higher than the i n i t i a l pore pressure, no flow from the borehole into the sand formation i s assumed. Since the f i n a l f l u i d support (4100 kPa) i s higher than the i n i t i a l pore pressure, flow into the borehole i s not considered herein. The stress solutions computed by the programme are i n good agreement with the closed form solutions mentioned i n Section 6.3, and shown i n Figure 6.4. The radius of the entire p l a s t i c zone i s small which shows that the borehole i s stable at the f i n a l f l u i d support pressure of 4100 kPa. 6.5.1 Undrained Response The undrained non-linear e l a s t i c - p l a s t i c model w i l l now be used to study the response of a wellbore on unloading. The f i n i t e element mesh and i n i t i a l conditions are shown i n Figure 6.1b. Only the long term undrained response w i l l be investigated i n t h i s thesis because this condition i s f e l t to be more r e a l i s t i c (t * 0) and i t was also shown that the long term undrained condition i s more c r i t i c a l than the immediate response (Sobkowicz, 1982) i n terms of s t a b i l i t y . The wellbore i s unloaded by decreasing the t o t a l stress at the wellbore wall, and the stress solutions and displacements are shown graphically i n Figures 6.5. A careful examination of these figures indicates some interesting r e s u l t s : 94 W <=> 10 20 30 40 50 Radii (r/r ) 0 E = 60 MPa v = 0.45 + = 35° i n i t i a l stress final stress inside radius outside radius Or O = 0 = Oz = = 600 k P a lb = R = r 8 4000 kPa 0.125 m 6.7 m (a) Fig.6.4-- Stresses around a wellbore in an e l a s t i c - p l a s t i c material 95 o o O 00 >< o o CD CL CO CO CD o o u CO o o (D CVJ > o Cd ° o o closed form program -oCO > o o CM O 0) V—I w ° J L 10 J 20 30 Radii (r/r ) 40 I 50 0 (b) Fig.6.4 - Stresses around a wellbore in an e l a s t i c - p l a s t i c material Fig.6.5 - Undrained response of a wellbore in o i l sand on unloading 97 Fig.6.5 - Undrained r e s p o n s e of a w e l l b o r e i n o i l sand on u n l o a d i n g 98 Fig.6.5 - Undrained response of a wellbore in o i l sand on unloading 99 1) The support pressure can be reduced below the i n i t i a l pore f l u i d pressure and the wellbore Is s t i l l stable; i n s t a b i l i t y i s defined when large displacements start to occur at the wellbore wall. 2) I n s t a b i l i t y occurs at a support pressure of approximately 2500 kPa. 3) The size of the p l a s t i c zone remains small as long as the support pressure i s higher than 2500 kPa (Figure 6.5 b). Once the support pressure drops below 2500 kPa, the size of the p l a s t i c zone increases rapidly (Figure 6.5 c ) . 4) Pore f l u i d pressure changes only occur i n the p l a s t i c zone. The evaluation of the f l u i d pressure response depends on volumetric s t r a i n Ae e l a s t i c zone, Ae e r v = Ae e P + Ae . v v But i n the = - A e^ and Ae = 0 so that Ae = 0, 6 z v e e and hence no change i n pore f l u i d pressure i s predicted. 5) Once i n s t a b i l i t y has been reached, a l i q u i d zone with zero e f f e c t i v e stress associated with a large p l a s t i c zone w i l l form adjacent to the wellbore. These zones w i l l extend into the sand formation rapidly upon further reduction of support pressure, leading to large displacements. 100 6.5.2 Drained Response For the drained condition, the pore f l u i d pressure i s assumed to be known. Dupuit's theory described i n Section 6.4.3 i s adopted to estimate the pore pressure p r o f i l e around the borehole i n this study. Two t y p i c a l pore pressure p r o f i l e s are shown i n Figure 6.7 with R = 100 and 150 m with the f l u i d support pressure f l u i d at 3200 kPa. There are some intermediate pore pressure p r o f i l e s between support pressure of 3500 kPa to 3200 kPa, depending on the number of increments on unloading, but they are not shown here. When the f l u i d support pressure i s above 3500 kPa, no flow from the wellbore into the sand formation i s assumed. The wellbore i s unloaded i n the same manner as f o r the undrained analysis. Figure 6.7. The results of the stress solutons are shown i n A careful examination of these r e s u l t s indicate some interesting points: 1) To maintain borehole s t a b i l i t y , the support pressure cannot be reduced to less than the i n i t i a l pore pressure; i n s t a b i l i t y i s defined as large displacements start to occur at the wellbore w a l l . 2) The stress solutions only d i f f e r by a few percent when the input pore pressure p r o f i l e s are generated by using R = 100 and 150 m. Therefore, only one set of stress solutons i s presented here i n Figure 3) 6.7. Once s t a b i l i t y has been reached, a l i q u i d zone with zero e f f e c t i v e stress associated with a large p l a s t i c zone w i l l extend into the sand formation rapidly upon further F i g . 6 . 6 - Pore pressure p r o f i l e around a borehole 102 200 300 400 500 600 700 Support Pressure (kPa) (X10 ) 1 (a) Fig.6.7 - Drained response of a wellbore in o i l sand on unloading 104 reduction of support pressure, leading to large displacements. 6.5.3 Implications of Undrained and Drained Analyses The analyses show that there are l i m i t s on the f l u i d support pressure reduction i n order to maintain borehole s t a b i l i t y . Comparisons of both analysis are made at support pressure of 3500 kPa and 3200 kPa. 3500 kPa i s the c r i t i c a l pressure below which i n s t a b i l i t y occurs i n drained condition. It i s noted that i n Figure 6.9 the p l a s t i c zone i n the undrained analysis i s much smaller than the one i n drained analysis. At a support pressure of 3200 kPa, the borehole i s obviously unstable under drained conditions (Figure 6.9), showing a large p l a s t i c zone and consequently large displacements. But the borehole only exhibits a small p l a s t i c zone at t h i s support pressure (3200 kPa) under undrained conditions (Figure 6.9). This i s because the pore pressure around the wellbore i s lower i n the undrained case, which results i n a higher e f f e c t i v e stress. Consolidation i s the process which bridges the f u l l y and drained conditions. undrained An interesting point i s that the pore pressure w i l l increase around the borehole during consolidation, leading to lower e f f e c t i v e stress. Hence, the long term drained condition i s less stable than the undrained condition. 6.6 Application to O i l Recovery For o i l production, the f i n a l support pressure must be reduced below the i n - s i t u pore f l u i d pressure. Based on the undrained with the wellbore supported only by f l u i d pressure, the support analyses g.6.8 - Comparisons of undrained and drained response of a wellbore in o i l sand at a support pressure of 3500kPa 106 Fig.6.9 - Comparisons of undrained and drained response of a wellbore in o i l sand at a support pressure of 3200kPa 107 pressure can be reduced below the i n - s i t u pore f l u i d pressure and the wellbore i s s t i l l stable. This allows the construction of the wellbore and i n i t i a l reduction of f l u i d support pressure below i n - s i t u pore f l u i d pressure. However, for the drained condition, the wellbore becomes unstable which causes collapse of the well and hence no o i l production. I n s t a b i l i t y r e s u l t s i n the formation of large l i q u i d and p l a s t i c zones (Figure 6.9) around the wellbore. Since the permeability i n the l i q u i d and p l a s t i c zones are higher due to the expansion of sand skeleton, i t i s desirable to have these zones around the o i l production well. This e f f e c t i v e l y increases the diameter of the w e l l . To enhance o i l production and maintain wellbore s t a b i l i t y , a screen may be i n s t a l l e d to provide e f f e c t i v e support pressure a f t e r the l i q u i d and p l a s t i c zones have formed. The three dimensional and viscous e f f e c t s have not been considered i n the analysis, however they may wellbore. help i n s t a b i l i z i n g the 108 CHAPTER 7 : SUMMARY AND CONCLUSIONS A new stress-stress r e l a t i o n s h i p f o r modelling the undrained response of o i l sand has been presented. s o i l skeleton and pore f l u i d s i s used. An analysis which couples the The pore pressure changes are computed from the constraint of volume compatibility. Separate s t r e s s - s t r a i n models are required f o r both s o i l skeleton and pore f l u i d s i n this analysis. The conventional hyperbolic s t r e s s - s t r a i n model described by Duncan et a l i s adopted f o r the s o i l skeleton. The pore f l u i d s s t r e s s - s t r a i n relationships are formulated on the basis of i d e a l gas laws. The developed model i s incorporated into a f i n i t e element programme for analysing the deformation behaviour of gassy s o i l s (e.g. o i l sand). Upon v a l i d a t i o n i n Chapter 5, i t i s shown that the undrained model i s capable of predicting the response of unsaturated to gassy s o i l s . Naylor has shown that t h i s model can be used to predict the response of saturated s o i l s . For non-rectangular QM-6 elements, equilibrium cannot be achieved for Poisson's r a t i o values greater than 0.4, but higher order element can remedy t h i s . In the study of the response of wellbore i n o i l sand upon unloading, the f l u i d support pressure can be reduced below i n - s i t u pore f l u i d pressure under undrained condition and the wellbore i s s t i l l stable. However, for the drained conditin, the f l u i d support pressure cannot be reduced below the i n - s i t u pore f l u i d pressure i n order to maintain wellbore s t a b i l i t y . For o i l production, the f l u i d support must be reduced below i n - s i t u pore pressure which results i n formation of large l i q u i d and 109 p l a s t i c zones. It i s desirable to have these zones around the wellbore because the diameter of the well i s e f f e c t i v e l y larger. To enhance o i l production and maintain wellbore s t a b i l i t y , a screen may be i n s t a l l e d to provide e f f e c t i v e l i q u i d and p l a s t i c zones have formed. support pressure after the BIBLIOGRAPHY Atukorala, U.D. F i n i t e Element Analysis of Fluid Induced Fracture Behaviour i n Oilsand. Biot, M.A. M.A.Sc. Thesis, U.B.C. 1983. General Theory of Three-Dimensional Consodilation. Journal of Applied Physics, V o l . 12, February 1941. Bishop, A.W. The Influence of an Undrained Change i n Stress on the Pore Pressure i n Porous media of Low Compressibility. Geotechnique, V.23, N.3. 1973. Bishop, A.W. The Influence of System Compressibility on the Observed Pore Pressure Response to an Undrained Change i n Stress i n Saturated Rock. Brooker, E.W. Geotechnique, V.26, 1976. Tar Sand Mechanics and Slopes Evaluation. 10th Canadian Rock Mechanics Symposium, Department of Mining Engineering, Queen's University, 1975. Byrne, P.M. and Cheung, H. Sand Masses. S o i l Parameters f o r Deformation Analysis of S o i l Mechanics Series No. 81, Department of C i v i l Engineering, U.B.C. 1984. Byrne, P.M. and Eldridge, T. A Three parameter Dilatant E l a s t i c Stress-Strain Model f o r Sand. S o i l Mechanics Series No. 57, Department of C i v i l Engineering, U.B.C. 1982. Byrne, P.M. and Janzen. W. Soilstress. S o i l Mechanics Series No. 52, Department of C i v i l Engineering, U.B.C. 1981. Byrne, P.M. and Janzen, W. INCOIL. S o i l Mechanics Series No. 80, Department of C i v i l Engineering, U.B.C. 1984. C h a t t e r j i , P.K., Smith, L.B., Insley, A.E. and Sharma, L. of Saline Creek Tunnel i n Athabasca O i l Sand. Geotechnical Journal, 1, 197 9. Construction Canadian Christian, J.T. Undrained Stress D i s t r i b u t i o n by Numerical Methods. J.S.M.F.D., A.S.C.E., S.M.6, 1968. Cook, R.D. Concepts and Applications of F i n i t e Element Analysis. 2nd Edition, John Wiley & Sons, 1981. Duncan, J.M. and Chang, C.Y. Nonlinear Analysis of Stress and S t r a i n in S o i l s . J.S.M.F.D., A.S.C.E., S.M.5, 1970. Duncan, J.M., Byrne, P.M., Wong, K.S. and Mabry, P. Strength, Stress- Strain and Bulk Modulus Parameters for F i n i t e Element Analysis of STresses and Movements i n S o i l Masses. Report No. UCB/GT/80-01, August 1980. Dusseault, M.B. Sample Disturbance i n Athabasca O i l Sand. Journal of Canadian Petroleum Technology, V.19, N.2, Dusseault, M.B. Undrained Volume and Stress Change Behaviour of Unsaturated Very Dense Sands. Canadian Geotechnical Journal, Volume 16, 1979. Dusseault, M.B. and Morgenstern, N.R. Sands. Shear Strength of Athabasca O i l Canadian Geotechnical Journal, V.15, N.2, 1978. Dusseault, M.B. and Morgenstern, N.R. Characteristics of Natural Slopes i n the Athabasca O i l Sands, V.15, N.2, 1978. Florence, A.L. and Schwer, L.E. Axisymmetric Compression of a Mohr-Coulomb Medium Around or Circular Hole. International Journal f o r Numerical and A n a l y t i c a l Methods i n Geomechanics, Vol. 2, 1978. Fredlund, D.G. Density and Compressibility of Air-water Mixture. Canadian Geotechnical Journal, V.3, 1976. Geertsma, J . Some Rock-Mechanical Aspects of O i l and Gas Well Completions. Paper EUR 38 presented at the 1978 SPE European Offshore Petroleum Conference and Exhibition, London, October 24-26. Hardy, R.M. and Hemstock, R.A. Athabasca O i l Sands. Shear Strength Characteristics of K.A. Clark Volume, Alberta Research Council, Information Series No. 45, 1963. Harr, Groundwater and Seepage. McGraw-Hill Harris, M.C., Poppen, S. and Morgenstern N.R. Tunnels i n O i l Sand. Journal of Canadian Petroleum Technology, 1979. Harris, M.C. and Sobkowicz, J.C. Engineering Behaviour of O i l Sand. The O i l Sands of Canada-Venezuela, 1977. CIM Special Volume 17. Hughes, J.M.O., Wroth, C P . and Windle, D. Sands. Pressuremeter Tests i n Geotechnique 27, 1977. Naylor, D.J. Discussion. Proceedings of the Symposium on the Role of P l a s t i c i t y i n S o i l Mechanics, Cambridge, 13-15, p. 291-294. September 1973. Pasley, P.R. and Cheatham, J.B. Rock Stresses Induced by Flow of Fluids Into Boreholes. Soc. Pet. Eng., J . , p. 85-94, March 1963. Risnes, R., B r a t l i , R.K. and Horsrud, P. Wellbore. Sand Stresses Around a Society of Petroleum Engineering Journal, Col. 22, No. 6, 1982. Smith, L.B. and Burn, P.M. Convergence-Confinement for Shafts and Tunnels i n Oilsands. Method of Design Conference on Applied Oilsands Geoscience, Edmonton, Alberta, 1980. 113 Sobkowicz, J.C. The Mechanics of Gassy Sediments. Ph.D. Thesis, University of Alberta, C i v i l Engineering Department, 1982. Timoshenki, S. Strength of Materials. V o l . 2, 1941, Van Nostrand, New Yord. Todd, . V a z i r i , H. Ground Water Hydrology. John Wiley and Sons, Inc. Forthcoming Ph.D. Thesis, Department of C i v i l Engineering, U.B.C. 1985. APPENDIX A The constitutive relationship may be written: (Aa') = [ D ' ] (Ae) i n which A.l (Ao ) i s the incremental stress vector, (Ae) i s the 1 incremental s t r a i n vector and [D*] i s the incrmental e f f e c t i v e stress s t r a i n matrix. The incremental s t r a i n vector i s related to the nodal displacements by: (Ae) = [B] (S) A.2 i n which [B] i s a matrix that depends on element geometry. By the p r i n c i p l e of v i r t u a l work, the external work done by the v i r t u a l displacement i s equal to the i n t e r n a l work done by the increment of v i r t u a l s t r a i n s : (5) i n which T (f) ( f ) - /A U e ) = T (Ao') dA + J . (Ae") (}) Au dA T A.3 the element force vector (Aa') = the element incremental e f f e c t i v e stress vector (Au) the element incremental pore pressure vector = Substituting for (Ae) and (Aa') from Equations A . l and A.2, (6) (f) = (6) [ B ] [D'] [B] (6) A + (8) [ ] (*) T T T o T B Au A A.4 e i n which A e i s the area of the element with unit thickness, Rearrangement of Equation A.4 y i e l d s (6) = ( f ) - ( K j (Au) [K] i n which [K] = [ B ] T [K T w ] = L [B] 1 A.5 [ D ] [ B ] A^ ( 1 )e ^o' A where [K] i s the element s t i f f n e s s matrix and (& ) i s a load vector w associated with the pore pressure. APPENDIX B Assumptions: 1) the volume^ of solids i s 1 unit, then the volume of voids i s e units by the d e f i n i t i o n of void r a t i o , 2) the solids are incompressible. The undrained response of the element under a change on external pressure w i l l be Skeleton Aa ' (Ae ) v'SK B m S K Fluid i n which (Ae )_ = v f B.l B.2 (Ae ) , (Ae ) , = volumetric s t r a i n : Skeleton, F l u i d v SK v f ' Bg^, = bulk modulus: Skeleton, f l u i d Aa ' m = mean e f f e c t i v e stress change Au = pore pressure change Since the s o i l skeleton and the pore f l u i d s deform together when the conditions are undrained, i n additions to Equations B . l and 117 Compatibility e (A£ ) y or - f (1 + e) t^)^ B.3a (Ae ) = - (Ae ) v f n v SK B.3b i n which e i s the void ratio and n i s the porosity. Substituting B.3b into B . l , the pore pressure change, Au, may be written as: A u = K < v>SK a A e f = —, a n B K i n which K the apparent bulk f l u i d modulus In f i n i t e element analysis Au = K = K a (Ae ) v a (Ae n + Ae 2 2 + Ae ) 33 ' 4 118 APENDIX C Poisson's r a t i o i s equal to 0.5 i n undrained analysis theoretically. This means that bulk modulus i s i n f i n i t e and numerical i n s t a b i l i t y w i l l arise (this sentence does not sound r i g h t ) . Therefore the default of Poisson's r a t i o , v, i s always l e s s than 0.5 (e.g. 0.495) to maiontain s t a b i l i t y and accuracy. In the t o t a l stress model developed i n Section 3.2, i t i s the combined Poisson's r a t i o , for matrix [D] that controls the o v e r a l l numerical s t a b i l i t y , not just v. The e l a s t i c moduli i n matrix [D'] are related as G = E 2(l+v) C.I B = E 3(l-2v) C.2 Assuming the e l a s t i c moduli i n matrix [D] = [ D ' ] + [ D ] are related as F G cb C.3 B C.4 cb 119 i n which s u f f i c cb means combined, G = G , as pore f l u i d does not cb transmit shear. Just consider the d i r e c t stress terms i n the constitutive relationship Aa Aa L +K n 1 2 a L*M+K a L*M+K ss Ao-33 L a L*M+K a L+K a L*M+K a L*M+K a L*M+K a L +K a Ae n Ae 2 2 A E 33 C.5 4 i M n w h i U. c T h E(l-V) (l+v)(l-2v) = 1-v Substituting C.2 into C.5, y i e l d s , Aff Ao 2 2 0- 3 3 A P+K H = a Q+K a Q+K a Q+K a P+K a Q+K a Q+K a Q+K a P+K a Ae u Ae 2 2 A E 33 C.6 120 v i n which „ 3B(l-v) P = / Q = 3 ( 1 B V + v 1+v Adding the d i r e c t stress i n Equation C.6 and rearranging T ( A q " + Ag 22 (Ae^i + Ae 2 2 + A g + 33) = B(l-v) ^£33) + 2BV l-v +K a 1+v = B cb Eliminating E cb from Equations C.3 and C.4, y i e l d s 3B V cb C.7 = c b - 2G 6B°. + 2G cb Substituting C.7 into C.8, rearranging 0 , 8 121 APPENDIX D In the incremental e l a s t i c method, two i t e r a t i o n s are performed to obtain the tangent moduli of the s o i l skeleton. Hence two i t e r a t i o n s are also employed to evaluate the tangent bulk f l u i d modulus i n order to make the procedure compatible. The parameters at the end of previous increment w i l l be used i n the analysis of the present increment, and then updated at the end of t h i s increment. The tangent bulk f l u i d modulus and the procedure to update the parameters are shown herein: Equation 3.24e and 3.26b are programmed as: e-n (1+e) + H n w ww L l 0 r 3 i - f + 3 N ww J D.l e-n (1+e) - n (1+e) + H,n, (1+e) + H n (1+e) o w b b ww L f + (1+e)] r „ _ 1 82 - e + (1+e) + H „n (1+e) co2 w p W 1 + e > ] Vw The parameters i n Equation D.l and D.2 may be updated according to the following f o a u l a e based on the assumptions: D.2 volume of solids i s 1 unit bitumen, water and s o i l solids are incompressible Ae = e^ = S e f f c £ (1+e) Ae v e = i + A e S.e, 1i D.3 ^"^ D.5 n (1+e) + Ae -81+e + Ae D.6 n (1+e) < - rzzn b D. 7 1+e+Ae n' "w = n (1+e) w r£—n 1+e+Ae P = u + Au D.8 D.9 123 APPENDIX E Refer to Figure 4.1b °1 9 = an °12 cosec 29 E.l + a2 2 03 = x + ^ ~ °12 cosec 28 E.2 = 0 1 3 E.3 The overstress i n the element i s removed by reducing AOj, A a , 3 AT 1 2 as follows AO! = Oi - a Ao-c, = 0 E.5 = 0 E.6 Ax 1 2 tan 3 2 (45° + |-) E.4 These change i n p r i n c i p a l stresses can be expressed i n terms o f stresses i n x-y space: A o ll AO} + A a = 0 3 AO} - A a + * 3 cos 29 E.7 124 Aoi + Aa-i Ao"i - Aero cos 29 Aa 2 = 2 - AOi Aa 1 2 AO: cos 29 = E.8 E.9 Equations E.7, E.8 and E.9 may be expressed i n matrix form Ac r Aoi n Ao 2 [T] 2 Aa Aa E.10 3 Ati3 1 2 [ T ] i s the transformation matrix [T] = 1_ 2 cos 29 2 1_ 2 cos 29 2 i 1_ cos 29 2 1_ 2 cos 29 2 2 s i n 20 s i n 20 APPENDIX F Back Calculation of S o i l Parameters Given the compressibility of s o i l 0 g = 9 E-b k P a e f f e c t i v e stress o§ = 751 kPa, the bulk modulus K - 1 and the number can be B backcalculated by assuming a value of m. From Equation 3.12 B = *B a p P ( *' ) m a s 1 a Substituting 3 , and P S and assuming m = 0.25 (Byrne and Cheung) cl F.2, = 665 Adopting the relationship K B K E = 1108 Other parameters are depicted Cheung reports. written as: = 0.6 K (Byrne and Cheung) E from Byrne and Eldridge on Byrne and A complete set of s o i l parameters f o r Test 11 may b \ = 1108 n = 0.5 h " 665 m = 0.25 = 0.8 R f <t>' = 42° A<f>' = 8° n = 0.3228 (e = 0.4767) S = 0.9975 APPENDIX G E l a s t i c Stress Solution If the f l u i d pressure i s included, the displacement u of an e l a s t i c material may be written as (X + 2G) ^dr (dr T + ~) r ' + 3 dr i where A = G = 3 C - 1 - -p-— b 0 G.l (1+v) (l-2v) E 2 (1+v) which i s assumed equal to 1 C = sand matrix compressibility C, .= sand bulk compressibility b The pressure may be expressed by Darcy's law i n r a d i a l form iP _ dr "q 2irhKr , The stresses are written as o r = X e e v + 2G e + P r 6 G.3 a. = X e e v 6 a = X e Z e e Where e , e r 9 + e V and e 0 e + 2G + P G.4 + 2G e + P z G.5 0 e e 6 6 are the e l a s t i c s t r a i n components and e = e + e z v r v G.6 e z Assuming the i n i t i a l loading cause a deformation only i n the v e r t r i c a l d i r e c t i o n , but no displacement i n the horizontal directions (e r = e n 9 = 0), the i n i t i a l v e r t i c a l s t r a i n e , i s given by G.5 as zo a : zo -P zo o X + 2G = G , / Assuming only r a d i a l displacement a f t e r i n i t i a l loading, the strains are e r e e = — dr u 9 = 7 e G.8 e = e z zo G ' y By solving Equation G.l with the boundary conditions a a and combining r r = a . when r = R. ri i = a ro when r = R, i the results with Equations G.8, G.9 and G.10, the stres solutions (Equations 6.2, 6.3, 6.4) can be found by i n s e r t i n g the result i n Equations G.3 through G.5 APPENDIX H P l a s t i c Stress Solutions (a < a < a„) r z 0 The conditions must be s a t i s f i e d within the p l a s t i c zone 1) Equilibrium da 2) a - o„ Mohr Coulomb f a i l u r e f = a. - a t a n 0 r The 2 criteria a + (tan 2 a-1) P - 2S tan a = 0 H.2 flow rule associated with y i e l d condition i s e P , 9f ^ 2 = X = - X tan a r da r z el 0 e - X ~ - - X 9a H.3a H.3b Q P-\||--0 z da H.3c From Equations H.3a and H.3b, i t follows that e£ + tan 2 a = 0 H.4 The t o t a l s t r a i n components may be written as e r e e = e r + r E H ' = Eg + e^ g e = e z H.5b p z + e 5 a = e zo z H.5c = 0 because i t i s assumed that there i s only r a d i a l displacement a f t e r i n i t i a l loading. Combining Equations H.5a and H.5b, i n s e r t i n g into H.4, gives e o + e. tan^ a = 0 r E e r + e o tan^ a 0 E q H.6 Applying Hooke's law of e l a s t i c i t y for porous material, y i e l d s Ee® = a - v (a + a ) - (l-2v)P r r o z H.7a E E Q = a . - v ( o - + a ) - (l-2v)P 0 9 r z H.7b fl v Ee e z = a z - v (a + a ) - (l-2v)P r 9 Q H.7c Substituting Equation H.5c a = Ee z into H.7c + v (a zo r + a.) + (l-2v)P 8 Combining Equations H.7 c r i t e r i o n Equation H.2 [tan * a + 1 - v ( t a n 1 yields and H.8, i n s e r t i n g into the y i e l d with the s t r a i n r e l a t i o n i n H.6, a + 1) 1 2 2 a J a - = 2G £ r + 2G e„ t a n r 9 + [tan 2 a (tan + (tan 2 a + 1) (1 - 2 v)]p - [ t a n + v (tan 2 a H.8 ' 2 1) - v (tan * a - a 1) 1 2 2 i t gives a ( l - v ) - v] 2 Stan a 1) 2G E + H.9 zo and [tan * a + 1 - v ( t a n 1 a + l) ] o 2 2 Q = 2Ge o + [v (tan * a - i ) - ( t a n 1 + tan 2 a (tan + v tan 2 2 2 tan 2 a + 2Ge a - 2 a + 1) 2Ge tan* a 1) a + 1) (l-2v)]p - [v ( t a n a (tan n 6 r 2 a + 1) - l ] 2Stan a H.10 zo Substituting Equations H.9 and H.10 i n t o equilibrium Equation H.l, together with the strain-displacement relations G.8 and G.9, the displacement equation may be written as r 2 2 — j+ r dr - utan + (tan 2 a = 7^7 (- [ t a n 2 a + 1) (l-2v)] + v (tan 2 j^hk + a - 1) 2G e a (tan 2 ( t a f l 2 a + 1 2 ) a - 1) - v(tan < 1 _ 2 v ) 2 S 1 t a n l+ - 1) a H.ll zo J The displacement solution of Equation H . l l i s 2Gu = , A.^ tan a , . - tan a , _ + A r + Br 2 „ .. H.12 2 2 and the corresponding strains are 2Ge r = tan 2 .0 ^ r 1 ^ - 1 - tan 2 a A 2 r _ t a n 2 a + B H.13 00 A tan a 2Ge = A ^ 2 Q Where Aj and A conditions. 2 . , + A 2 -tan a r 2 ._ + B „ ,, H.14 are constants of integrations which depend on i n i t i a l 134 a v ,-i [f2? tan'•'a+l ^ l-2v -i tan^a-l z B = L - v 2G e Substituting + - — 2 — r JJ uq o ui 2irhk - l-2v „ tan ct-l 0 7—2 — r 2 Stana z H.15 zo Equations H.13 and H.14 into Equations H.9 and H.10 together with Darcy's law f o r r a d i a l flow, P = P * + i -Hi- £n — 2TThk X R n H.16 i gives a r = P, + -^rr In-Ii 2irhk R, oz + 2tan A l p u y t a ^— r a = - + o^T * 9 I 2irhk a - -\- (2Stana - ^ 3 - ) t 2TThk A l + 2tan*a — T n t r I H.17 7" (2Stana - t a n t 2 a -^tr") 2trhk H.18 a = P. + z 0 i .* Jin ZTrhk ^S_] + 2irhk ( 1 J A 2tan where — 2 R 7 ) ( 1 1-v 4Stan a - ( t a n t " 2 V (a ) z a + 1), L -P ) - + v ( t a n a + 1), zo o 2 l a - r ' t = tan 2 a - 1 T = tan 4 a + 1 - v (tan H.19 2 a +l ) 2 The constant of Integration Aj can be found by i n s e r t i n g the boundary condition, o r = P. when r = R., into H.17. i i ' 136 APPENDIX J The quantity of steady-state seepage at any distance r from the centre of the w e l l i s Q = k 2irrh ^ dr J«l x on setting the l i m i t of integration, _ R r , h yields 2 hi w 1 and Q = IT k (h2 - h') A n I J.2 w Equating J . l and J.2, 2rh§ yields = (4 -h?) j^sbr-'C J.3 w integrating J.3 h 2 = C Jin r + D J.4 137 where D i s the constant of integration Substituting the boundary conditions r = and h = h^ i n J.4, yields, 2 D = hi - C Jin r J.5 XJ * Substituting J.5 into J.4, y i e l d s 2 h l ~ ? = hi + - — 5 - 7 — In R/r w 2 1 h h r An — r J.6 w
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A stress-strain model for the undrained response of...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A stress-strain model for the undrained response of oil sand Cheung, Ka Fai Henry 1985
pdf
Page Metadata
Item Metadata
Title | A stress-strain model for the undrained response of oil sand |
Creator |
Cheung, Ka Fai Henry |
Publisher | University of British Columbia |
Date Issued | 1985 |
Description | An efficient undrained model for the deformations analysis of oil sand masses upon undrained loading is presented in this thesis. An analysis which couples the soil skeleton and pore fluids is used. The soil skeleton is modeled as a non-linear elastic-plastic isotropic material. In undrained conditions, the constitutive relationships for the pore fluids are formulated based on the ideal gas laws. The coupling between the soil skeleton and the pore fluids is based upon volume compatibility. The undrained model was verified with the experimental results and one dimensional expansion of soil sand cores. Comparisons between computed and measured responses are in good agreement and suggest that this model may prove useful as a tool in evaluating undrained response of oil sand. The response of a wellbore in oil sand upon unloading was analyzed using the developed model. Such analysis are important in the rational design of oil recovery systems in oil sand. |
Subject |
Oil sands |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-05-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062954 |
URI | http://hdl.handle.net/2429/25084 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1985_A7 C43.pdf [ 4.23MB ]
- Metadata
- JSON: 831-1.0062954.json
- JSON-LD: 831-1.0062954-ld.json
- RDF/XML (Pretty): 831-1.0062954-rdf.xml
- RDF/JSON: 831-1.0062954-rdf.json
- Turtle: 831-1.0062954-turtle.txt
- N-Triples: 831-1.0062954-rdf-ntriples.txt
- Original Record: 831-1.0062954-source.json
- Full Text
- 831-1.0062954-fulltext.txt
- Citation
- 831-1.0062954.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062954/manifest