STRESS, DEFORMATION A N D FLOW ANALYSIS OF OIL SAND MASSES UNDER APPLIED LOAD AND T E M P E R A T U R E CHANGES by T H I L L A I K A N A G A S A B A I S R I T H A R B. Sc (Engineering), University of Peradeniya, Sri Lanka A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA Apr i l , 1989 © T H I L L A I K A N A G A S A B A I S R I T H A R , 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Civil Engineering The University of British Columbia Vancouver, Canada Date A-f»vve 1959. DE-6 (2/88) Abstract Stress and temperature changes in oil sand masses associated with oil recovery can give rise to stability and deformation problems. A n analytical formulation is devel-oped to analyze such kind of problems and incorporated in the finite element program ' C O N O I L IF. O i l sand is a multiphase material, which can be analyzed as an equivalent two phase material; solid grains saturated with a single fluid phase, having an equivalent compressibility. The equivalent compressibility is obtained by combining the concepts of unsaturated soils with the gas laws and the gas-liquid interaction effects. Classical approach, as used in thermal elasticity is not appropriate to analyze the temperature induced effects, unless the time steps are made very small. A better formulation is presented, which involves introducing additional terms in the stress-strain relation and in the flow continuity equation, to account for temperature induced effects. The new approach reduces the oscillations predicted by the classical approach and leads to more realistic results. The analytical equations for the coupled stress, deformation and flow problems have been solved by finite element analysis. The finite element formulation involves nonlinear variations of stress-strain response, compressibility and flow, and is per-formed in an incremental manner. The developed finite element program has been verified by comparing its results with the closed form solutions and laboratory data. Then, the program has been applied to predict the responses associated with an oil recovery scheme. Such kind of analyses are important in the rational design of oil recovery schemes in oil sands. ii Table of Contents Abstract ii List of Tables vi List of Figures viii Acknowledgement ix 1 General Introduction and Scope 1 1.1 Introduction 1 1.2 Characteristics of O i l Sand 2 1.3 The Scope and Organization of Thesis 6 2 Review of Previous Work 7 2.1 Introduction 7 2.2 Derivation of Equivalent F l u i d Compressibility for O i l Sands 8 2.3 Stress-Strain Mode l Used 15 2.3.1 Shear Volume Coupling 19 2.3.2 Stress Transfer Concept 20 2.4 Analysis of Drained and Undrained Conditions 22 2.5 Analysis of Transient Conditions 25 2.6 Predictions of Pore Pressure and Volume Changes due to Temperature Variations 25 2.6.1 Temperature Effects on O i l Sand 26 2.6.2 Theoretical Analysis 30 i i i 2.6.3 Discussion on Application 33 3 Proposed Approach for Thermal Consolidation 34 3.1 Introduction 34 3.2 Modifications to Basic Equations 35 3.2.1 Equilibrium Equation . 35 3.2.2 Stress-Strain Relation 36 3.2.3 Strain Displacement Relation 37 3.2.4 Continuity of Flow 37 3.2.5 Thermal Energy Balance 40 3.2.6 Resulting Equations . 41 3.2.7 Summary and Discussion 42 4 Finite Element Consolidation Analysis 43 4.1 Introduction 43 4.2 Development of Analytical Model 44 4.2.1 Representation of Permeability . 45 4.2.2 Finite Element Formulation 47 4.3 Finite Elements and the Procedure Adopted 55 5 Verification and Application of The Program 58 5.1 Introduction 58 5.2 The Aspects Checked by Previous Researchers . , 58 5.3 Validation of Thermal Analysis 64 5.3.1 Drained Thermal Analysis 64 5.3.2 Undrained Thermal Analysis 68 5.3.3 Thermal Consolidation Analysis 70 5.3.4 Comparison of C O N O I L and C O N O I L - I I 74 5.4 Application to Well-Bore Problem 77 iv 6 Summary and Conclusions 86 6.1 Summary and Concluding Remarks 86 6.2 Recommendations for Further Research 87 Bibliography 89 appendices 97 A Load Shedding Techniques 97 A . l Constant Principal Stress Directions ; . 97 A.2 Constant Mean Stress 101 A.3 Constant Minor principal stress 102 A . 4 Method for Transferring to The Adjacent Elements 103 B Representation of Permeability 105 B. l Computation of Variation in e and S 106 B.2 Factor Fe 107 B.3 Factor F, 109 C Closed Form Solution for 1-D Thermal Consolidation 110 D Back Calculation of Coefficients of Thermal Expansion 114 v List of Tables 2.1 Parameters for typical quartz sand 18 2.2 Parameters for oil sand 18 5.1 Values of the parameters used 70 5.2 Values of the parameters used for consolidation analysis 74 5.3 Hyperbolic soil parameters for Athabasca clearwater formation oil sand 82 5.4 Moduli values for elastic and nonlinear elastic analyses 84 vi List of Figures 1.1 In-situ Structure of oil sand (after Dusseault,1980) . . . 4 1.2 Schematic spring analogy (after Dusseault,1979) 5 2.1 Conditions of stress state violations . . ; 21 2.2 Drained thermal expansion tests, (after Agar et al., 1986) 27 2.3 Undrained thermal expansion tests, (after Agar et al.,1986) 29 2.4 'Modelling of an element for temperature effects 30 4.1 Experimental and predicted values of viscosity (after Puttagunta et al., 1988) 48 4.2 Element types used 56 5.1 Stresses and displacements around a circular opening for an elastic material (after Cheung, 1985) 60 5.2 Stresses and displacements around a circular opening for an elastic-plastic material (after Cheung, 1985) 61 5.3 Comparison of observed and predicted pore pressures (after Cheung, 1985) 62 5.4 Comparison of observed and predicted strains (after Cheung, 1985) . 63 5.5 Pore pressure dissipation (after Vaziri , 1986) 65 5.6 Results for a circular footing on a finite layer (after Vaziri , 1986) . . . 66 5.7 Finite element mesh and the applied temperature field 67 5.8 Stresses and displacement in a circular cylinder 69 5.9 Finite element mesh for the soil sample 71 5.10 Undrained volumetric expansion 72 vii 5.11 Finite element mesh and the applied temperature for a column of soil 73 5.12 Pore pressure variation with time 75 5.13 Comparison of effective stress and pore pressure . . . 76 5.14 Finite element mesh for a layer of oil sand 77 5.15 Temperature profile with time 78 5.16 Effective stresses and pore pressure around the well 80 5.16 Effective stresses and pore pressure around the well 81 5.17 Variation of radial displacement 83 5.18 Comparison of effective stresses and pore pressures for incremental elastic and linear elastic skeleton behaviour 85 A . l Correction of stresses normal to the yield surface 98 A.2 Correct and incorrect Mohr circles for constant principal stress direction 100 A.3 Correct and incorrect Mohr circles for constant mean stress 101 A . 4 Correct and incorrect Mohr circles for constant minor principal stress 102 B. l Variation of permeability with void ratio (after Lambe and Whitman, 1969) 108 D . l Thermal volume change of U T F oil sands (after Kosar, 1989) . . . . . . 116 viii Aknowledgement The author is greatly indebted to his supervisor Professor P. M . Byrne for his guid-ance, valuable suggestions and the encouragement throughout this research. The author would also like' to express his gratitude to Professor Y . P. Yaid for reviewing the manuscript and making constructive criticisms. Appreciation is extended to Dr. M . Yogendrakumar for the helpful suggestions and discussions. Finally, the research assistantship awarded by Dept. of C iv i l Engineering. Uni-versity of British Columbia and the research grant provided by Alberta Oil Sand Technology and Research Authority ( A O S T R A ) are gratefully acknowledged. IX Chapter 1 General Introduction and Scope 1.1 Introduction The oil reserves in northern Alberta are one of the major resources in Canada. A p -proximately 5% of these oil sand deposits are found at depths less than 50 m and the rest are encountered at depths of 200-700 m. O i l recovery schemes involve open pit mining in shallow oil sand formations, and in-situ extraction techniques such as tun-nels and shafts in deep formations. In the in-situ extraction procedures, some form of heating is often required as the very high viscosity of the crude bitumen makes conventional recovery by pumping impractical. Analytical treatments are needed to design such oil recovery schemes, and to check their responses during operation so that modifications and improvements can be made. In the process of bitumen recovery by means of steam injection, temperature ef-fects are of prime concern. Therefore, in addition to the changes in stresses (loads), the changes in temperature have to be taken into account to predict the realistic re-sponses. The changes in temperature induce changes in volume, pore fluid pressures and consequently some of the engineering properties of soil mass such as strength, compressibility and permeability. The changes in volume comprise of change in vol-ume of solid particles, pore fluid and gas. The changes in pore pressure will occur due to the difference in the volume changes of the soil constituents. If the volume change of the pore fluid is greater than that of the voids in the soil skeleton, there will be an increase in pore pressure and a consequent reduction in effective stresses. The effective stresses may become zero and liquefaction may occur, if the oil sand is 1 Chapter 1. General Introduction and Scope 2 subjected to a rapid increase in temperature and, if an undrained condition prevails. Mathematical models have been developed by Harris and Sobkowicz (1977), Dusseault (1979), Byrne and Grigg (1980), Byrne and Janzen (1984) and Byrne and Vaziri (1986) to investigate these problems. A n analytical formulation that involves a fully coupled stress-deformation and fluid flow phenomena was developed by Byrne and Vaziri (1986). They further developed a finite element program namely C O N O I L to perform this analytical treatment. The program has been checked against closed form solutions and laboratory results where possible and gave excellent agreements. More recent application of the program to problems involving changes in tem-perature indicate that the method of including temperature changes is not stable when the temperature changes are rapid over a short period of time. This is because the commonly used classical 'Staircase' method of including temperature effects was used in the program. A better method of solving this problem involves applying the temperature as a rate of change or a gradient over the step. This thesis involves the development of such a method and incorporation of that in the finite element code C O N O I L II. 1.2 Characteristics of Oil Sand Oil sand is considered as a multiphase geologic material. Analytical methods dealing with the stress-strain and volume change behaviour of these materials, therefore, should include all constituent phases. The two dominant physical characteristics of oil sand are the quartzose mineralogy and the large quantity of interstitial bitumen. The quartz grains of the oil sand are 99% water-wet,which means the bitumen occupies interstitial space and the water phase forms a continuous film around the sand grains. The in-situ condition of oil sand is very dense and large number of contacts per grain exist. The pore space is filled with water, bitumen and part of that may be occupied by occluded gases. A n illustration of oil sand structure is shown in figure 1.1 (Dusseault, Chapter 1. General Introduction and Scope 3 1980). The simple one dimensional response of the oil sand may be modelled by a set of springs as shown in figure 1.2. The oil sand is split into two load carrying com-ponents; soil skeleton and pore fluids. The soil skeleton spring with compressibility 0T, characterizes the stress-deformation response of soil skeleton, which results in a change in effective stress. The pore fluid spring with compressibility 3j, character-izes the deformation of the pore fluid, including free gas, which result in changes in pore pressure. The relative stiffness of these springs control the sharing of the load between soil skeleton and pore fluid which is directed to the stiffer phase. Because of the relatively low permeability of bitumen-rich sands, the rapid unload-ing in oil sand associated with oil recovery methods occurs in an undrained condition and is particularly important in the analysis of oil sand. This is because, for a de-crease in confining stress, there will be a decrease in pore pressure which results in pore gas expansion and exsolution of dissolved gas. The increase in volume of pore gas will result in an increase of pore fluid compressibility, which in turn keeps the pore fluid pressure close to the liquid-gas saturation pressure. Therefore, the drop in the stresses due to unloading will be carried by the soil skeleton and be reflected as a decrease in effective stresses. This will continue until the effective stresses reduce to zero and instability conditions can occur. The coefficient of thermal expansion is also an important parameter for the pre-diction of thermally induced stress changes and deformations. Excess pore pressures are generated when thermal expansion is constrained and also owing to the fact that pore fluids expand more than the solids (Agar et. al, 1986). Therefore, for a multi-phase material like oil sand, thermal expansion also has to be considered for all its constituents. Figure 1.1: In-situ Structure of oil sand (after Dusseault,1980) Chapter 1. General Introduction and Scope o + Ac i i 1 4' i i i 1 •o \ g Soil I — § Skeleton ) p. <o 1 , o o Water o o § Oil (bitumen) § Stress = o' + Ad <=> ) o § Gas o Pore Fluid Pressure = u + A u Rigid Frame, Mo Lateral Yield o = a' + u Ao = Aa'+ Au AV « /(Aa, Au) Figure 1.2: Schematic spring analogy (after Dusseault,1979) Chapter 1. General Introduction and Scope 6 1.3 The Scope and Organization of Thesis Review of previous work is given in chapter 2. It describes the work carried out by Vaziri (1986) on this topic and the basic theories involved in developing an analytical formulation. The separate treatments for drained, undrained and transient conditions are also described briefly. The main purpose of this study is to present an efficient formulation for the deformation and flow analysis for oil sands due to temperature changes. The new for-mulation is to directly include the temperature induced stresses and volume changes in the stress strain relation and in the volume constraint equations. This results in a new set of equations to be solved for displacement and pore pressure unknowns. Details of this formulation is described in chapter 3. The finite element formulation for consolidation analysis, including the new ap-proach to analyze temperature effects is described in chapter 4. A n empirical cor-relation for changes in viscosity due to temperature and pressure changes, which is incorporated in the program, is also explained. Some important aspects in the adopted finite element procedure are also discussed in chapter 4. Chapter 5 is concerned with verification of the aspects considered in the analytical formulation. The finite element program results have been compared with closed form solutions and laboratory results where possible. Another major objective of this study is to validate the numerical procedure for simulating the oil recovery process. A n example of bitumen recovery process involving steam injection is considered and the program predictions are also described in chapter 5. Summary and conclusions are given in chapter 6. Possible future research in this area, and possible modifications to the program are also suggested in chapter 6. Chapter 2 Review of Previous Work 2.1 Introduction Theoretical expressions for the behaviour of oil sand are mainly derived from the concepts for the behaviour of unsaturated soils. In the analysis of a multiphase material, the key parameters which govern the response for applied loads are the compressibilities of the phase components. In the general analysis of unsaturated soil, it is considered as two load carrying components; soil skeleton and pore liquid with gas. The equation for the compressibility of the pore fluid components includes liquid pressure(pj), gas pressure(p f l) and for some cases surface tension between gas and liquid. For the analysis of oil sand, this equation is extended to its gassy nature by considering the gas laws and gas-liquid interaction effects. For a material with two phase fluid, both gas and liquid, Bishop (1957) derived a theoretical expression for the equivalent compressibility without considering surface tension effects. Sparks (1963) suggests that the surface tension effects should also be included in the effective stress equation for unsaturated soils. Fredlund (1973) claims that partly saturated soils should be regarded as four phase materials, with solids, water and air, . air water interface (referred as contractile skin) also being considered as an independent phase. Appropriate stress state parameters for this four phase material have been derived by Fredlund and Morgenstern (1977). The conceptual basis which can be used to develop an appropriate engineering practice for unsaturated soil is given by Fredlund (1979). However, the number of parameters required in describing the behaviour of unsaturated soil makes it unreasonable to 7 Chapter 2. Review of Previous Work 8 . accept that a simple treatment will apply widely. Harris and Sobkowicz (1977) first used the concepts of unsaturated soil and developed an analytical model to explain the behaviour of oil sand. Byrne and Janzen (1984) extended this for a nonlinear stress strain relation. A finite element program I N C O I L was developed to handle both undrained and drained analysis of oil sand. Cheung (1985) further modified the analytical model to get realistic response of expansion of oil sand core during unloading and to avoid the iterative procedure for obtaining the pore pressure. Vaziri (1986) claims that it is not possible to derive a unique relationship between stress and pore pressures for unsaturated soils which can be broadly applied to all soils and'under all conditions, and a successful application of any formulation requires certain idealization. He developed an analytical formulation, based on Cheung (1985) for oil sands, to treat most of the conditions that arises in practice. This formulation was incorporated in a finite element program C O N O I L . The following sections in this chapter briefly describe Vaziri's formulations and theoretical foundations for his formulations. 2.2 Derivation of Equivalent Fluid Compressibility for Oil Sands A general approach has been established in this section, to treat the soil regardless of its saturation state. The gas and liquid components of the pore, in partially saturated soils, have been replaced by a single fluid phase having an equivalent compressibility and thereby, it can be assumed that the soil is saturated with the equivalent fluid phase. Under isothermal conditions, it can be stated that (2.1) Chapter 2. Review of Previous Work 9 where, Bf - fluid bulk modulus ( A e v ) / - change in volumetric strain of the fluid Apj - change in fluid pressure Therefore, if we derive a equivalent fluid compressibility, stresses and strains for multiphase fluid can be related through the above equation. If the gas and liquid are immiscible the physical laws governing compression of gas will determine the fluid compressibility, while if the gas is soluble to some extent in the liquid, the fluid compressibility wil l also be influenced by the solubility relationship. The basic gas laws governing the volume and pressure relationships are Boyle's law and Henry's law. According to Boyle's law, under constant temperature conditions peg = wg RT where, p - absolute pressure T - absolute temperature R - universal gas constant eg - volume of gas uig - weight of gas For undrained conditions, the weight of gas does not change and therefore, equa-tion 2.2 can be written as peg = K (2.3) where, (2.2) Chapter 2. Review of Previous Work 10 K - a constant The presence of gas can be in both dissolved and free states. According to Henry's at constant temperature, is directly proportional to the absolute pressure of the gas above the solution. Mathematically. Wdg - weight of dissolved gas superscripts 0 & 1 - refer to initial & final conditions In other words, Henry's law implies that the volume of dissolved gas in a fixed quantity of liquid is constant at constant temperature and at a confining pressure p , when the volume is measured at p. Thus, law (Sisler et al., 1953); the weight of gas dissolved in a fixed quantity of a liquid, (2.4) where edg = H e, (2.5) where, H - Henry's constant, which is temperature dependent and over a wide range of pressure is also pressure dependent e; - volume of liquid Since the volume of dissolved gas is constant, free and dissolved gas components can be combined. Then application of Boyle's law to entire volume yields (Fredlund, 1976), Chapter 2. Review of Previous Work 11 (edg ^ e°fg) p° = (edg + e]g) p1 (2.6) where, ejg - volume of free gas The liquid and gas pressures can be related by considering the effect of surface tension. Kelvin's equation for a single capillary tube, relates the difference in the gas and liquid pressure to the surface tension and effective radius of the gas-liquid interface. Schuurman (1966) has used such a relation to develop an expression for fluid compressibility of unsaturated soils. As a result of the surface tension, the gas and liquid pressures in a system differ by an amount equal to the capillary pressure (p c). Pl=Pg-Pc (2.7) Let the average capillary radius in the system be rc, then, Pc = 2-7 1 (2.8) where, TS - surface tension Mitchell et al. (1965) determined an empirical equation for the average capillary radius. They claim that reasonable agreement -with the test results can be obtained if the following expression for the capillary radius is used. r . = r M . ^ 4 4 (2.9) (i - sf) where, Chapter 2. Review of Previous Work 12 rcs - capillary radius at saturation Sf - the lowest degree of saturation at which the liquid begins to flow freely 5 - current degree of saturation The value of Sj may be determined from permeability tests for samples with different degrees of saturation. Now by applying Boyle's law to the total volume of gas in the system, (e 0fg+edg)(Pa+P°g) = (e)g + edg)(pa+p lg) (2.10) where, pa - atmospheric pressure subscripts fg & dg - refer to free and dissolved components superscripts 0 &: 1 - refer to initial and final conditions pg - gas pressure Combination of equations 2.10, 2.7 and 2.8 leads to, ' (e)g + e<ig) rc Sparks (1963) states that for the case where gases exists as occluded bubbles, the assumption that fluid pressures and liquid pressures are identical, is reasonably valid. Then the fluid compressibility can be expressed as, „ 1 / de fa d e< \ , Since, P = Pf+Pa and Pj =Pi Chapter 2. Review of Previous Work 13 equation 2.12 can be written as, In case of oil sands, the liquid has two components, namely water and bitumen. Therefore, (dei)/(dpi) in equation 2.13 has to be considered as, dei _ de^ de^ dpi dpi dpi where de^ is the change in volume of bitumen. Differentiation of 2.11 gives, dp} (e% + edg)pa + (e% + edg)p°g 27\ dr^ by using the relationships, efg = e ( l ~~ s) a n < ^ edg — S.e.H Equation 2.15 can be written in the following form. dp} = (e0 - S0e0 + HSoe0){p0g + pa) 2T, drc Sx de)g (e1-Slel^HS1elf rldS.e, l ' b } where, e0 and ej - initial and final void ratios So and Si - initial and final degrees of saturation Now differentiation of equation 2.9 gives, drr TV, dSx 1 - Sf substitution of equation 2.17 in 2.16 results in, (2.17) dp} = (eo-gpeo-f HS0e0){p°g+pa) 2T, rcl S1 de)g (ei - Sle{ + HS1ex)* + r\ (1 - S,) e, 1 ' * j Chapter 2. Review of Previous Work 14 Equation 2.14 can be further written as, dei dew de^ dpi dpi dpi — ew.Cw e0.Cb = -S.e[{l-b).Cw + b.Cb} = -S.e.Ci (2.19) w here, b - percentage of bitumen volume in liquid C, =(l-b)Cw + bCb Cw, C[, - compressibilities of water and bitumen substituting equations 2.15 and 2.19 in 2.13 gives, where, L _ (e0 - S0eo + HS0e0)(p° + pa) (ei - Siei + HS-iei) 2 2TS rcs Sl M =~Tn c ^ - e° In most practical situations, the liquid phase has a very small compressibility (C( ~ 0) and the surface tension effects are negligible (T, ~ 0). Further, if we assume the load increments are small and there are no significant changes in void ratio and degree of saturation: i.e. Si ~ S0 and ea ~ e 0, equation 2.20 can be simplified to (l-S + SH) °' ~ (* + » ) ( 2 ' 2 1 ) Chapter 2. Review of Previous Work 15 The expression for fluid compressibility stated by Bishop and Henkel (1957) is the same as equation 2.21. In the normal course of analysis, it is expected that the stress changes are imposed in small increments so that we can assume S\ ~ So and e-[ ~ eo- Then equation 2.20 can be reduced to Cf= . O . N ^—+SCl (2.22) l-S + HS V c 2 ( l - 5 / ) 2.3 Stress-Strain Model Used The prediction of deformations of soil for applied load and temperature changes is mainly dependent on the stress-strain model used. However, it generally appears that no real soil can be accurately represented by a unique stress-strain model. Various stress-strain models have been proposed to represent the behaviour of soil; from linear elastic to highly sophisticated elastic plastic models; and each of the proposed model can only represent some particular constitutive feature of real soil. Because of the complexity in more sophisticated models and the need for more elaborate complex laboratory tests to obtain parameters for that model, an incremental elastic hyperbolic model (Duncan et al., 1980) is adopted. In the incremental elastic model, the stress-strain relationship over each load in-crement can be described in terms of any pair of elastic parameters; Young's modulus E and Poisson ratio v or bulk modulus B and shear modulus G. The later pair B and G are more fundamental as they represent the volume and distortion components. The shear modulus is somewhat difficult to measure accurately, therefore, in practice G is obtained from E through the relationship, Chapter 2. Review of Previous Work 16 where, 2 V W) The the stress-strain relation can be written in the incremental form as, Ao^ = (Bt - ^G t) Aejtfc^ - + 2GtAetj (2.24) where, Et - tangent Young's modulus for the increment Bt - tangent bulk modulus for the increment 6ij - Kronecker delta; = 1 when i = j = 0 when i ^ j A hj :perbolic stress-strain relation for soils, was first proposed by Konder (1963). Duncan and Chang (1970) extended Konder's work and Duncan et al. (1980) further studied in more detail and proposed the following expressions for the tangent moduli. Rf(l - sin^)(o-1 - o-3) 2c. cos cf> + 2a'3 sin <f> kK P« 11 Pa (2.25) Bt = kB Pa Pa (2.26) where, <f> = <f>\ (2.27) kE & kB • Young's and bulk modulus numbers n & m - Youngs and bulk modulus exponents Rf - failure ratio; maximum deviatoric stress to the asymptotic value predicted from the hyperbola Chapter 2. Review of Previous Work 17 c &i <b - Mohr-Columb strength parameters <f>i & A(b - value of friction angle at 1 atm. and reduction in that for ten fold increase in confining stress (This represents a. nonlinear failure envelope) <r'm - effective mean normal stress Pa - atmospheric pressure a'3 - effective confining stress Determination of these parameters from laboratory experiments or from S P T blow counts is well described by Byrne et al. (1982), Byrne (1983) and Byrne and Cheung (1984). Some typical values for typical quartz sands are given in table 2.1. The typical values for some oil sands obtained from laboratory tests conducted b}-Kosar (1989) are given in table 2.2. The nonlinear incremental elastic model has some limitations. Therefore, it is important to discuss these limitations and how some of them have been relaxed. The major limitations are listed below. 1. It does not adequately account for dilatancy of real soils when they are sheared. 2. It assumes the principal axis of strain increment always coincide with the prin-cipal axis of stress increment. This is perfectly valid in the region considerably below the failure surface, but not when approaching failure. 3. Inability to incorporate anisotropy, effects of principal stress rotations, stress path dependency etc. The approaches described in the subsections below circumvent some of the limita-tions on nonlinear incremental elastic model and make the soil to respond in realistic manner. Agar et al. (1987) carried out experiments to investigate the influence of Chapter 2. Review of Previous Work Table 2.1: Parameters for typical quartz sand Dr Ni kE n kB m fa A<f) <t>cv Rj 25 4 250-500 0.5 150-300 0.25 33 0 33 0.9 50 10 450-900 0.5 270-540 0.25 36 2 33 0.8 75 23 580-1500 0.5 450-900 0.25 41 4 33 0.7 100 >50 1500-3000 0.5 900-1800 0.25 50 9 33 0.6 where, Dr represents the relative density, Ni represents the S P T value normalized to 1 atmospheric pressure and' <j)cv represents the constant volume friction angle Table 2.2: Parameters for oil sand Type of Oi l Sand kE n kB m Rf fa A<f> fav Athabasca McMurray formation 1300 0.52 1050 0.34 0.7 60 18 31 Athabasca clearwater formation 1300 0.59 1000 0.49 0.7 65 21 31 Athabasca McMurray formation, interbedded 1200 0.5 1000 0.47 0.6 50 13 29 Chapter 2. Review of Previous Work 19 such factors as stress path dependency, microfabric disturbance and heating to ele-vated temperatures, on the shear strength and stress-strain relation of oil sand. In triaxial tests, since the sample orientation is kept as same as in-situ, the strain in-crement direction coincides with the stress increment direction. Therefore, it may be not possible to check the effects when the oil sand is approaching failure with the triaxial test results. However, curve fitting of their experimental data indicated that the hyperbolic model is a useful empirical technique for modelling the nonlinear stress-strain behaviour of oil sand. The effects of anisotropy which is most probable in nature is not studied, and a sophisticated model may be needed to analyze those conditions. 2.3.1 Shear Volume Coupling The dilatancy effects are incorporated in the hyperbolic model as suggested by Byrne and Eldridge (1982). A n additional term has been included to account for shear induced volume change and then equation 2.24 becomes, Aai:j = (Bt - ^ G e ) Aekk6tj + 2GtAe{j + BtAe'Jij (2.28) where, Ae^ - Volume change due to dilatancy Ae^ can be related to shear strain by means of dilation angle v, as given below. Ae-• - s i n t ; = - ^ (2.29) Hughes et al. (1977) obtained an expression to evaluate v, using modified form of Rowe's stress dilatancy theory. Chapter 2. Review of Previous Work 20 where, <f>m - mobilized friction angle (f>cv - constant volume friction angle It should be noted down that the formulation give accurate results for (f>m > (f>cv, but overestimates the volume contractions for (f)m < (j)cv. Therefore, as recommended by Byrne and Janzen (1984), no dilation is assumed for <bm < (j)cv. 2 . 3 . 2 Stress Transfer Concept In the incremental elastic stress-strain model adopted, the possibility of violating the failure criteria is quite likely. Some of such situations are schematically presented in figure 2.1. A n element having a stress condition represented by point A is considered elastic and capable of carrying additional load. But it is possible that an additional load increment will take the stress state into the inadmissible zone as shown by point B in the figure. One way of rectifying this problem is to control the magnitude of the load increment, but this involves too many load increment steps and it may not be easy to find the optimum magnitude in advance. . Another situation of stress violation may occur when an element whose stress state lies on the failure envelope as indicated by point C in the figure is subjected to unloading of normal stress. The stress state wil l move to point D and even reducing the magnitude of the load increment as for the earlier case will not solve the problem. For this kind of situations, the shear modulus is set to zero (or very close to zero) which means that the shear stress in that element will not change. Since the imposed load change is a decrease in normal stress, unless the strength envelope is horizontal, the predicted stress will violate the failure criterion. Under these situations, the resulting stress state have to be brought back to the Chapter 2. Review of Previous Work 21 Figure 2.1: Conditions of stress state violations Chapter 2. Review of Previous Work 22 failure surface. This is achieved by converting the resulting overstress as an equiv-alent load vector and transferring to the adjacent elements. The various techniques employed in such load shedding procedure are discussed in appendix A . 2.4 Analysis of Drained and Undrained Conditions In the finite element analysis, the undrained response of soil is usually analyzed with the total stress parameters. Therefore the total stress parameters have to be obtained from undrained tests. If pore pressures are desired, they are computed from Skempton equation relating total stress changes to pore pressure parameters. Byrne and Vaziri (1986) adopted an approach similar to the one proposed by Nay-lor (1973) to model the undrained behaviour of oil sand. This approach is developed from the principle of effective stress and the stiffness matrix for a total stress analysis is obtained from the effective stress parameters and the compressibility of the fluid components as shown in this section. Then the solution procedure is carried out in the usual manner to obtain deformations. The pore pressures can be evaluated from the computed deformations using the relative contribution of the fluid and the skeleton, without the use of the Skempton equation. It is also possible to analyze the drained conditions with the same formulation as explained later in this section. The incremental effective stresses are related to the incremental strains by the following relationship. {Aa'} = [D'}{Ae} (2.31) where, {e} - strain vector {cr'} - effective stress vector [D'\ - matrix relating effective stress and strain Chapter 2. Review of Previous Work 23 The volumetric strain can be expressed as, (Ae t,) = {m}T{Ae} (2.32) where, Ae„ - volumetric strain in the skeleton {m}T = { 1 1 1 0 0 0} - a vector indicating only direct strains wil l be involved For undrained conditions volume compatibility requires that the volume change in the skeleton equals the volume change in the fluid, i.e, (Aev)f = ~{Aev) (2.33) n where, (Ae„ ) / - volume change in the pore fluid n - porosity In the previous section, all the pore fluid components in the oil sand have been con-sidered as a single fluid phase and the equivalent compressibility has been obtained. The changes in pore pressure can therefore be expressed as, Ap = Bj(Aev)f (2.34) where, Ap - change in pore pressure Bf - inverse of the equivalent compressibility Chapter 2. Review of Previous Work 24 Substitution of equations 2.33 and 2.32 in 2.34 gives, A p = -±{m}T{m} (2.35). n From the definition of effective stress, in an incremental form, {Ao-} = {ACT'} + { m } A p (2.36) where, {cr} - total stress vector By substituting equations 2.35 and 2.31 in 2.36 { A a } = \{D'} + ^{m}{m}T]{Ae} (2.37) Equation 2.37 express the relation in terms of total stress where matrix [D] for the total stress is given by, [D] = \D'} + ^{m}{m}T (2.38) n Equation 2.37 is used in the finite element formulation for undrained conditions. The pore pressure is not an unknown in the resulting system of equations, but is obtained from equation 2.35, once the deformations are computed. It is also possible to perform a drained analysis with the same data base by simply setting the value of Bf to zero. Then equation 2.37 will become as an effective stress equation. Byrne and Vaziri (1986) claim that this method has definite advantages such as, adaptability for saturated or unsaturated soils and for any stress-strain relation. In particular, this method gives stable solutions when the effective stresses go to zero and all of the load is carried by the pore fluid. For an incompressible fluid Bj becomes infinity and the above formulation becomes i l l conditioned. However, this can be overcome by setting the value off?/ to a suitably high but finite value (Naylor, 1973). Chapter 2. Review of Previous Work 25 2.5 Analysis of Transient Conditions Extension of Biot's three dimensional consolidation theory is used to analyze the transient conditions. The theory has been extended to allow for nonlinear stress-strain behaviour of soil, nonlinear flow of pore fluid and for the effects of temperature changes. A coupled stress-flow finite element formulation is developed on this basis which is also applicable to unsaturated soils. The nonlinear stress-strain relationship can be modelled by deriving the governing equations in a general incremental form. The nonlinear flow of the pore fluid is dealt with expressing the permeability as a function of void ratio, degree of saturation and coefficient of viscosity. Changes in void ratio are directly related to the changes in volumetric strain, and for unsaturated soils the variation in the degree of saturation can be determined from the changes in void ratio and quantity of flow in or out of the soil mass. As explained in earlier sections of this chapter, saturated soil with three phase fluid can be replaced by a equivalent compressible phase, which is phenomenolog-ically treated as a saturated soil. It is therefore possible to apply Biot's theory to unsaturated soil provided the effect of fluid compressibility is accounted for in the for-mulation. Details of this stress-flow finite element formulation with the new method of including temperature effects are described in chapter 4. 2.6 Predictions of Pore Pressure and Volume Changes due to Temperature Variations In this section, the behaviour of oil sand and the formulation used in the program earlier to predict the excess pore pressure and the volume changes due to temperature changes are described. Even though this formulation is not adopted in the analysis any more, it will be useful to describe it here, so that the basic concepts of the oil Chapter 2. Review of Previous Work 26 sand behaviour due to temperature changes and the shortcomings of this approach can be understood. 2.6.1 Temperature Effects on Oil Sand A complete examination of the effect of temperature on oil sand is beyond the scope of the present study. However, some of the engineering parameters which are signifi-cantly influenced by the temperature variations are described. As explained in the previous sections, drained thermal expansion can be consid-ered to represent a lower bound to thermal expansion behaviour while the upper bound occurs under undrained conditions. The change in volume is composed of the expansion of solid particles and the pore fluid components and the structural change of the soil skeleton. Under undrained conditions, the expansion of the fluid is more dominant than the expansion of the solid particles and the volume changes are com-parably large. Furthermore, if gas exsolution occurs, there wil l be additional large volume changes. Under fully drained conditions, the overall volume change of the soil skeleton is the volume change of immobile solid phase and is very small. Agar et al. (1986) conducted a number of drained and undrained tests in the high temperature consolidometer to study the thermal expansion behaviour of oil sand. The results from their drained tests are shown in figures 2.2. A small volume decrease has been observed during the-initial stages of first time heating. This may be attributed to rearrangement of solid grains. Increase in temperature causes changes in interparticle contact forces and leads to a particle rearrangement towards the denser side. This phenomenon has also been identified hy Campanella and Mitchel (1968) for ilhtic clay and by Scott and Kosar (1984) for undisturbed oil samples. Campanella and Mitchel quantified the magnitude of the volume change in terms of a parameter a>t known as the coefficient of structural reorientation. Scott and Kosar analyzed this parameter and concluded that ast is not a soil property, but varies with void Chapter 2. Review of Previous Work 27 1.35 105 ;075 ,0.45 • 0.15 -0.15 POKE PRE SSURE = 5 MPo i 1 f G f N D I X P A N S I O N OF S A N D GRAINS EFFECTIVE C O N F I N I N G MRESS = 6 0 0 0 k P o EFFECTIVE C O N F I N I N G STRESS = S O l P o • PORE PRESSURE =15MPo SOIL STRUCTURE COMPRESSION 40 80 120 160 200 T E M P E R A T U R E C O -240 280 a) Drained volume expansion T E M P E R A T U R E l " C ) b) Drained coefficient of cumulative thermal expansion Figure 2.2: Drained thermal expansion tests, (after Agar et al., 1986) Chapter 2. Review of Previous Work 28 ratio and temperature and dependent on the soil type. Since the greatest structural reorientation occurs in soils having their loosest possible state, one might deduce that such an effect is not a dominant feature of oil sands which are generally very tightly packed under high confining stresses. However, this aspect has been considered in the analysis to provide a complete modelling. Furthermore, figure 2.2(a) shows that thermal expansion followed by the initial structural rearrangement occurs at approximately same rate as that of the individual sand grains. The sudden volume decrease at 200° C is due to the increase in the back pressure from 5 M P a to 15 M P a . Variation of drained coefficient of cumulative thermal expansion is given in figure 2.2(b). The drained coefficient of thermal expansion at initial stages during first time heating deviates by a large amount because of the structural collapse and because cumulative volume change has been considered for the estimation. Figure 2.3 shows the results from undrained thermal expansion tests. The pore pressure and the effective confining stress were maintained constant by heating very slowly and allowing continuous vertical expansion of the sample in order to isolate thermally induced volume change. The beginning stages of gas exsolution were also observed as a stage at which a drastic increase in volume occurs. The drained thermal expansion is also given in the figure. It can be observed that the volume change in undrained condition is much more compared to the drained condition, even without gas exsolution. Figure 2.3(b) shows the undrained coefficient of cumulative thermal expansion and its variation with the temperature, due to the increase in the coefficient of thermal expansion of the pore fluids at high temperatures. Another important effect of the temperature increase is the change in viscosity. At in-situ temperatures, the viscosity of bitumen is very high and therefore, the permeability of oil sand is very low. Increase in temperature of oil sand will increase the permeability, which allows high oil recovery in practice. Experiments have been Chapter 2. Review of Previous Work 29 25 20 60 KX) 140 180 220 260 300 T E M P E R A T U R E 1'C) a) Thermal expansion of oil sand i 7 8 10 - ' T ' I • 1 —, , ( GAS EXSOLUTION LIMIT ^ ^ ^ i S M P o ^ • PORE FLUID PRESSURE 0 . 0 5 M P o s ^ - - ^ ^ ^ M P t a jT 2MPb ^ • 1 i L. J — i i „, — i • 1 l 1 1 . 1 i 1 1 1 1 1 1 I , t 20 60 100 140 180 220 260 300 T E M P E R A T U R E CO b) Undrained coefficient of cumulative thermal expansion Figure 2.3: Undrained thermal expansion tests, (after Agar et al.,1986) Chapter 2. Review of Previous Work 30 done by several researchers to study the variation of viscosity with temperature and pressure. Recently Puttagunta et al. (1988) proposed a empirical correlation for the changes in viscosity due to the changes in temperature. This has been incorporated in. CONOIL- I I and the details can be found in chapter 4.2.2. 2.6.2 Theoretical Analysis Consider a general element of soil as shown in figure 2.4(a). A change in temperature will cause thermal expansion of soil constituents and that will result in changes in boundary stress, pore pressure and volume changes. The analysis of this problem is carried out by following the conventional thermal elastic analysis approach as in solid mechanics. In this approach the response for a temperature increment is obtained as follows. Ac9 (a) A e 0 = 0 A p 0 Acr' A / 00 Nodes Fixed + - A / (c) Fixity according to Boundary Conditions Figure 2.4: Modelling of an element for temperature effects Chapter 2. Review of Previous Work 31 1. Ad is applied to each element under constrained conditions (zero strain and flow) and the stresses and pore pressure changes Aa'0 and Apo are computed. 2. The forces Af to restrain the element as mentioned above, are then evaluated (figure 2.4(b)). 3. The global force system A F arising from the element forces Af is assembles and applied to the domain in the opposite direction. The stresses Acr[ , pore pressure Api and strains A e i are computed for this condition (figure 2.4(c)). 4. The stresses and deformation changes are then computed as follows, ACT' = Aa'0 + Aa[ (2.39) Ap = A p 0 + A p a (2.40) Ae = Aeo + Aex (2.41) For the situation depicted by figure 2.4(b) the changes in effective stress and pore pressure can be computed, using the criterion that the sum of separate volume changes of the soil constituents must be equal to the change in volume of the soil skeleton. By invoking such strain compatibility, theoretical expressions can be obtained for the pore pressure and stress changes. If we consider the changes in volume of soil constituents and soil skeleton due to stresses and the changes in volume due to thermal expansion and structural reorien-tation, volume compatibility for undrained conditions requires that, A e . + Ae/ + Ae f l + Ae,t = 0 (2.42) A e . , = 0 • (2.43) where, Chapter 2. Review of Previous Work 32 A e , - change in volume of solids due to stresses and thermal expansion Ae ; - change in volume of liquid Aeg - change in volume of gas Aetk - change in volume of soil skeleton Aett - change in volume due to structural reorientation By substituting the relevant components, the following expressions can be ob-tained (Vaziri (1986)) for stresses and pore pressures, _ A9(escts + e,a, + egag - etka,t - eska, B, E>i Bg Bt B, Aa'0 = Bsk a , A T - ^ A P o (2.45) where, e - volume a - coefficient of thermal expansion B - bulk modulus (inverse of compressibility) subscripts s,l,g,sk - solid, l iquid, gas, skeleton respectively a3t - coefficient of structural reorientation To represent drained condition, where there is no change in pore pressure, effective stresses can be explicitly determined from equation 2.45, as a function of change in temperature. Once A p 0 a n d Aa'0 are known, the equivalent nodal forces ( A / ) to restrain the element can be evaluated in the normal fashion, by applying the principle of virtual work. Chapter 2. Review of Previous Work 33 {AF} = [ [B}T {APo + Aa'0}dv (2.46) where, [B] - the matrix of shape function derivatives { A F } - change in global force vector 2.6.3 Discussion on Application The main disadvantage in this approach is that the temperature changes are applied as discrete increments. Although this is the most widely used approach to include temperature effects, it may predict unrealistic results for a coupled consolidation analysis, unless the temperature increments are kept very small. This effect is partic-ularly important where the rates of pore pressure dissipations are high, and the pore pressures predicted by this approach shows high oscillations. When we analyze the temperature changes by this approach, the pore pressures are evaluated under constrained conditions (figure 2.4(b)), and then the equivalent nodal forces are applied and analyzed for the conditions explained in figure 2.4(c). If we use such an approach to analyze transient conditions, it is first analyzed outside the time frame for constraint conditions and then analyzed inside and the final results are obtained by adding the results from these two conditions. There is an uncoupled time factor involved and to reduce its influence the time steps have to be made very small. Therefore it is better if we include the effects of temperature changes in the governing equations and solve for both pore pressure and displacement unknowns. The formulations according to this new approach are described in chapter 3. Chapter 3 Proposed Approach for Thermal Consolidation 3.1 Introduction Temperature variations can significantly influence the behaviour of soil and give rise to stability and deformation problems. Temperature variations in soil can arise in a number of energy and mineral resource development areas in addition to ther-mally enhanced oil recovery form oil sand reservoirs. Underground storage of nuclear waste, power plants, underground power cables, liquefied natural gas storage and underground pipe lines are some examples: Chapter 2.4 describes the importance of temperature changes in the analysis of oil sand and the disadvantages of the classical thermal elastic method used earlier. A better method has been developed to analyze the transient condition with temper-ature changes, which involves additional terms in the governing equations. Booker and Savvidu (1985) included temperature effects in the stress-strain and continuity of flow equations, combined with thermal energy balance equation and derived a closed form solution for consolidation around a point heat source. Aboustit et al. (1985) also used the same method of having additional terms in the basic equations to account for temperature effect and solved the equations by finite element formulation using variational principle. But in both of their works and in most of the other research works, it has been assumed that the soil is a linear elastic material, having constant permeability and fully saturated and incompressible fluid. The above mentioned assumptions are not valid for the case of problems associ-ated with oil sands. In this chapter, the formulations are extended to overcome these 34 Chapter 3. Proposed Approach for Thermal Consolidation 35 assumptions and solved by a finite element scheme using Galerkin's weighted resid-ual method. The fluid components occupying the pores are assumed compressible, and are modelled as a single component with an equivalent compressibility. Details of the skeleton and the pore fluid characteristics and the derivation of equivalent compressibility were described in chapter 2. 3.2 Modifications to Basic Equations The modifications in the basic governing equations to account for the effect of temper-ature variations are described below. The generalized three dimensional consolidation equations derived by Biot (1941) are used here to develop the formulation. However, in the present formulation the equations have been extended for unsaturated soils, nonlinear stress-strain relation, nonlinear permeability variation and for the effects of temperature changes. The governing equations to be solved for a consolidation problem with temperature effects are, 1. Equihbrium equation. 2. Stress-strain relation. 3. Strain-displacement relation. 4. volume constraint or continuity of flow. 5. Thermal energy balance. 6. Boundary conditions. 3.2.1 Equilibrium Equation Using the conventional Cartesian tenser notation, for the equihbrium of stress, Chapter 3. Proposed Approach for Thermal Consolidation 36 a^-F^O (3.1) where, (Tij - total stress tenser subscript ,j = 9 F{ - body forces a_ 3.2.2 Stress-Strain Relation In tenser notation the effective stress can be written as, P 6ij (3.2) w here, a'- - effective stress tenser p - pore pressure 6ij - Kronecker delta The effective stresses are related to the strains by a tenser containing bulk and shear modulus terms ( D ^ ; ) . In a nonlinear stress-strain relation this tenser will vary with stress level, but since we represent the nonlinearity as incremental elastic this is constant for an increment. Thus, Aa'{j = Dijkl Aekl (3.3) where. - tenser relating effective stress and strain - strain tenser Chapter 3. Proposed Approach for Thermal Consolidation 37 If we include temperature induced stresses, ACT',. = D i j k l A e « + V A 0 (3.4) where, b' — B,k-cts B,k - bulk modulus of the skeleton At9 - change in temperature (in °C or in °K) a, - coefficient of thermal expansion of solids The notation 6 is used hereafter to denote the temperature to avoid confusion with t, which represent time. 3.2.3 Strain Displacement Relation Assuming infinitesimal strains, ^i^\{U^+Ujti) (3.5) where, U{ - displacement vector 3.2.4 Continuity of Flow For a completely incompressible pore fluid, the change in volume of the soil skeleton must balance the volume of liquid entering or leaving the system. However, in the case of a compressible fluid, the decrease in the volume of the remaining fluid should also be taken into account. Then the volume constraint equation can be written as, Chapter 3. Proposed Approach for Thermal Consolidation 38 where velocity vector volumetric strain superscript at n 1 — n Bf ' Bs bulk modulus of the homogenized compressible phase bulk modulus of the equivalent fluid phase B. bulk modulus of the solids The term Bcp in equation 3.6 also includes the compressibihty of the solids, thus all the soil constituents are considered as a homogenized compressible phase. The term pjB^ takes account of the volume change due to the compressible components of the soil. By assuming the flow through the soil is governed by Darcy's law, the velocity vector is given by, Now by considering the volume change due to thermal expansion, equation 3.6 can be written as, Vi = — ^ pj If (3.7) where, 7/ - unit weight of the permeant k^ - coefficients of permeability in generalized Darcy's law - Vi,i + uiti JL Bc cp + cxeq 9 = 0 (3.8) where, Chapter 3. Proposed Approach for Thermal Consolidation 39 cteq - equivalent coefficient of thermal expansion The equivalent coefficient of thermal expansion can be obtained by considering the volume and coefficients of thermal expansion of all soil constituents, i.e., cteq = ctt(l — n) + aw nw + ab nb + ag ng + cx,t (3-9) where, subscripts s,w,b,g - denoting solid, water, bitumen and gas respectively n - porosity a,t - coefficient of structural volume change (expansion is assumed positive here) n with subscripts - fraction of volume of components in pore The coefficient of thermal expansion of gas can be obtained as follows, According to gas laws, = (3.10) Only the volume change due to the temperature change has to be considered to evaluate the coefficient of thermal expansion. Thus, by assuming the gas pressure remains the same, Hi - ^ Vo #1 v0 d0 Chapter 3. Proposed Approach for Thermal Consolidation 40 Adopting the usual notations, (3.12) Thus, 1 (3.13) Where 6 is the temperature in K. Now equation 3.9 can be written as oceq — 0:5(1 - n) + aw nw + ab nb + - ng + a,t (3.14) 3.2.5 Thermal Energy Balance In many applications mechanical contributions to energy balance may be neglected when compared to thermal contribution. Then, the net heat inflow into an element of soil should be equal to increase in the internal energy. By neglecting the convection terms, The equivalent specific heat can be obtained in a similar way as for obtaining equivalent thermal expansion. - hiti - ceq 6 . (3.15) where K - heat flux vector c e q - equivalent specific heat ceq = (1 - n)psct + nwpwcw + nbpbCb + ngpgcg (3.16) where, Chapter 3. Proposed Approach for Thermal Consohdation 41 p - density c - specific heat By assuming that the heat flow is governed by Fourier's law of heat conduction, The thermal energy balance will give the temperature profile over the domain and its variation with time. In fact, this part is not included in the analytical formulation. It has been solved separately with the heat flow boundary conditions, by a separate program which can handle a sophisticated heat flow equation including convection, and the temperature profile is evaluated and considered as a known condition for the analysis. 3.2.6 Resulting Equations The fundamental unknowns to be solved in a consohdation analysis are the displace-ment and the pore pressure. This can be done by combining the basic equations together as follows, Substituting equations 3.5 and 3.2 in 3.1, hi = - (3.17) where, Kn - coefficient of heat conduction (3.18) Similarly substituting equation 3.7 in 3.6, (3.19) Chapter 3. Proposed Approach for Thermal ConsoHdation 42 These equations are solved for a given boundary conditions by a finite element method using Galerkin's weighted residual scheme as explained in chapter 4.2.3. 3.2.7 Summary and Discussion The above mentioned approach to treat the temperature changes has been successfully applied in thermo-elastic analysis of solids by several researchers. It has been adopted to analyze thermo-elastic consolidation analysis for soils by Aboustit et al. (1982), Aboustit et al. (1985) and Booker and Savvidu (1985). Here this has been further extended to handle nonlinear stress-strain relation, nonlinear flow, compressible fluid and unsaturated soils. In this method, the temperature changes and the time are coupled and analyzed within the governing equations. The predicted results shows smooth variations with time while the earlier approach adopted in C O N O I L predicts highly oscillating results with time. In fact, this approach does not need any additional parameter but only changes the method of analysis. To validate this technique a number of examples are examined and the computed results are checked against the analytical solutions and experimental and an excellent agreements are found. Details of these are described in chapter 5. Chapter 4 Finite Element Consolidation Analysis 4.1 Introduction Consolidation is a time dependent process in which an applied total stress is trans-ferred as effective stress through a gradual dissipation of pore pressure to an equi-hbrium value. The change in temperature wi l l also cause changes in stress and pore pressure and lead to consolidation processes, if the permeability is low. It is a cou-pled phenomenon of fluid flow and soil deformation; the fluid flow depends on the soil permeability and the induced pore pressure gradient, whereas the soil deformation depends on stress-strain properties of the soil. This chapter deals with the finite element formulation to solve the equations ob-tained in chapter 3 for consolidation in oil sands. The equations were derived by considering Biot's three dimensional theory and extending it for unsaturated soils, nonlinear stress-strain relation, nonlinear permeability variation with void ratio, sat-uration and for the effects of temperature. The derived equations are solved by finite elements using Galerkin's weighted residual scheme. In the finite element formulation the nonlinear stress-strain relation is simulated in an incremental elastic manner. Nonlinear flow is modelled by treating the perme-ability varying with void ratio, degree of saturation and coefficient of viscosity. As described in chapter 2.2 the three phase fluid is considered as a single phase with an equivalent compressibility which is further considered as varying with saturation. 43 Chapter 4. Finite Element Consolidation Analysis 44 4.2 Development of Analytical Model The analytical models used in consolidation analysis are mainly based on the theo-ries developed by Terzaghi (1923) or Biot (1941). Terzaghi's theory involves a one dimensional flow and settlement solution using diffusion theory while Biot's theory involves a three dimensional case and both using the conventional elastic approach. Closed form solutions have been derived by a number of researchers for consoli-dation problems based on the above theories. For instance, Biot (1941 (b)) derived solution for consolidation under a rectangular load distribution. De Jong (1957) obtained solution for consolidation under uniformly loaded circular area on a semi-infinite soil. MacNamee and Gibson (1960) obtained solutions to plane strain and axisymmetric problems of strip and circular footings on a consolidating half space. Booker (1974) derived solutions for square, circular and strip footings on a finite soil layer. Solution for consolidation around point heat source in a saturated soil mass was derived by Booker and Savvidov (1985). However, the closed form solutions are possible only for simple boundary conditions. It is very difficult to obtain closed form solutions for complex boundary conditions and for varying load and soil properties with time. Computer aided techniques such as finite element methods has made possible the analysis of the consolidation phenomenon for more involved boundary conditions in a realistic manner. Sandhu (1968) developed the first finite formulation for two di-mensional consolidation using a variational principle. Christian and Boehmer (1970) and Sandhu (1972) advanced the finite element formulation towards stable solutions. Hawang et al. (1972) proposed a formulation using Galerkin's weighted residual method, however, in their procedure the pore pressures are obtained from the total stresses through the pore pressure parameters. Aboshi et al. (1985) developed finite element formulation for thermo elastic consolidation using a variational principle. A l l the finite element formulation described so far are based on Biot's theory. Thus, Chapter 4. Finite Element Consolidation Analysis 45 the solutions are valid for an incompressible fluid and for fully saturated soils having an isotropic and linear elastic skeleton behaviour with no variation in permeability during the course of consolidation. Lewis et al. (1976) assumed a incremental elastic stress-strain model for the soil skeleton and use a nonlinear permeability of soil in the consolidation analysis. A method of analysis for consolidation problems considering compressible fluid was proposed by Ghaboussi and Wilson (1973). Ghaboussi and K i m (1982) analyzed consohdation in saturated and unsaturated soils with nonlinear skeleton behaviour and nonlinear fluid compressibility Byrne and Vaziri (1986) developed a finite ele-ment formulation which includes nonlinear stress-strain relation; nonlinear variation of permeability with void ratio, saturation and viscosity; nonlinear variation of fluid compressibility with saturation and temperature induced effects. This formulation has been modified here to handle the temperature induced effects on deformation and flow analysis in a more efficient and realistic manner. In consohdation problems associated with oil sand, the relative movements of liquid and gas phases have to be considered. The gas may move with the liquid, it may not move if the bubbles become large or it may flow separately from liquid. The formulation developed here only considers the movement of gas with liquid which is often encountered in oil recovery. However, it should be noted that the formulation is limited to a saturation level of around 80%, below this the treatment of gas and liquid as a single phase may not be strictly valid. 4.2.1 Representation of Permeability For the case of problems considered in this study, variations in the coefficient of permeability is an important concern. The coefficient of permeability depends on the characteristics of the fluid and of the soil skeleton. The particle size distribution, shape of the particle, irregularity of the pore space and changes in viscosity of the Chapter 4. Finite Element ConsoHdation Analysis 46 fluid influence the permeability of the soil. Therefore, a semi-empirical approach may be appropriate to deal with the practical situations. Hence by adopting an approach similar to that proposed by Chang and Duncan (1983), and extending it for changes in viscosity, the following equation can be obtained. k = ^A (4.1) where, k - permeability of unsaturated soil k, - permeability of fully saturated soil = e 3 / ( l + e ) ' 6 eo7(l + e0) e, e 0 - current and initial void ratios Sf - threshold degree of saturation at which liquid begin flow freely F M - factor to account for changes in viscosity of bitumen due to temperature and pressure A detail description and determination of the factors F, and Fe can be found in appendix B. The justification for using the above expression is also given in that appendix. Since there was no successful theoretical or empirical techniques available to account for the change in viscosity of bitumen due to temperature and pressure changes, the change in viscosity had to be experimentally determined and considered as an input factor F^ at desired stages. But recently Puttagunta et al. (1988) proposed an empirical correlation for the viscosity variation due to temperature and pressure changes, which is given below. Chapter 4. Finite Element Consolidation Analysis 47 m/Z(gf>) = 2.3026 6 3.0020 + Bo P exp(d 6) (4.2) e-3o 303.15 where, 9 temperature in degree Celsius V pressure in M P a gauge b log/*(3o,o) + 3.0020 a 0.0066940.6 + 3.5364 Bo = 0.0047424.6 + 0.0081709 d -0.0015646.6 + 0.0061814 This empirical correlation needs only a single viscosity value at a reference temper-ature of 30 °C and at a reference pressure of one atmosphere, to predict the variation of viscosity with temperature and pressure. The authors claimed that their correla-tion is valid for most of the bitumens, and the error using this formula is less than 6% when compared to the experimental results . Figure 4.1 shows a very good agreement between experimental and predicted values of viscosity for two different bitumens. This correlation has been incorporated in the program C O N O I L II. The viscosity factor, Fp is determined by employing equation 4.2, instead of inputting F M . 4.2.2 Finite Element Formulation The method of solving the set of equations derived in chapter 3.2.6 is described in this section. The finite element formulation to solve these equations are developed adopting the Galerkin's weighted residuals method. Basically there are a number of ways available to adopt a finite element procedure, such as, variational principle, principle of virtual work and minimizing the residuals. The choice of the method depends on the type of the problem and boundary conditions on one. hand, and the knowledge of the mathematics involved on the other hand. In the Galerkin's scheme Chapter 4. Finite Element Consolidation Analysis 48 50000 10000 10 —- empirical equation * experimental 20 40 60 80 TEMPERATURE, *C a) Wabasca bitumen too 120 50000 •-. 10 — empirical equation * experimental _ L 23 40 SO 80 100 120 140 TEMPERATURE, *C b) Cold Lake bitumen Figure 4.1: Experimental and predicted values of viscosity (after Puttagunta et al. , 1988) Chapter 4. Finite Element Consolidation Analysis 49 a single application of Green's theorem is needed and that will give a set of integral equations for an element. This can be easily formulated in matrix form and solved. From chapter 3.2.6, the equations to be solved are, Diju- \{Uk,i + Ui,k) + pj + b'dj - F i = 0 (4.3) I J ,j \--kijPj) + U i t i - - ^ - + a J = 0 (4.4) The above differential equations have at most second order derivatives of displace-ments and pore pressures, which are the unknowns to be solved for. But by applying Green's theorem it can be brought down to first order. Therefore, to solve the re-sulting integral equations the assumed functions for displacement and pore pressure should be continuous. Let us assume that the approximate solutions for displacement and pore pressure variations U" and p" within an element are given by, U" = Nu8e (4.5) Np qe (4.6) where, S'T ={81,62,...,Sn} - nodal displacements qeT = {<7l ,92 , • • • ,Qn} - nodal pore pressures Now if we substitute these approximate solutions to the differential equations, they will not exactly satisfy, but will give some residual errors r\ and r 2 as given below. Chapter 4. Finite Element Consohdation Analysis 50 + p-. + b'6j - Fi = n (4.7) (-— k{j p") + U-iti -£--raJ=r2 (4.8) In the weighted residual method these residual errors are minimized in some fash-ion to give the best approximate solution. Thus, for best solution, I w r dv = 0 (4.9) where. r - residual error w - weighting function In Galerkin's method the weighting function is chosen to be the same as the assumed shape function. Therefore, J V r,dv = 0 (4.10) Furthermore, / Np r2dv = 0 Jv (4.11) e - Bu8' Ui,i = m e (4.12) (4.13) (4.14) where, Chapter 4. Finite Element Consolidation Analysis 51 mT = {1 1 1 0 0 0 } Bu &: Bp - shape function derivatives The possible boundary conditions in a consolidation problem can be specified as follows. 1. Applied tractions crij nj = fi for t>0 (4.15) 2. Specified or Known displacements Ui = Ui for t>0 (4.16) 3. Specified pore pressures (Drained boundary can be expressed as setting the specified value to a known value) p = p for t>0 (4.17) 4. Undrained boundary Vi n{ = 0 for t>0 (4.18) where, - normal vector to the boundary bar symbol - indicates a known quantity These boundary conditions will cause additional residuals, if the assumed solution does not satisfy them. Except the applied tractions, other boundary residuals can Chapter 4. Finite Element Consohdation Analysis 52 be eliminated by choosing appropriate weighting functions, while applying Green's theorem. Green's theorem for integration of two functions £ and £ over the subdomain fie of an element can be written as, Jn vc v£ dne = -J^ c v 2 1 dn + ^ / | | <*re (4.19) where, Te is the boundary around Qe and 7; is the normal to the boundary. By substituting equations 4.12 to 4.14 and 4.7 to 4.8 in equations 4.10 and 4.11 and applying Green's theorem, I Bl D BJdv + / Blm Npqdv = f / V j Tds + / JVj Fdv - [ B* m b'6dv (4.20) Jv J V J 5 */ U Jv ( — Bl k Bpqdv + I Nl mT Bjdv - I - j - A ' p T Npqdv = - [ i V j a u ^ « (4.21) </« 7/ *'t' t't' •'-'cp •'f For a time increment A£ the above set of equations can be written in matrix form as follow. [K]{AS} + [L){Aq}={AA} (4.22) where, [Lf {AS} - At [E] {q} - [G] {Aq} = -{A<?}_ (4.23) [K] =IvBlDBudv [L] =JvBTumNpdv \E] =fv±BlkBp l°] =SvirpNlNpdv {-4} = S.NfTds-r IvNlFdv - IvBTumb'6dv {C} =JvNpaJ Chapter 4. Finite Element Consolidation Analysis 53 [K], [E] and [G] matrices vary with the values of stress, void ratio, saturation and permeability. The equation 4.23 is considered over a time increment At, and therefore, qin that equation has to be expressed as, q = (1 - a)qt + aqt+At (4.24) Where a is a parameter corresponding to some integration rule. For example, a = 1/2 implies trapezoidal rule, a = 0 implies a fully explicit method and a = 1 gives a fully implicit method. But Booker (1974) showed that for an unconditionally stable numerical integration | < ot < 1. Then equations 4.21 and 4.22 can be written as, [K]{A6} + [L]{qt+±t-qt} = {AA} (4.25) [L]T {A8} - At [E] {(1 - a)qt + aqt+At} - [G] {qt+At - qt} = - { A C } (4.26) By rearranging the terms, [K] {A8} + [L]{ql+At} = {AA} + [L] {<?«} (4.27) [L}T {AS} - [[E] a At- [G]} {qt+At} = [[E] (1 - a) At - [G]} {qt} - {AC} (4.28) This can be further written in matrix form as, \K] [L] [L)T [~[E] a At - [G]} A8 < > = A-4 + [L]{qt} [{E} ( l - a ) At-[G}){qt}-{AC} J (4.29) Chapter 4. Finite Element Consolidation Analysis 54 It has been said that an unconditionally stable solution can be obtained if | < Q < 1. Therefore, to simplify the formulation expressed by equation 4.29, a value of 1 is assumed for a . Hence, \L)T. [-[£?] At-[G\] AA + [L}{gt} -{G}{qt}-{AC} (4.30) This can be further rearranged as, f AS AA \ [L)T [- [E] At - [G]} 1 ^ . _ {E]At{qt}-{AC} t (4.31) By changing the notation, equation 4.31 can be written in the usual matrix form as. ' [K] r < [E']_ where, {AC} = At{E]{qt} - {AC} \E>] =-At[E]-[G] (4.32) Equation 4.32 gives the matrix equation to be solved for an element. From the element matrix equations a global matrix equation is formed and solved for displace-ment and pore pressure unknowns. Stresses and strains are then evaluated from the displacements. It should be noted that it may be not possible to use the above consolidation routine to get the initial condition results; i.e. at t — 0. This is because for At = 0 {AC} = 0 and [E'\ = [G\. Chapter 4. Finite Element Consolidation Analysis 55 If the fluid is incompressible \G] wil l become zero and equation 4.31 will become i l l conditioned. For this situation an appropriate solution can be obtained by assuming a small value for At. This will circumvent the i l l conditioning, however, a better way to get around this problem is to use the undrained routine (chapter 2) to get the initial condition and then use the consolidation routine. 4.3 Finite Elements and the Procedure Adopted The type of finite elements considered, the way of forming the stiffness matrices and on the way of solving the simultaneous equations in the program are identical to that described by Vaziri (1986), and are summarized in this section. Lower order elements are incapable of predicting reasonable results for the prob-lems in which the failure state is approached or in highly constrained situations. Hence, higher order elements were used in the program. Furthermore, to obtain accurate results, the pore pressure and the stresses should have the same order of variation. Since stresses are related to the derivatives of displacements, pore pressure variation should be one degree lower than the displacement variation. On these con-siderations two types of triangular elements have been used in the program and they are shown in figure 4.2. The technique suggested by Taylor (1977) has been adopted for the generation of the stiffness matrix to save computer time and cost. This technique avoids unneces-sary operations associated with zero terms in the stiffness matrix. The higher order elements used cause a considerable increase in the band width and this requires much more time to solve the resulting simultaneous equations. Therefore, instead of a traditional band solver, the frontal solution technique (Irons, 1979) is used in the program. In the frontal solution technique, the equations to be solved for nodal variables are assembled and eliminated on an element by element basis. The frontal width (analogous to band width) is dependent on the number of Chapter 4. Finite Element Consohdation Analysis A Displacement nodes (2 d.o.f) O Pore pressure nodes (1 d.o.f) Linear strain triangle 6 displacement nodes 3 pore pressure nodes 6 nodes and 15 d.o.f. Cubic strain triangle 15 displacement nodes 10 pore pressure nodes 22 nodes and 40 d.o.f. Figure 4.2: Element types used Chapter 4. Finite Element Consohdation Analysis 57 nodes per element and the arrangement of the elements within the domain. There-fore, element numbering is important here to reduce the size of the front width. For a relatively complex finite element mesh, this wil l be a difficult task to be imposed on the user. This inconvenience is circumvented by incorporating an element order optimizer (Solan and Randolph, 1981) in the program and the user need not to be concerned about the numbering of nodes and elements. The improved modified Eular method is used to perform nonlinear finite element analysis. In this method, any incremental changes to the loading and boundary conditions is analyzed twice. In the first cycle of analysis, parameters based on the initial conditions of the increment are used. In the second cycle, parameters based on the mean conditions of the increment are used. In the modified Eular method this process is carried out until the difference between the successive results satisfy the specified tolerance. Since such an iterative procedure can increase the computer cost drastically, the improved modified form of Eular method is adopted. In this method, the final results are obtained at the end of second cycle. However, to maintain the equilibrium and to follow the nonlinear path closely the imbalanced loads at the end of second cycle is estimated and added to the next increment. Chapter 5 Verification and Application of The Program 5.1 Introduction The validity of the finite element program CONOIL- I I is examined in this chapter. The program has an excellent capability to handle some special aspects such as; partial saturation of soiL stress transfer concept, gas exsolution phenomenon and thermal consolidation analysis. Since the program deals with a number of aspects, the best way of validating the program is by considering each aspect separately. Examples for which closed form solutions or experimental results are available have been used to check the performance of the program on the important aspects. Once validated, the program is used to predict the responses for a real oil recovery problem. 5.2 The Aspects Checked by Previous Researchers The following aspects have been checked by Cheung (1985). 1. The general performance of the program in prediction of stress and strain. 2. Stress transfer treatment of the program. 3. Gas exsolution phenomenon. Since these aspects have been kept in CONOIL- I I without any changes the veri-fications are still valid. A brief description and the results of the aspects verified by 58 Chapter 5. Verification and Apphcation of The Program 59 Cheung are given in this section in order to provide a complete documentation on verification of CONOIL-II . The general performance has been verified by considering a thick wall cylinder under plane, strain conditions, with a linear elastic stress-strain model. Closed form solutions for this problem were obtained from Timoshenko (1941) and the computed and the closed form results are found to be in good agreement as shown in figure 5.1 taken from Cheung (1985). The stress transfer operation of the program has been examined by considering a problem of a tunnel subjected to unloading, under undrained condition. In a tun-nel, reduction of support pressure causes yielding and creates a plastic annular zone extending from the inside face. Closed form solutions for this kind of problems, for a elastic-plastic material have been developed by Hughes et al. (1977). Figure 5.2 taken from Cheung (1985) shows the comparison of program results and the closed form solutions and are in good agreement. Laboratory triaxial test results for gassy soil samples taken from Sobkowicz (1982) have been compared with the program results to check its validity on gas exsolution phenomenon. Sobkowicz (1982) carried out triaxial tests to predict the short term undrained response, i.e. no gas exsolution, and the long term undrained response, i.e. with complete gas exsolution. The comparison of these test results with the program results is shown in figures 5.3 and 5.4. The measured and computed results are in remarkably good agreement. Further to the aspects mentioned in the earlier paragraphs, the formulation for consolidation analysis, for soils with compressible fluids, has been validated by Vaziri (1986). Analytical solutions were derived by modifying Terzaghi's one dimen-sional consolidation solution for linear elastic material. The coefficient of consolida-tion has been modified as, Chapter 5. Verification and Apphcation of The Program 60 Radii (r/r 0 ) E = 3000 MPa v = 1/3 i n i t i a l s t r e s s : o r <= o. = 6000 kPa f i n a l s t r e s s : o r = 2500 kPa i n s i d e r a d i u s : ry >= 1 m Figure 5.1: Stresses and displacements around a circular opening for an elastic ma-terial (after Cheung, 1985) Chapter 5. Verification and Apphcation of The Program 61 2 ° CO 00 W o o Tf CO > 1H o 0) JCL c l o s e d form o program 6 E co C £ cn ' O " a CO 4 Radii u -E = 0 3000 1/3 40° (r /r . ) B 10 kPa i n i t i a l f i n a l i n s i d e s t r e s s s t r e s s r a d i u s Or Or To oe •= 1500 1m 6000 kPa kPa Figure 5.2: Stresses and displacements around a circular opening for an elastic-plastic material (after Cheung, 1985) Chapter 5. Verification and Application of The Program 62 Figure 5.3: Comparison of observed and predicted pore pressures (after Cheung, 1985) Chapter 5. Verification and AppHcation of The Program 63 0 20 40 60 Effective SigmaP (kPa) (X101 ) Figure 5.4: Comparison of observed and predicted strains (after Cheung, 1985) Chapter 5. Verification and AppHcation of The Program 64 C" = T T where, c'v - modified coefficient of consolidation R = D'/Bc D' - ID constrained modulus Bcp - equivalent fluid bulk modulus Closed form solutions and computed results agree ver}r well and are shown in figure 5.5. The consolidation routine was also validated by Vaziri (1986) for the two dimen-sional case. Closed form solutions developed by Gibson et al. (1976) for a circular footing resting on a layer of finite thickness, and for a fully saturated, elastic material was considered. Comparison of computed results and closed form solutions are shown in figure 5.6 and are in good agreement. 5.3 Validation of Thermal Analysis Validation of the stress-strain model, stress transfer technique, gas exsolution phe-nomenon, drained, undrained and transient conditions have been given in the previous section. Verification of the remaining aspect and the major undertaking of this thesis namely temperature analysis is described in this section. 5.3.1 Drained Thermal Analysis In order to check the implementation of the proposed formulation, the axisymmetric, plane strain closed form solutions obtained by Timoshenko and Goodier (1951) for a long elastic circular cylinder subjected to temperature changes has been considered. (5.1) Chapter 5. Verification and Apphcation of The Program 65 e* t T U i F.etor T* - — — * »I Figure 5.5: Pore pressure dissipation (after Vaziri , 1986) Chapter 5. Verification and Apphcation of The Program 66 T 1 r a)Amount of settlement i 1 I A n a l y t i c a l S o l u t i o n 0 F i n i t e Element A n a l y s i s x/B - 0 y/B - 0 v - 0.3 \ \ v • 0.0 D/B - 1 -1 • I irr4 io-3 10-2 1 0 - i i.o 10 c t b) Degree of settlement Figure 5.6: Results for a circular footing on a finite layer (after Vaziri, 1986) Chapter 5. Verification and Application of The Program 67 The finite element solution is obtained by considering a disk within the long cylinder as shown in figure 5.7. 100 c Applied Temperature Profile 1 100 m M e s h c o m p o s e d o f 4 0 l i nea r s t r a i n t r i a n g l e s A B & CD - v e r t i c a l l y r e s t r a i n e d Figure 5.7: Finite element mesh and the applied temperature field The boundaries A B and C D are allowed to move only in lateral direction and a drained analysis is carried out. A non uniform temperature field is applied which is also shown in figure 5.7. The closed form solution for this problem at a steady state condition is as follow 1 - p, r2 Ja 1 + ft \ 1 E I C\ C2 (5.2) ag = — 6.r.dr + - 7r + ~ \ ( 5-*> \ - ur 2 Ja l - / x I - p\\-2p r 2) Chapter 5. Verification and Application of The Program 68 q.E.9 2.p.E.C1 1 ^ 7 " (!+#*)(! -2 / * ) ( 5 - 4 ) 5 = I* d.r.dr ~dr+ — (5.5) 1 - / i r J a r where, crr - radial stress at radius r erg - hoop stress at radius r crz - vertical stress at radius r 8 - radial displacement at radius r 9 - temperature at radius r C\,C2 - constants can be evaluated according to the boundary conditions a - internal radius b - external radius Ci and C-2. can be determined by considering the boundary conditions on the stress free boundaries, i.e., crr = 0 at r = a <rr = 0 at r = b The values used for the analysis and the corresponding values of C\ and C 2 are listed in table 5.1. The closed form solutions and the finite element results are shown in figure 5.8 and are in excellent agreement. 5.3.2 Undrained Thermal Analysis The program is validated for the undrained thermal analysis by comparing the exper-imental results obtained from Kosar (1989) with the program results. The test was Chapter 5. Verification and Apphcation of The Program 69 3000 -1000 b 0 25 50 75 100 Radial Distance(m) Figure 5.8: Stresses and displacement in a circular cylinder Chapter 5. Verification and Application of The Program 70 Table 5.1: Values of the parameters used E r1 a b 0a eb C2 MPa m/°C m m °C °c m 2 180 0.25 l x l O - 4 0.5 100 100 0 1.389 x l0~ 3 6.945 XIO" 4 carried out on undisturbed oil sand samples in a high temperature consolidometer. The pore fluid pressure and the total confining stress were kept constant by heating slowly to prevent mechanically induced volume changes. The results volume changes were measured- at steady state conditions, i.e. there is no temperature gradients within the sample. The oil sand sample is modelled by four elements as shown in figure 5.9. The coefficients of thermal expansion were back calculated from test results as described in appendix D. The coefficient of thermal expansion of the fluid varies from 5.7 x 1 0 _ 4 / ° C at 20 °C to 7 x 10~4/ o C at 200 °C and then increase up to 1.4 xl0 3/°C at 250 °C . The coefficient of thermal expansion of solids is constant and have a value 4 x 1 0 ~ 5 / ° C . Computed and measured results agree very well and are shown in figure 5.10. 5 . 3 . 3 Thermal Consolidation Analysis The proposed formulation has been verified for a drained analysis in the previous section. Validity of this formulation for a consolidation routine has to be checked next. Closed form solution for one. dimensional thermal consolidation is derived from the solution for constant rate of loading obtained by Aboshi et al. (1970) for a one dimensional case. The details of the derivation are given in appendix C. A column of soil subjected to a uniform temperature rise is examined. The column of soil is modelled by 10 linear strain triangular elements and drainage is allowed from Chapter 5. Verification and Application of The Program 71 in C\l d B 0.75 AB — totally restrained BC.AD - laterally restrained Figure 5.9: Finite element mesh for the soil sample the top end only, as shown in figure 5.11. The temperature is assumed to be applied at a constant rate, 100° C / h r throughout the sample. The closed form solution for one dimensional thermal consolidation is given by the equation below. 16 M n A0 (a , - a.) ^ 1 . muz v = zr* 2^ ~i sin —vr \ 1 - e x P 7T3 T m" 2H ' 2 2\ T m TT where, p - pore pressure at distance z at time t T - time factor |jr! n - porosity Ad - change in temperature at time t M - constrained modulus (5.0) Chapter 5. Verification and Application of The Program 8 Temperature(° C) Figure 5.10: Undrained volumetric expansion Chapter 5. Verification and Application of The Program 73 Free Surface Uniform Temperature Change applied at a rate of 100°/hr < > 1 m AD, BC - Laterally Restrained AB — Totally Restrained Figure 5.11: Finite element mesh and the applied temperature for a column of soil Chapter 5. Verification and Application of The Program 74 a; - coefficient of thermal expansion of liquid Q, - coefficient of thermal expansion of solids The soil properties used for this analysis are given in table 5.2. The pore pressures have been obtained at two different depths, z/H = 0.875 and zjH = 0.5 denoted by P I and P2 in figure 5.11. Closed form solutions and the computed results of pore pressure agree very well and shown in figure 5.12. Table 5.2: Values of the parameters used for consolidation analysis E V Oil a, k M H z MPa m 3/°C m 3/°C m/s MPa m m 15 0.25 l x l O - 5 l x l O " 3 2 x l 0 ~ 6 18.3 2 1.75 5.3.4 Comparison of CONOIL and CONOIL-II The same example as in the previous section is considered for comparison. C O N O I L has been used to predict the response for the column of soil subjected to uniform temperature changes. The temperature changes are taken in steps as shown in figure 5.13(a). The predicted pore pressures and the effective stresses for point P I are shown in figure 5.13(b) and 5.13(c). The pore pressures predicted by C O N O I L shows high oscillations while C O N O I L -II predict smooth results which are in excellent agreement with the closed form results. The oscillations also reflected in the effective stresses and the results from C O N O I L implies negative stresses are occurring. Such negative stresses, if real, could cause fracture in the oil sand. Consequently the oscillations predicted by C O N O I L could result in misleading predictions concerning fracturing, amount of oil recovery, stability of the well etc. Chapter 5. Verification and Application of The Program z / H = 0 .875 -cr o o o o o CONOIL-II closed form solutions Figure 5.12: Pore pressure variation with time Chapter 5. Verification and AppUcation of The Program 76 75 o CO to o D_ 50 H 25 ^ I T •*0—c -t—e-o 50 CO (A in ^ 0 o o > - 5 0 -CONOIL o o o o o CONOIL-II Closed Form -t-rQ- T I / I / I / J-t 4 I I I / I I I / I / I I 1/ -r—r i / 11 i/ A i i ' O T-r i / i/ i/ 0 i I i r 1000 —I r—r 2000 Time(s) I i i i [ — 3000 4000 Figure 5.13: Comparison of effective stress and pore pressure Chapter 5. Verification and Application of The Program 77 5.4 Application to Well-Bore Problem The developed finite element program has been applied to predict the response of a steam-injection oil recovery process. Figure 5.14 represents a layer of oil sand at a depth of 200 m. Typical in-situ stress and pore pressure values for this depth are assumed; av = 3500 kPa , uT = 2000 kPa, cre = 2000 kPa and p = 500 kPa . Finite element mesh of the wellbore through a layer of soil Well Center Line Mesh composed of 44 linear strain triangles 650 m AD - Outflow boundary, laterally restrained BC - Inflow boundary, totally restrained AB & C D - No flow boundaries, laterally free Figure 5.14: Finite element mesh for a layer of oil sand The oil sand layer is heated by injecting steam through the well bore. The steam is injected at a temperature of 200°C and the resulting temperature profiles at various times during the steaming process are shown in figure 5.15. These temperature profiles were obtained from a separate analysis and are considered as known temperature loading condition. The fluid pressure at the well bore is raised and maintained at 1100 kPa (600 kPa above the in-situ pore pressure), to avoid fractures in the formation Chapter 5. Verification and Application of The Program 78 Figure 5.15: Temperature profile with time Chapter 5. Verification and Apphcation of The Program 79 during heating. Sufficient time is allowed (10 days is assumed here), to obtain a heated zone to the desired radial distance. Then the steam injection is stopped and the well bore pressure is allowed to return to its in-situ value. Material behaviour is initially assumed to be linear elastic. However, the problem was also analyzed with nonlinear stress-strain law using the hyperbolic incremental elastic representation to check their influence as described later in this section. The properties of the oil sand used are; E = 1000 M P a , v = 0.25, Bf = 2.7 x 106 kPa , a, = 3 x 1 0 ~ 5 m / ° C , aj = 3 x 1 0 ~ 4 m / ° C , the in-situ permeability, k = I0~9m/s and the dynamic viscosity at 3 0 ° C and zero gauge pressure, p — 20 P a s. Computed effective stresses and pore pressure are shown in figure 5.16. Figure 5.16(a) shows the immediate effects on the stress conditions obtained 1 minute after steam injection started. Only the elements adjacent to the well bore experience temperature changes. Pore pressure in these elements has increased to a value of 1900 kPa and the radial and tangential effective stresses have decreased, because of the rise in pore pressure. Figure 5.16(b) shows the results one day after the steam injection started. At this time the soil elements up to 7 m from the face of the well bore have the same temperature changes. Since the permeability of the elements close to the well bore have increased due to previous temperature rises, and since they are close to the drained boundary, the pore pressure in these elements have dropped and now have a value close to the well bore pressure (1100 kPa). The most severe condition occurs at this time a distance of about 6 m from the well bore wall, where the effective hoop stresses decreased to a value of 900 kPa from their initial in-situ value of 2000 kPa. In some oil sand formations, depending on their properties and the applied temperature, the effective hoop stress may go negative and would cause a fracture in the oil sand. Figure 5.16(c) shows the stress and pore pressure profiles, 8 days after the steam injection process started. A t this stage, the temperature rises have caused a large Chapter 5. Verification and Application of The Program 80 Figure 5.16: Effective stresses and pore pressure around the well Chapter 5. Verification and Application of The Program 81 6000 D CL °-. 4000 D_ c o co in to 2000 UJ 20 After 8 day; 0"r 30 40 6000 o C L -X '—s 0. a! 4000 c o CO CD CO CO £> -t-< c/i 2000 Just after s team injection stopped IA. 0"r a Q P 1 1 1 r 10 20 Radial Distance(m) 30 40 Figure 5.16: Effective stresses and pore pressure around the well Chapter 5. Verification and Application of The Program 82 reduction in viscosity and a subsequent increase in permeability. This result in a complete dissipation of pore pressure. The pore pressures are equal to the well bore pressure of 1100 kPa in the heated zone which is up to 18 m from the face of the well bore. Further injection of steam for the rest two days just extends the heated zone further about one meter which is not shown here. Figure 5.16(d) gives the stress and pore pressure profiles just after steam injection stopped (1 minute and 10 days in the total time scale). The pore pressure in elements next to the well bore have returned to their in-situ values of 500 kPa. Figure 5.17 shows the variation of radial displacement. The displacements reach their maximum at completion of the steam injection process (after 10 days). Max-imum displacement of 16mm is computed at a distance of 15m from the well bore. Once the steam injection is stopped, the well bore pressure returned to its in-situ value of 500 kPa , and the displacements start to decrease. The problem is also analyzed with nonlinear elastic soil parameters obtained from triaxial test results given in Kosar (1989) for Athabasca clearwater formation oil sand. The values for the parameters used are listed in table 5.3. Table 5.3: Hyperbolic soil parameters for Athabasca clearwater formation oil sand kE n kB m Rf <f>i 0cv 1300 0.59 1000 0.49 0.7 65 21 32 The moduli values used in the earlier linear elastic analysis and in the nonlinear elastic analysis are listed in table 5.4. The values given for the nonlinear elastic analysis are the initial moduli values corresponding to the parameters listed in table 5.3 for a effective confining stress of 2000 kPa . The effective stresses and pore pressures are examined for the most severe con-dition predicted earlier, i.e. 1 day after the steam injection processes started. The Chapter 5. Verification and Apphcation of The Program 83 0.02 c £ o _o CL CO 0 . 0 1 H 0.00 •- After 1 day After 1 0 days (Steam injection completed) After 1 0 days •+- 1 min (1 min after s team injection stopped) i 1 1 r 25 50 Radial Distance(m) Figure 5.17: Variation, of radial displacement Chapter 5. Verification and Apphcation of The Program 84 Table 5.4: Moduli values for elastic and nonlinear elastic analyses E B V MPa MPa Elastic 1000 667 0.25 Nonlinear elastic 976 535 0.2 stresses and pore pressures from both analysis are shown in figure 5.18. The stresses and pore pressures from the nonlinear elastic analysis are slightly lower, because the moduli values used in the nonlinear analysis are low compare to those in the linear elastic analysis. However, the stresses and pore pressures and the shapes of varia-tions are comparable because the developed strains are small in the example problem considered. Therefore, it can be concluded that for this kind of problems, analysis with linear elastic parameters is quite adequate to make any conclusion regarding the variations of stresses and displacements, if the temperature profiles resulting in steaming process and the material properties are comparable to those considered in the example. Chapter 5. Verification and Apphcation of The Program 85 symbols — nonlinear elastic lines — linear elastic Figure 5.18: Comparison of effective stresses and pore pressures for incremental elastic and linear elastic skeleton behaviour Chapter 6 Summary and Conclusions 6.1 Summary and Concluding Remarks Oil sand is a four phase material composed of solids, water, bitumen and gas. The behaviour of oil sand differs from the general soil behaviour in a number of aspects. Gas exsolution upon unloading, very low in-situ permeability and effects of tempera-ture changes are some important aspects that have to be considered wrhen analyzing problems associated with oil sands. Effects of temperature changes are important in the analysis of oil sand because, heating by of steam injection is often used in practice to enhance oil recovery. Changes in temperature will cause changes in viscosity, stresses and pore pressures and con-sequent^ in some of the engineering properties such as strength, compressibility and permeability. A finite element program C O N O I L has been developed by Byrne and Vaziri (1986) to handle the problems associated with oil sand including all important aspects. The fluid components have been treated as a homogenized fluid with an equivalent com-pressibility. The compressibility is considered as varying with void ratio, saturation and the amount of gas exsolved. The fluid flow has been assumed nonlinear by treat-ing the permeability as changing with viscosity, void ratio and saturation. Hyperbolic stress-strain relation has been considered for sand skeleton. In addition, shear-volume coupling and stress transfer technique are considered to obtain a better stress-strain response. Temperature effects have been analyzed by a method similar to that commonly 86 Chapter 6. Summary and Conclusions 87 used in thermal elasticity. Application of this method for thermal consolidation pre-dicts highly oscillating results, unless the time steps are kept very small. In this study a different method which avoids such oscillations is formulated to analyze the stress-deformation and flow problems related to temperature changes in oil sands. In this formulation additional terms are included in the stress-strain relation and in the volume constraint equation to account for temperature induced effects. The new method is incorporated in the finite element program CONOIL-I I . A n empirical correlation for changes in viscosity with the temperature and pressure proposed by Puttagunta et al. (1988) is also included in the program, instead of inputting experimentally determined factors. The important aspects of the developed model have been checked by comparing the program results with closed form solutions and laboratory results. The new for-mulation for thermal consohdation predicts very smooth results, which are in excellent agreement with closed form solutions. Finally to show the capability of the program, it has been applied to predict the responses of a oil recovery project. The results give insights into the likely behaviour and indicate that for the problem analyzed the linear and nonlinear responses are similar. 6.2 Recommendations for Further Research The following recommendations are made towards future research and possible mod-ifications to the program. 1. The analysis procedure should be applied to various recovery schemes, such as the underground test facility, where field observations are possible and predicted results can be compared with the measured response. Chapter 6. Summary and Conclusions 88 2. Incorporating different stress-strain models in the program will be useful. This will give the user a choice and allow the influence of different stress-strain models to be determined. 3. The possibilities of gas movement should be considered in the program. (a) Gas moves with the liquid; presently included. (b) Gas does not move, and (c) Gas flows separately from the liquid 4. The heat conduction theory may also be included in the program. Presently a separate program is used to obtain the temperature profiles during steam injection process. 5. In practice, a number of steam injection and production wells are employed in a oil recovery project. In such cases, specifying the boundary conditions become difficult in a two dimensional analysis. Therefore to predict realistic results the analysis technique needs to be extended to three dimensions. Bibliography [1] Aboshi, H . , Yoshikuni, H . and Maruyama, S. (1970), " Contant Loading Rate Consolidation Test", Soils and Foundation, vol. X, No. 1, pp 43-56. [2] Aboustit, B.L. , Advani, S .H. , Lee, J .K. and Sandu, R . S . (1982), " Finite Element Evaluation of Thermo-Elastic Consolidation", Issues on Rock Mech., 23rd Sym. on Rock Mech., Univ. of California, Berkely, pp 587-595. [3] Aboustit, B.L. , Advani, S .H. and Lee, J .K. (1985), " Variational Princi-ples and Finite Element Simulations for Thermo-Elastic Consolidation", Int. J. for Num. and Anal. Mtds. in Geomech., Vol. 9, pp 49-69. [4] Agar, J .G . , Mogenstern, N . R . and Scott, J .D. (1987), "Shear Strength and Stress Strain Behaviour of Athabasca Oil Sand at Elevated Temperatures and Pressures", Can. Geotech. J., No. 24, pp 1-10. [5] Agar, J .G . , Mogenstern, N . R . and Scott, J .D. (1986), " Thermal Expan-sion and Pore Pressure Generation in Oi l Sands", Can. Geotech. J., Vol. 23, pp 327-333. [6] Biot, M . A . (1941), "General Theory of Three Dimensional Consohdation", J. of App. Physics, Vol. 12, pp 155-164-[7] Biot, M . A . (1941(b)), "Consolidation Under a Rectangular Load Distribu-tion" , J. of App. Physics, Vol. 12, pp 426-430. [8] Bishop, A . W . (1973), " The Influence of an Undrained Change in Stress on The Pore Pressure in Porous Media of Low Compressibility", Geotechnique, Vol. 89 Bibliography 90 23, No. 3, pp 435-44%-[9] B i shop,A.W. and Henkel, D. J . (1957), " The Measurements of Soil Prop-erties in The Triaxial Test", 2nd Ed., William Arnold Publ., U.K. [10] Booker, J .R . (1974), " The Consolidation of a Finite Layer Subject to Surface Loading", Int. J. of Solids Structures, Vol. 10, pp 1053-65. [11] Booker, J .R . and Savvidov, C . (1985), " Consolidation Around a Point Heat Source", Int. J. for Num. and Anal. Mtds. in Geomech., Vol. 9, pp 173-184-[12] Byrne, P . M . (1983), " Static Finite Element Analysis of Soil Structure Sys-tems", Soil Mech. Series, No. 71, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. [13] Byrne, P . M . and Cheung, H . (1984), " Soil Parameters for Deformation Analysis of Sand Masses", Soil Mech. Series, No. 81, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. [14] Byrne, P . M . and Grigg, R . G . (1980), " O I L S T R E S S - A Computer Program for Analysis of Stresses and Deformations in Oi l Sand " , Soil Mech. Series, No. 42, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. [15] Byrne, P . M . and Eldridge, T . L . (1982), " A Three Parameter Dilatent Elastic Stress-Strain Model for Sand", Int. Sym. on Num. Mtds. in Geomech., Switzerland, pp 73-79. [16] Byrne, P . M . and Janzen, W . (1984), " INCOIL- A Computer Program for Nonlinear Analysis of Stress and Deformation of O i l Sand Masses", Soil Mech. Series, No. 80, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. Bibliography 91 [17] Byrne, P . M . , Vaid, Y . P . and Samerasekera, L. (1982), " Undrained De-formation Analysis Using Path Dependent Material Properties", Int. Sym. on Num. Mtds. in Geomech., Switzerland, pp 294-302. [18] Byrne, P . M . and Vaziri, H . H . (1986), " C O N O I L : A Computer Program for Nonlinear Analysis of Stress Deformation and Flow in Oi l Sand Masses", Soil Mech. Series, No. 103, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. [19] Campanella, R . G . and Mitchell, K . J . (1968). " Influence of Temperature Variations on Soil Behaviour", J. of SMFD, ASCE, Vol. 94, SM 3, pp 709-734. [20] Chang, C.S. and Duncan, J . M . (1983), "Consoliation Analysis for Partly Saturated Clay Using A n Elasto-Plastic Effective Stress-Strain Model " , Int. J. for Num. and Anal. Mtds. in Geomech., Vol. 17, pp 39-55. [21] Cheung, H . K . F . (1985), " A Stress-Strain Model for The Undrained Response of Oi l Sand", M.A.Sc Thesis, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. [22] Christian, J .T . and Boehmer,J .W. (1970), "Plane Strain Consolidation by Finite Elements", J. of SMFD, ASCE, Vol. 96, SM4, pp 1435-57. [23] Duncan, J . M . and Chang, C Y . (1970), " Nonlinear Analysis of Stress and Strain in Soils", J. of SMFD, ASCE, Vol. 96, SM 5, pp 1629-1651. [24] Duncan, J . M . , Byrne, P . M . , Wong, K.S. and Mabry, P. (1980) " Strength, Stress-Strain and Bulk Modulus Parameters for Finite Element Anal-yses of Stresses and Movements in Soil Masses", Report No. UCB/GR/78-02, Dept. of Civil Engr., Univ. of California, Berkeley. Bibliography 92 [25] Dusseault, M . B . (1980), " Sample Disturbance in Athabasca Oil Sand", J. of Can. Petroleum Tech., Vol. 19, pp 85-92. [26] Dusseault, M . B . (1979), " Undrained Volume ans Stress Change Behaviour of Unsaturated Very Dense Sands", Can. Geotech. J., Vol. 16, pp 627-640. [27] Fredlund, D . G . (1973), " Discussion: Flow and Shear Strength", Proc. of 3rd Int. Conf. on Expansive Soils, Haifa, Vol. 2, pp 71-76. [28] Fredlund, D . G . (1976), " Density and Compressibility Characteristics of A i r Water Mixtures", Can. Geotech. J., Vol. 13, pp 386-396. [29] Fredlund, D . G . (1979), " Appropriate Concepts and Technology for Unsatu-rated Soils", Can. Geotech. J., Vol. 16, pp 121-139. [30] Fredlund, D . G . and Mogenstern, N .R . (1977), " Stress State Variables for Unsaturated Soils", J. of Geotech. Eng. Div., ASCE, Vol.103, pp 447-466. [31] Ghaboussi, J. and K i m , K . J . (1982), " Analysis of Saturated and Partly Saturated Soils", Int. Sym. on Num. Mtds. in Geomech., Zurich, pp 377-390. [32] Ghaboussi, J. and Wilson, E .L . (1973), "Flow of Compressible Fluid in Porous Elastic Media" , Int. J. for Num. Mtds. in Engr., Vol. 5, pp 419-442. [33] Gibson, R . E . , Schiffmann, R .L . and P u , S.L. (1970), " Plane Strain and Axisymmetric Consohdation of a Clay Layer on a Smooth Impervious Base", Quart. J. of Mech. and Applied Math., Vol. XXIII, Part 4, pp 505-519. [34] Harris, M . C . and Sobkowicz, J .C . (1977), " Engineering Behaviour of O i l Sand", The Oil Sands of Canada Venezula, The Canadian Inst, of Mining and Metallurgy, Special Vol. 17, pp 270-281. Bibliography 93 [35] Hwang, C . T . , Mogenstern, N .R . and Murray, D . W . (1970), " Applica-tion of Finite Element Method to Consohdation Problems", Application of FEM in Geotech. Engr., Proc. of Symp., Vicksburg, Mississippi, Vol. II, pp 739-765. [36] Hughes, J . M . O . , Worth, C P . and Windle, D . (1977), " Pressuremeter Tests in Sands", Geotechnique, Vol. 27, No. 4, PP 455-477. [37] Irons, B . M . (1979), " A Frontal Solution Program for Finite Element Analy-sis", Int. J. for Num. Mtds. in Engr., Vol. 2, pp 5-32. [38] Konder, R . L . (1963), " Hyperbolic Stress-Strain Response: Cohesive Soils", J. ofSMFD, ASCE, Vol. 89, SM 1, pp 115-143. [39] Kosar, K . M . (1989)' " Geotechnical Properties of Oi l Sands and Related Strata", Ph.D Thesis, Dept. of Civil Eng., Univ. of Alberta, Edmonton, Canada. [40] Lambe, T . W . and Whitman, R . V . (1969), " Soil Mechanics " , John Wiley & Sons [41] Lewis, R . W . , Roberts, G . K . and Zienkiewicz, O . C . (1976), " A Nonlinear Flow and Deformation Analysis of Cosolidation Problems", 2nd Int. Conf. on Num. Mtds. in Geomech., Blacksburg, Virginia, Vol.2, pp 1106-1118. [42] McNamee, J. and Gibson, R . E . (1960), " Plane Strain and Axially Sym-metric Problems of Consolidation of A Semi-Infinite Clay Stratum", Quart. J. of Mech. and Appl. Maths., Vol. 13, pp 210-227. [43] Mitchel, J .K. , Hooper, O.R. and Campanella, R.J . (1965), " Permeabil-ity of Compacted Clay" , J. of SMFD, ASCE, Vol. 91, SM 4, pp 41-65. [44] Naylor, D . J . (1973), " Disscussion", Proc. of The Symp. on The Role of The Plasticity in Soil Mech., Cambridge, pp 291-294-Bibliography 94 [45] Puttagunta, V . R . , Singh, B. and Cooper, E . (1988), " A Generalized Viscosity Correlation for Alberta Heavy Oils and Bitumen", UNITAR/UNDP Int. Conf. on Heavy Crude and Tar Sands, Edmonton, Alberta, Canada. [46] Risnes, R., Bratli, R . K . and Horsrud, P. (1982), " Sand Stresses Around a Well Bore", Society of Petroleum Engineers J., pp 883-898. [47] Sandhu, R.S. (1972), " Finite Element Analysis of Consolidation and Creep" State of the Art Led., Application of FEM in Geotech. Engr., Proc. of Symp., Vicksburg, Mississippi, Vol. II, pp 697-738. [48] Schuurman, L E . (1966), " The Compressibility of an Air-Water Mixture and a Theoretical Relation Between The A i r and The Water Pressures", Geotechnique, Vol. 16, pp 269-281. [49] Scott, J .D. and Kosar, K . M . (1985), " Foundation Movements Beneath Hot Structures", Proc. of 11th Int. Conf. on Soil Mech. and Foun. Eng., San Francisco, California. [50] Scott, J .D. and Kosar, K . M . (1984), "Geotechnical Properties of Oi l Sands", . WRI-DOE Tar Sand Symp., Vail, Clorado. [51] Scott, J .D. and Kosar, K . M . (1982), " Thermal Expansion of Oi l Sands" Proc. of Forum on Subsidence Due to Fluid Withdrawal, U.S. Dept. of Energy and The Rep. of Venezula Ministry of Energy and Mines, Oklahoma, pp 46-57. [52] Singh, R. (1965), " Unsteady and Saturated Flow in Soils in Two Dimensions", Tech. Report No.54, Dept. of Civil Engr., Stanford Univ. [53] Sisler, H . H . , Vanderwarf, C . A . and Davidson, A . W . (1953), "General Chemistry - A Systematic Approach", McMillan Co., New york. Bibliography 95 [54] Sobkowicz, J. (1982), " Mechanics of Gassy Sediments", Ph.D Thesis, Dept. of Civil Eng., Univ. of Alberta, Edmonton, Canada. [55] Sobkowicz, J. and Mogenstern, N .R . (1984), " T h e Undrained Equlibrium Behaviour of Gassy Sediments", Can. Geotech. J., Vol. 21, No. 3, pp 439-448. [56] Sloan, S .W. and Randolph, M . F . (1982), <; Numerical Prediction of Col-lapse Loads Using Finite Element Methods", Int. J. for Num. and Anal. Mtds. in Geomech., Vol. 6, pp 47-76. [57] Sparks, A . D . W . (1963), " Theoretical Consideration of StressEquations for Tartly Saturated Soils", Proc. of 3rd Regional Conf. for Africa on Soil Mech. and Foun. Eng., Salisburg, Southern Rhodesia, Vol. 1, pp 215-218. [58] Taylor, R .L . (1977), Finite Element Method, by Zienkiewicz, O.C., McGraw Hill, pp 677-757. [59] Terzaghi, K . and Peck, R .B . (1948), "Soil Mechanics in Enginerring Prac-tice", John Wiley & Sons [60] Timoshenko, S. (1941), Strength of Materials, Vol. 2, Van Nostrand, New York. [61] Timoshenko, S. and Goodier, J . N . (1951), Theory of Elasticity, McGraw Hill, New York. [62] Vaziri, H . H . (1986), "Nonlinear Temperature and Consolidation Analysis of Gassy Soils", Ph.D Thesis, Dept. of Civil Eng., Univ. of British Columbia, Van-couver, Canada. [63] Wallace, M.I . (1948), "Experimental Investigation of The Effect of Degree of Saturation on The Permeability of Sands", S.M. Thesis, Dept. of Civil Engr., M.I.T., Cambridge. Bibliography 96 [64] Yan, L. (1986), " Numerical Studies of Some Aspects W i t h Pressuremeter Tests and Laterally Loaded Piles", M.A.Sc Thesis, Dept. of Civil Eng., Univ. of British Columbia, Vancouver, Canada. Appendix A Load Shedding Techniques Three types of load shedding technique are described here; assuming constant prin-cipal stress directions, assuming constant mean stress and assuming constant minor principal stress. The Mohr-Columb failure criterion is assumed here and the yield criterion can be expressed as, cr ' - i V ; - 2 c i v j (A.l) where, l - s i n ^ , * 1 + sin <f> <f) - friction angle c - cohesion intercept A . l Constant Principal Stress Directions In this approach, the correcting path ( P"P in figure A . l ) is assumed to be normal to the yield surface. The principal effective stresses cr\ and a'3 for point P can be obtained by simple vector algebra. The vector n can be expressed as, n = {k,-kN+) (A.2) where, k - an arbitrary constant 97 Appendix A. Load Shedding Techniques 98 Figure A . l : Correction of stresses normal to the yield surface Appendix A. Load Shedding Techniques 99 From figure A . l , / = rh + n K X ) = (*i">03") + (*» ~kN<i>) = K t * ) , ( « J - ^ ) ] (A.3) ( t r ^ c T j ) lies on the failure surface and therefore, should satisfy equation A . l . Then, « + k) = N+ia'; - kNJ - 2cNl (A.4) By solving this, T h en, a\ = a'-'-N^k (A.6) a'3 = a'3-+k (A.7) The cartesian components of the corrected principal stresses a[ and o~'z can be obtained as follows, a'x = S + Rcose (A.8) a'y = S - R cos 6 (A.9) T. Rsm9 (A.10) Where R and 5 are the radious and centre of the corrected Mohr circle as shown in figure A.2 and given by the following equations. Appendix A. Load Shedding Techniques 100 S = ^{a[+a'3) 9 = tan" 1 2T. ( A . l l ) (A.12) (A.13) Figure A.2 shows the Mohr circles for the stress conditions before and after cor-rections, for a plane strain case. For a three dimensional problem, the correction may-be interpolated as a correction normal to the failure plane keeping the intermediate principal stress constant. Yield Surface Figure A.2: Correct and incorrect Mohr circles for constant principal stress direction Appendix A. Load Shedding Techniques 101 A.2 Constant Mean Stress The Mohr circles for this approach before and after correction are shown in figure A.3. Since the mean stress is constant for the correcting procedure the centre of the Mohr circles will coincide. The corrected stress can be easily determined from the consideration of the geometry. Yield Surface Figure. A.3: Correct and incorrect Mohr circles for constant mean stress < = <-(T - ; y - S s i n < / > ) c o s 0 (A.14) < •= ( r ; y - S s i n e * ) c o s 0 (A.15) Txy = r£ (1 — sin#) + 5 s i n ^ s i n 0 (A.16) where, S and 6 are as defined previously by equations A.11 and A.13. Appendix A. Load Shedding Techniques 102 A.3 Constant Minor principal stress Figure A.4 represents the corrected and uncorrected Mohr circles for this approach where the minor principal stress <T'3 kept constant. Since a'3~ and a'3 are equal, the yield condition can be expressed as, Figure A.4: Correct and incorrect Mohr circles for constant minor principal stress a\ = N^+2cNl (A. 17) From the consideration of the geometry the following equations can be derived. R = sin 0(1 - N+) (A. 18) 2 Appendix A. Load Shedding Techniques 103 i 2 -Rcos$ (A. 19) a'y = S" + S"P cos 6 a'x = S~ — S"P cos 6 — <r'{ - f a[ Txy — S'Psh\9 A.4 Method for Transferring to The Adjacent Elements The previous sections explained how the corrections to a violated stress state are esti-mated with different considerations. This corrections of stresses have to be transferred to the adjacent elements to maintain overall equihbrium. Byrne and Janzen (1984) devised a suitable scheme for this purpose. W i t h this scheme, further iterations may be required in order to ensure that all the elements adequately satisfy the yield cri-terion. The procedure in the finite element analysis is described below. At the termination of every load increment the stress states which have vio-lated the yield criterion at each integration point (or optionally at each element) are recorded. The magnitude of the overstress is then determined in accordance with one of the techniques described in the previous sections for every offending stress state. Denoting these corrections by Aa'x, Aa'y and A r ^ , the next step is to calculate the equivalent nodal loads which balance the stress changes in the normal fashion. AF= f BTAa'dv Jv where, S'P - 0", ( c o s - l ) 2 + .R2 (A.20) (A.21) (A.22) (A.23) AF - equivalent nodal force vector for load transfer AcrT = { A < , Aa'y> Arxy} Appendix A. Load Shedding Techniques 104 The strains and stresses throughout the body due to these correcting factors are then computed. Additional iterations may be required in problems where the elements undergo a change in normal stress due to these correcting forces, because the computed stresses may again exceed the failure criterion. Appendix B Representation of Permeability For the problems associated with oil sands, the factors which have a significant role on permeability are the coefficient of viscosity(/z), void ratio(e) and the degree of saturation( S). Chang and Duncan (1983) proposed a semi empirical formulation for the variation of permeability with void ratio and saturation. In this study the formu-lation is extended to include the changes in viscosity as described in chapter 4.2.1. This appendix describes how the factors in the equation for permeability variation for changes in void ratio and degree of saturation are evaluated and how they are justified. The equation for changes in permeability is given by, k = k. Fe F. (B.l) where. k permeability of unsaturated soil permeability of fully saturated soil Fe factor whose value depends on void ratio F, factor whose value depends on degree of saturation factor whose value depends on viscosity e, e 0 current and initial void ratios threshold degree of saturation at which liquid begin flow freely 105 Appendix B. Representation of Permeability 106 B.l Computation of Variation in e and S Changes in void ratio(de) are related to the changes in volumetric strain(e t,) as given below. de = , v , (B.2) ( l + e 0 ) v ' ei ~ eo + de (B-3) Where subscripts 0 and 1 refer to initial and final conditions of the increment. In the derivation which follows, it is shown that changes in the degree of saturation can be obtained from a knowledge of changes in void ratio and quantity of n'ow(dQ). From definition, Vi • V Where \\ and Vv are volume of liquid and voids in an element of soil of volume V. Differentiating equation B.4 d S - g - ^ + S. (B.5) By denoting Si = S0 + dS, Equation B.5 can be written as, but, dVt dQ Vv ( V . e i ) / ( l + e0) and (B.7) Appendix B. Representation of Permeability 107 f = (B.8) Substitution of equations B.7 and B.8 in B.6 yields 5 i = rf(?(l + eo) + 5^o For a two dimensional problem V becomes the area of the element. Since in the normal course of finite element analysis both the quantity of flow and volumetric strains are evaluated, no extra effort is needed to calculate the changes in void ratio and degree of saturation. B.2 Factor Fe Lambe and Whitman (1969) have collected considerable experimental data to illus-trate the variation of permeability with void ratio. Although there is a considerable scatter in the data, they found that a linear relationship passing through the origin exists between k and e 3 / ( l + e) for quite a wide range of granular materials as shown in figure B . l . It can be argued that various other relationships could be established for the variation of k with e. However, without the need for much specific details about the soil, the relationship obtained form Lambe and Whitman (1969) is quite reasonable for most engineering purposes. In the present study Fe is therefore considered to be given by the following expression. F e-e | 7 ( l T ^ ) ( B - 1 0 ) where. e - current void ratio eo - initial void ratio Appendix B. Representation of Permeability 108 1000 Void ratio function Figure B . l : Variation of permeability with void ratio (after Lambe and Whitman, 1969) Appendix B. Representation of Permeability 109 B.3 Factor Fs There is no unique relationship observed for the variation of permeability with satu-ration. The experimental results (Wallace, 1948 ; Lambe and Whitman, 1969) shows that the permeability increase with saturation but greatly influenced by the soil fab-ric. Although the variation in dependent on the soil type, Singh (1965) came up with the following relationship for clayey soils. 0.0 and 0.5. A value of 0.0 for Sf confirms to the variation F, with 5 , implicit in Poiseiulle's equation. In the absence of other relationships, the above formulation is used in the present study as a reasonable relation. However, if there is no experimental data available to define a suitable value for F,, it may be justifiable to use Poiseiulle's equation by specifying Sf = 0. ( B . l l ) where S - current degree of saturation Sf - threshold degree of saturation at which the liquid in the pores begin to flow freely The value of Sf depends on the soil type but is generally found to be between Appendix C Closed Form Solution for 1-D Thermal Consolidation In this appendix closed form solution for one dimensional consolidation in a soil mass subjected to uniform temperature changes is derived, from the solution for constant rate of loading given by Aboshi et al. (1970). The soil is assumed to be fully saturated, the solid particles and the liquid are considered incompressible and the permeability is assumed to be constant. The temperature in the soil mass is assumed to increase linearly with time. Changes in temperature wil l cause changes in pore pressure and effective stresses. These changes can be determined by considering the volume compatibility between the soil skeleton and its components. The change in pore pressure can be considered as an equivalent effect of a change in applied load, hence the equivalent change in applied load can be determined in terms of the soil parameters and the applied temperature change. If the soil is undrained and free from constraints, the change in effective stress (ACT') caused by the temperature change (A0) will be equal and opposite to the change in pore pressure (Ap). Then, Acr' = - A p ( C . l ) The change in the pore pressure can be determined by equating the volume change in the soil skeleton to the total volume change of the soil components. Volume change in the soil skeleton due to change in change in effective stress,-110 Appendix C. Closed Form Solution for 1-D Thermal Consolidation 111 Ac' 4 ^ = — (C.2) where the volume of the soil skeleton is equal to unity. Substitution of equation C . l yields, Ae,fe = =— (C3) Since the liquid and solid particles are considered as incompressible, there will be no volume change due their compressibilities. Changes in volume due to thermal expansion are, Ae ; = -n ot\. A6 (C.4) A e , = - ( l - n ) a . A 0 ( C o ) A e a f c = - a , AO (C.6) where n is the porosity and a with subscripts denote the coefficients of thermal expansion. Volume compatibility between soil skeleton and its components requires that, Aesk = Ae , + A e , (C.7) By substituting the relevant terms, - a,A6 = -n a, A 0 - (1 - n) a . Ad (C.8) B,k By solving the equation, A p = B . f c n ( a / - a . ) A 0 (C.9) Appendix C. Closed Form Solution for 1-D Thermal ConsoHdation 112 If the soil deformation is constrained to one dimensional the constrained modulus will govern the equation C.9. Then, Ap = M n(cn - a,)A6 ( C I O ) where M is the constrained modulus. For a saturated soil application of a total stress increment will cause the equal amount of pore pressure increase under undrained conditions. Therefore the equivalent change in load ( A P ) to cause the pore pressure increase given by equation C.10 will be, A P = M n(ai - a,)&6 ( C . l l ) The closed form solution for one dimensional consolidation subjected to constant rate of loading is given by (Aboshi et al., 1970), where, L P 16 oo Tr-3 T ^ 7 i 1 m = l,3,. 1 . mit: —r sm 2H exp 2 2 \ T (C.12) p - pore pressure at distance z at time t T - time factor ^ | P - external load at time t M - constrained modulus on - coefficient of thermal expansion of liquid a , - coefficient of thermal expansion of solids By considering equations C . l l and C.12, for a soil mass subjected to uniform temperature increasing linearly with time, under one dimensional flow conditions, the pore pressure at distance z from the drained boundary at time t wil l be given by, Appendix C. Closed Form Solution for 1-D Thermal Consolidation 113 16 M n A0(a, - a.) ~ P = zx^ 1, TT3 T 1 . rrnrz sin < 1 — exp m = l,3,. 2H ' 7. 2 \ T 771 7T \ (C.13) where, Ad - change in temperature at time t n - porosity Hence, a closed form solution for one dimensional consolidation of a soil mass subjected to uniform temperature increasing linearly with time is derived as given by equation C.13. Appendix D Back Calculation of Coefficients of Thermal Expansion The change in volume of a soil mass subjected to temperature, changes under drained condition is given by Avak = v.k a, Ad ( D . l ) where, vak - volume of the soil skeleton a3 - coefficient of thermal expansion of solids Then, the volumetric strain will be 'Aev = ^ = a , A 6 (D.2) The drained test results shows the volumetric strain increases almost linearly with temperature. By considering the results from 2 5 ° C to 175°C (figure D.l(a)) and substituting the values in equation D.2, 0.006 = as x 150 a, = 4 x 10~h/°C (D.3) The total volume change of a saturated soil mass subjected to temperature changes under undrained conditions with constant effective stress is given by 114 Appendix D. Back Calculation of Coefficients of Thermal Expansion 115 Avlk = on v,A6 + a v,A6 (D.4) where, vi - volume of liquid v, - volume of solids a; - coefficient of thermal expansion of liquid The volumetric strain will be given by, Aev = [n cti + (1 - n)aa] A6 (D.5) Where n is the porosity. B y considering the undrained thermal expansion from 25°C to 150°C (figure D.l(b)) and substituting the corresponding values in the above equation, 0.025 = [0.375 a, + (1 - 0.375)4 x 10 - 5 ] 100 a, - 5.9 x 1 0 ~ 4 / ° C (D.6) It should be noted down that the value of a; is not constant with temperature but increasing slightly up to 200°C and then increasing considerably thereafter. endix D. Back Calculation of Coefficients of Thermal Expansion cn c o sz ( J o > > '-*— o E o I-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1-• UFC0S3 g . ^ 2.3 MPo (Cycled) o UFC0S4 q . ' = 2.5 MPa « UFC0S5 ai '= 2.5 MPa Quartz Grains 18 hrs —I <—I 1—1 100 125 150 175 200 Temperature (°C) 300 a) Drained test results 16-3 o 14 12 10 0) cn c o -C o V 8-> 6-(D > 3 2-Drained Tests ° UFC0S1 ai'= 2.5 MPa (Cycled) " UFC0S3 o,'= 2.3 MPa (Cycled) • UFC0S4 o>'= 2.5 MPa « UFC0S5 oC= 2.5 MPa Undrained Test —1 1 1 1 1 1 ) > 1 . 1 1 1 r-25 50 75 100 125 150 175 Temperature (°C) i 200 — i — 1 — i — • — i — i — r 225 250 275 300 b) Undrained test results Figure D . l : Thermal volume change of U T F oil sands (after Kosar, 1989)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Stress, deformation and flow analysis of oil sand masses...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Stress, deformation and flow analysis of oil sand masses under applied load and temperature changes Srithar, Thillaikanagasabai 1989
pdf
Page Metadata
Item Metadata
Title | Stress, deformation and flow analysis of oil sand masses under applied load and temperature changes |
Creator |
Srithar, Thillaikanagasabai |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | Stress and temperature changes in oil sand masses associated with oil recovery can give rise to stability and deformation problems. An analytical formulation is developed to analyze such kind of problems and incorporated in the finite element program 'CONOIL II’. Oil sand is a multiphase material, which can be analyzed as an equivalent two phase material; solid grains saturated with a single fluid phase, having an equivalent compressibility. The equivalent compressibility is obtained by combining the concepts of unsaturated soils with the gas laws and the gas-liquid interaction effects. Classical approach, as used in thermal elasticity is not appropriate to analyze the temperature induced effects, unless the time steps are made very small. A better formulation is presented, which involves introducing additional terms in the stress-strain relation and in the flow continuity equation, to account for temperature induced effects. The new approach reduces the oscillations predicted by the classical approach and leads to more realistic results. The analytical equations for the coupled stress, deformation and flow problems have been solved by finite element analysis. The finite element formulation involves nonlinear variations of stress-strain response, compressibility and flow, and is performed in an incremental manner. The developed finite element program has been verified by comparing its results with the closed form solutions and laboratory data. Then, the program has been applied to predict the responses associated with an oil recovery scheme. Such kind of analyses are important in the rational design of oil recovery schemes in oil sands. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062814 |
URI | http://hdl.handle.net/2429/28067 |
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 |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1989_A7 S64.pdf [ 4.56MB ]
- Metadata
- JSON: 831-1.0062814.json
- JSON-LD: 831-1.0062814-ld.json
- RDF/XML (Pretty): 831-1.0062814-rdf.xml
- RDF/JSON: 831-1.0062814-rdf.json
- Turtle: 831-1.0062814-turtle.txt
- N-Triples: 831-1.0062814-rdf-ntriples.txt
- Original Record: 831-1.0062814-source.json
- Full Text
- 831-1.0062814-fulltext.txt
- Citation
- 831-1.0062814.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-0062814/manifest