Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Finite element analysis of fluid induced fracture behaviour in oilsand Atukorala, Upul Dhananath 1983

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

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1983_A7 A88.pdf [ 6.09MB ]
Metadata
JSON: 831-1.0062766.json
JSON-LD: 831-1.0062766-ld.json
RDF/XML (Pretty): 831-1.0062766-rdf.xml
RDF/JSON: 831-1.0062766-rdf.json
Turtle: 831-1.0062766-turtle.txt
N-Triples: 831-1.0062766-rdf-ntriples.txt
Original Record: 831-1.0062766-source.json
Full Text
831-1.0062766-fulltext.txt
Citation
831-1.0062766.ris

Full Text

FINITE ELEMENT ANALYSIS OF FLUID INDUCED FRACTURE BEHAVIOUR IN OILSAND By UPUL DHANANATH ATUKORALA B.Sc. (Engineering), University of Peradeniya, Sri Lanka, 1981 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of C i v i l Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA JULY 1983 © UPUL DHANANATH ATUKORALA, 1983 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission for extensive copying of t h i s t h e s i s for scholarly purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of C i v i l E n g i n e e r i n g The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 I5th J u l y 1983 DE-6 (3/81) ii ABSTRACT In-situ recovery of o i l from Oilsand deposits at depth, rely on fl u i d induced fracturing of Oilsand. A method of analysis for predicting the i n i t i a t i o n and propagation of fractures in Oilsand has been developed. The prediction of fracture i n i t i a t i o n involves a stress analysis of the domain. Analysis of fracture propagation requires a coupled stress and flow analysis'. In the method proposed herein, the stress and flow analyses are f i r s t considered as two separate analyses that are later coupled through volume compatibility. For the f l u i d flow analysis, the fractures are appro-ximated as parabolic in shape. Fracture propagation has been analyzed incrementally, the direction of propagation being perpendicular to the computed minor principal stress direction at the tip of the fracture. The results are in good agreement with both laboratory and f i e l d exper-ience. The results indicate that the presence of a propagating fracture causes significant changes in both the magnitude and orientation of the minor principal stress, at the tip of the fracture. As a result, a path of propagation that is significantly different from the i n i t i a l in-situ minor principal stress plane direction has been observed. The rate of pumping of f l u i d , i s an important factor in predicting the propagation behaviour of fractures in Oilsand. In fracture operations carried out at shallow depths, that involve high rates of pumping, the fractures must be dominently horizontal in orientation. This is a result of the increase in f l u i d carrying capacity of the horizontal fractures, observed with propagation. TABLE OF CONTENTS Abstract i i Table of Contents i i i List of Tables v List of Figures vi Acknowledgements v i i i CHAPTER-1 INTRODUCTION 1 .1 Introduction 1 1.2 The concept of fluid fracture 2 1.3 The interpretation of flu i d fracture data 3 1 .4 The scope 6 1.5 The organization of the thesis 7 CHAPTER-2 REVIEW OF PREVIOUS WORK 2.1 Introduction 9 2.2 Fracture initiation 9 2.3 Fracture propagation 2.3.1 Theories of propagation 17 2.3.2 Path of propagation 23 2.3.3 Fluid flow aspects 26 2.4 Ground displacement predictions 29 CHAPTER-3 METHOD OF ANALYSIS 3.1 Introduction 33 3.2 Criterion for fracture initiation 38 3.3 The fini t e element representation of a fracture 41 3.4 Criterion for fracture propagation 42 3.5 The proposed method of stress analysis 45 3.6 The proposed method of fl u i d flow analysis . 3.6.1 Introduction 50 3.6.2 Fluid flow formulations 56 CHAPTER-4 FINITE ELEMENT FORMULATIONS 4.1 Introduction 66 4.2 The axisymmetric formulation 4.2.1 The Constitutive matrix[D] 66 4.2.2 The Strain displacement matrix [B] 68 4.2.3 The Stiffness matrix [K] 70 4.3 The plane strain formulation 4.3.1 The Constitutive matrix [D] 72 4.3.2 The Strain displacement matrix [B] 74 4.3.3 The Stiffness matrix [K] 75 CHAPTER-5 RESULTS AND DISCUSSIONS 5.1 RUSULTS 5.1.1 Fracture initiation 5.1.1(a) A comparison with theoretical results 77 5.1.1(b) A comparison with f i e l d data 79 5.1.2 Fracture propagation 5.1.2(a) A comparison with laboratory observations.83 5.1.3 The results from the fluid flow analysis 89 5.1.4 A comparison of displacements with Sun's closed form solution 93 i v 5.2 DISCUSSIONS 5.2.1 Fracture initiation 101 5.2.2 Fracture propagation 102 5.2.3 The flu i d flow aspects 103 5.2.4 Shapes of the fractures 112 CHAPTER-6 SUMMARY AND CONCLUSIONS -. 116 CHAPTER-7 RECOMMENDATIONS FOR FURTHER RESEARCH 120 CHAPTER-8 REFERENCES 122 Appendix-A 126 Appendix-B 1 30 Appendix-C 142 Appendix-D 149 V LIST OF TABLES Table Page 2.1 The c r i t i c a l [p -po ]/av values for different soils 16 5.1 A comparison of initiation pressure/overburden pressure ratio 77 5.2 Data used in the analyses i 94 5.3 The t i p fluid pressure data with length 104 LIST OF FIGURES Figure Page 1.1 Idealized injection pressure rate histories 4 2.1 The different modes of fracturing 18 2.2 The cross section of the fracture considered for Sun's solution 30 3.1 The plane strain domain considered for vertical fractures 34 3.2 The axisymmetric domain considered for horizontal fractures 35 3.3 A macroscopic view of the region in the vicinity of the pressurized zone 39 3.4 The variation of the stiffness of a fracture segment(in the direction of the minor axis) with the pore f l u i d pressure 43 3.5 The f i n i t e element mesh in the vicinity of the pressurized zone 44 3.6 Mesh reorientation 47 3.7 The displacements of nodes of the fracture face 49 3.8 The parabolic approximation of the shape of a fracture 54 3.9 The approximated parabolic shape vs f i e l d observations of shape 55 3.10 The quasi-static fracture 58 3.11 A typical cross section of a fracture tip 63 3.12 The quasi-static fracture with a blunt tip 64 5.1 The injection histories for Oilsand at Coldlake area 81 5.2 The injection history for Oilsand at Athabasca region - 82 5.3 Laboratory fracture test conditions 85 5.4 The cross section of the slice-6 indicating a vertical fracture 85 5.5 The pressure flow rate variations for different fracture fluids 86 5.6 The predicted path of propagation for a horizontal fracture 88 5.7 The flow rate-length-propagation velocity variations 91 5.8 A comparison of the predicted and assumed pressure profiles 95 5.9 A comparison of the vertical ground displacements predicted with Sun's(43) closed form solution 96 5.10 A comparison of the horizontal ground displacements predicted with Sun's(43) closed form solution 97 5.11 The finite xelement mesh in the vicinity of a horizontal fracture 99 5.12 A comparison of the predicted fracture shape with Sun's closed form solution 100 5.13 Results of the flow analysis with identical B-L curves for both fractures 107 5.14 The assumed B-L curves 108 5.15 The velocity-length predictions for the B-L curves shown in fog(5.14) ....108 5.16 The variation of aperture with length 110 v i i 5.17 The variation of wellhead pressure with rate of flow and time 113 5.18 The shape of the fracture with uniform and non-uniform pressure profiles 114 ACKNOWLEDGEMENTS The author would like to express his gratitude to his supervisor Dr. Peter M. Byrne for his guidance, valuable suggestions and the encouragement throughout this research. The author is indebted to Dr. Donald L. Anderson for his continued interest and advice particularly in the developement of the Flow Model. The valuable information provided by Dr. 0. Hungr are appreciated. The author is thankful to Mr. Thomas Vernon for the suggestions and help in presenting this thesis. Finally, the Research Assistantship awarded by the Department of C i v i l Engineering, University of British Columbia during the years 1981-1983 is gratefully acknowledged. 1 CHAPTER 1 INTRODUCTION 1 .1 Introduction Fluid fracture is an important concept in in-situ recovery of o i l from Oilsand and has been considered in detail by a number of researchers including Byrne et al(6), Dusseault(16,17,18), Dusseault & Simmons(l9) and Settari & Raisbeck(37). Oilsand is comprised of a dense sand matrix with its pore spaces f i l l e d with bitumen, water and free or dissolved gasses(7). The presence of bitumen reduces the effective permeability of the material and the presence of dissolved gasses makes sampling and testing problems unique to this material(3). Oilsand, at zero stresses, possess negligible tensile strength(17). Hence, virtually no energy is expended in creating an open mode fracture in such a material. The material behaviour of Oilsand is complicated being non-linear, inelastic and stress level dependent(6). Fracture Mechanics analyses traditionally deal with materials of considerable tensile strength and the stress-strain behaviour is often idealized as linear-elastic. Thus, the applicability of fracture initiation and propagation c r i t e r i a based on Linear Elastic Fracture Mechanics(LEFM) theories, to the fluid fracture initiation and propagation in Oilsand, is questionable. Fluid induced fracture propagation is associated with the rapid transmission of flu i d inside the open fracture. Flow of fl u i d is accompanied by a drop in fl u i d pressure due to viscous dissipation of energy. The resulting fluid pressure distribution 2 is dependent on the size of the opening created and vice versa. Hence, the prediction of f l u i d fracture behaviour involves a coupled stress and fluid flow analysis. The purpose of this thesis is to develope a method for predicting initiation and propagation of fluid fracturing in Oilsand, where fracture propagation is a coupled problem of stress and f l u i d flow analyses. 1.2 The concept of fluid fracturing and i t s objectives The concept of generating fluid fractures in s o i l or rock by injecting fluids at high pressures and rates of flow, often referred to as "Hydraulic Fracturing" in the literature, has been recognized by the Petroleum Industry since 1947(30). In i t s simplest form f l u i d fracturing consists of sealing off a section of a wellbore and injecting fluid, at sufficiently high pressure .and rate, to overcome the in-situ stresses and the tensile strength of the formation until an open fracture is created at the wellbore. This fracture is then extended by continued injection of the fracturing flu i d . The prime objective of such a process is to either enhance the effective reservoir permeability or create paths for introducing steam or air for heating the viscous hydrocarbons. In addition to the in-situ o i l recovery processes, fl u i d fracturing has also been employed in underground waste disposal and geothermal energy recovery programs(26). 3 1.3 The interpretation of fluid induced fracture data. In general, the variation of injection pressure measured at the ground surface and the volume of fl u i d pumped in, with time, are recorded throughout the entire fracture operation. In open hole fluid fracture operations, sealing off a section of the wellbore is achieved by using a straddle packer configuration. This consists of two rubber sealing elements, at a set distance apart, which are inflated to isolate the zone of pressurization from the rest of the hole. Perforated steel casings which contain compression packers to permit sealing off a section are also being used. It is commonly assumed that the pressure loss due to flow of fluid and leakages inside the wellbore is balanced by the gain in pressure head with depth(37). Based on this assumption, the injection pressure measured at the ground surface is essentially the same pressure at the sealed off section, at some depth. In some . fracture operations ground displacements are also measured to aid in the estimation of the extent and the orientation of a created fracture. The variation of injection pressure with time provides a preliminary means of interpretation of the initiat i o n , propagation and the orientation of a fracture. Two such idealized curves obtained for br i t t l e rock and Oilsand are shown in fig(1.1). The i n i t i a l abrupt drop in the injection pressure concomitant with an increase in the rate of flow defines the Breakdown pressure(BDP) of the formation, or the fracture initiation pressure. It is believed that the pressure at which a Fig (1.1) Idealized injection pressure - rate histories for (a) brittle rock (b) oilsand 5 fl u i d f i l l e d fracture closes, under a decreasing state of pressure, nearly defines the minor principal total stress across the fracture faces. The injection pressure observed after an instantaneous shut down of pumps(i.e.zero flow rate), is defined as the Instantaneous Shut-in pressure(ISP). Under above conditions, the fluid pressure in the fracture just balances the stresses forcing the closure of the fracture and hence is a measure of the magnitude of the minor principal total stress. The difference in magnitude of the Breakdown pressure and the Instantaneous Shut-in pressure, under low rates of flow and with low viscosity fluids, is considered to be the tensile strength of the formation. Any pressure higher than the Instantaneous Shut-in pressure causes a fl u i d pressure greater than the stresses forcing the closure, and fracture propagation takes place. This is defined. as the Fracture Propagation pressure(FPP). The various pressures defined above are illustrated in f i g ( 1 . 1 ) together with the flow rate variations with time to provide a better understanding. It has been shown both by laboratory and f i e l d experiments that f l u i d induced fractures propagate in planes perpendicular to the in-situ minor principal stress direction(4,11,31,36). Other criteria(28,39) based on energy considerations have also been proposed . for the direction of propagation and will be discussed in Chapter 2. Assuming that fractures propagate in planes perpendicular to the in-situ minor principal stress direction, and considering 6 an idealized vertical-horizontal principal stress orientation, the fracture orientation can be identified by making a comparison ' of the Instantaneous Shut-in pressure with the overburden pressure corresponding to the depth of the fracture. Hence, i f , ISP £ oy fractures are horizontal i f , ISP < av fractures are vertical where, av is the overburden pressure corresponding to the depth of fracturing. In the absence of ISP data, FPP data are used for the above purpose. 1.4 The scope of the thesis The modelling of fluid induced fracture behaviour, as mentioned in Section 1.1, requires a coupled stress and fluid flow analysis. The analysis can be simplified, by f i r s t separating the fluid flow analysis from the stress analysis. The propagation of a fluid f i l l e d fracture is associated with flow of fluid to f i l l the newly created volume. Thus, at any given time, with appropriate assumptions regarding the leakage of fluid to the surrounding region, the volume of fluid pumped in is seen to be directly related to the volume of the cavity. By imposing this compatibility condition the stress and the flow analyses are coupled together. The stress analysis is performed using a non-linear, incremental elastic finite element program. Fracture propagation 7 is analysed incrementally, taking account of the changes in orientation of the fractures during propagation. The path of propagation of a fracture, in a given domain with given i n i t i a l stress conditions, is not known beforehand. For this reason, the analyses are performed with a mesh reorientation process that allow for the changes in orientation of fractures. The mesh reorientation process is carried out at the end of each load increment. For the fl u i d flow analysis the quasi-static fracture is assumed to be represented by a series of parallel plate segments. Between each segment the flow is described by parallel plate flow equations. The fl u i d flow analysis predicts the variation of velocity of propagation of the fracture with the rate of flow of fl u i d inside the fracture and the resulting fluid pressure distribution, at quasi-static(i.e. equilibrium) fracture lengths. 1.5 Organization of the thesis. The thesis consists of seven chapters each of which serves a selected purpose. In Chapter 2, a review of previous work is given with special consideration of the fracture c r i t e r i a , problem idealizations, numerical modelling techniques and the effects of flu i d flow on propagation behaviour. The proposed methods of Fracture Mechanics analysis and the Fluid flow analysis are presented in Chapter 3. 8 Chapter 4 summarizes the mathematical formulations used in the development of the finite element program. The predictions of the fin i t e element fracture model are compared to existing theoretical solutions and laboratory and f i e l d observations in Chapter 5. A summary of the work together with the major conclusions is presented in Chapter 6. Recommendations for further research work are outlined in Chapter 7. 9 CHAPTER 2 REVIEW OF PREVIOUS WORK 2.1 Introduction Since the introduction of the concept of fluid fracture in 1947, a considerable amount of research has been done on this subject. Classical fracture mechanics consider materials that are strong in tension. Rock is strong in tension whereas soil is not. As a result, a greater portion of the research work on flu i d induced fracture behaviour is appropriate to rock rather than s o i l . Nevertheless, i t is of interest to examine the basic fracture c r i t e r i a , problem idealizations and the numerical techniques that are employed in the analysis of fl u i d fracture in rock. Fluid fracture propagation modelled using coupled stress and f l u i d flow analyses are briefly reviewed. The change in orientation of fractures, which has lead to injection pressure time records that are d i f f i c u l t to interpret, are examined. Further, the theoretical models developed for predicting the ground displacements are also reviewed. 2.2 Fracture initiation An extensive analysis on fracture initiation pressures including theoretical formulations as well as experimental observations has been reported by Haimson and Fairhurst(20). Rock is idealized as a b r i t t l e , elastic, homogeneous, isotropic linear and porous material. The theoretical formulations are based on the analogy that exists between the elasticity of a porous material and thermoelasticity. Fluid flow through the 10 pores is assumed to take place according to Darcy's law. A distinction between the penetrating and non-penetrating fracture fluids has been incorporated into the theoretical formulations. In the case of penetrating fluids, in the vicinity of the wellbore wall, a pore fluid pressure that is equal in magnitude to the pressure at the wellbore face, p^, is assumed. For non-penetrating fluids, a pore fluid pressure equal to the i n i t i a l in-situ pore pressure, p 0, which is different from the f l u i d pressure at the wellbore face, is assumed. The effective normal stress (compressive taken positive), for the above two cases, is expressed as; for non-penetrating fluids; a'= a - po 2.1 for penetrating fluids; CT'= a - p b 2.2 where, a'= effective normal stress a = total normal stress If the tensile stress at rupture is denoted by -o|, fracture initiation occurs when; a'<c -at' 2.3 Therefore, i f aQ'$ a vertical fracture results and i f -o{ a horizontal fracture results. where, o'Q is the effective circumferential stress and oz' is the effective a x i a K i . e . vertical) stress in the vicinity of the wellbore face. 11 The c r i t i c a l fracture initiation pressures derived from the stress conditions for a cylindrical hole in an infinite medium, and for the fluid conditions described by equations 2.1 & 2.2 are as follows; (a) For vertical fracturing with penetrating fluids; pP = at + 3 a2 ~ ai -p0 2.4 [2 - a{(l-2*)/(1 -*v )} ] and with non-penetrating fluids; p"P= at + 3 az - a, -p 0 2.5 Note that a l l symbols are defined after equation 2.8. (b) For horizontal fracturing away from the ends of the hole, with penetrating fluids; Pc = gt + P 3 ~Po 2.6 [\-a{(\-2v)/{\-v)}] For non-penetrating fluids the pore fluid pressure everywhere is equal to the i n i t i a l in-situ pore fluid pressure. Thus, the application of a wellbore pressure does not influence the magnitude of the axial effective stress. As reported by Haimson & Fairhurst(20), fracture initiation for this case is possible only i f the wellbore wall is prefractured or notched. (c) For horizontal fracture initiation near the ends of a pressurized hole with solid packers and penetrating fluids; D V „ _ Pc = °\ + " 3 "Po 2.7 [ l . 9 4 - o { ( l - 2 t v ) / ( l - r v ) } ] and with non-penetrating fluids; 12 c = "Po 2.8 .94 where, Po = i n i t i a l pore fluid pressure = wellbore pressure c r i t i c a l initiation pressure with penetrating fluids . P C K = c r i t i c a l initiation pressure with non-penetrating fluids = maximum horizontal principal total stresses(comp +ve) a2 = minimum horizontal principal total stresses(comp +ve) o3 = vertical principal total stress(comp.+ve) = the fl u i d pressure p nP in equation 2.5, of a thick walled hollow cylinder under no lateral loads(i.e. fractured vertically) aj = the fluid pressure pcnP in equation 2.8, of a thick walled hollow cylinder under horizontal lateral loads and a negligible vertical load(i.e. fractured horizontally). a = Biot's constant for a poro-elastic material. This is defined as [1-material matrix compressibility/ material bulk compressibility]. v = Poisson's ratio The c r i t i c a l fracture initiation pressures obtained for penetrating fluids are reported to be lower than those obtained for non-penetrating fluids, in theoretical predictions as well as experimental observations. In wellbores sealed with rubber packers the initiated fractures are found to be vertical 13 irrespective of the stress conditions. Horizontal fracture initiation is reported to be possible only with rigid packers. The f l u i d pressure that acts on the packer faces results in an axial load which is transmitted to rock as a boundary load at the wellbore face. With rigid packers that are in contact with the wellbore wall, this axial load is transmitted to rock as a shear load. Rubber is a liquid-solid elastomer. Rubber packers, under the above axial load would apply a similar radial load to rock. The application of a radial load would have a negative effect on the stress concentrations in the vertical direction and the chances of horizontal fracture initiation would thus be minimized(20). The placement of a flexible rubber liner at the wellbore face prevents the penetration of f l u i d into the material. Fluid fracture experiments carried out on rock samples with and without rubber liners are reported by Schonfeldt & Fairhurst(36). The magnitude of the fluid pressure required for fracturing of the lined samples are about 100% higher than those for the unlined samples. Furthermore, the magnitude of the initiation pressure for the unlined samples increased with increasing rate of pressurization. Based on the above findings, the authors suggest an initiation mechanism that is strongly dependent on fl u i d intrusion into the medium. Fluid intrusion into rock through a multitude of fine cracks eventually leading to a macroscopic fracture on planes perpendicular to the minor principal stress direction is proposed. Bjerrum et al(5) report the results of the f i e l d 14 permeability tests(i.e.falling head) carried out in s i l t y sands found near the Dead Sea area in Israel. The results indicate a rapid increase in the in-situ coefficient of permeability once the water pressure exceed a certain c r i t i c a l value. These high in-situ coefficients of permeability were found to be inconsistent with the type of material found in the area. The above inconsistency lead Bjerrum et al(5) to believe that the apparent increase in permeability is due to hydraulic fracture taking place in s o i l . As a result, theoretical studies were undertaken to predict safe water pressures to be used in in-situ permeability tests. The theoretical study carried out by Bjerrum et al(5) is based on the effective stress approach, the pore water pressure varying logarithmetically with the radial distance from the wellbore wall. The c r i t i c a l water pressures causing fracture have been derived for a radial plane strain domain consisting of a linear elastic material. The f i e l d permeability test involves the installation(by driving) of a piezometer which alters the state of stress around i t , the magnitude of which is material dependent. Bjerrum et al(5) take account of these stress changes by introducing two material dependent coefficients a & ^  , in equations 2.9 & 2.10. Furthermore, when the piezometer is pressurized, s o i l could deform outwards and separate from i t . Such a condition is referred to as a 'Blow-off condition. Fracturing could take place prior to, or, after 'Blow-off. If fracturing takes place after 'Blow-off, the resulting changes encountered in the state of stress and the geometrical boundary conditions of the problem would significantly influence the c r i t i c a l pressure for 15 fracturing.' For this reason, the c r i t i c a l water pressure for hydraulic fracture is considered to be the lesser of the c r i t i c a l pressures for fracturing, computed prior to and after 'Blow-off. The c r i t i c a l excess water pressures for fracturing, for the various stress conditions are as follows; (a) for vertical fracturing(i.e. minor principal stress horizontal) prior to blow-off; [p " PolA; = U\-v)/v}[(\-a)K0 + ot'/av'] 2.9 subjected to the condition {(l-»')/i>}[(1-a)Ko + o[/a^] < (1 + ^8 )K0 and for vertical fracturing after blow-off; [p " PcJAj = (1-a)[(2+ (4 -a)K0 + o{ /aj ] 2.10 subjected to the condition K 0O+p) < (1-a) [ (2-a+ f, )K0 + o[/o^] (b) for horizontal fracturing(i.e. minor principal stress vertical); [p " p 0 ] / * v =1 2.11 where, a,p = coefficients to account for stress changes due to driving of a piezometer, dependent on the so i l type. v = Poisson's ratio K0 = earth pressure coefficient(at rest) = effective overburden pressure (compressive positive) af' = effective tensile stress of so i l at rupture 16 p - p0= excess pressure for fracturing The c r i t i c a l excess pressures based on the above equations, for different s o i l types are presented in Table 2.1. The values presented in Table.2.1 are based on a Poisson's ratio of 1/3. The observed c r i t i c a l pressures from in-situ permeability tests are reported to be in the range of values predicted by the above equations for the respective s o i l types. «. & ^ in Table 2.1, correspond to the same in equations 2.9 & 2.10. Table 2.1 The c r i t i c a l [p - pc ]/<^  values for different soils Method of installation Soil type Ratio of horizontal/vertical principal effective stress, K 0-3 0-5 0-7 1-0 Ideal installation (no disturbance) —0, j8«0 All soils C4 b 0-7b b 09 1-0C High compressibility 0-4-0-5^ 0-6-0 8 ° o - s ^ - i o 0 1-0C Driven Medium compressibility 0-5-07a 0-8 - ! O c l « c l-Oc Low compressibility 0-7a-10c 1-0C 1-0° a before blow-off I b after blow-off j vertical fracture C horizontal fracture In analysing fluid induced fracture in Oilsand, Byrne et al(6) consider an initiation criterion in the light of Terzaghi's principle of effective stress. Accordingly, in a porous medium, fracture initiation occurs when the boundary fl u i d .pressure f i r s t exceeds the total normal stress on a potential fracture plane, which outcrops on the boundary, by more than the tensile strength of the formation. The analysis considers penetrating fracture fluids. Accordingly, the stress 17 conditions for fracture initiation are identical to those given by equations 2.2 & 2.3. In summary, i t can be concluded that the magnitude of the flui d induced fracture initiation pressure, in a porous material, is strongly dependent on the penetration of fluid into the medium, being significantly greater for non-penetrating media. When either high viscosity fluids or low viscosity fluids at high rates of pumping are used, causing a partial non-penetrating condition, significantly higher fracture initiation pressures are observed. The orientation of the fracture is dependent on the orientation of the minor principal stress direction. 2.3 Fracture propagation Classical fracture mechanics considers three fracture modes; a tension mode and two shear modes which are shown in fi g ( 2 . l ) . Depending on the loading conditions, a fracture could propagate as a single mode fracture or a mixed mode fracture. This behaviour results in single mode as well as mixed mode Fracture Mechanics theories. These theories have found extensive application in predicting the conditions under which fl u i d fracture propagations occur and their possible paths of propagation, in materials strong in tension. 2.3.1 Theories of propagation Fracture mechanics analyses assess the stability of a fracture, in single mode propagation, by comparing the stress 18 mode-1 mode-2 mode-3 Fig{ 2.1) The different modes of fracturing 19 intensity factors in the vicinity of the tip of a fracture with the c r i t i c a l stress intensity factor for the material. The stress intensity factors are a function of the fracture length and the distribution and magnitude of the boundary loads(34). If the stress intensity factor(for the mode) computed at the tip of the fracture is greater than or equal to the c r i t i c a l stress intensity factor for the material, the fracture propagates. If not, the fracture would not propagate. The methods of assessment of the stability of a fracture, in mixed mode propagation, are briefly examined in the following sections (a) to (c). (a) The maximum hoop tensile stress theory(28) This theory investigates the stress state near the ti p of a fracture in plane stress or plane strain conditions under mode-1 and mode-2 loading. Fracture extension of a b r i t t l e material under slowly applied plane loads are considered. Fracture propagation starts at the tip of the fracture in the plane perpendicular to the greatest tension when the hoop tensile stress reaches a c r i t i c a l material value. Accordingly, ae =j^p Cos(t9/2)[K1Cos2(t9/2)-1 .5 K 2Sin ( 6) ] + 0 [ r 2 , r 2 , etc ] 2.12 r r e = 7^7Cos(t9/2)[K1Sin(c9)+K2(3Cos(t9)-l)] + t 9 [ r 2 , r 2 , e t c 3 2.13 K! = a[ 7TL ] 2 2.14 K2 = r U L ] 2 20 where, oQ = the hoop tesile stress Tre= the shear stress 0 = fracture branch angle w.r.t the i n i t i a l orientation K,2 = mode-1 & 2 stress intensity factors r = the radial distance measured from the fracture t i p L = fracture half length The direction of propagation is obtained by setting rre =0 From equation 2.13, i t can be seen that for a pure mode-1 f racture (i .e. K2 = 0), the fracture branch angle, 0 = 0 (0=±7r is a t r i v i a l solution for the equations). This indicates propagation in the direction of the i n i t i a l orientation. For any combination of mode-1 & 2 (i.e. K;1 ,K2 ^ 0) or for pure mode-2 (i.e. K,=0), 65 is different from 0. Further, i t is seen that at the tip of the fracture(r=0), equations 2.12 and 2.13 are singular. In this case, the maximum strain energy release rate is considered as the governing parameter. Fractures propagate in a direction where the strain energy release rate is a maximum, when this maximum reaches a c r i t i c a l material value. (b) The maximum strain energy release rate theory(27) 4 9/% 1-0A [(i+3Cos 2e) + G(0) = E[3+Cos20] 1+0/7T 8Cos0Sin0 K, K2+(9-5Cos20) K 2] 2.15 K, = a[ 7rL] 2 2.16 21 K2 = T [ T T L ] 2 where, G(0) = the strain energy release rate 0 = propagation branch angle E = Young's modulus of the elastic medium K = mode-1 & 2 stress intensity factors L = fracture half length a = tensile normal stress T = shear stress Again, i t can be shown that for pure mode-1(i.e. K2=0), fracture propagation occurs without change of i n i t i a l orientation. The mathematical singularity observed in the previous theory, at the tip of the fracture, is not present in this theory. The strain energy density near the point of propagation is the governing parameter in this theory. Fracture extension occurs in the direction of minimum strain energy density, when the strain energy density factor reaches a c r i t i c a l material value. (c) The minimum strain energy density factor theory(39) S(0) = r [dW/dA] = a,,K* +2 a 1 2 K, K 2+a 2 2K 2 2.17 where, a n = [ (1+cos0) (X -cost?) ]/l6G a 1 2= sin0[2 cos0-(# -1)]/16G a 2 2= [(* +1)(1-cos0)+(1+cos0)(3 COS0-1)]/16G where, 22 S(t?) = strain energy density factor dW = the incremental strain energy dA = the area of the element r = radial distance from the fracture tip G = shear modulus ^ = (3-4f) for plane strain <& = (3-*>)/(\+v) for plane stress where v is the Poisson's ratio K 1 2 = stress intensity factors for modes-1 & 2 As in the previous theories for the special case of pure mode-1, fracture propagation is predicted in the direction of the i n i t i a l orientation. The strain energy density dW/dA, given by equation 2.17 above, becomes unbounded at r=0. From the examination of the above theories of mixed mode fracture propagation it follows that; (1) Initiation of fracture propagation and the direction of propagation are controlled by the stresses and strains in the vicinity of the tip of the fracture as opposed to the far f i e l d conditions. (2) Fracture extension without change of orientation is possible only if the loading conditions prefer a pure mode-1 fracture. 23 2.3.2 Path of propagation The path of propagation of a flu i d induced fracture has been assumed to be perpendicular to the in-situ minor principal stress direction(23,25,26,36,37,46). Laboratory (4,5) and f i e l d (11,31,36) experimental data obtained from small scale f l u i d fracture tests justify such an assumption. Dusseault(17) proposes the consideration of mixed mode fractures for the analysis of fluid fracture propagation in Oilsand. Wiles and Roegiers(45), incorporate the c r i t i c a l strain energy release rate theory for predicting the paths of propagation of fractures in rock. Pollard and Holzhausen(34) consider a two dimensional linear elastic model, and using the method of successive approximation (discussed in section 2.4) predict horizontal as well as vertical fractures, initiated at shallow depths, propagating towards the ground surface. These predictions are for materials that are strong in tension and are based on the significant increase in the normalized mode-2 stress intensity factors(i.e. K2/Ko0= rJltL/ph J~TfL= r/p b where, T=shear stress, L=fracture length/2 and p^ =wellbore pressure) observed at the tip of the fracture. For horizontal fractures a significant increase in the normalized mode-2 stress intensity factors is observed when the fracture propagates to a distance greater than about the magnitude of the depth of fracturing below the ground surface. During fracturing, with high temperature fluids(i.e. steam), heat is transferred to the material in the vicinity of the created fracture. Therefore, an analysis of the heating 24 patterns of the region of fracturing provides information on the location of the fracture in the ground. Based on the heating patterns observed by Jenkins & Kirkpatrick(1978), Dusseault(17) concludes that vertical and horizontal fractures climb towards the ground surface during propagation. Holzhausen et al(24) report a f l u i d fracturing test carried out in Athabasca Oilsand with ground displacement measurements. The ground response indicates that steam injection into the Oilsand is not a continuous process but rather is characterized by numerous episodic events. During the events the wellhead pressures dropped, injection rates increased and ground surface rose. The i n i t i a l propagation pressure based on the criterion described in Section 1.3 corresponds to that for a vertical fracture but in contrast the ground surface displacements have a symmetrical dome shaped distribution indictive of a horizontal fracture. After two weeks of steam stimulation the propagation pressure increased to that required for a horizontal fracture. The thermal expansion of the medium leading to the closure of the i n i t i a l l y created vertical fracture and an increase in the horizontal state of stress has been suggested as reasons for such observations. Dusseault(17) concludes that in massive hydraulic fracturing schemes in Oilsand, where fracture initiation favours vertical fractures, horizontal fractures would eventually prevail. The proposed hypothesis for such a behaviour is seen to be different from that suggested by Holzhausen et al(24). In Dusseault's view, the i n i t i a l fractures created are vertical. He 25 suggests that the high fluid pressures acting on the vertical fracture faces subsequently lead to an increase in the horizontal state of stress in the vicin i t y of the fracture. Once the horizontal stresses are increased beyond the vertical overburden total pressure, a horizontal fracture is initiated at the wellbore face and the injection pressure now corresponds to the overburden pressure which is higher than the i n i t i a l horizontal stress. If the rate of pumping is high enough to exceed the f l u i d carrying capacity of a fracture, an increase in the wellbore pressure results. The result would be an increase in the stresses around the fracture. Therefore, the effects of rate of pumping of f l u i d are important. However, to Dusseault, the variations in the rate of pumping are less dominant for the stress changes. The thermal strains induced by heating the reservoir, in Dusseault's view, are negligible. Dusseault and Simmons(l9) justify the above hypothesis by carrying out a two dimensional f i n i t e element analysis. The analysis considers three vertical and parallel fractures and predicts the altered state of stress due to the excess fl u i d pressures acting on the fracture faces. In summary, i t can be seen that the theoretical predictions as well as f i e l d observations indicate that, at shallow depths, horizontal and vertical fractures propagate towards the ground surface. Results from small scale laboratory experiments indicate propagation perpendicular to the in-situ minor principal stress direction. Therefore, the extrapolation of results from laboratory experiments to large scale f i e l d operations, is questionable. Furthermore, the magnitudes of the 26 injection induced stresses are seen to have a significant influence in deciding the final orientation of the fractures. 2.3.3 Fluid flow aspects An integral part of propagation of fractures is the flow of fluid within them. The resulting aperture of the fracture and the pressure distribution due to flow of fluid inside the fracture are highly inter-dependent. The stress analysis can be simplified by specifying arbitrary distributions of pressure thereby ignoring the flow of fl u i d inside the fracture. Zoback and Pollard(46) report that many authors including Secor(l975), Pollard(1976) and HSU(1975) prescribe arbitrary pressure distributions in their two dimensional models. Furthermore, Abed976), Perkins(1961), Nordgren(1972) and Kern(1969) consider the coupled fluid-flow stress-analysis problem but assume uniform pressure in the entire fracture length to obtain simplified expressions for the extent and the width of the fracture. Hagoort et al(22) and Settari and Raisbeck(37), in their two dimensional linear elastic models, incorporate mode-1 fractures e l l i p t i c in shape. The former authors analyse fracture in rock whereas the latter in Oilsand. Both models assume single phase compressible flow of fl u i d with a uniform pressure distribution inside the fracture. The propagation has been coupled with the mass balance of fluid within the fracture. In effect, this means that the fracture propagation occurs only i f there is fl u i d available to f i l l up the fracture. The model 27 developed by Settariand Raisbeck(37), predicts fracture lengths of several miles for an injection period of a few days. Such predictions are reported to be rea l i s t i c for air or cold water injection fracturing. Zoback and Pollard(46) investigate the coupled fluid flow-stress analysis problem in detail using a two dimensional(plane strain) fracture model in an infinite medium of linear elastic, homogeneous and isotropic material. The predictions are based on the mode-1 stress intensity factors computed at the tips of the fractures. The analytical solution for uniform loading of a portion of the wall of a fracture in an infinite region is employed in the prediction of the stresses, stress intensity factors and the displacements. For this purpose, the non-uniform pressure distribution is approximated as a series of uniform pressure segments and superimposed to obtain the final stresses and displacements. The fluid flow and elasticity equations are solved iteratively to obtain a consistent fracture shape-pressure variation for a prescribed wellbore pressure, viscosity and fracture length. Steady, constant property flow of a Newtonian fl u i d is assumed. For given rate of flow of fluid and viscosity, the analysis predicts the variation of the mode-1 stress intensity factors computed at the tip of the fracture with fracture length. The stress intensity factors indicated a drop in magnitude with length up to a c r i t i c a l length. Beyond this c r i t i c a l length the stress intensity factors increased rapidly. For higher fl u i d viscosities, larger c r i t i c a l lengths are observed. Based on these results the authors conclude the existence of a c r i t i c a l fracture length of stable propagation. 28 This conclusion indicates that stable fracture growth can precede breakdown when high viscosity fluids are used for fracturing and hence can lead to erroneously low estimates of o, predicted by equation 2.5, in regions where vertical fractures occur. Wiles and Roegiers(45) consider a two dimensional (plane strain) linear elastic fracture model that is appropriate for rock, using a displacement discontinuity approach. The coupled f l u i d flow-stress analysis problem is solved iteratively to obtain a consistent fracture shape-pressure distribution. Non-linear stress-strain variations are assigned to the discontinuities. Hanson et al(23) consider a numerical model that simulate the process of flu i d fracturing, in porous media, including the porous fl u i d flow and pore pressure effects in the surrounding region due to f l u i d penetration. The model consists essentially of a solution (in fi n i t e difference form) of the two dimensional, time dependent porous flow equations which numerically overlays the static two dimensional description of the elastic continuum. The analysis considers mode-1 fractures that are e l l i p t i c a l in shape. The mode-1 stress intensity factors computed at the tip regions of the fracture indicate a reduction with time when the porous flow effects are accounted for. Hence, if a fracture becomes stationary for a period of time the depressurization due to leakage of fl u i d to the surrounding low pressure regions, would reduce i t s tendency to propagate further. 29 The above review highlights the fact that the true solution for the coupled fluid flow-stress analysis problem can be attempted by an iterative type stress and flow analysis. Such analyses have been restricted to two dimensions and linear elastic materials. Only mode-1 or open mode fractures allow the necessary transmission of flu i d inside the fracture. The viscosity of the fracture fluid and the leakage of flu i d to the surrounding region are seen to be important for the propagation of a fracture. The f i l l i n g of the fracture with fluid during propagation offers a volume compatibility condition through which the stress and flow analyses are linked together. 2.4 Ground displacement predictions The mathematical closed form solution developed by Sun(42) predicts the vertical and horizontal ground surface displacements due to a horizontal , uniformly pressurized and penny shaped fracture embedded in a semi-infinite elastic continuum. The original derivations for a fracture in an infinite continuum are later modified to account for the presence of a free surface. The formulations consider the fracture length to consist of two separate regions, an inner region and an outer or an edge region as illustrated in fig(2.2). The essential difference in the two types of regions is that the fracture fluid does not penetrate into the edge region and is present only in the inner region. In effect this defines two different lengths, the fluid pressurized length of 3 0 Fig (2.2) The cross-section of the fracture considered for Sun's solution 31 the fracture and the true fracture length, the latter being always the larger. Accordingly, for a fracture with a fl u i d pressurized length, L^, and a true fracture length, L; W = [/k Sin(0/2) - d Cos(6/2)/{hfk }] A 0 2.18 U = [{L+L/k" Sin(6?/2)/A1 + {/)T L Sin(0/2)-d/T Cos(0/2) - L k Cos(0)/A 2]A o r d/L 2.19 where, the true length and the f l u i d pressurized length are related by, [1-(L^/L ) 2 ] 2 = (Ujj-ar.cD/Up 2.20 in which, j_ 1 k2 = [(r 2+d 2-L 2)/L 2+(2d/L) 2]« 6 = Cof 1[(r 2+d 2-L 2)/2 d L] A 0 = the maximum aperture= 8 (1 -v2 ) (U^-tf.djL/TrE A, = (d + L/k Cos 2(0/2)) 2 + ( L + L/k Sin(0/2)) 2 A 2 = ( d/k Cos(65/2) - L/k Sin(65/2) + L k Cos (65 ) ) 2 + ( L/k Cos(65/2) + d/k Sin(£5/2) + L k Sin ( t 5 ) ) 2 where, W = the vertical ground displacement U = the horizontal ground displacement d = the depth of the fracture from the ground surface L = the true fracture length Lp= the fluid pressurized fracture length ~X = the unit weight of the material Uj3= the flu i d pressure inside the fracture r = the radial distance. 32 E,v= the elastic parameters of the medium Pollard and Holzhausen(34), using the method of successive approximation, predict ground displacements for arbitrary fracture orientations and pressure distributions. The method of successive approximation is useful for solving problems in regions bounded by several contours ( for the present problem the free surface and the fracture boundary are considered as two contours). The analytical solutions for uniform loading of a planar surface of a semi-infinte region and for uniform loading of a portion of the fracture wall in an infinite medium are required. The two solutions are successively(iteratively) superimposed. Each superposition results in the correct boundary conditions on one of the contours but alters the boundary conditions on a l l other contours. The analysis is performed in a two dimensional(plane-strain), elastic, homogeneous and isotropic medium. The numerical model formulated by Hungr and Morgernstern(26) is related to both the Boundary Integral Equation Method and the Displacement Discontinuity Method. The approach ut i l i z e s the simple classical solutions for point loads in an elastic half-space. The fl u i d pressure is considerd as a series of point loads. The analysis is restricted to linear elastic, homogeneous and isotropic materials. The advantage of such an approach is that planar three dimensional penny shaped fractures of any arbitrary orientation subjected to arbitrary flu i d pressure variations can be analysed. 3 3 CHAPTER 3 METHOD OF ANALYSIS  3 . 1 Introduct ion The fl u i d induced fracture behaviour, in reality, is essentially a three dimensional problem. However, in a medium with an idealized vertical-horizontal principal stress orientation and for fractures initiating from the circular periphery of a vertical wellbore on planes perpendicular to its axis, the problem simplifies to permit an axisymmetric analysis. On the other hand, if fractures are to occur on planes parallel to the axis of the wellbore, a plane strain analysis can be performed on the assumption that a sufficient length of the wellbore is pressurized to permit such an analysis. The analyses carried out in this study have been restricted to the above two types of fractures. The infinite domain in which fractures occur in reality has been assumed to be represented by the finite domains shown in figs 3 . 1 & ' 3 . 2 . For the plane strain domain considered(i.e.for the vertical fracture analysis), the effect of infinity is incorporated by the inclusion of springs which represent the elastic stiffness of the medium from the outer boundary to in f i n i t y . Theoretically, i t can be shown that in a radial plane strain domain, vertical fracture initiation occurs simultaneously at diametrically opposite p o i n t s ( 1 l ) . Further, i t can be expected that the symmetry of the problem would lead to radial propagation of such a fracture. Hence, only half the domain needs to be considered. The lateral extent of the 0 0 f t . 4 0 5 5 7 o So ^ 10 o IIS 130 (£•0 l i l l l l l l l l | Kofe Fig(3-2) The axisymmetr ic domain cons ide red for ho r i zon t a l fractures 7> 36 axisymmetric domain(i.e. for the horizontal fracture analysis) is selected such that a fracture with a prescribed maximum length causes negligible changes in the state of stress of the elements in the extreme right boundary. The lower boundary represents the rockbed. The wellbore walls are assumed to be supported by smooth steel casings, the pressurized zone being perforated. The loss of pressure head due to flow of fluid through the perforations are ignored. The stress changes caused by the installation of these, casings are assumed to be negligible. Oilsand is basically plastic in shear and b r i t t l e only in tension. As discussed in Section 2.3.3, for fl u i d induced fracture propagation i t suffices to consider only an open mode fracture. Only an open mode fracture permits the rapid transmission of fluid along the fracture which is essential for the propagation of the fracture. The analysis of fluid fracture propagation requires a coupled stress and flow analysis. The above analysis can be simplified by f i r s t uncoupling the fl u i d flow analysis from the stress analysis. If the flu i d pressure distribution inside the fracture is known or assumed, the size and shape of the fracture that correspond to the known or assumed flu i d pressure distribution can be obtained from a finite element analysis. The fluid flow analysis is then coupled to the problem through volume compatibility based on the fracture shape predicted from the finite element analysis(i.e. from the stress analysis). The fractures are assumed to be symmetric about their major axes. 37 The fracture shape is approximated by a series of parallel plate segments and the flow at a point is assumed to be described by parallel plate flow equations. The fluid flow analysis in return predicts a pressure distribution which can be compared with the assumed pressure profile. The iterative type of analysis thus formulated, has been restricted to a single iteration in the present method of analysis. In a porous medium, i t is the effective stresses that govern the response of the material. The effective stress-strain relations for Oilsand are highly complex being non-linear, inelastic and stress level dependent. Further, the presence of dissolved gasses makes the estimation of effective stresses complicated. Because of the inherent d i f f i c u l t i e s associated with an effective stress approach a simpler approach based on total stresses has been taken for the present study. Linear elastic and linear elastic-perfectly plastic stress-strain relations are assigned to the unfractured and the fractured material, respectively. However, the non-linear incremental elastic f i n i t e element program developed in this thesis is capable of handling non-linear stress-strain variations. The simplification mentioned above is only the f i r s t step towards the understanding of the coupled stress-fluid flow analysis problem. The problem is further simplified by ignoring the presence of bedding features, pre-existing fractures, dilation behaviour of dense Oilsand and the thermal aspects of the problem. 38 3.2 Criterion for fracture initiation A fluid induced fracture is initiated . by sealing off a section of the wellbore, at the depth where a fracture is required, and pressurizing this sealed section with a flui d . The magnitude of the fracture initiation pressure, as concluded in Section 2.2, is highly dependent on the degree of penetration of the fracture fluid into the medium. A macroscopic view of the region in the vicinity of the wellbore face is shown in fig(3.3). The compressive normal effective stress, a', in the vicinity of the wellbore face, can be expressed as; a' = a - u 3.1 where, a is the total normal stress and u is the pore f l u i d pressure. For fully penetrating fluids; u = u p 3.1a where u^ is the fluid pressure at the wellbore face. For fully non-penetrating fluids; u = u' 3.1b where, u' is the pore f l u i d pressure in the material when subjected to the boundary f l u i d pressure,u p, under undrained conditions. For partially penetrating fluids; u' < u < U p 3.1c The state of stress of an element of s o i l , adjacent to the wellbore face, can be represented by three mutually perpendicular principal stresses, which are in general different in magnitude. If a', in equation 3.1, represents any one of these stresses, material separation occurs when, 39 boundary fluid pressure = pore fluid pressure=u - sand grains nominal boundary Fig(3-3) A macroscopic view of the region in the vicinity of the pressurized zone o 40 o' ^  -o\ 3.2 where, -a\ = effective tensile stress at rupture (comp. ,+ve). and a fracture is said to be initiated. For a given material, this f i r s t occurs when a' in equation 3.2 corresponds to the minor principal effective stress. Oilsand is a porous material exhibiting a porosity of approximately 30%. In such a medium, most of the fracture fluids of common use can be considered as fully penetrating fluids. Oilsand, at zero stresses, possess essentially zero tensile strength. Accordingly, u = u^and -a\ - 0. Substituting for u and in equations 3.1 & 3.2, fracture initiation f i r s t occurs when; or, u b ~ < W > 0 3.3 where, a m i n = total minor principal stress (comp.+ve) Fracture propagation, which wi l l be discussed later in this chapter, requires the flow of fluid into the fracture. Therefore, the opening should outcrop on the flu i d boundary. From an effective stress point of view, fracture initiation occurs when the minor principal effective stress f i r s t becomes zero on a plane outcropping the fluid boundary. From a total stress point of view, fracture initiation occurs when the boundary fl u i d pressure f i r s t exceeds the minor principal total stress on a plane outcropping the fluid boundary. 41 3.3 The Finite Element representation of a fracture Physically, fractures can be visualized as thin line discontinuities and can be discretized into a series of short line segments. Prior to fracturing, the region enclosed by a segment offers a near infinite stiffness in the direction of its minor axis. This is a direct result of the negligible thickness of the segment in comparison to its length. Once fractured, the region enclosed by the same segment is f i l l e d with the fracture f l u i d . The stiffness contribution of the fracture f l u i d in comparison to the surrounding material is negligible. Hence, after fracturing the stiffness of the segment in the direction of the minor axis is nearly zero. The f i r s t step in the Finite Element Method of analysis is the discretization of the domain into f i n i t e elements with acceptable aspect ratios (i.e. length/width ratio). In such a process the presence of thin line discontinuities, such as fractures, generate a large number of elements. Such an approach becomes cumbersome in the analysis of fracture propagation since the path of propagation is not known beforehand. Furthermore, the generation of a large number of elements increases the computation time required for the analysis. For the above reason, a fracture segment is represented by a f a i r l y thick element with an aspect ratio < 2.5. The fracture initiation criterion, as discussed in Section 3.2, can be incorporated into the f i n i t e element analysis by specifying a stiffness (in the direction of the minor axis of a fracture susceptible element) pore fluid pressure variation as 42 shown in fig (3.4) . The fluid pressure in excess of the minor principal total stress of the fracture plane leads to additional deformations in the process of fracturing. This is quantified in the finite element analysis by converting the f l u i d pressure on the fracture faces into equivalent nodal forces, and treating them as additional nodal loads for the increment. The above process is easily carried out by modifying the nodal force vector obtained from the f i n i t e element formulations(refer to Chapter 4 for details). 3.4 Criterion for fracture propagation The propagation behaviour of a fracture, which is essentially a new initiation, has been analysed incrementally. Once a fracture is initiated, the domain is analysed with the new boundary fluid pressures. Thereafter, the propagation is analysed as a new initiation but with the current stresses and the current boundary fluid pressures. The prediction of fracture propagation is accomplished by defining a quantity called the the Fracture Potential. The Fracture Potential(F.P) for the element i+1, shown in fig (3.5), is defined as; [U l > r " <W ] 3.4 the fluid pressure at the right boundary of element i the computed minor principal total stress in element i+1 F.P where, U j r a = miiA. I+I c a> E CD l/> 3 O c 9999 ( G,B) ^•000001( G ;B) prior to • f racturi ng after fracturing pore fluid pressure Fig(34) The variation of stiffness of a fracture segment (in the direction of the minor axis) with the pore fluid pressure i Ul J k FigC3.5) The f i n i t e element mesh in the vic i n i t y  of the p r e s s u r i z e d zone 45 If F.P < 0 stable fracture - does not propagate If F.P ^ 0 unstable fracture - propagates When F.P ^  0 an increment in fracture length occurs in the element which has the highest F.P, which by definition is the most fracture susceptible element. The direction of propagation is assumed to be perpendicular to the minor principal stress direction in the region of the fracture t i p . A mesh reorientation process is built into the program to trace the path of a fracture and to include the effects of previous fracturing on subsequent fracturing. The mesh reorientation process is useful only for the special case of axisymmetric fractures initiating in a domain where the in-situ minor principal stress direction is vertical. 3.5 The proposed method of stress analysis The analyses are performed using a version of the computer program NLSSIP (8), modified to incorporate axisymmetric formulations and mesh reorientation. The program uses isoparametric elements and a tangent modulus approach (see Chapter 4). The load is applied in increments starting from known i n i t i a l in-situ stress conditions. Each increment in fracture length is treated as a new load increment, and each increment is analysed twice. The stresses, strains and displacements are computed for moduli corresponding to the average stresses, at the beginning and end of an increment. The moduli corresponding to the stresses at the end of the increment are stored to be used in the beginning of the next increment of 46 loading. The stresses and strains are computed at the centre of an element. The method of stress analysis is described below. The f i r s t step in the stress analysis is to specify the in-situ state of stress for the domain under consideration. The boundary pressure in the region where a fracture is required is increased until i t exceeds the minor principal total stress of any one of the elements in the vicinity of the fl u i d boundary that leads to a fracture which outcrops on the f l u i d boundary. At this value of the boundary fluid pressure the orientation of the above element(i.e. the one with the lowest minor principal total stress) is changed to that corresponding to the direction of the current minor principal stress plane, as illustrated in fig 3.6(a). The moduli of the element are reduced to near zero values and the nodal loads that correspond to the pressure that is in excess of the minor principal total stress of the element are applied at the nodes of the element. The nodal loads are equivalent of the excess fluid pressure acting on the two fracture faces. The domain is then analysed for the displacements and the stresses and the strains of the elements are computed. Once fracture is initiated the next fracture susceptible element is obtained by computing the Fracture Potential of the elements in the vicinity of the tip of the fracture, as described earlier in Section 3.4. If the boundary flu i d pressure is high enough to cause fracture (i.e. F.P ^0), the second element is oriented in the direction of its current minor principal stress plane, as shown in fig 3.6(b). Thereafter, as for the f i r s t element, the moduli are reduced to near zero values and the excess pore fluid pressure is converted (a) min i < i 1 i 1 2 • 1 i 1 < i (b) 1 f x • -€ j ( > V 2 / i / t • nodes changed due to reorientation of element 1 o nodes changed due to reorientation of element 2 Fi g( 3.6) Mesh reorientation (a) for element 1 (b) for element 2 48 into equivalent nodal loads and added to the nodal load vector. The domain is analysed for the displacements and the new stresses and strains in the elements are computed. For subsequent elements the procedure is identical to the one described for the second element. This process is repeated until the fracture propagates to a prescribed length. As described above, fracturing of each element is considered as a separate load increment. Therefore, i f a fracture extends for 'n' number of elements there are 'n' number of load increments. The shape of the fracture is computed by summing the incremental nodal displacements, for each increment of loading, of the nodes that f a l l on the fracture boundary. A fracture extending to four elements together with the relevant nodes is shown in fig(3.7). The incremental nodal displacements mentioned above are computed after reducing the moduli of the elements 1,2,3 & 4 to near zero values and applying the equivalent nodal loads. These nodal loads correspond to a fluid pressure in excess of the current minor principal total stress, acting on the fracture faces. The nodal displacements of the two nodes that correspond to the tip of the fracture are assumed to be zero in shape calculations. The ground surface displacements are computed in a similar manner, the relevant nodes being those in the upper boundary of the domain in this case. The details of the method of analysis are presented in Appendix A in step form. 4 9 1' 2< -© . L /fracture face •—e—-1' 2' 1 2 3 U ^ ^ ^ \ . !—fracture ti • initial location of nodes o final location of nodes note: d i sp l a cemen t s of node 5 restricted to zero in shape calculations Fig (3-7) The displace ments of nodes of the fracture face 50 3.6 The proposed method of Fluid flow analysis  3.6.1 Introduction The objectives of the fluid flow analysis, as described in Section 1.4, are to predict the variations of, (a) the rate of flow of fl u i d within the fracture, (b) the fluid pressure distribution within the fracture and (c) the velocity of propagation of the fracture with the length of fracture. In this analysis the wellbore pressure and the fl u i d viscosity, for both vertical and horizontal fracture analyses, are assumed to be constant and are specified. For a given wellbore pressure and fluid viscosity, the pressure distribution within a fracture of given length and the resulting fracture shape are highly inter-dependent. The solution for this coupled problem, as concluded in Section 2.3.3, can be obtained by an iterative type stress and flow analysis. Such an approach requires assumptions regarding the i n i t i a l shape of the fracture or the i n i t i a l pressure distribution inside the fracture. The design of fluid fracture operations are often based on prescribed f l u i d pressure distributions and are referred to as geometric pressure profiles. Wiles and Roegiers(45) compare the fluid pressure distribution predicted from a coupled fluid flow-stress analysis with the existing geometric pressure profiles (in their figures 4,6 & 8). This comparison corresponds to a viscosity that is representative of water. A l l the geometric pressure profiles except the one proposed by Daneshy(l2) are seen to be significantly different from the predicted 51 distribution. Daneshy's geometric fl u i d pressure profile proposed for a 2D case, is mathematically expressed as; i Px = P,, + t P b " H M x / L f >"]T 3.5 where, P b = the wellbore pressure p_ = the far f i e l d stress p x = the pressure at a distance x measured from the wellbore wall L_p = the fl u i d pressurized half length of the fracture n = an exponent 2,4 or 8 For no leakage of fluid, the fl u i d pressure distribution given by equation 3.5 with n=4 indicates the best f i t for the pressure distribution obtained from the coupled fluid flow-stress analysis. The far f i e l d stress, p^ , in equation 3.5, is assigned the in-situ minor principal stress which inturn is the tip fluid pressure. According to the fracture c r i t e r i a used in the fin i t e element fracture model, the t i p f l u i d pressure changes with the length of propagation. Furthermore, the fluid pressurized fracture length and the actual fracture length are considered to be the same. Hence, the pressure distribution given by equation 3.5, modified for the tip fluid pressure and length of fracture, provides an i n i t i a l pressure distribution that can be assumed for the iterative analysis. Accordingly, the following modifications are performed on 52 equation 3.5. (a) replace p^ by p ^ where, p^ . is the fluid pressure at the t i p of the fracture. (b) replace by L where L is the true fracture length as discussed in Section 2.4. P> - Rap + [Pb - P 4 i p][l-(x/L)"]"^ 3.6 For a given pressure distribution, the size of the opening is dependent on the boundary conditions of the domain in which the fracture is present, the location of the fracture with respect to the boundaries, and the material properties of the domain. The size and shape of a fracture that corresponds to the pressure distribution given, by equation 3.6, at an equilibrium length, is obtained by inputting the pressure distribution into the respective fractured domains.. For a given length(L) and a constant wellbore pressure(p b), a pressure profile leading to an equilibrium fracture is obtained by varying p such that the Fracture Potential in the next, fracture susceptible element records a value just below zero. This is a t r i a l and error procedure and is repeated for various equilibrium fracture lengths to obtain the variations in size and shape with length. The shapes obtained from the fin i t e element analysis after one iteration are observed to be approximately parabolic. If the shape is approximated by a parabola the f l u i d flow formulations can be simplified considerably. Such an approach, however, limits the iterative stress and flow analysis to a single 53 iteration. A comparison of the shape predicted for a vertical fracture, with the approximated parabolic shape, is shown in fig(3.8). The shape is seen to f i t well to the parabolic variation. It was found that the fracture shape observed for short horizontal fractures(those with depth/length ratio 3-4) does not correspond very well with the parabolic shape. At larger fracture lengths, the f i t is significantly better. Fluid induced fracture operations are usually carried out at great depth below ground surface. Hence, the visual examination of the created fractures is an impossibility. For this reason, Pollard(33) suggests an analogy between the fluid fractures and sheet intrusions caused by the intrusion of liquid magma into the earth, which could be visually examined at outcrops. The thickness measurements of vertical dikes as a function of length, as reported by Pollard(33), are shown in fig(3.9). It can be seen that the f i e l d data closely approximates a parabolic shape. In the f l u i d flow formulations described in Section 3.6.2, both horizontal and vertical fractures are described as parabolic in shape. The analysis is restricted to a single iteration. Hence, the size of the opening of the fracture and the t i p f l u i d pressure correspond to the results obtained from the i n i t i a l pressure distribution given by equation 3.6. The fracture shape is approximated by a series of parallel plates and hence the flow at a point is governed by the parallel plate flow equations. Steady state laminar flow of a constant property Newtonian flu i d between each parallel segment is assumed. The parabolic approximation finite element results 0 0-4 1 2 3 distance from the wellbore axis (ft) Fig(3-8 ) The parabolic approximation for the shape of a fracture Fig( 3.9) The approximated parabolic shape vs field observations of shape obtained from Pollard ( 33 ) 56 leakage of fluid into the surrounding medium is ignored. Also, the vertical fractures created are assumed to propagate at constant height. 3.6.2 The fl u i d flow formulations In this section, the mathematical formulations involved with the fl u i d flow model are described. As outlined in Section 3.1, the symmetrical fractures considered herein are either vertical or horizontal in orientation. A propagating fracture is analysed as a quasi-static fracture. The aperture of a fracture, is dependent on a number of factors. They are, (a) material properties of the domain (b) boundary conditions of the domain (c) location of the fracture w.r.t the boundaries of the domain (d) wellbore pressure, p b (e) fracture half length, L (f) distance from the wellbore wall to the location within the fracture, where the aperture is required(i.e. x) The two domains(see figs 3.1 & 3.2) analysed are assumed to consist of linear elastic and isotropic material with similar properties. Wellbore pressures that are equal in magnitude are considered for the analysis of the two fracture types. Therefore, in a given domain, the aperture, B(see f i g . 3.10), of a fracture depends only on (e) & (f) and can be expressed as B=B[x,L]. As outlined in the latter part of Section 3.6.1, the 57 shapes obtained for the vertical and horizontal fractures, from the finite element analyses, are approximated as parabolic in shape. Therefore, as shown in Appendix C, B[x,L] can be written as, B[x,L] = B 0 [ L ] B'[x/L] 3.7 where, B 0 [ L ] , represents the maximum half aperture of the fracture and B'[x/L] is a shape function which is quadratic in x/L, representing the parabolic variation used in this analysis. It should be noted that B 0 is also dependent on the factors (a) to (e) and that i t can be written as B 0 [ L ] only with the simplifications mentioned above. Consider a fracture exhibiting a length L at time T and a length L+8L at time T+5T , as shown in figure 3.10. The equation for volume balance across section A-A, at a distance x from the wellbore axis, can be written as; for a vertical fracture, 5V = J~B 0[L+5L]B' [x/L] dx- _f B 0 [ L ] B ' [x/L] dx 3.8a x * for a horizontal fracture, L + S L L SV = 2it J"B 0[L+5L]B'[x/L]x dx-27r J B 0 [ L ] B' [ x/L ] x dx 3.8b * x where, 5V is the net volume of fluid flown across the section at A-A. As shown in Appendix C, equations 3.8 can be rearranged to take the form, 58 Fig(3.10) The quasi-static fracture 59 for a vertical fracture, 6V = 6L f,[L,x] 3.9a for a horizontal fracture, SV = SL f'jL,x] 3.9b Dividing equations 3.9 by ST and taking the limit as ST 0, we have; for a vertical fracture, cW dL f j L , x ] 3.10a 5T dT for a horizontal fracture, 2>V dL f', [L,x] 3.10b $T dT where, 3V volume flow rate across A-A 3T and dL fracture propagation velocity dT Therefore, for volume compatibility, the flow rate is linear in dL/dT. The flow of fluid between the fracture faces is accompanied by a drop in fluid pressure. Let 'p' represent the resulting pressure at a point in the parabolic fracture. Since the flow at a point within the fracture is assumed to be described by parallel plate flow equations, the gradient of the resulting pressure distribution can be written as; 60 for a vertical fracture, 3p 3V [-12 M/B 3] ax B T 3V f 2 [ B , M ] 3.11a aT for a horizontal fracture, ap av [-6 M/{TTB 3X}] 3x &T SV f'2[B,/x,x] 3.11b aT Li = absolute viscosity B = separation of the parallel plates x = distance from the wellbore axis From equations 3.11 i t can be seen [ap/dx], is linear in [ 5 V/aT]. Substituting for [aV/^T] and respect ively, for a vertical fracture, a p dL f 3 [ /x,x,L] 2»x dT for a horizontal fracture, 3>p dL f 3 [ i i , x , L ] 5x dT Integrating equations 3.12 with that the pressure gradient, from equations 3.10 & 3.7 3.12a 3.12b respect to x, the f l u i d pressure 61 distribution can be expressed as; for a vertical fracture, p(x,L,p b,ju)=P b+ d L [ A - f 0 ( x , L , M ) ] 3.13a dT for a vertical fracture, p(x,L,pb ,M) = p b + dL [A' - f'„ (x,L,M) ] 3.13b dT P b = wellbore pressure A = f,(r0,L,|x) and A' = fn(r0,L,u) are constants for a given fracture r 0 = wellbore radius From equations 3.13 i t can be seen that inorder to compute dL/dT, the following should be known. (?) . wellbore pressure (h) u, f l u i d viscosity (i) L, fracture half length (j) p , fluid pressure at a distance x(r0<x<L). Out of the above quantities, (g) & (h) are externally controlled and hence are knowns. In the following section, the procedure for estimating the quantities listed under (i) & (j) will be described. In Section 3.4, i t was outlined that fracture propagation is analysed incrementally. Within each increment a static stress analysis is carried out with the fl u i d pressure distribution given by equation 3.6. Criteria for stable fracture, as 62 described by equation 3.4, are imposed by varying the fluid pressure at the tip of the fracture until i t becomes approximately equal to the minor principal total stress at the t i p of the fracture. This involves a t r i a l and error procedure. The above process is repeated for each increment and the magnitude of the f l u i d pressure at the tip of the fracture at varying fracture lengths are obtained. From the above results and equations 3.13, dL/dT can be computed at the end of each increment in fracture length. From equations 3.10 and the dL/dT values computed at different fracture lengths, av/aT that correspond to the same fracture lengths can be computed. When x=L or X=x/L=1, equations 3.13 denote the fluid pressure at the tip of the fracture. But as shown in Appendix-C, f 4(x,L,/i) and f„(x,L,m) contain terms such as 1/(L 2-x 2), 1/(L2-x 2 ) 2 and log(L-x). Thus, for a non-trivial solution (i.e. dL/dT^O), the tip fluid pressure should be negative in f i n i t y . However, i f the fracture tip is defined differently, the above singularity can be eliminated. A typical cross section of a fracture tip is shown in fig(3.11). The material in the shaded area may be highly non-linear and visco-elastic-plastic and may not form a continuum(35). The classical fracture mechanics theories indicate stress singularities at the t i p of the fracture. However, in reality, these stresses will get redistributed and the shape of the tip may adjust i t s e l f to form a finite state of stress. The blunt fracture tip shape assumed in this study, to 63 Fig(3.11) A typical cross section of a fracture tip 64 FigC3.12) The quasi-static fracture wi th o blunt-1ip 65 eliminate the pressure singularity, is shown in fig(3.12). The formulations described above are modified to account for the new shape of the fracture. For the new shape, the tip of the fracture is defined at x = a 0L. The limits of integration in equations 3.8 would change to x - a0(L+5L) and x - a 0L. The resulting changes are as follows; Equations 3.8 5V = J B 0 [ L + 8 L ] B ' [ X / L ] dx- JB 0[L]B'[x/L] dx 3.14a SV = 27r / B 0 [L+5L]B' [ X / L ] X dx-27r JB0 [L]B' [x/L] x dx 3.14b x x Equations 3.9 6V = 6L f,[L,x,a 0] 3.15a SV = SL f,[L,x,a 0] 3.15b Equations 3.12 ^p dL f 3 [ ( i,x,L,a 0] 3.16a dx. dT dp dL £ 3[n,x,L,a 0] 3.16b Sx dT Equations 3.13 p = p, + d L [ A - f ( , ( x , a 0 , L , M ) ] 3.17a D dT p = p. + dL [ A ' - f ( x , a 0 ,L,M) ] 3.17b b dT where,the fracture tip is given by X = a0 and the fracture tip velocity is given by [d(a 0L)/dT]. Once dL/dT is known at any given length of a fracture, the pressure distribution within the fracture is given by equations 3.17. 66 CHAPTER 4  THE FINITE ELEMENT FORMULATIONS 4.1 Introduction For the two fracture types mentioned in Section 3.1, the analyses simplify to the following two classes; (a) Axisymmetric - horizontal fractures (b) Plane strain - vertical fractures The s o i l is modelled by isoparametric quadrilateral or triangular elements. The finite element formulations for the above two classes are outlined briefly in the following sections. A complete formulation is given in Appendix-B. 4.2 The Axisymmetric Formulation  4.2.1 The constitutive matrix [D] An axisymmetric problem is characterized by its symmetry with respect to both geometry and loading. Because of the symmetry, the stress components are independent of the angular (0) coordinate. Hence, a l l the derivatives with respect to 6 vanish. Accordingly, 6e 2 3 = 6e 3 2 = 0 4.1a 8e,3 = 5e3, = 0 4.1b v = 0 4.1c Subscripts 1,2 & 3 refer to r, z and 6 directions with corresponding displacements u,v and w, respectively. From the generalized Hooke's law of incremental elasticity, i t 67 follows that; 5a 2 3 = 5o3 2 = 0 4.2a 5a,3 = 6a 3 1 = 0 4.2b With the strains and stresses given in equations 4.1 & 4.2, the incremental stress vector {6a}, and the incremental strain vector {5e} reduce to; {5af = [ 6 a n 6a 22 § " 3 3 Sa, j] 4.3 T {Se} = [ 6c,, Se 2 2 oe-33 6e 1 2 ] 4.4 where, T denotes the transpose of the matrix. The constitutive relations for such a case are given by; 6a, , \-v V V 0 6e , , ha 22 E V \-v V 0 § E 2 2 6a3 3 (1+y)(1-2?) V V \-v 0 5e3 3 6a, 2 0 0 0 (1-2l0/2 8e , 2 where, E = tangent Young's modulus v = tangent Poisson's ratio Alternately the constitutive relations can be written in terms of the shear and bulk moduli as; 68 8a, , B' +G' B' -G' B' -G' 0 8e , , SO 2 2 B' -G' B' +G' B' -G' 0 5e2 2 8a3 3 B' -G' B' -G' B' +G' 0 8a, 2 0 0- 0 G' 5e , 2 where, G' = E/{2(^+^>)} and B' = 3B/{2(l+v)} in which E,B & v are the tangent Young's modulus, bulk modulus and the Poisson's ratio respectively. In matrix notation equation 4.6 is written as; {6a} =[ D ] { 8 e } 4.7 where, [D] = the constitutive matrix. 4.2.2 The strain displacement matrix [B] In compatible isoparametric elements both the geometry and the displacements are expressed by the same shape functions. The displacement u,v and the coordinates r,z at a point are approximated by; u N, 0 N2 0 N3 0 N, 0 {8} 4.8 V 0 N, 0 N2 0 N-3 0 N, r N, 0 N2 0 N3 0 Nfl 0 15'} 4.9 z 0 N, 0 N2 0 N3 0 where, {8} is the incremental nodal displacement vector g.iven by; {8} = [u, v, u 2 v 2 u 3 v 3 u, v a] and {6'} is the nodal coordinate vector for the element given 69 by; {6' }T = [r, z, r 2 z 2 r 3 z 3 r 8 z, ] and N, =(1-s)(1-t)/4 N2 =(1-s)(!+t)/4 N3 «(1+s)(1+t)/4 N4 =(1+S) d-t)/4 U j , V ; being the nodal displacements and rj ,Zj being the nodal coordinates of the element in r & z directions, respectively. s,t are the coordinates of a point in the local frame. The incremental strains for an axisymmetric problem, in terms of the displacements are expressed as; Ue-} = Se , , 6e 2 2 6e 3 3 6e , 2 au/dr dv/az u/r dv/dr + 3u/dz 4.10 substituting for u & v from equation 4.8 into equation 4.10, 3NX 3r~ 0 3N2 3r~ 0 3N3 3r 0 3N4 3r~ 0 u l v l 0 3~Z~ 0 3N2 3Z~ 0 3N3 3Z~ 0 3N4 3Z~ "2 v2 r 0 N2 r 0 ! l r 0 ^4 r 0 "3 v3 3NX 3Nj 3N2 3N2 3N3 3N3 3N4 3N4 u4 3Z 3r 3Z 3r 3Z 3r 3Z 3r v4 4-11 70 or, in matrix notation; {6e} = [B] {8 } 4.12 where [B] is called the strain displacement matrix and contains the f i r s t derivatives of N(with respect to r and z. Since the shape functions, N,', are functions of s & t, [B] is not a constant for the element. For this reason, i t is evaluated numerically at selected points within an element. The details of the selection of the Gauss quadrature for numerical integration are given in Appendix-B. 4.2.3 The stiffness matrix [K] The force displacement relationship which results in the formulation of the stiffness matrix for an element is based on the principle of virtual work. For a virtual nodal displacement vector {8e}, the resulting virtual strain increment vector, as given by equation 4.12 i s ; {Se} = [B]{5e} 4.13 The internal work done, W(int), is expressed as; W(int) = j, {5e} {So} dv 4.14a where, dv is the volume of the element. Substituting from equation 4.13, J{5e} [B] [D] [B] {8} dv 4.14b 71 The external work done, W(ext), by the external load vector {f} due to the virtual increment in displacements is written as; W(ext) = {Se} {f} 4.1 5 From the principle of virtual work, W(ext) = W(int) hence, 4. 16 (f) = J"[B]T[D] [B] {5} dv {f} = [K] {6} 4.17 4.18 where [K] =J[B][D][B]dv, which is called the stiffness matrix V for the element. If a pore f l u i d pressure of magnitude 6U is present in the element, {6a} = {6a'} + 6U 4. 19 where {8a'} is the matrix containing the element effective stresses. Substitution of equation 4.19 in 4.14 gives; W(int) = J { oeHB]~[D' 1 [B]{5}dv v ° 8U dv 4.20 where, [D'] is a new constitutive matrix in terms of effective 72 stresses. From equations 4.15,4.16 & 4.20, {f} = [K'] {6} + 6U{k'} 4.21 {f'} = {f} -6U{k'} = [K']{6} 4.22 where [K' ] = J [BHD'] [B]dv {f} = the modified nodal load vector which includes the presence of 5U. - j | 1 |dv which is a matrix associated with {k'} = I [B] T v the pore fluid pressure. 4.3 Plane strain formulation  4.3.1 The constitutive matrix [D] Plane strain problems are characterized by the following conditions; (1) no deflection in z direction (2) f i r s t derivative of the deflections in x & y directions with respect to z are zero. Accordingly, 5e 3 3 = 6e 1 3 = 6e 23 =0 4.23a w = 0 4.23b Subscripts 1,2 & 3 refer to x,y & z directions respectively, with corresponding displacements u,v,w. From the generalized Hooke's law, i t follows that; 73 6aT 3 = 6a 3 1 =0 5a 2 3 - oa 3 2 =0 4.24a 4.24b where Str-is the incremental stress with a similar notation as before. With the strains and stresses given in equations 4.23 & 4.24, the incremental stress vector {6a}, and the incremental strain vector {5e} reduce to; {6a} = [5a,, 6a2 6a 1 2] {Sef = [6e,, 6e2 2 5e , 2 ] 4.25 4.26 as; 5a, , 5a2 2 6a, 2 E ( ( 1-2v) for a plane v 0 V \-v 0 0 0 (1 - 5e , , 8£ 2 2 5e , 2 4.27 where E = tangent Young's modulus v = tangent Poisson's ratio Alternately the constitutive relations can be written in terms of the shear and bulk moduli as; {5a} = B' +G' B'-G' B'-G' B' +G' 0 0 4.28 {6} 74 where, G'= E/{2(1+j>)} and B'= 3B/{2(1+»v)} in which E,B & v are the tangent Young's modulus, bulk modulus and the Poisson's ratio respectively. In matrix notation equation 4.28 is written as; {6a} = [D]{6e} 4.29 where [D] is the Constitutive matrix. 4.3.2 The strain displacement matrix [B] For plane strain formulation, the geometry of the element(i.e. x,y) and displacements(i.e. u,v) are approximated as, x = [x1N1+x2N2+x3N3+xllNfl ] 4.30 y = [yiN1+y2N2+y3N3+y(tNfl ] u = r u 1 N 1 + u 2 N 2 + u 3 N 3 + u , N , ] 4.31 v = [v1N,+v2N2+v3N3+vl|Nl) ] where, Xj ,y- = nodal coordinates in x & y directions. Uj,vj = incremental nodal displacements in x & y directions. N i , 2 , 3 , « are identical to those given in eqn. 4.8. 75 T h e i n c r e m e n t a l s t r a i n v e c t o r {Be}, i n t e r m s o f d i s p l a c e m e n t s i s e x p r e s s e d a s ; {6e} = 3 u / d X 4.32 6 e 2 2 — 8e , 2 3u/ay + dv/dx a n d t h e r e s u l t i n g [B] m a t r i x c a n b e w r i t t e n a s , 3Nj 3N-, 3Y 0 3N2 3X~ 0 3N3 3X~ 0 3N4 3X~ 0 3Y~ 0 3N2 3~Y~ 0 3N3 3Y~ 0 3N4 3Y~ 3NX 3N2 3N2 3N3 3N 3 3N4 3N4 3X 3Y 3X 3Y 3X 3Y 3X u 3 v 3 u 4 4.33 4 . 3 . 3 T h e s t i f f n e s s m a t r i x [ K ] T h e f o r m u l a t i o n f o l l o w s e s s e n t i a l l y t h e same p r o c e d u r e a s o u t l i n e d i n S e c t i o n 4 . 2 . 3 . E q u a t i o n s 4.21 & 4.22 t a k e t h e f o r m ; { f } = [ K 1 ] {6}+ 6U{k'} { f } = { £ } - 5 U { k ' } = [ K ' J {6} 4.34 4.35 w h e r e [ K ' ] = 1 [ B H D ' ] [ B ] t d A Ck'} = [ [ B f d U t d A w h e r e , t = t h i c k n e s s o f t h e e l e m e n t i n t h e z d i r e c t i o n a n d dA = t h e a r e a o f t h e e l e m e n t S i n c e t h e s t r a i n d i s p l a c e m e n t m a t r i x i s 3 x 8 , t h e r e s u l t i n g 76 s t i f f n e s s m a t r i x w i l l be 8 x 8 . The d e t a i l s of the f o r m u l a t i o n of {f} f o r the p r e s s u r e l o a d s f o r both p l a n e s t r a i n and a x i s y m m e t r i c f o r m u l a t i o n s a r e g i v e n i n Appendix-B. The user manual and the f l o w c h a r t f o r the program a r e g i v e n i n Appendix-D. 77 CHAPTER 5 RESULTS & DISCUSSIONS  5.1 RESULTS 5.1.1 Fracture initiaton 5.1.1(a) A comparison with theoretical results The fracture initiation pressures predicted from the finite element analyses together with the theoretical fracture initiation pressures obtained from equations described in Section 2.2, are presented herein for comparison. Finite element analyses are performed for two different principal stress orientations. A state of stress which represents a minor principal stress that is horizontal in orientation(i.e. K0<1) and a state of stress which represents a minor principal stress that is vertical in orientation(i.e. K0>1) are considered. The domains considered are shown in figs(3.1) and (3.2), respectively. The results, expressed as ratios of the initiation pressure to the overburden pressure, are shown in Table 5.1. Table5.1 A c ompari son of initiation pressure/over burden pressure ratio Orientation of the f racture Results from Haimson& Fairhurst(20) Results from Bjerrum etal(5) Results from the finite element analyses away from the ends closer to the ends penetr. non-penetr penetr. non-penetr ver t ic al K = 1/2 0 0-67 (2-4) 1-00 (2-5 ) — — 0-7 0 (2-9 & 2-10) 0-50 horizontal 2-00 (2-6) — 0-6 9 (2-7) 1-06 (2-8) 1-00 (2-11) 1-0 0 note: the numbers in parentheses refer to equation numbers from which they are computed. 78 where, K0 = coefficient of lateral earth pressure. The theoretical pressures correspond to a Poisson's ratio of 1/3. The tensile strength of the material, has been assumed to be zero. For sands, the material matrix compressibility is negligible when compared to the material bulk compressibility. As a result, the Biot's constant(see Section 2.2) for sands is close to 1. Hence, in computing the c r i t i c a l pressures from equations 2.4,2.6 & 2.7, a Biot's constant of 1 has been assumed. Furthermore, the i n i t i a l in-situ pore fl u i d pressure has been assumed to be zero. For K0<1, the finite element analysis predicts the initiation of a vertical fracture, and for K0>1, a horizontal fracture. The magnitude of the initiation pressure is equal to the magnitude of the in-situ minor principal total stress that correspond to the location of the fracture in ground. For the vertical and horizontal fractures, the theoretical initiation pressures, computed for penetrating fluids, are dependent of the Poisson's ratio for the material. However, for noh-penetrating fluids, the above pressures are independent of the Poisson's ratio. It is important to note that the initiation pressures predicted in the f i n i t e element analyses, for both vertical and horizontal fractures, are also independent of the Poisson's ratio. The theoretical pressures computed for vertical fracture ini t i a t i o n , using penetrating fluids, are about 40% higher than the results obtained from the fin i t e element analysis. For horizontal fracture initiation away from the ends of the 79 wellbore and with penetrating fluids, the theoretical pressure is about 100% higher than the pressure predicted in the finite element analysis. For horizontal fracture initiation closer to the ends of the wellbore and with similar f l u i d conditions as above, the theoretical pressure is 30% lower than the finite element results. The horizontal fracture initiation pressure that corresponds to the equation given by Bjerrum et al(5) is identical to the pressure predicted in the f i n i t e element analysis. The theoretical fracture initiation pressures obtained with non-penetrating fluids are seen to be higher than those obtained with penetrating fluids. 5.1.1(b) A comparison with f i e l d data Data on actual fracture operations carried out in Oilsand, are extremely limited. Settari & Raisbeck(37) have presented results of the injection pressure, rate and time histories obtained from fracture operations in the Coldlake area in Alberta. The paper by Holzhausen et al(24) presents injection pressure-time-rate data on a fracture operation carried out in the Athabasca Oilsand region in Alberta. The fracture initiation pressures obtained from the above fracture operations are compared with the predictions of the finite element fracture model. Because of the different magnitudes of the stress conditions involved, the comparison is based on the ratio of the init i a t i o n pressure to the minor principal total stress across the fracture faces. 80 The injection histories reported in the paper by Settari & Raisbeck(37) are shown in f i g ( 5 . l ) . From these results, Dusseault(17) deduces a fracture in i t i a t i o n pressure ratio that is approximately one. The results presented by Holzhausen et al(24), shown in fig(5.2), indicate a fracture initiation pressure ratio of about 1.5, which is 50% higher than the estimated ratio. The f i n i t e element analyses predict a fracture initiation pressure ratio of one. The magnitude of the fracture in i t i a t i o n pressure for a cohesionless.material, as discussed in Section 2.2, is dependent on the degree of penetration of the f l u i d into the medium. The degree of penetration depends on the rate of pumping of f l u i d , the fluid viscosity and the porosity of the material. Thus, a proper comparison is feasible only i f the degree of penetration of the fracture fluid is similar for the two cases. For the aforementioned data, both the material and the injection fluids are almost the same whereas the pumping rates of the two studies differ. On the other hand, the material may possess some tensile strength which would result in an increased fracture initiation pressure. The finite element analyses performed assume zero tensile strength for the material and hence the initiation pressures predicted are a lower bound. It can be concluded that the magnitudes of the fracture initiation pressures predicted from the finite element analyses are only in partial agreement with the theoretical predictions as well as the limited f i e l d data. The differences between the finite element results and the field data could be due to the 82 Fig {5-2) The injection history for Oilsand at Athabasca area (adapted from Holzhausen (2 A )) 83 varying degrees of penetration of fracture fluid and the tensile strength of the formation that are not taken account of in the analyses. 5.1.2 Fracture propagation In this section, the fracture propagation results obtained from laboratory experiments are compared with the results from the f i n i t e element analyses. The f i n i t e element analyses are carried out assuming; (a) fully penetrating fracture fluids, and (b) a material of negligible tensile strength. 5.1.2(a) A comparison with laboratory observations Hydraulic fracture tests carried out on laboratory t r i a x i a l samples of a s i l t y clay are reported by Bjerrum and Anderson(4). A small piezometer (shown in fig 5.3) is installed in the centre of the sample and consolidated to a prescribed anisotropic state of stress that leaves a horizontal minor principal stress. The piezometer is then pressurized at a constant rate of flow of water mixed with fluorecene. The important observations made in these tests are; 1. The pressure in the piezometer increased to a c r i t i c a l value, remained constant for a short time, and started to decrease. In a short time the rubber membrane started bulging out and the pore pressure in the top of the specimen started to increase. These observations indicate that a fracture was formed a l l the way out to the rubber membrane. 84 2. The horizontal and vertical slices, when visually examined, did not indicate the presence of a fracture. When illuminated by ultraviolet light, the horizontal slices indicated large concentrations of fluorecene on a diametrical plane, as shown in fig (5.4). These observations indicate that the created fracture was vertical and upon removal of the excess pressure the fracture had closed completely. 3. The fracture closure pressure provides a reliable estimate of the minor principal total stress and is independent of the type of f l u i d used in fracturing. (both paraffin and water indicated identical closure pressures). These results are illustrated in f i g (5.5) as obtained from Bjerrum & Anderson(4). The finite element fracture model developed in this study has been employed to model the propagation behaviour observed in the hydraulic fracture tests in t r i a x i a l samples. For the above case, a horizontal section of the sample subjected to an identical vertical-horizontal stress ratio of 1/2 (i.e. K0=1/2) has been analysed. The analysis assumes plane strain conditions. The domain selected is shown in fig(3.1). The radial symmetry of the domain considered indicates equal probability for a fracture to initiate at any diametrically opposite pair of points. However, if the existence of a weaker plane is assumed, fracture initiation could take place in that plane. With this assumption, the finite element fracture model predicts fracture initiation on a diametrical plane. Because of the small dimensions of the specimen mentioned 85 soil sample C z = -s cr, Fig( 5-3) The laboratory f racture test condi tions (adapted from Bjer rum!- A nde rson^Aj) Fig(5-4) The cross-section of the slice-6 indicating a vertical f racture (adapted from Bjerrum <rAnrierson(4 )) — initial pore pressure •us 3 ' i . flow rate (cm / mm) Fig( 5-5 ) The pressure-flow rate  variation for different fracture fluids (adapted from Bjerrum h Anderson(4 )) 87 in the above experiment, i t could be assumed that the viscous pressure drops would have been negligible. With such an assumption the finite element analysis predicts the vertical fracture propagating to the outer boundary of the domain. Experimental results reported by Haimson & Fairhurst(20) on small rock samples and Schonfeldt & Fairhursts(36) on large f i e l d rock samples, with a minor principal stress horizontal in orientation, indicate behaviour similar to that predicted by the analysis. Since the vertical fracture analysis is restricted to a plane strain analysis, predictions regarding the climbing behaviour of vertical fractures cannot be made. The propagation behaviour of a horizontal fracture has been modelled next, using the f i n i t e element fracture model developed. The domain shown in fig(3.2) has been analysed for a vertical to horizontal stress ratio of 2 (i.e.K 0=2). The analysis assumes axisymmetric conditions. Laboratory results obtained from samples subjected to a minor principal stress that is vertical, indicate horizontal fractures propagating to the outer boundaries. As shown in fig(5.6), at larger fracture lengths the predictions of the f i n i t e element analysis deviate from the above experimental observations indicating a fracture climbing towards the ground surface. This occurs once the fracture half length exceeds about the depth of the fracture from the ground surface. The presence of a free surface makes the geometry of the problem asymmetric about a horizontal axis. As a result, the stress distribution caused by the pressure inside a horizontal opening is also ground surface elevation = 00 for a uniform pressure distribution & a weightless fracture fluid U 6 distance from the wellbore axis(ft) 10 Fi g(56 ) The predicted path of propagation for a horizontal fracture 89 asymmetric about the axis of the opening. The significant shear stresses induced in the vicinity of the tip of the fracture, due to the cause mentioned above, reorient the minor principal stress plane towards the free surface. For this reason, the analysis predicts fractures climbing towards the free surface. The theoretical work of Pollard & Holzhausen(34) and Dusseault's(17) conclusions regarding the climbing behaviour of fractures support the above predictions. In conclusion, i t can be seen that the finite element fracture model satisfactorily predicts the fluid induced fracture propagation behaviour. Field observations as well as laboratory observations are modelled correctly. 5.1.3 The results from the fluid flow analysis The results presented herein are for a specified wellbore pressure and fluid viscosity. Domains shown in figs(3.l) and (3.2) are analysed for the same isotropic stress state with similar material properties. The mesh reorientation process has been bypassed. Therefore, only fractures propagating on vertical or horizontal planes are considered. For the above idealized conditions, the analysis predicts the rate of flow of f l u i d between the fracture faces, the velocity of propagation of the t i p of the fracture and the f l u i d pressure within the parabolic fracture at equilibrium fracture lengths. These predictions are based on the tip f l u i d pressure and the magnitude of the aperture predicted from the f i n i t e element analyses and correspond to the pressure distribution at the start of the 90 f i r s t i t e r a t i o n . The v e r t i c a l f r a c t u r e s a r e assumed t o propagate at c o n s t a n t h e i g h t . The maximum r a t e of f l o w of f l u i d , Q , a c r o s s the f r a c t u r e f a c e s , o c c u r s a t the w e l l b o r e f a c e . The v a r i a t i o n of Q w i t h l e n g t h f o r t h e two f r a c t u r e s ( i . e . v e r t i c a l & h o r i z o n t a l c a s e s ) ar e shown i n f i g ( 5 . 7 ) . The v a r i a t i o n of the v e l o c i t y of p r o p a g a t i o n w i t h l e n g t h , i s a l s o shown i n the same f i g u r e . The v a r i a t i o n s of the r a t e of f l o w w i t h l e n g t h o b t a i n e d f o r the two f r a c t u r e s i n d i c a t e t o t a l l y d i f f e r e n t t r e n d s . For the h o r i z o n t a l f r a c t u r e , Q n i n c r e a s e s w i t h L w i t h an i n c r e a s i n g g r a d i e n t dQ n/dL. For the v e r t i c a l f r a c t u r e , Q v i n c r e a s e s w i t h l e n g t h but l e s s r a p i d l y . Q v, i s d i r e c t l y p r o p o r t i o n a l t o the h e i g h t assumed f o r the v e r t i c a l f r a c t u r e and the r e s u l t s shown i n f i g 5.7 c o r r e s p o n d t o a h e i g h t e q u a l t o the l e n g t h of the p r e s s u r i z e d zone of the w e l l b o r e . The Q-L c u r v e s f o r the two f r a c t u r e s i n t e r s e c t each o t h e r , the f r a c t u r e l e n g t h c o r r e s p o n d i n g t o the p o i n t of i n t e r s e c t i o n b e i n g dependent on the h e i g h t assumed f o r the v e r t i c a l f r a c t u r e . P r i o r t o i n t e r s e c t i o n , the v e r t i c a l f r a c t u r e i n d i c a t e s a h i g h e r r a t e of f l o w whereas a f t e r i n t e r s e c t i o n the r a t e of f l o w i s h i g h e r f o r the h o r i z o n t a l f r a c t u r e . The v e l o c i t y of p r o p a g a t i o n of the h o r i z o n t a l f r a c t u r e i n c r e a s e s w i t h an i n c r e a s i n g g r a d i e n t dL^/dL, whereas f o r the v e r t i c a l f r a c t u r e i t i n c r e a s e s w i t h a d e c r e a s i n g g r a d i e n t . U n l i k e the Q-L c u r v e s , the f r a c t u r e l e n g t h t h a t c o r r e s p o n d s t o the p o i n t of i n t e r s e c t i o n i s independent of t h e h e i g h t assumed f o r t h e v e r t i c a l f r a c t u r e . 91 v -ver tical fracture h-horizontal fracture 1 2 3 4 5 distance from wellbore a x i s ( L ) ft 1 a u a o 5 O o CJ o Fig (5*7 ) The flow rate-length-propagation veloci t y variations 92 The analyses carried out in this thesis correspond to a constant wellbore pressure and assume no leakage of fluid to the medium surrounding the fracture. Further, the rate of flow of fluid was allowed to vary with no fixed upper bound. However, in reality, no geotechnical material is fully impermeable and often the maximum rate of pumping is limited to the capacity of the pumps available. In fracturing of permeable materials, significant leakage of fluid takes place and the resulting fluid pressure inside the fracture can be significantly different in distribution and magnitude from that for no leakage of f l u i d . The effects of the above factors, discussed later in Section 5.2.3, are to reduce the velocity of propagation of fractures. In an isotropic stress f i e l d , the pressurization of a section of a wellbore would induce vertical and horizontal fractures simultaneously. The results shown in f i g 5.7 indicate that at i n i t i a l stages of propagation, the vertical fractures grow faster. This is a result of the larger fluid carrying capacity of vertical fractures, observed during the i n i t i a l stages of propagation. The horizontal fractures grow, but at a slower rate. However, after the vertical fractures propagate a certain distance, the increase in velocity of propagation as well as the fluid carrying capacity reduce significantly. A horizontal fracture of longer length of propagation exhibits a larger f l u i d carrying capacity than a vertical fracture. High rates of pumping of fluid force the fractures to accommodate large quantities of fl u i d . Therefore, i f only the rate of flow is considered, in fracture operations with high rates of pumping of fluid, fractures should be dominently horizontal in 93 orientation. The pressure distributions predicted for the parabolic fractures are compared with the pressure distribution assumed at the start of the f i r s t iteration, in fig(5.8). The predicted pressure distribution is seen to be in good agreement with the assumed distribution except for the region close to the tip of the fracture, where the maximum deviation is about 10% of the wellbore pressure. 5.1.4 A comparison of displacements with Sun's closed form  solution The closed form solution presented for ground surface, displacements by Sun(43), is based on highly idealized conditions. In order to permit a comparison of the numerical solution with the closed form solution, similar idealizations are made to the model developed in the present study. Thus, only uniformly pressurized penny shaped horizontal fractures embedded in a linear elastic, homogeneous and isotropic medium are considered. The radius of the wellbore is reduced to zero and no self weight stresses are considered. Further, the two fracture lengths defined in the closed form solution, and L, discussed in section 2.4, are assumed to be the same for the analysis. Figures (5.9) & (5.10) illustrate the comparisons of the horizontal and vertical ground displacement predictions from the fi n i t e element analysis with the closed form solution. The data used for these predictions are shown in Table 5.2. The closed 94 Table 5-2 data used in the analyses Young's modulus 950.0 lbf/ft* Poisson s r at io 0.20 Cohesion 0-0 0 Friction angle 40-0 depth of fracture from free surface 7.5 ft. well bore rad'i us 0-4 f t. excess fluid pressure ( uniform ) 1000.0 I b f / f t2 radius of the fracture 4-9 ft 95 ro O 3 73 l / l CU (a) (b) predict ed a s s u m e d \ length (ft) predicted assumed \\ l e n g t h ( f t ) Fiq(5-8 ) A comparison of the predicted  and a s s u m e d p r e s s u r e p r o f i l e s (a) horiz ontal fracture (b) v e r t i c a l fracture o x 0I = 0 5 10 15 radial distance (ft) Fig(5-9)A comparison of the vertical ground displacements predicted  with Sun's(43) closed form solution closed form for section 1 closed form for section 2 • finite element analysis 5 10 radial distance (ft) 15 Fig(5-10) A comparison of the horizontal ground displacements predicted with Sun's(43) closed form solution 98 form solution requires the depth of the centre of the fracture from the ground surface. As shown in the fig (5.11) because a fracture segment is represented by a fa i r l y thick element in the fini t e element analysis, the exact location of the centre of the fracture that suits the closed form solution is debatable. For this reason, the displacements are computed for the two possible locations of the centre of a fracture, as shown in fig (5.11) by sections 1 & 2. It is seen that the results are in close agreement with the closed form solution. The results agree extremely well with the closed form solution that corresponds to a depth as given by section 1. The displacements predicted are at most 13% higher than those obtained from the closed form solution when a depth that corresponds to section 2 is selected. The shape of the fracture predicted in the analysis is compared with that given by Sun's solution, in f i g (5.12). The fi n i t e element results are in excellent agreement with the closed form solution. 99 ground surface Up -' fracture face' •section 1 _^ section 2 fractured elements shaded Fig(5-11) The finite element mesh in the vicinity of a horizontal fracture 100 Fig(5-12) A comparison of the predic ted f racture shape with Sun's closed form solution 101 5.2 DISCUSSIONS 5.2.1 Fracture Initiation The state of stress of a normally consolidated (i.e.K 0<l) s o i l consists of a major principal stress which is vertical(i.e.overburden pressure) and intermediate and minor principal stresses which are horizontal in orientation. From the f i n i t e element analysis i t was seen that when K0<1, the fractures created were vertical. As pointed out in Section 1.3, the magnitude of the fluid pressure at which a fl u i d f i l l e d vertical fracture closes, corresponds to the minor principal total stress(i.e. horizontal) at the location of the fracture. If the minor principal total s t r e s s ( a 2 ) , the tensile strength of the soil(o ) and the fracture initiation pressure obtained with a non-penetrating fl u i d are known, the magnitude of the intermediate principal total s t r e s s C a , ) can be estimated from the theoretical relation given in equation 2.5. However, if a penetrating f l u i d is used, CT, can be estimated from equation 2.4, provided the Poisson's ratio and the Biot's constant for the material are also known. If overly viscous fluids or high rates of pumping of even low viscosity fluids are used, due to partial penetration of the fracture fluid , a ^ computed from the above equations would be in error. In over-consolidated soils(i.e.K 0>1), the overburden pressure corresponds to the mi.nor principal total stress. For K0>1, the fractures created are horizontal in orientation. As a result, the f l u i d induced fracture data are useful only for estimating the overburden pessure and the tensile strength in 102 over-consolidated soils. 5.2.2 Fracture Propagation The common assumption that the path of propagation is perpendicular to the in-situ minor principal stress direction is relaxed in the present analysis. Instead, the new state of stress is estimated after each increment of the fracture length. The direction of propagation is assumed to be perpendicular to the current minor principal stress direction in the tip of the fracture. Such a process reflects the effects of previous fractures on the process of subsequent fracturing. Even though Pollard & Holzhausen(34) and Dusseault(17) predict climbing fractures, their hypotheses for the climbing behaviour are different. Pollard & Holzhausen(34) predict climbing fractures based on the significant increase in the normalized mode-2 stress intensity factors(i.e. K2/K00=r/pb in which r is the shear stress and p^ is the wellbore pressure) computed at the tip of the fracture. Dusseault proposes a work minimization criterion for the path of propagation and concludes that i t is a reduction in stress and hence the lesser work required to part the material that causes the fractures to climb towards the ground surface. The finite element analysis indicates that i t is because of the reorientation of the minor principal stress plane, in the vicinity of the tip of the fracture, that causes the shallow, horizontal fractures to propagate towards the free surface. The larger normalized mode-2 stress intensity factors, 103 predicted by Pollard & Holzhausen(34), are a result of larger shear stresses induced near the tip of the fracture. These shear stresses reorient the minor principal stress plane towards the ground surface. Hence, i t is seen that the results of the f i n i t e element analyses as well as the results obtained from Pollard & Holzhausens' analysis, predict climbing fractures for the same reason. However, because of the qualitative nature of Pollard & Holzhausens' predictions regarding climbing, a quantitative comparison of the change in orientation of the principal stress plane cannot be done. 5.2.3 The fluid flow aspects The objectives of the fluid flow analyses are to predict the variations in velocity of propagation, rate of flow and the fluid pressure distribution inside a fracture, with the length of propagation. As outlined in Section 3.6.2, inorder to achieve the above objectives, i t is necessary to know the variation of the maximum aperture of fracture and the fluid pressure at the tip of the fracture, at different fracture lengths. These data were obtained by carrying out finite element analyses in the domains shown in figs 3.1 & 3.2. The tip fluid pressures obtained from the finite element analyses are shown in Table 5.3. These results correspond to the same wellbore pressure of 6000 l b s / f t 2 . It should be noted that for the horizontal fracture, P^ .ip indicates a reduction with length whereas for the vertical fracture i t increases with length. 104 Table 5*3 The tip fluid pressure data with length horizontal fracture vertical fracture 1 ength* (.ft) P. dbf/ft2) tip length*(ft) P. (Ibf/ft1) tip 1-90 4 700 0-90 41 00 3-40 4500 1-80 4100 4- 90 4 300 3-50 4200 6-40 4 100 5-00 4220 * measured from wellbore axis 105 The fi n i t e element program used for the analyses computes the stresses at the centres of the elements. However, the t i p fluid pressure corresponds to the fluid pressure at the boundary of an element and is compared with the minor principal total stress at the centre of the element adjoining the boundary, to predict the stability of the fracture(i.e. propagation). When the element gets bigger in size, the minor principal stress at the centre of the element could be significantly higher than the same near the edge, where the t i p of the fracture is located. The radial domain in which the vertical fracture analysis is carried out, consists of elements that grow bigger in size with the radial distance. Hence, the increase in the tip fl u i d pressure could be a result of the change in element size with length of propagation. For the horizontal fracture analysis, the domain consists of rectangular elements that are approximately equal in size. Therefore, even though the stresses computed are higher than the stresses at the boundary(closer to the tip) of the element, the effect of the size of the elements will not be as significant as for the radial domain. When a free surface is introduced near a fracture, the resistance to deformation provided by a l l the material above that surface is removed. As a result, a significant portion of the strains are relieved as displacements perpendicular to the free surface. Also, a significant portion of the load caused by the excess fl u i d pressure, acting on the fracture face, is transferred to the material in the tip region. The load mentioned above, increases with length of propagation. Hence, the minor principal total stress at the tip of the fracture 106 decreases with length, .leading to lower tip f l u i d pressures( since they are identical for an equilibrium fracture). From equations 3.10 & 3.17, i t can be seen that the velocity of propagation and the rate of flow are directly proportional to Pb -P^p . As shown in Table 5.3, Pb-Ptip varies between 0.23Pb and O . 4 O P 5 , the influence of which, on velocity of propagation and rate of flow, is significant. In the following sections, the influence of the aperture-length variations on velocity of propagation and rate of flow are discussed, on the assumption that the tip fl u i d pressure is constant with length. Inorder to investigate the effect of aperture(B)-length(L) variations on velocity and rate of flow, i t was f i r s t decided to assume identical B-L curves for both types of fractures. Results obtained with a linear B-L relationship are shown in fig 5.13. The variation of velocity(L) with L, for both the fractures, are seen to.be identical. Hence, i t is seen that the velocity-length curves are useful in studying the influence of various aperture-length relations. The rate of flow, however, should be different because of the axisymmetric nature of the horizontal fracture and the constant height assumed for the vertical fracture. Next, the effect of two different shapes of B-L curves on L-L are examined. The B-L curves assumed are shown in fig 5.14 and the results that correspond to them are shown in fig 5.15. From f i g 5.15, the following could be deduced; (a) i f aperture is linear in length, L is also linear in L (b) i f the magnitude of the aperture is constant with L, L decreases with L (c) for a given length of fracture, for higher apertures the 107 B = (.533* 1.333L )«10' P = 4500 lbf/ft2 p,= 6000 distance from wellbore axis (ft) Fig(5-13) Results of the flow analysis with identical B-L curves for both fractures 108 BrL/IOOCK / *B=-004 / ^^C% - L / 2000 0 4 , 8 12 length (L) -Fi. Fig(5-H) The assumed B-L variations 4 length(L) -f-t. 8 Fig (5-15) The velocity-length predictions for B-L var iat ions in fig (5*1 109 L c o m p u t e d a r e a l s o h i g h e r The v a r i a t i o n s o f t h e a p e r t u r e w i t h l e n g t h o b t a i n e d f r o m t h e f i n i t e e l e m e n t a n a l y s e s , f o r t h e d o m a i n s shown i n f i g s 3.1 & 3.2, a r e shown i n f i g 5.16. From ( a ) t o ( c ) shown a b o v e , t h e s t e a d y i n c r e a s e i n t h e v e l o c i t y o f p r o p a g a t i o n o f t h e h o r i z o n t a l f r a c t u r e , shown i n f i g 5.7, c a n be a t t r i b u t e d t o t h e n e a r l i n e a r v a r i a t i o n o f i t s a p e r t u r e w i t h l e n g t h . The g r a d u a l f l a t t e n i n g o f t h e B-L c u r v e o f t h e v e r t i c a l f r a c t u r e h a s l e a d t o v e l o c i t i e s t h a t i n c r e a s e a t a d e c r e a s i n g r a t e ( i . e . d L / d L ) . A g r a d u a l i n c r e a s e i n t h e t i p f l u i d p r e s s u r e r e d u c e s t h e v e l o c i t y o f p r o p a g a t i o n w h e r e a s a g r a d u a l d e c r e a s e i n t i p f l u i d p r e s s u r e i n c r e a s e s i t ( a s i n t h e v e r t i c a l a n d h o r i z o n t a l f r a c t u r e s , r e s p e c t i v e l y ) . The v e l o c i t y o f p r o p a g a t i o n w i t h l e n g t h , o b s e r v e d f o r t h e two f r a c t u r e s , w o u l d be s i g n i f i c a n t l y d i f f e r e n t i f , (1) t h e a c t u a l p r e s s u r e d i s t r i b u t i o n i n s i d e t h e f r a c t u r e i s c o n s i d e r a b l y d i f f e r e n t f r o m t h e d i s t r i b u t i o n u s e d i n t h e a n a l y s i s t o o b t a i n t h e B-L c u r v e s , o r , ( 2 ) an u p p e r b o u n d f o r r a t e o f p u m p i n g i s i m p o s e d . From t h e c l o s e d f o r m s o l u t i o n ( 4 2 ) d e v e l o p e d f o r a penny s h a p e d h o r i z o n t a l f r a c t u r e , embedded i n a l i n e a r e l a s t i c medium an d s u b j e c t e d t o a u n i f o r m p r e s s u r e d i s t r i b u t i o n , B oC L ( s e e s e c t i o n 2.4) H e n c e , i t c a n be s e e n t h a t f o r a n e a r l y u n i f o r m f l u i d p r e s s u r e d i s t r i b u t i o n , a r e d u c t i o n i n t h e v e l o c i t y o f p r o p a g a t i o n o f a h o r i z o n t a l f r a c t u r e c a n be a c h i e v e d by l i m i t i n g t h e r a t e o f pumping o f f l u i d . I f l e a k a g e o f f l u i d i s a l l o w e d , 110 F i g (5 -16 ) The var iat ion of aperture with lengt h 111 which is the case in reality,' the fluid pressure inside the fracture could be significantly lower than those predicted from the flow analysis witg no leakage. With increasing length of propagation, the pressure losses are likely to be higher since the area of material exposed is high. Hence, at larger fracture lengths, a considerable flattening of the B-L curves could be expected, which would eventually lead to lower velocities of propagation. The rate of flow inside a fracture varies with B 3. Therefore, the change in magnitude of B has a significant influence on the rate of flow. The horizontal fracture is axisymmetric in nature and the vertical fracture has been assumed to propagate with a prescribed height. The aperture(B) obtained for the horizontal fracture is seen to be higher than that for the vertical fracture at larger fracture lengths. As a result, the rate of flow for the horizontal fracture is seen to increase rapidly with length of propagation. If the rate of pumping of f l u i d exceeds the rate of flow of fluid that a fracture can accommodate, a build up of wellbore pressure results. If this increased wellbore pressure exceeds the secondary or the major principal total stress on a plane outcropping the fluid boundaries, a secondary fracture w i l l form in that plane. The resulting increase in the capacity will automatically reduce the wellbore pressure and w i l l be followed by an increase in the rate of flow to satisfy volume compatibility. Therefore, i f the rate of pumping changes significantly during fluid fracture operations, injection 1 12 pressure time variations that are significantly different from those corresponding to primary fracture propagation could result. The injection pressure-time data reported by Holzhausen et al(24) have been obtained under irregular rates of pumping of steam and may have been a result of such secondary fracturing, the rate effects probably playing a dominant role. The inter-dependency of the rate of pumping and the wellbore pressure, as reported by Settari & Raisbeck(37), is illustrated in fig(5.17). The consistent pattern of increase of the injection pressure with the rate of pumping supports the above reasoning. The flow of fl u i d inside a fracture is accompanied by a drop in f l u i d pressure with distance. If the pressure drops are ignored, the finite element analyses predict fractures that are unstable at a l l lengths of propagation. Hence, for proper modelling of fluid induced fracture behaviour this effect must be included. 5.2.4 Shape of the fractures It is of interest to examine the shape and the volume of the fractures that correspond to the different pressure distributions. The shapes obtained with uniform and non-uniform (given by equation 3.6) pressure distributions that correspond to the same maximum pressure, are compared in figs 5.18(a) and 5.18(b). The shape of the vertical fracture obtained for a uniform pressure profile is seen to be close to the theoretical e l l i p t i c a l shape and is shown in fig 5.18(a). The non-uniform 113 114 (a) uniform pressure non-uniform pressure (b) 1 2 3 4 fracture length(ft) uniform pressure non-uniform pressure 2 3 fracture length(ft) Fig(5-18) The shape of the fracture with uniform and non uniform pressure profiles (a) ve r t i ca l fracture (b) horizontal fracture 1 15 pressure distribution has lead to a flattening of the tip region. The cross-sectional area of the fracture for this case is seen to be somewhat lower than for the uniform presure distribution. The shape of the horizontal fracture obtained in the case of uniform pressure, digresses from the e l l i p t i c a l shape. Such a variation is thought to be a result of a major portion of the strains being relieved in the vertical direction due to the presence of a free surface. For the non-uniform pressure distribution, as in the vertical fractures, the tip area indicates a flattening together with a slight change of curvature. The cross sectional area of the fracture is seen to be smaller for the non-uniform distribution. The magnitude of the maximum opening obtained for the horizontal fracture is seen to be higher than that of the vertical fracture for both uniform and non-uniform pressure distributions. Even though the cross-sectional areas of the horizontal and vertical fractures, for a given pressure distribution, are approximately equal(for the case shown in figs 5.18(a) & 5.18(b)),the volume is seen to be higher for the horizontal fracture because of i t s axisymmetric nature and the height assumed for the vertical fracture. 1 16 CHAPTER 6 SUMMARY AND CONCLUSIONS A method of analysis for the initiation and propagation of fluid induced fracture in Oilsand has been developed. The propagation of a fluid induced fracture involves a coupled stress and fluid flow analysis. This is solved herein by f i r s t considering the stress and the flow analyses separately. The results are then coupled through volume compatibility between the quantity of fluid pumped in and the volume of the fracture. For the stress analysis, linear-elastic and elastic-plastic stress strain relations have been assumed for the unfractured and fractured material, respectively. The stress analysis is based on a total stress approach. For the fluid flow analysis, the fractures are approximated as parabolic in shape. The flow of f l u i d between the fracture faces is assumed to be described' by parallel plate flow equations. The flu i d flow analysis has been restricted to steady, incompressible, laminar flow of a Newtonian f l u i d . The analyses performed assume a vertical-horizontal in-situ principal stress orientation. For horizontal fractures, three dimensional axisymmetric analyses have been performed. For vertical fractures, the analyses have been restricted to two dimensions. The magnitude of the fracture initiation pressure is strongly dependent on the intrusion of fluid into the medium. The use of high viscosity fluids or low viscosity f l u i d with high rates of pumping, causing partial fluid penetration, leads 1 17 to a fracture initiation pressure that is significantly higher than for a low viscosity fluid with low rates of pumping. However, the estimation of the instantaneous shut-in pressure is independent of the type of flu i d used. For f l u i d fracturing in Oilsand, with penetrating fluids, the analysis predicts a fracture initiation pressure that is equal in magnitude to the in-situ minor principal total stress that correspond to the location of the fracture. Fracture propagation can be considered as a new initiation subjected to the current stresses of the domain and the current boundary fl u i d pressures. Fluid induced fracture propagation requires the transmission of fluid inside the opening, which is possible only with an open mode fracture. The common assumption that the path of propagation is perpendicular to the in-situ minor principal stress direction has been removed. Instead, a direction perpendicular to the minor principal stress direction at the tip of the fracture is considered. The fracture c r i t e r i a predict the path of propagation based on the computed stresses at the end of each increment of fracture length. The presence of a flu i d pressurized fracture changes the state of stress of the domain in the vic i n i t y of the fracture and thus the orientation of the principal stresses. Hence, the true behaviour of a fluid induced fracture can only be simulated by analysing the propagation of the fracture incrementally, taking account of the stress changes caused in its course of propagation. For fractures initiating in a medium where the i n i t i a l in-1 18 situ minor principal stress direction is vertical, a path of propagation inclined towards the ground surface is predicted. This occurs once the fracture length propagates to a distance about the depth of the fracture from the ground surface. These results are in agreement with f i e l d observations of actual fracture operations and theoretical predictions. Propagation of a vertical fracture, under high rates of pumping of fl u i d , is limited in extent. For a selected wellbore pressure and fluid viscosity, a horizontal fracture exhibits a rapidly increasing fluid carrying capacity with length, than a vertical fracture. The velocity of propagation of a horizontal fracture increases with length with an increasing gradient, whereas for a vertical fracture i t increases with a decreasing gradient. Loss of fluid pressure resulting from leakage of fluid to the surrounding medium and the limiting values of rates of pumping experienced in practice, would reduce the velocity of propagation of the fractures. In fluid induced fracturing with rates of pumping higher than those corresponding to the maximum possible rate of flow in a vertical fracture, fractures must be dominently horizontal in orientation. The effects of rate of pumping of fl u i d on the propagation behaviour of fractures, which has been overlooked by a number of researchers in the literature, is seen to be important. The injection pressure time records obtained under high rates of pumping are considerably different from those obtained for primary fracture operations under low rates of pumping. The predictions of ground displacements and fracture shape 119 from the fin i t e element analyses are in good agreement with the theoretical closed form solution results. 120 CHAPTER 7 RECOMMENDATIONS FOR FURTHER RESEARCH As described in Chapter-3 and discussed in Chapter-5, the analysis approximates the shape of the fracture with a parabola. It was also seen that such an approximation simplified the fluid flow formulations considerably. The shape of the fracture is important in determining the fluid pressure distribution inside the fracture. The correct pressures can be obtained only if the correct shapes are used. Such an analysis requires the fracture to be represented by its real shape. The above can be achieved by incorporating a numerical integration scheme to perform the volume calculations. The analysis can then be extended for a number of iterations to obtain a consistent shape-pressure distribution. The pressure profile described by Daneshy(l2) can be used as a f i r s t approximation in such an analysis. The flow model can further be developed by including multiphase, non-Newtonian fluids. The boundary traction forces exerted by the flow of fluid and the effects of leakage on the resulting pressure distribution could also be examined. Further, the temperature aspects of the problem can also be included. The behaviour of a porous material is governed by the effective stresses. By incorporating a pore pressure model to represent the pore fl u i d behaviour of Oilsand, the analysis can be extended to an effective stress approach. The complexity of the pore fluid of Oilsand requires the consideration of ideal gas laws. The actual f l u i d induced fracture problem is a truly 3D 121 problem. The model developed in this thesis cannot incorporate vertical and horizontal fractures simultaneously. For this reason, the observations such as fracture roll-over from horizontal to vertical or vice versa cannot be examined. Such behaviour can be examined only by a true 3D modelling program. 122 REFERENCES 1 . Advani S.H "Finite element model simulations associated with hydraulic fracturing", Paper SPE 9260, presented at SPE/DOE, Symposium on unconventional gas recovery, Pennsylvania, May 1980. 2. Amodt R.L, Potter R.M "Anomalous fracture extension pressures in granitic rocks", 19th U.S Symposium of Rock Mechanics, Reno, Nevada, 1978, pplO-13. 3. Bawden W.F "Hydraulic fracturing in Alberta Tar sand formations - A unique material for in-situ stress measurements", workshop on hydraulic fracturing and stress measurements, Dec 3-5, 1981, Montery, California. 4. Bjerrum .L, Anderson K.H "In-situ measurement of lateral pressures in clay", Proceedings, 5th European Conference in Soil Mechanics, 1972, Vol-1, pp11-98. 5. Bjerrum L, Nash J.K.T.L, Kennard R.M, Gibson R.E "Hydraulic fracturing in f i e l d permeability testing", Geotechnique 22, No2, 1972, pp 319-322. 6. Byrne P.M, Atukorala U.D, Charlwood R.G "Analysis of fl u i d fracture in Oilsands", Proceedings, 33rd annual technical meeting of the Petroleum Society of CIM, June 6-9, Calgary, 1982. 7. Byrne P.M, Grig R.F "OILSTRESS- A computer program for nonlinear analysis of stress and deformation in Oilsands", Soil Mechanics series No:42, Deptartment of C i v i l Engineering, UBC, Vancouver, Canada, July, 1980. 8. Byrne P.M, Duncan J.M "NLSSIP - A computer program for non-linear analysis of so i l structure interaction problems", Soil Mechanics series No:41, Department of C i v i l Engineering, UBC, Vancouver, Canada, July, 1 979. 9. Cook R.D "Concepts and Applications of Finite Element Analysis" 2nd Edition, John Wiley & Sons, 1981. 10. Cruickshank D.J, Curran J.H "A BIE approach to modelling fl u i d flow in discontinuous rock masses", Proceedings, 7th Canadian Congress of Applied Mechanics, Scherbrooke, May 27- June 1, 1979, pp 887-888. 11. Daneshy A.A "A study of inclined hydraulic fractures", Society of Petroleum Engineers Journal, April 1973, pp 83-97. 12. Daneshy A.A 123 "On the design of vertical hydraulic fractures", Society of Petroleum Engineers Journal, April 1973, pp61-68. 13. Desai C, Abel J "Introduction to the Finite Elment Method", VNR Company, 1972. 14. Desai C, Christian J "Numerical Methods in Geotechnical Engineering", McGraw-Hill Book Company, 1977. 15. Dieterich J.H, Decker R.W "Finite element modelling of surface deformations associated with volcanism", Journal of Geophysical Research, Vol 80, No:29, pp 4094-4102. 16. Dusseault M.B "Stress state and Hydraulic fracturing in Athabasca Oilsands", Journal of Canadian Petroleum Technology, Vol:16, 1977, pp 19-27. 17. " A conceptual geomechanical model for Hydraulic fracture in Oilsand", presented at the AOSTRA, University Industry seminar, Fort McMurray, Alberta, 1979. 18. ' "The behaviour of hydraulically induced fractures in Oilsands", 13th Canadian Rock Mechanics Symposium, University of Toronto, Ontario, May 28-29, 1980. Published by the CIM Montreal, Quebec, 1980, pp 36-41 . 19. Dusseault M.B, Simmons J.V "Injection induced stress and fracture orientation changes", Canadian Geotechnical Journal, Vol:19, November 1982, pp 483-493. 20. Haimson B, Fairhurst C "In-situ stress determination at great depth by means of Hydraulic Fracturing", 11th Symposium on Rock Mechanics, University of California, Berkeley, June 16-19,1969, pp 559-589. 21 . "Hydraulic Fracturing in Porous permeable materials", Journal of Petroleum Technology, July 1969, pp 811-817. 22. Hagoort J, Weatherill ,Settari.A "Modelling the propagation of water flood-induced Hydraulic fractures", Journal Society of Petroleum Engineers, August 1980, pp 293-303. 23. Hanson M.E, Anderson G.D, Schaffer R.J "Theoretical and Experimental research on Hydraulic fracturing", Journal of Energy Resources Technology, Vol:102, June 1980, pp 92-98. 24. Holzhausen G.R, Wood M.D, Raisbeck J.M, Card C C 124 "Results of deformation monitoring during steam stimulation in a single-well test", Applied Oilsands Geoscience, Department of Mineral Engineering, University of Alberta, Edmonton, 1980. 25. Howard G.C, Fast G.R "Hydraulic Fracturing", Monograph vol:2, 1970. 26. Hungr 0, Morgernstern N.R "A numerical approach to predicting stresses and displacements around a 3D pressurized fracture", Department of C i v i l Engineering, University of Alberta, Edmonton, Canada, April 1980. 27. Hussain M.A, Pu S.L, Underwood J "Strain energy release rate for a crack under combined Mode-1 & Mode-2", Fracture Analysis, Proceedings of the 1973, National Symposium on Fracture Mechanics, Part-2, STP 560, ASTM, pp 1-28. 28. Irwin G.R "Analysis of stresses and strains near the end of a crack transversing a plate", Transactions of the American Society of Mechanical Engineers, Vol:79, 1957, pp 361-364. 29. Leach E.R "Hydraulic fracturing of Oilsands ; a literature review", Miscallaneous paper s-77-6, U.S Army Engineers waterway experiment station, March 1977. 30. Massarch K.R * "New aspects of Soil fracturing in Clay", Journal of Soil Mechanics & Foundation Division, ASCE, August, 1978, pp 1109-1123. 31. Medlin W.L "Laboratory investigation of fracture initiation pressure and orientation", Society of Petroleum Engineers Journal, April, 1979, pp 129-144. 32. Ozawa Y, Duncan T.M "ISBILD - A computer program for analysis of stresses and movements in embankments", Report No TE-73-4 to National Science foundation, 1973. 33. Pollard D.D "Forms of Hydraulic fracture as deduced from f i e l d studies of sheet intrusions", Proceedings, 19th Symposium on Rock Mechanics, Reno, Nevada ,1978, pp 7-9. 34. Pollard D.D, Holzhausen G "On the mechanical interaction between a fluid f i l l e d fracture and the earth's surcface", Tectonophysics:53, 1979, pp 27-57. 35. Schapery R.A " A theory of crack initiation and growth in viscoelastic media-Part-I", International Journal of Fracture, Vol:11, No-1, February 1975, pp 141-159. 125 36. Schonfeldt H.V, Fairhurst C "Field experiments on Hydraulic Fracturing", Society of Petroleum Engineers Journal, February, 1972, pp 69-77. 37. Settari A, Raisbeck J.M "Fracture mechanics analysis in in-situ Oilsands recovery", Journal Canadian Petroleum Technology, 1979, Vol:18, No:2, pp 85-94. 38. Shaw P "Coupling Boundary Integral Equation Method to "Other" Numerical Techniques", Recent advances in Boundary Element Method, by Brebbia CA, pp 137-1 47, 1978. 39. Sih G.C "Strain energy density factor applied to mixed mode crack problems", International Journal of Fracture, Vol:10, 1974, pp 305-321 . 40. Skempton A.W "Effective stresses in Soil, Concrete & Rocks", Conference on Pore pressure and suction in soils , London, 1960, pp 4-16. 41. Sobkowicz J "A state of the art bibliography on Hydraulic Fracturing", Vol:1, Alberta environment earth sciences division, May, 1979. 42. Streeter V.L "Fluid Mechanics", 6th Edition McGraw-Hill. 43. Sun R.H "Theoretical size of hydraulically induced horizontal fractures and corresponding surface uplift in an idealized medium", Journal of Geophysical research, Vol:74, No:25, November 15, 1969, pp 5996-6011. 44. Terzaghi K, Peck R.B "Soil mechanics in engineering practice, New York:Wiley. 45. Wiles T.D, Roegiers J.C "Modelling of hydraulic fractures under in-situ conditions using a Displacement Discontinuity Approach", Proceedings, 33rd Annual technical meeting of the Petroleum Society of CIM, June 6-9, 1982, Paper No:82-33-70 46. Zoback M.D, Pollard D.D "Hydraulic fracture propagation and the interpretation of pressure time records for in-situ stress determinations", Proceedings, 19th Symposium on Rock Mechanics, Reno, Nevada, 1978, pp 14-22. 1 26 APPENDIX A  The proposed method of stress analysis STEP-1 Specify the in-situ state of stress in the domain under consideration STEP-2 Begin increasing the boundary fl u i d pressure(Up) until i t exceeds the minor principal total stress(ow) of any one of the elements adjacent to the fluid boundary, (i.e. elements H,I,J shown in fig A.1(a)) When Up f i r s t exceeds omA of an element on a plane the outcropping the fl u i d boundary, a fracture is initiated in the element that corresponds to the above amn.( say i t occurs in element I) STEP-3 Determine the current direction of the plane of crmua and reorient element I in that direction. The mesh after reorientation is shown in f i g A.1(b). If nodes 3 or 2 of element I (ref f i g A.1(b)) exceeds nodes 3 or 2 of elements H or J respectively, in the process of reorientation, a prescribed tolerence of .1(z 3 -z 2 ) or .1(z 3 -z 2 ) is assigned to AB or A,B2 to avoid overlapping. A .Ka) H I J K A.1(b) A H 3 2 B — I -—~~3 2 -—~~3 Al 1 J 2 K F ig (A-1) M e s h pr ior to and after reor ientat ion 128 2 denotes the vertical coordinate of node i in element j . NOTE; Together with the reorientation process are the changes in elevation of the relevant elements. This variation, however, is restricted at most to the length AB,. Hence, by selecting AB, sufficiently small, the errors can be minimized. STEP-4 Reduce the moduli G' & B'(refer to Chapter 4 for notation) of element I to .00001(G'&B') and insert a fluid pressure of magnitude [U p-a m i^]. Modify the nodal load vector to eliminate the forces resulting from the boundaries 1-4 & 2-3 in element I. STEP-5 Analyse the domain for the displacements and compute the stresses and strains of the elements. STEP-6 Compute the Fracture Potentials of the elements H+1,I+1 & J+1. To compute the Fracture Potential, consider the current boundary fl u i d pressure i f a non-uniform pressure distribution is considered. For a uniform pressure distribution this will be the wellbore pressure. The next fracture susceptible element is the element with the highest Fracture Potential. (say i t occurs in element 1+1). STEP-7 Return to STEP-3. Substitute 1+1 for I and repeat the analysis. Continue until the fracture extends the required length. NOTE: Orientation of the elements already fractured are preserved in subsequent reorientations. 130 APPENDIX-B  B.1 Axisymmetric formulation  B.1.1 The Constitutive matrix [D] The generalized Hooke's law for incremental elasticity can be stated as; 5*11 = = [6a, ,-»>( 6a2  + 6a3 3 ]/E B.1 5*2 2 = = [ 6 a 2 2 ~ " ( 5 a 3 3 + 5a,,]/E B.2 5*33 = = [5a 3 3 -p (5a,,+6a 2 23/E B.3 5e i 3 = = 6a,3/G B.4 5e , 2 = = 6a,2/G B.5 5*3 2 = = 6a 3 2/G B.6 For axisymmetric problems; 5e 1 3 = 0 5*3 2 = 0 B. 1 to B . 6 can be represented 5*11 1 - V - V 0 5*2 2 1 - V 1 - V 0 6*3 3 = E -v - V 1 0 5*12 0 0 0 2(l+v) B.7 B.8 5a, , 6a2 2 6a 3 3 6a, 2 B.9 Let, B'= 3B/{2(1+»»)} and G' = E / { 2 ( 1 + J 0 } Inverting the matrix given in equation B.9 and substituting for E & v from the above expressions; 131 6a, , 6a2 2 5a3 3 6a, 2 B' +G' B'-G* B' -G' 0 or, in matrix notation; {5a} = [D]{5e} B'-G' B' +G' B'-G' 0 B'-G' B'-G' B' +G' 0 0 0 0 G' 5e 1 1 6e2 2 6e 3 3 5e , 2 B. 10 B. 1 1 B.1.2 The strain displacement matrix [B] In compatible isoparametric elements, the coordinates and the incremental displacements at a point are approximated as; and where r z u V [N] [N] {5'} B. 1 2 [N] {6} N, 0 N2 0 N3 0 N« 0 0 N, 0 N2 0 N3 0 N„ B. 13 B. 14 It can be shown easily that equation B.13 satisfies compatibility as well as the completeness c r i t e r i a and hence gurantees convergence (from the Potential Energy Theorem). For an axisymmetric problem, the incremental strains are expressed as; 132 5e 1 1 3u/3r B. 1 5 be2 2 5*3 3 u/r -8e, 2 du/dz+dv/dr 3NX 3r~ 0 3N2 3r~ 0 3N3 3r~ 0 3N4 3r~ 0 u l v l 0 3Z~ 0 3N2 3Z~ 0 3N3 3Z~ 0 3N4 3Z~ °2 v2 ^1 r 0 r 0 N3 r 0 r 0 u3 v3 3NX 3N2 3N2 3N3 3N3 3N4 3NA u4 a z 3r 3Z 3r 3Z 3r 3Z 3r V4 B.16 or , in matrix notat ion; {6e} = [B] {6} B.17 where, {B} i s the s t r a i n displacement matrix The shape functions Nj are functions of the l o c a l coordinates s & t . Therefore, a d i r e c t evaluation of the f i r s t der iva t ives of N,-with respect to r & z i s not p o s s i b l e . From the chain rule i t follows that; 3_ ds 3_ 3t 3jr 3S 3r 3t 3S 3Z 3t 3_ 3r 3_ 3Z B.18 [J] 133 The partial derivatives with respect to r & z are written as; 3_ 3r 3_ 3Z [J] -1 as 3 at B.19 The strain displacement matrix formed, contains s,t terms and hence is not a constant for the element. Therefore, i t is evaluated numerically at the Gaussian points with a 2x2 Gauss quadrature. B . I .3 The stiffness matrix [K] Let {5e} represent a virtual nodal displacement vector. From equation B.17 i t follows that; {Se } = [B] {6%} B.20 The internal work done , W(int), w i l l be; W(int) v T Se}{So) dv J"{SeHB]"[D][B]{6} dv V For axisymmetric problems, B.21 dv = 27rr dr dz B.22 Since [B] is in terms of s & t, i t is easy to carry out the integration interms of the local coordinates. 134 dr dz =lJ|ds dt B.23 = {6^}[K]{5} T B.24 where [K] = H 2ir r [B][D][B] |J| ds dt which is called the - 1 -1 stiffness matrix of the element. For rectangular elements|J| is a constant. Hence, the product of the matrices [B][D][B] .IJ| contains terms such as s 2 , t 2 , s t . Therefore, a Gauss quadrature of 2x2 is sufficient for the numerical evaluation of the integral given in equation B.24 (i.e. 2*2 is good enough for a cubic polynomial). Accordingly, For n=4, w = .1 and hence, equation B.25 can be written as; where, f(s,t) represents the integrand of equation B.24 and f ( s j , t f ) represents the value of the integrand at the Gaussian points and w is a weight attached to the Gauss quadrature. B.1.4 The force displacement relationship T 1 I ia J f f ( S , t ) d S d t = i : W f ( S ; , t ; ) B.25 B.26 If the nodal load vector of the element is defined by {f}, for a virtual nodal displacement increment {5 e }, the external work done is given by; 135 W(ext) = { 5 e H f } From the principle of virtual work; W(ext) = W(int) Hence, from equations B.24, B.27 and B.28; {f} = [K]{5} B.27 B.28 B.29 If a pore fluid pressure of magnitude 5U is present, the stress vector for the element can be written as; {6a} = {5a'} + Su-Substituting B.30 in B.21, B.30 W(int) v 5{6?}{5o} dv I{8gjf 6U dv i i • Lf{5eHB]LD][B] [J | ds dt 27rr{6} II T T + ([ {6e}[B] -I-I T SU 27r lJ | r ds dt B.31a B.31b B .31C = {6e}[K']{5} + {5 }{k'} 6U From equation B.28; {f} = [K*]{6} + 5U{k'} {£'}= {f} - 6U {k'} = [K'j {5} {k'} is evaluated in the same manner as [K'] using the 2 2 Gauss quadrature. B.32 B.33 136 B.1.5 The consistent load vector for a distributed load  (perpendicular to the edge) As shown in fig B.I, for an axisymmetric ring element; p' « = 27rr4p f l p' 3 = 27rr 3p 3 B.34 Assume that p' vary linearly from p'„ to p' 3 between the nodes 4 & 3 of the element. On the edge 3-4, t=1. Hence, for a pressure load normal to the edge, the normal displacements are given by; u = {(1+s)/2} u 3 + +{(1-s)/2} u, B.35 The variation of pressure can be represented as; p = p'»{(1-s)/2} + p'3{(l+s)/2} B.36 Let [u 3 uft] represent the nodal virtual displacements. From equation B.35; T u = [u 3 u„] 11+s/2 1-S/2 The external work done can be written as; T W(ext) = I u p ds - 1 r V T T = {uffi} (1+s)/2 (1-S)/2 2/3 2/6 2/6 2/3 (1-S)/2 (1+s)/2 P'« P'a P'« P'a ds 137 Fig(B.1) The axisymmetric isoparametric element 138 2/3p\ + 1/3 p' 3 l/3p'« + 2/3p'3 B.37 Define, {f} = 2/3p'„ 1/3p'3 f 1 B.38 1/3p'« 2/3p'3 f 2 where, {f} is called the consistent load vector. These results are independent of the orientation of the edge in the global system. The r,z global components of {f} can be obtained from trignometry(9). f a x " fi6y/2 " f a y f,8x/2 f HK f 26y/2 f «Y f 26x/2 B.38 in which f,j denotes the component of force in j direction at node i . B.2 Plane strain formulation  B.2.1 The Constitutive matrix [D] From the generalized Hooke's law for incremental elasticity, 6e,! = [6a,, - v{6a22 + Sa 3 3)]/E B.40 6e 2 2 = [6a 2 2 ~ v(5a,,+6a33]/E B.41 6e 3 3 = [6a 3 3 - v(5a,,+8a22]/E B.42 6e,2 = 6a,2/G B.43 139 6e 2 3 = 5a23/G B.44 6e 3 1 = 6a31/G B.45 strain problems; 5e 2 3 = 0 B.46 6e 3 1 = 0 B.47 6e 3 3 = 0 B.48 In matrix form eqns B.40 to B.48 can be written as; 6e, ! 1 -v 0 5a, , 6e2 2 1 -v 1 0 SO 22 5*1 2 ' ~ E 0 0 2(1+p) 6a, 2 Inverting B.49,and substituting for E & v as B.49 5a, 1 B' +G' B'-G' 0 6e , , So2 2 — B'-G' B' +G' 0 6e2 2 Soy 2 0 0 G' Se , 2 B.50 or, in matrix notation; {So} = [D] {Se} B.51 B.2.2 The strain displacement matrix [B] The displacement f i e l d and the geometry of the element are approximated as, r z N, 0 N2 0 N3 0 Na 0 0 N, 0 N2 0 N3 0 N« {sn B.52 1 40 U = [ U ^ ^ U z N z + U a N a + U a N , ] v = [ v ^ + V j N j + V a N a + V / j N a ] B.53a B.53b The element strain vector, {e} for a plane strain case i s , fie, , du/dX B.54 6*22 = av/sy 6*1 2 au/ay + 3v/ax Substituting for u & v in B.54 from B.53, 0 9N 2 0 3N 3 0 3N 4 0 3X~ 3X~ 3X~ 3X~ 0 3N X 0 3N 2 0 3N 3 0 3N 4 3Y~ 3Y~ 3Y 3Y~ 3N X 3N 2 3N 2 3N 3 3N 3 3N 4 3N 4 3Y~ 3X~ 3Y 3X 3Y~ 3X~ 3Y 3X U / B.55 B.2.3 The stiffness matrix [K] The formulation of the stiffness matrix wi l l be similar to that for the axisymmetric case. Hence, 1 1 T [K] = |/[B] [D][B] |JI ds dt t -1-1 B.56 A similar Gauss quadrature as for the axisymmetric element is employed to evaluate the 8x8 stiffness matrix. B.2.4 The force displacement relationship Equation B.31b wil l take the form, 11 T T W(int) =jj{6eJ[B][D' ][B] U l t ds dt {6} -i-i 11 T T " +j({6e}[B] 6U |J| t ds dt B.57 141 Equation B.32 takes the form, {f} = [K'] {5} + 5U{k'} i 1 B.58 where, [K'] = J"J [B]T[D'] [Bjilt ds dt -i - i i i and U ' l = \\ [B] -I -1 SU |J| t ds dt B.2.5 The consistent load vector for a distributed load  (perpendicular to the edge) Only equation B.34 is to be changed for this case. Accordingly, p \ = p« B.59 P'3 = P3 Therefore, the consistent load vector is written as; {f} -|2/3 p'ft l/3p' 3 l/3p'fl 2/3p'3 and the global components are obtained from eqn. B.38. B.60 1 42 APPENDIX-C C.1 Flow in vertical fractures Consider the quasi-static parabolic fracture shown in figure(3.10), exhibiting a length L at time T and L+6L at time T+ST. Due to reasons mentioned in Section 3.6 the modified lengths of the fracture are a 0(L) and a0(L+5L), where ao=0.99 At time T; where, B 0[L] is the maximum half width of the fracture at the wellbore axis. At time T+6T; B[x,L] = B 0[L][1-(x/L) 2] C.1 B[x ,L+6L] = B 0 [ L + 5 L ] [ 1 - { X / ( L + 6 L ) } 2 ] C.2 With reference to f i g 5.14(a), B 0[L] = [ a,+a 2L+a 3L 2+a«L 3] C.3 hence, B0[L+6L] = [a1+a2(L+6L)+a3(L+5L)2+all(L+6L)3] C.4 Let a = a,+a2L+a3L2+aaL3 and therefore, B 0[L] = a C.5 and for f i r s t order terms in 5L, 1 43 B0[L+5L] = a+6.8L C.6 where 6 = [a 2+2a 3L+3a 4L 2] Writing the equation for volume balance across any section, a distance x away from the wellbore axis, SV = |2B0[L+5L][1-{x/(L+5L)}2]dx- J2B 0[L][1 -(x/L) 2]dx C.7 x x. For small 5L, [L+5L]-2 = [1 - 2.5L/L]/L2 C.8 Substituting equations C.5, C.6 & C.8 in C.7, 6V = \2[a+6.8L][1-x2{1-2.5L/L}/L2]dx- J2[a][1 -(x/L) 2]dx * x and defining X = x/L, = 26L[ {ao-a3/3}{a+enL} + {-0.L}X+{0.L-2a}X3/3] C.9 Dividing equation C.5 by 5T and taking the limit as 6T-«-0, 5V a V srZo8T dT 2 dL [ {a0-a3/3}{a+69.L} + {-r3.L}X+{fJ.L-2a}X3/3] C.10 dT Let, K, = [a0-a3/3][a+B.L] C.1 1 144 K2 = [-0.L] C.12 K3 = [0.L-2a]/3 C.13 Substituting C.11 to C.13 in C.10, 9V 2 dL [ K, + K2.X + K3.X3 ] C.14 2>T dT From parallel plate flow equations for steady laminar flow of a Newtonian fluid, av -[ B3 ] r>p C.15 dT [ 12M ] Sx Substituting for B from C.2 in C.15, Sp -12M 3V C.16 Sx 8[B(L+5L)]3[1-{x/(L+5L}2]3 9T Combining equations C.6,C.8,C.14 & C.16; SP dL -3M [ K, + K2X + K 3X 3 ] C.17 dx dT U+0.6L]3 [ 1-(1-26L/L)X2 ] 3 where X = x/L. dx = L dX C. 1 8 Substituting C . 1 8 in C . 1 7 , X P = P B + [ dL [-3ML] [ K , + K 2 X + K 3 X 3 ] dX C . 1 9 x 0dT [a+0.6L]3 [1-(1-26L/L)X 2] 3 as 5L-— 0, 145 x p = p b + j t " 3 * x L ^ d L f R i + K2 X + K3 X3 ] dX C.20 [a 3] dT [1-X 2] 3 P = P b + [e][K 1[ l/{4( l-X 2) 2}+3X/{8( l-X 2)} + (3/8){log(l+X 0)/(l-2g}] + [K 2+K 3 ni/{4(1-X 2 0) 2}3 +K3[-1/{2(1-X20)}] -[e][K,[1/{4(1-X2)2}+3X/{8(1-X2)}+(3/8){log(1+X)/(1-X)}] ~[K 2+K 3][1/{4(1-X 2) 2}] -K3[-1/{2(1-X2)}] C.21 where [e] = 3/zL[dL/dT]/a3 and X0 = r 0/L where r 0 is the wellbore radius p, = flu i d pressure at the wellbore face b For the two halves of the fracture, opening simultaneously, the total SV/dT at the wellbore is given by 2[9V/£T]. The fracture tip velocity is given by [d(a 0L)/dT], where, the tip of the fracture is defined at x=a0L or X=a 0(<l). 146 C.2 Flow in horizontal fractures The shape of the fractures are approximated by similar variations as expressed by equations C.1 and C.2. Since the variation of maximum half width B 0[L] with length is different for this case, the coefficients a,,a 2fa 3,a 4 are different. B 0[L] = [a ,+a2L+a3L2+aItL3 ] C.22 B0[L+6L] = [a,+a2(L+6L)+a3(L+5L)2+a«(L+6L)3] C.23 where, as before, B 0 [ L ] represents the maximum half width of the fracture at the wellbore axis. In writing the equation of volume balance across a section a distance x from the wellbore axis, for the axisymmetric fracture, 6V = 4TT/B 0[L+6L][1-{X/(L+5L)} 2] X dx - 4 7TJB 0[L] [ 1-(x/L) 2 ] x dx C.24 x let a = [a 1+a 2L+a 3L 2+a I |L 3 ] hence, B 0[L] = a C.25 as in C.5, B 0[L+6L] = a+0.SL C.26 where, 6 = [a 2+2a 3L+3a„L 3] Substituting C.25 and C.26 in C.24, SV = 4TT I [a+0.6L][x-x 3/(L+SL) 2] dx X <*oL - 47T I a[x-x 3/L 2 ] dx 1 47 defining X = x/L, 5V = TT 5L[a 2L(2-a 2) (2a+0.L)+(-2.L2.0)X2 Dividing equation C.27 by 5T and taking the limit as Let, K, = [a 2L(2-a 2)(2a+0.L)] K2 = [-2.L2.0] K3 = [L20-2.L.a] Substituting C.28, C.29 & C.30 in C.27, dV ir dL [K, + K2.X2 + K3. X* ] dT dT From parallel plate flow equations for radial flow, -[irB3] x Sp 3T [6M3 9x Combining C.6,C.8,C,31 & C.32; 3 p -6M [ K, + K 2X 2 + K3X" ] dX 8[a+0.5L]3 X[1-X 2(1-25L/L)] 3 where, X = x/L. C.28 C.29 C.30 C.31 C.32 dL C.33 dT +L(L.0-2.a)Xa] C.27 As Sir* 0, 1 48 p =. pb+ 3M 4[o] X K, + K 2X 2 + K3Xa x 0 X[1-X 2] 3 dL dT dX C.34 P = P b + [ e][K 1 { l/{4 ( l - X 2 ) 2 } + l/{2 ( l - X 2 ) } + ( l/2)log { X 2 / ( l - X 2 ) } + K 2[1/{4(1-X 2 ) 2 ] + K 3[-1/{2(1-X 2)}+1/{4(1-X 2) 2] - [ e][K 1 { l/{4 ( l - X 2 ) 2 } + l/{2 ( l - X 2 ) } + ( l/2)log { X 2 / ( l - X 2 ) } -K 2[1/{4(1-X 2 ) 2 ] -K 3[-1/{2(1-X 2)}+1/{4(1-X 2) 2] C.35 where, [ e ] = 3MIdL/dT]/4[a]3 xo = r 0/L where r 0 is wellbore radius p^ = wellbore pressure The fracture tip is defined at X = a0 = .99. For the axisymmetric fracture [ 2>V/ ^ T] in equation C.35 gives the total rate of flow for the f u l l fracture. 149 APPENDIX-D PROGRAM ORGANIZATION D.1 An outline of the modifications The development of the present program involved the following changes to the 1982 version of the Non-Linear Soil Structure Interaction Program; (1) The inclusion of axisymmetric formulations (2) A mesh reorientation process (designed specially for the fracture problem) Accordingly, the following subroutines are affected; DERIVE FORMST DILAT FVECT POREF I SQUAD ISRSLT LAYOUT ELAW In the process, the subroutine MODNOD was added to the program. In order to simulate the forces on the fracture faces, due to the presence of the excess fluid pressure, the pore pressure nodal load vector obtained for a quadrilateral element has been modified in the subroutine POREF. The Flow chart for the present 150 program is shown in fig D.1. D.2 User manual  D.2.1 Added variables KS = -1 minor principal stress plane inclined downwards KS = +1 minor principal stress plane inclined upwards KTA = 1 fracture in r-0 plane w i t h °m^°z > ° ^ = a e KTA = 2 fracture in r-0 plane w i t h <W = a z ' °m^ar KTA = 3 fracture in r-z plane W i t h ° w = ° 9 ' am^°Z KTA = 4 fracture in r-z plane KTA = 5 fracture in z-0 plane with am. = aY , om^=a& KTA = 6 fracture in z-0 plane w i t h °m«=0r ' <V« = a z 151 T A — _ / INPTNP FOR MST L N = 1 yes yes no CAL BAN BEAM Al CALBAN LAYOUT ELAW yes DERIVE yes F 0 R M S T MX = 1 DERIVE L N= L N < contd. 152 NELW: 1 yes DERIVE FO RM S T MODNOD CALBLK MX^ MX .1 I S QU ADDSTF SYM B A N TAPER L J M O D D U J Al SRSLT A1 ELAW c ontd. 2 J y e s 4 c ontd. 1 54 Fig C D-1 ) The Flow-chart for the PROGRAM 155 MPASS = 0 mesh reorientation is bypassed MPASS = 1 mesh reorientation is carried out MTA = +1 reorientation of an element leading to overlapping of upper node (i.e. node 3 as shown in fig D.2) MTA =0 no overlapping MTA = -1 reorientation of the element leading to overlapping of lower node (i.e. node 2 as shown in f ig D.2) MX counter on number of elements to be fractured NOS counter on number of elements fractured D.2.2 Control cards  3rd control card 1 - 5 IELEM IELEM=4, four node isoparametric axisymmetric element. 21-25 MN1 number of elements in each layer (for mesh reorientation this has to be costant) 26-30 NELW number of elements to be fractured. i I * 1 o lower node overlapping with node 2 • upper node overlapping with node 3 Fig (D-2 ) Possible overlapping of nodes due to reorientat ion Force cards Refer NLSSIP manual for this. The pressure cards and the load cards are input in an identical manner). Fracture cards Card-1 1 -10 NSTRT element to be fractured Card-2 1 -5 6 -10 NUMFC 0 nodal point force cards NUMPC 0 nodal pressure force cards Card-3 1 -5 NPW =1 (only a single element should be read for the fracture problem) Card-4 1 -5 6 -10 NMEL above element number PW(NMEL) pore pressure in NMEL Card-5 Repeat cards 3 & 4 for as much elements as required (presently dimensioned for 100 elements) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062766/manifest

Comment

Related Items