A S U R V E Y OF B O U N D A R Y E L E M E N T M E T H O D S FOR A TWO-DIMENSIONAL F R A C T U R E M O D E L By Douglas L. James B. Sc. (Applied Math) University of Western Ontario, 1995 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES INSTITUTE OF APPLIED MATHEMATICS DEPARTMENT OF MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 1997 © Douglas L. James, 1997 In p resen t ing this thesis in partial fu l f i lment of the requ i rements for an advanced deg ree at the Univers i ty of Brit ish C o l u m b i a , I agree that the Library shall make it freely avai lable fo r re ference and study. I further agree that pe rm i s s i on fo r ex tens ive c o py i n g of this thesis fo r scholar ly pu rpose s may be granted by the head of m y depa r tmen t or by his o r her representat ives. It is u n d e r s t o o d that c o p y i n g o r pub l i ca t i on of this thesis for f inancial ga in shall not be a l l owed w i thou t my wr i t ten pe rm i ss i on . Depa r tmen t of The Univers i ty of Brit ish C o l u m b i a Vancouve r , C anada Date ' S f p T . 11, ffl7 DE-6 (2/88) Abstract This thesis surveys boundary element methods for solving a two-dimensional pressurized line crack problem in a homogeneous infinite elastic medium. Starting with an overview of linear elastic fracture mechanics, the exact solution to the pressurized line crack problem is stated and results for several pressure distributions are provided. Various numerical methods for solving the crack problem are then introduced, such as the displacement discontinuity method ( D D M ) using point collocation and a Galerkin method. Crack tip elements and higher-order D D M are discussed. Self-effect correction methods, stated and developed for the piecewise constant D D M on a uniform grid, are shown to significantly improve the numerical solution. A n accurate D D M correction for modelling crack tip element extensions is also presented. Finally, numerical results for the various methods are given and stress intensity factors are also presented for comparison. ii Table of Contents Abstract ii List of Tables vi List of Figures vii Acknowledgements viii 1 Introduction 1 2 Linear Elastic Fracture Mechanics (LEFM) 4 Crack Tip Loading and Modes of Displacement 5 Asymptotic Crack Tip Solutions and Stress Intensity Factors 5 Pressurized Line Crack 7 Displacement and Stress Solution 8 Stress Singularities and the Stress Intensity Factor 10 Constant Pressure Distribution Solution 10 Polynomial Pressure Distribution 11 Pressure Distributions for Smoothly Closing Cracks 12 Strain Energy of the Line Crack 13 3 The Displacement Discontinuity Method (DDM) 14 The Displacement Discontinuity Element (DDE) 15 D D M Solution of the Pressurized Line Crack Problem 17 Calculation of Stress Intensity Factors 19 Crack Tip Elements 20 iii CONTENTS iv "Squareroot Displacement" Crack Tip Element 21 Piecewise Linear Collocation and Higher-order DDMs 22 Self-effect Tip Corrections for Piecewise Constant D D E 24 The Quarter-grid Correction Scheme 24 Numerically Determined Crack Tip Self-effect Correction Schemes . 26 Far-field Collocation Correction Strategy 26 Tip Collocation Correction Strategy 27 A Scaled Quarter-grid Correction 27 A Fractionally Mined Crack Tip Element 28 4 A Galerkin Finite Element Method 32 Piecewise Linear (PWL) Basis Functions 33 PWL Basis Functions with a Crack-tip Correction 34 Matrix Equations 34 Strain Energy Calculation Method 35 5 Comparison of Numerical Methods 37 Relative Error Notation 37 Crack Width Results 38 Quarter-grid Methods and Fractional Mining 39 Fractional Mining Tip Correction 42 Stress Intensity Factor Estimation 42 Asymptotic Method 43 Strain Energy Method 44 6 Conclusions 52 Bibliography 55 CONTENTS v Appendices 59 A . Stresses and Displacements for a DD of Finite Length 59 B. Influence Matrix Simplification for the Symmetric Problem 60 C. P W L C Matrices 61 D. Galerkin Matrices using PWL Basis Functions 63 Stress Matrix . 63 Weighting Matrix 65 E. Galerkin Matrices using Tip-corrected P W L Basis Functions 66 Stress Matrix 66 Weighting Matrix 69 List of Tables 3.1 Alpha values for the far-field correction 27 3.2 Alpha values for the tip collocation correction 28 5.1 Naming conventions adopted for crack width methods 38 5.2 Mid-crack relative error for width calculation methods 39 5.3 Crack tip relative error for width calculation methods 40 5.4 Average relative error for width calculation methods 41 5.5 Stress intensity factors calculated using the asymptotic method. . . . 43 5.6 Stress intensity factors calculated using the strain energy method. . 45 vi List of Figures 2.1 Three principle modes of fracture with accompanying forces F. . . . 5 2.2 Pressurized line crack 8 3.1 The fractional mining tip 30 5.1 Plots of constant DD and related methods (N=10) 46 5.2 Plots of constant DD and related methods (N=20) 47 5.3 Plots of Galerkin and PWCL methods (N=10). 48 5.4 Plots of Galerkin and PWCL methods (N=20) 49 5.5 Quarter-grid and scaled variants (N=10) 50 5.6 Fractional mining tip results (N=10) 51 v i i Acknowledgements I would first like to thank Dr. Anthony Peirce for being such an enthusiastic and insightful supervisor. I have enjoyed my many discussions with him. This work could not have been completed without his encouragement. Next I would like to thank my friends and family for their love and support during the last year. Finally, I would like to thank NSERC for its financial support of this and other work. viii Chapter 1 Introduction Fracture mechanics has undergone considerable development in the last 20 years. In partiuclar the field has seen the development of sophisticated nonlinear and dynamic fracture theories for various media. Nevertheless, the classical theory of linear elastic fracture mechanics ( L E F M ) continues to provide a useful approximation for the modelling of fractures in such materials as brittle rock. The equations of linear elasticity that arise in the context of L E F M may be numerically approximated using boundary element methods ( B E M ) . The discretized B E M equations result in a dense system of linear equations that has dimensions proportional to the number of boundary elements used. There-fore considerable computational effort is required to solve systems involving many boundary elements. The use of simple boundary elements with error correction techniques and fast-solution techniques or higher-order boundary elements are com-mon ways of achieving more accurate and efficient solutions. The following thesis will provide an overview of a few B E M solution techniques by way of their application to the analytically solvable problem of a pressurized line crack in an infinite elastic medium. The simple fracture problem will serve as a model problem to illustrate and compare the numerical methods and also to provide a test case for which numerical correction strategies may be developed. It is hoped that the results may provide guidance for numerical methods used in more 1 2 complicated problems, e.g., problems with more complicated crack geometries or inhomogeneous media. In Chapter 2 L E F M will be introduced and its relation to the pressurized line crack problem discussed. The exact solution to the pressurized line crack problem along with asymptotic relations near the crack tip and their relationship to stress intensity factors are given. Crack width solutions are provided for various commonly used pressure distributions. Expressions for strain energy and strain energy release rates are stated. In Chapter 3 the B E M referred to as the displacement discontinuity method (DDM) is described and its application to the line crack problem is discussed. Constant and higher-order displacement discontinuity elements (DDE) are intro-duced along with crack tip elements. Strain energy calculations and approaches for calculating the stress intensity factors are discussed. Tip correction strategies are introduced for constant D D E on a constant grid. In particular, tip self-effect corrections which allow high accuracy at the crack tip or the middle of the crack are numerically obtained. Also, a nearest neighbour crack tip influence coefficient correction for describing cracks that are longer than an integral multiple of the uniform grid spacing is constructed. In Chapter 4 a Galerkin B E M formulation of the line crack problem is intro-duced. The Galerkin method will be provided for comparison with the D D M since one method is often preferred over the other [28]. A crack tip correction for the Galerkin method as well as strain energy calculations are also provided. In Chapter 5 numerical results and comparisons for the methods contained in Chapters 3 and 4 are presented. The D D M and Galerkin results are compared and the accuracy of the various constant D D E correction strategies are presented. Stress intensity factors are calculated using many of the methods. 3 This thesis also contains several Appendices. The first of these, Appendix A, provides expressions for the stresses and displacements that accompany a general higher-order functional variation D D E . Appendix B contains a D D M matrix sim-plification for the symmetric pressurized crack problem. Appendix C details the construction of the matrix equations for the piecewise linear collocation D D M . The matrices constructed for the Galerkin method in Chapter 4 are calculated in Appendices D and E. Chapter 2 Linear Elastic Fracture Mechanics (LEFM) L E F M has proven to be a successful theory for many fracture problems in elastic solids and brittle materials such as rock. L E F M analysis is an idealization of the crack problem in a number of ways: the crack is assumed to lie in an infinite plane; the fracturing material is assumed to form a continuum with no distinct local material microstructure; nonlinear or plastic elasticity effects occur on such a small scale so as to be negligible. The last point is a significant idealization since, in the neighbourhood of a crack tip, linear elastic theory predicts the existence of a singular stress field. The singular stress field is accompanied by the presence of a plastic process zone of nonlinear behaviour ahead of the crack tip. While a stress singularity implies that the material will be forced to yield in a nonlinear fashion, the linear solution is useful in rock mechanics as long as the criterion for small scale yield (SSY) is satisfied (see [37], pl57). SSY essentially requires that the size of the plastic zone is sufficiently small with respect to the characteristic dimensions of the crack, e.g., length, thickness, specimen size. In brittle rock, the size of the process zone is related to microcracking at the crack tip and is best estimated using by a maximum normal tensile stress criterion than more classical yield criterion, such as the Von Mises criterion ([29]). The following sections provide an overview 4 CRACK TIP L'OADING AND MODES OF DISPLACEMENT 5 Figure 2.1: Three principle modes of fracture with accompanying forces F. of L E F M sufficient for the development of the numer ica l methods that follow. Crack Tip Loading and Modes of Displacement Cons ider an infinite plate containing a crack and also having far-field tract ions. T h r e e different ways of loading the plate, designated by modes I, II and III, are i l lustrated in F i g . 2.1. M o d e I is the "opening m o d e " that results f rom tensile stress perpendicu lar to the plane of the crack. M o d e II is the "sl iding m o d e " that results f rom in-plane shear forces. M o d e III is cal led the "tearing m o d e " and is accompanied by out-of-plane shear forces. T h e crack is assumed to be slit-like, of length 2L, and to lie on the line {{x,y)\ — L < x < L,y = 0}. Fur thermore , the crack surfaces are assumed to be traction-free. B y considering un i fo rm far field stresses which correspond to each of modes I, II and III, a plane stra in elasticity p rob lem has been formulated. Asymptotic Crack Tip Solutions and Stress Intensity Factors T h e plane stra in elasticity p rob lem of the previous subsection may be solved using Westergaard functions (see [30], [22]) or Goursat -Muskhe l i shv i l i potent ia l functions [20]. T h e mode I asymptot ic stress solutions for b o t h plane stress a n d plane stra in are ASYMPTOTIC CRACK TIP SOLUTIONS AND STRESS INTENSITY FACTORS 6 ~ ^ f e cosCf) { 1 - ain ( f ) Bm ( ^ ) } ^~^ C 0 S(f) sMI) s i n(¥) as r —> 0, w h e r e if/ i s c a l l e d t h e m o d e I s t r e s s i n t e n s i t y f a c t o r . F o r a r e m o t e l y a p p l i e d u n i f o r m t e n s i l e t r a c t i o n C Q O , t h e s t r e s s i n t e n s i t y f a c t o r i s if/ = c T o o v ^ - (2 .2) T h e a s y m p t o t i c d i s p l a c e m e n t f i e l d s i n t h e n e i g h b o u r h o o d o f t h e c r a c k t i p a r e f o u n d t o b e «« ~ % s c o s (I) iK~1 + 2 s i n 2 (I)} ,„ (2 .o) %~^v / ^ s i n ( l ) { K + 1 - 2 c o s 2 ( l ) } as r->0, w h e r e 2—^ p l a n e s t r e s s 1 + " " (2-4) 3 — 4J/, p l a n e s t r a i n a n d v i s P o i s s o n ' s r a t i o [18]. I n s i d e t h e c r a c k , a l o n g t h e 9 = TT l i n e , t h e a s y m p t o t i c f o r m f o r t h e n o r m a l d i s p l a c e m e n t is ^ { 4 ( 1 - * , ) } = ^ ^ ( p l a n e s t r a i n ) . (2 .5) T h e s o l u t i o n s f o r m o d e s I I a n d I I I h a v e t h e s a m e a s y m p t o t i c r d e p e n d e n c e (see [35] f o r d e t a i l s ) a n d c a n b e d e r i v e d i n a s i m i l a r w a y . I n f a c t , m a n y o f t h e p r o b l e m s o f L E F M c a n b e s o l v e d u s i n g s i m i l a r t e c h n i q u e s s i n c e t h e y a l s o i n v o l v e p l a n e e l a s t i c i t y p r o b l e m s . N . I . M u s k h e l i s h v i l i [20] h a s d e v e l o p e d a g e n e r a l m e t h o d f o r s o l v i n g p l a n e e l a s t i c i t y p r o b l e m s t h r o u g h a c o m p l e x v a r i a b l e m a p p i n g s c h e m e a n d u s e o f t h e K o l o s o v e q u a t i o n s . B r i e f i n t r o d u c t i o n s t o t h e t e c h n i q u e c a n b e f o u n d i n [34] a n d [32]. PRESSURIZED LINE CRACK 7 A s (2.1) indicates, the strength of the stress singularity in the ne ighbourhood of the crack t ip is control led by the stress intensity factor, Kj. A s a result, the stress intensity factor plays an important role in L E F M for predict ing when a sol id wi l l y ie ld and crack elongation wil l occur. A s imple cr i ter ion for decid ing w h e n the crack wi l l propagate is whether or not the stress intensity factor is greater than a cr it ical stress intensity value, KJC, e.g., KL > KIC • (2.6) T h i s concept was first studied by Gr i f f i th ([16], [17]) and is often cal led Grif f i th's stabil ity cr i ter ion. T h i s and other fracture propagat ion cr i ter ia can be interpreted as condit ions on the energy release rate in order for the crack to open, e.g., see (2.29). M o r e detai led discussions of dynamic fracture are given by K a n n i n e n a n d Pope lar [18] or F r e u n d [13]. Pressurized Line Crack T h e two-dimensional line crack p rob lem for a solid in equ i l ibr ium is important for the s tudy of fracture mechanics. In this case, it is possible to obta in closed fo rm solutions wi th relative ease compared to fractures of a more compl icated geometry. W h i l e the growth of cracks is inherent ly a dynamic phenomena, an understanding of the equ i l ibr ium case is a natura l start ing point . A l so , the simple line crack p rob lem in an elastic m e d i u m provides a test bed for developing more sophist icated numerica l methods for fracture mechanics. F i rs t studied by Gr i f f i th [16] in 1920, the pressurized line crack p rob lem consists of a line crack (or hollow slit) of length 2L, s i tuated along the x axis and centred about the or ig in (see F i g . 2.2). T h e p rob lem is assumed to be symmetr i c about the x and y axes and s ituated in an infinite plane. T h e stresses in the m e d i u m are PRESSURIZED LINE CRACK 8 F i g u r e 2.2: Pressur ized l ine crack given by crxx(x,y), ayy(x,y) and crxy(x,y) — ayx(x,y). T h e crack interface located in the y > 0 region is called the upper crack face, and the displacement of the upper crack face is given by uv(x) = uy(x,0). T h e elasticity boundary condit ions for the x > 0 upper crack face are cryy(x, 0) = — p(x), 0 < x < L, axy(x,0) =0, x > 0, uv(x) = 0) x > L, (2.7) where p(x) is an even funct ion representing the net pressure responsible for ho ld ing open the crack. A n addit ional boundary condit ion is that al l stresses are zero a n infinite distance f rom the crack. T h i s p rob lem wi l l be referred to as the symmetr i c crack prob lem. Displacement and Stress Solution T h e symmetr ic crack p rob lem has been solved using the plane elasticity solut ion techniques developed by Muskhel ishv i l i [20], e.g., E n g l a n d and G r e e n [11], G r e e n and Zerna [15] or Sneddon [31]. T h e norma l displacement of the upper crack face i s Vx^=T2 0<x<L, (2.8) PRESSURIZED LINE CRACK and the stress ahead of the crack t ip is given by avy(x,0) = 7T xg{L) y/x2 - L2 x > L, where Jo p{x)dx 0 < x < L. (2.9) (2.10) T h e plane stra in modulus , E', can also be expressed in terms of Young's modulus , E, the shear modulus , G, and Poisson's ratio, v, as E 2G E' l - l / 2 l - l / (2.11) These crack results were der ived for the case in which g(£) is differentiable. N o t e that if p(x) is differentiable then <?(£) is a continuous funct ion. It turns out [31] that (2.8) and (2 .9) are val id even if p(x) has a finite number of j u m p discontinuities. T h e stress in the ^-direction is given by the equations of plane stra in as crxx(x,0) = ayv(x,0), x > L. (2.12) It can also be verif ied that the condi t ion axy(x,0) = 0 is satisfied for x > 0. A n alternate expression of uy(x) is given by Spence and Sharp [33, append ix A] as UV(X) ^E'So^ VL2 - x 2 - ^L2 - £2 (2.13) I VL2 - x 2 + ^/L2 - i2 | Alternat ively , the solut ion to the crack elasticity p rob lem m a y express the pres-sure in terms of the crack width . Spence and Sharp [33, see A p p e n d i x A] have shown that the crack solution may be wr i t ten as duy(s) ds L ds s — x' (2.14) where p.v. denotes the C a u c h y pr inc ip le value. For the symmetr i c prob lem, (2 .14) simplifies to ?' r L duv(s) s K JO ds s2 — x •As (2.15) PRESSURIZED LINE CRACK 10 where the pr inc ip le value is assumed. Integration by parts and the condi t ion uy(L) = 0 yields ™ JO (Sz — Xz) B o t h (2 .15) and (2 .16) wi l l be useful for developing numer ica l approx imat ion schemes to the crack stress boundary value prob lem. Stress Singularities and the Stress Intensity Factor F r o m (2 .9) it is clear that a stress singularity can exist at the crack t ip , i.e., at x = L. A n asymptot ic expansion about x — L, assuming <?(£) is nonsingular, identifies the square-root singularity typ ica l of L E F M : T h e stress intensity factor, Ki, can then be defined as x —> L+. (2.17) Ki= lira V2Trrayy(L + r,0) (2.18) Us ing ( 2 . 17 ) , Ki is found to be F r o m (2 .19) it is clear that g[V) has a direct effect on the strength of the singularity. A l so , note that Kj has no dependence on mater ia l propert ies. Constant Pressure Distribution Solution B y far the most popular pressure d is t r ibut ion in the l i terature is that in which the pressure is constant. In this case, the funct ion g(£) is also constant and given by PRESSURIZED LINE CRACK 11 Us ing (2.8), the displacement solut ion for the upper crack face is /(z>0) 4 J Jx IfE' Jx jx2 _ £2 E and therefore a m a x i m u m crack w id th of }4V 0 < x < L, (2.21) wo (2.22) occurs at the middle of the crack (x = 0). T h e stress ahead of the crack t ip is given by O-yy(X,0) =PQ 1 x > L. (2.23) _Vx2 - L 2 T h e stress intensity factor, determined f rom (2.23) using (2.18) or f rom (2.19), is found to be Ki = poV^L. (2.24) Polynomial Pressure Distribution If the pressure is given by a po lynomia l funct ion, p(x) =p0 +pix+p2x2 -\-p3x3 + ... Valk6 and Economides [36] have shown that rc\ nP0 , c , ^P2^ , 2P^ , 9(0 = ~Y +Pi£ + — 2 ~ + —3— + ,(*,0) = E< " ~ ~ ' nE> +^7 ( L 2 + 2x2)VV^ + L^/L2 -x2 + x2ln L + VL2 - x2 + 3E' P3 IxE' + Lx2^j sjL2 - x2 + x4 In 'L + VL2-X2' (2.25) + + and PRESSURIZED LINE CRACK 12 +P2 + 2P3 ^ - x3^j /y/x2 - L 2 - x 2 + ^-y- ~ Lx^ / V ^ 2 - L 2 - x3 arctan ^~ + .. ^ x 2 - L 2 , T h e mode I stress intensity factor, found by combin ing (2.19) and (2.25), is Pressure Distributions for Smoothly Closing Cracks Unders tand ing the crack t ip region is, by nature, fundamenta l to understand-ing the fracture process. Indeed, in hydraul ic f ractur ing the crack t ip region is the source of m u c h concern [9], [28]. A popular boundary condit ion, often cal led the Barenblat t t ip condit ion based on the work of Barenblat t [2], is that the displace-ment field, uy, satisfy the smooth closing condit ion duy(x = L,0) dx = 0 In order for this to occur, it is sufficient to have g(L) = 0 (2.27) which is the so-called "z ipper crack equat ion" [36]. F r o m (2.19) it is clear that the z ipper crack equation is equivalent to having if/ =0. In terms of the pressure funct ion, p(z), a region of "negative pressure" is required in order to force the sol id to close smoothly. A popular pressure d istr ibut ion used by Khr i s t ianov i t ch and Zhel tov ([19], [38]), Barenblatt [2] and G e e r t s m a and de K l e r k [14] is a two-level piecewise constant funct ion given by r Pi, x<x0 -p2, x0 < X < L p(x) = { STRAIN ENERGY OF THE LINE CRACK 13 where xn is the locat ion of the pressure discontinuity, and p\ a n d pi are posit ive constants. T h e locat ion of the discontinuity is determined by the z ipper crack equation. E q u a t i o n (2.10) implies 9(0 = | and so (2.27) implies [19] P I T T 2 ' -^f- + (pi +£>2)arctan X < XQ , X > £ 0 #o = I<sin ixp2 2(pi+P2)' Strain Energy of the Line Crack W h e n a crack is opened, energy is stored in the elastic m e d i u m in the form of stra in potent ia l energy. U s i n g the definit ion of work, the energy stored in the m e d i u m for x > 0, per unit thickness (in the z d irect ion), for the symmetr i c crack p rob lem is Wo = / p(x)uy(x)dx. Jo T h i s result assumes that the pressure varies l inearly w i th the displacement. Sned-don [31] showed that this may also be wr i t ten in terms of g(£) as W0 = ^ l (2.28) Tak ing the part ia l derivative of (2.28) wi th respect to L yields the stra in energy release rate per crack length d W0 dL ^ 2 | f (using (2.19)), (2.29) T h i s last result, given by R i ce in [26], wil l prove useful for est imating Kj numerical ly. Chapter 3 The Displacement Discontinuity Method (DDM) T h e displacement discontinuity m e t h o d ( D D M ) is a b o u n d a r y element m e t h o d ( B E M ) suitable for the numerica l descr ipt ion of cracks in a linear elastic m e d i u m . Or ig ina l ly developed by C r o u c h ([6], [7]), the D D M m a y be used to approximate the crack boundary by a sequence of b o u n d a r y elements called displacement discon-t inuity elements ( D D E ) , or D D elements. T h e D D E utilizes the analyt ic solut ion for a finite segment of a (constant) displacement discontinuity in an infinite sol id in order to approximate the result ing tractions at the locations of other D D E on the crack. A s a result, for a proper ly formulated elasticity p rob lem (see T i m o -shenko [34]), it is possible to construct a system of equations so that the unknown norma l and tangential stresses and/or displacements m a y be solved for in terms of the known ones [8]. T h e following s tudy wi l l be concerned w i th stress b o u n d -ary value problems, in which the tractions are specified on the boundary and the boundary displacements are to be determined. A n example of such a p rob lem is the pressurized crack for which the pressure is specified and the crack w id th is to be determined. A n excellent reference for the D D M , which contains m u c h of the following mater ia l , is the text by C r o u c h and Starf ield [8]. N u m e r i c a l results for the various methods presented in this chapter wi l l be delayed unt i l C h a p t e r 5. 14 THE DISPLACEMENT DISCONTINUITY ELEMENT (DDE) 15 The Displacement Discontinuity Element (DDE) Consider the p rob lem of a D D E in an infinite elastic sol id. In this case, the displacement field, Uj — (ux, uy), is everywhere continuous except along a finite line segment in the xy-plane. Denote the line segment by C = {(x,y)\ —a<x<a,y = 0}. Since the D D E wi l l be used to represent the displacement between two crack faces, it is useful to refer to the y = 0_ side of C as the negative side and the y = 0+ side of C as the positive side. W h i l e more general funct ional forms are possible (see p. 22), in what follows, it wil l be sufficient to consider the simplest and most popu lar case of a constant D D , i.e., the D D is constant along C. T h e D D E is descr ibed by the coordinates Dx and Dy, where T h e Dx and Dy coordinates describe the D D E ' s norma l and shear displacements, respectively, that occur u p o n crossing C. N o t e that the convention is for Dy to be negative when the crack faces are physical ly displaced f rom one another, i.e., the crack has been opened. T h e result ing displacement and stress fields associated w i th the constant D D E were given by C r o u c h as Dx = ux(x,0-) - ux(x,0+) Dy = uy(x,0-)-uy(x,0+). ux(x,y) = Dx [2(1 - !/)/,„ - yf>xx] + £>„[-(! - 2v)f,x - yf,xy] (3-1) uy(x,y) = Dx [(1 - 2v)f,x - yf,xy] + Dy [2(1 - u)fiV - yf,yy] (3-2) and = 2G{DX [2ftXy + yf,Xyy] + Dy [fiVy + yf,yyy]} (3.3) Vyy(x,y) = 2G {Dx [-yf,XVy] + Dy Hyy - yf,yyy}} (3.4) axy(x,y) = 2G {Dx [f}yy + yf,yVy} + Dy [-yf,xVy]} (3.5) THE DISPLACEMENT DISCONTINUITY ELEMENT (DDE) where -1 16 4-7r(l - v) i y y axctan arctan • y arctan \ x — a x + a —(x — a) In y/(x — a)2 + y2 + (x + a) In y/(x + a)2 + y2 (3.6) W h e n construct ing a l inear D D M approx imat ion of the solut ion to the pres-surized line crack, it wil l be useful to know the stress result ing from a single D D E at the posit ion of another D D E on the discretized crack. W i t h this in m i n d , the displacements and stresses along the x-axis can be the calculated f rom (3.1)-(3.6). For \x\ > a , y = 0 the displacement components are (1 - 2i/) ux(x,0) 47r(l — v\ In x — a U^0) = +11(1-1)* x + a x — a x + a Dv Dx while for |a;| < a they are ux(x,0±) = = F - i ? B _ i - ^ _ ^ l n / n x 1 ^ (1-2^) , uy(x,0±) = T ^ s + ^ x — a x + a x — a x + a Dv Dx (3.7) (3.8) 2 ~ t f ' 4TT (1- I / ) T h e constant discontinuities Dx and Dy are clearly visible from (3.7) and (3.8) for |x| < a. F r o m (3.3)-(3.6) the stresses along y = 0 are found to be -G ' crxx(x,0) O-yy(x,0) = Pxy(x,0) = 7T(1 -V) X2 - a2 G a 7T(1 -v) X2 - a2 G a 7r(l -v) X2 -a2 Dy Dy Dx (3.9) (3.10) (3.11) with axy(x,0) = ayx(x, 0). Not i ce that stresses are singular at the D D E ' s endpoints (x = ±a) but continuous everywhere else. A l so , the stresses are continuous across the y = 0 l ine, by construct ion. DDM SOLUTION OF THE PRESSURIZED LINE CRACK PROBLEM 17 D D M Solution of the Pressurized Line Crack Problem Once again, consider the pressurized line crack p rob lem in an infinite elastic body. T h e boundary condit ions are cryy(x,0) = -p(x), \x\<L, <JXy(x, 0) = 0, — O O < X < O O , uy(x, 0) = 0, |a;| > L, along w i th the addit ional condit ion that al l stresses and displacements are zero at infinity. B y d iv id ing the crack into JV l ine segments, each representing a D D E placed end-to-end wi th adjacent D D E , the w id th of a crack may be approx imated by a piecewise constant ( P W C ) funct ion. Le t the ./V D D E be denoted by (Dx,Dy), k = 1,..,N and be of length hk = 2afc. T h e midpoint of the kth D D E , called the noda l point , is therefore located at the x posi t ion fc xk = —L + 2~2 hi — ai • A l t h o u g h curved cracks w i th shapes more general than a straight line may be dis-cretely approx imated using the D D M [8], the line crack p rob lem wi l l be sufficient to i l lustrate the method . A l so , since the line crack p r o b l e m has been solved ana-lytically, the qual ity of the D D M numer ica l solut ion wi l l be evident. T h e stress relations for a single D D E , (3.9)-(3.11), in conjunct ion w i th the crack boundary condit ions, allow the construct ion of a l inear system of equations which may be solved for the result ing D D ampl itudes. Specifically, the no shear stress B C at the D D E midpoints , CTxy{Xk,0) = 0, k = l,..,N, along wi th (3.11) implies that £>* = 0, k = l,..,N. DDM SOLUTION OF THE PRESSURIZED LINE CRACK PROBLEM 18 T h e norma l stress B C appl ied at the D D E midpoints , o-vv(xk,Q) =-p(xk), k = l,...,N, along wi th the linear superposi t ion of the resultant D D E stresses, N O-yy{xk,0) =Y^<?yy{Xk,Q), k = l,...,N, j = l and (3.10) y ie ld the condit ion N ~V{Xk) = Y^O-yy{Xk,Q), k = l,...,N, 0=1 where cjy(xk,0) = -- ^ 52V, j = l,...,N, y y 7r ( l - v) [Xj - xk)2 -aj v is the normal stress at x — Xk result ing f rom the jth D D E . T h e result ing mat r ix equation to be solved for the unknown D D E norma l displacements, {£>^ }j^ =1, is N -Pk = Y^AkjDi> k = l,...,N, (3.12) where pk — p(xk) and G_ 7r(l - v) (XJ — xk)2 - a2-A-kj — ~i 7Kt„. .„?\2_„2> k,j— 1,...,N, (3.13) are called the D D E influence coefficients. T h e result ing D D E ampl i tudes, Dy, m a y be interpreted as the (negative) w i d t h of the crack at the midpo int of the j t h element located at the noda l points, Xj. T h i s method is also called a P W C col locat ion method since stress is col located at the nodes of the D D E . For the symmetr i c p rob lem where p(x)=p(-x), the result ing displacement solution is also symmetr ic : Dy = D»+1~j, j = l,...,N: CALCULATION OF STRESS INTENSITY FACTORS 19 In this case, it is unnecessary to solve the dense N-by-N system of equations (3.12) when an ^ -by-^ system (for N even) will produce an equivalent result. A formula for constructing the ^-by-^r influence matrix from A^j, for the positive x-axis, is given in Appendix B. Calculation of Stress Intensity Factors Since stress intensity factors play such a central role in characterizing fractures in L E F M , it is necessary to be able to numerically approximate them with sufficient accuracy. A direct approach that is commonly used is to use the asymptotic form of the normal displacement, uy(x,0+), at the crack tip. Using (2.5) the asymptotic relation is uy(x,o+) ~ — y - ^ — . (3-14) The stress intensity factor is then approximated as K<~T\[I~XZU;> ^ where U* « 2uy(x, 0+) is the width of the crack at a distance L — x* <C L from the crack tip. Typically U* is the magnitude of the D D E closest to the tip. In order to be effective, this approach clearly requires high resolution at the crack tip. Higher accuracy crack tip elements, and/or relatively more dense distributions of D D E near the tip are often used to improve the accuracy of the crack tip solution. For example, in addition to a special squareroot displacement crack tip element, Raveendra and Cruse [25] have used constant elements of width h, except for near the tip where four individual D D E of width | were used. An alternate approach for estimating Kj is to use equation (2.29). Since the strain energy, per unit thickness, of a half-crack of length L is (3.16) CRACK TIP ELEMENTS 20 a simple numerical approximation, using constant D D E of width hk, is fc=i Then, using a difference approximation such as dWp(L) _ W0(L + h)-W0(L-h) dL ~ 2h + 0(h2) (3 .17) the stress intensity factor may calculated from (2.29) using Kj = (3 .18) In cases where higher accuracy is required and/or more sophisticated DDMs are used, a more appropriate approximation of (3.16) and its derivative should be used. Examples of the asymptotic and the strain energy release rate methods for calculating Kj are presented in Chapter 5. A simple way to improve the accuracy of the D D M crack solution is to combine the constant D D E with special crack tip elements. These tip D D E incorporate the known asymptotic dependence of the displacement solution in the neighbourhood of the tip (see (2.3)) in order to calculate more accurate influence coefficients for stress BVPs. In general, for a mixed B V P (i.e., a B V P where a mixture of normal and shear component stresses and displacements are specified on the boundary), asymptotic crack tip elements may be used for both the stress and displacement fields. For example, Raveenda and Cruse [25] consider displacement, u(r), and traction, t(r), tip parametrizations of the form Crack Tip Elements CRACK TIP ELEMENTS 21 where {Uj}j=1, {Tj}j=1 are expressed in terms of the nodal element amplitudes, e.g., D D E amplitudes and the traction element amplitudes. However, in the following, only stress BVPs will be considered. Calculation of the normal stress along the a>axis resulting from a general D D E shape, ,(x) = ux(x, 0_) - ux(x, 0+) < 0, (3.19) when x 6 [x i , X2], will use the fact that <Tyy(x,0) = 77—" \ U m / %(£) (3.20) for x ^ Xi, x<2 [8]. It is easy to verify that avv(x,0), as given by (3.10), results from (3.20) when •ly(x) = l Dy, x G (—a,a), 0, otherwise. "Squareroot Displacement" Crack Tip Element Since the asymptotic crack tip width depends on the square-root of the distance from the tip (see (3.14)), consider a tip element of length h = 2a and (-ve) width given by ( A/vf. * e ( 0 , 2 a ) , 0, otherwise. Note that the crack is located in the x > 0 region. Evaluation of (3.20) yields uy(x) — ux(x, 0_) - ux(x, 0+) = (3.21) CTyy(x,0) = -G 2TT(1 - v) x—2a lyjax \/o.{-x) arctan y/x—\/2a 2a y/2a(-x) (—x) 2a—x si x > 0, x < 0. (3.22) PIECEWISE LINEAR COLLOCATION AND HIGHER-ORDER DDMS 22 T h e same argument used to derive (3.13) leads to the use of Akj — 2TT(1 - v) \Xkj - 2aj 2^Xl~ -G \ V2 1 In j = 1 or N where Xkj = aj + \xk — Xj\. Hence, the t ip element amounts to changing the influence coefficients for 27V elements of the influence matr ix , Akj. A l so , the long distance character of the t ip element influence coefficient is different f rom that of the constant t ip element, since while the constant t ip has the inverse square distance dependence (see (3.13)). W h i l e the squareroot t ip element does reduce overall crack error and improves the t ip accuracy, there is a l imit to the increase in accuracy that a single D D E t ip correct ion can provide. A l so , since the squareroot t ip correct ion only involves one element, the relative importance of the crack t ip element becomes negligible as N gets larger. O n e disadvantage of the crack t ip element is that O(N) influence coefficients must be modif ied. T h e use of higher-order D D E wi th more degrees of f reedom t h a n the constant D D E is a natura l choice for achieving a more accurate approx imat ion of the crack width . However, higher-order D D M s are more difficult to formulate and implement , especially in higher dimensions. A l so , the result ing linear system of equations has different propert ies which can make it more difficult to solve approx imate ly in the case when N is very large. Nevertheless, higher-order methods can be rel ied u p o n to y ie ld significantly more accurate results than the ord inary P W C D D M . In the following, the case of l inear elements wi l l be considered. T h e P W L D D M wil l be referred to as the piecewise l inear col location ( P W L C ) method . -G 1 \/27r(l - v) \xk - Xj Piecewise Linear Collocation and Higher-order D D M s PIECEWISE LINEAR COLLOCATION AND HIGHER-ORDER DDMS 23 Cons ider a P W L D D E of length h = 2a a n d ly ing along the a>axis f rom x — —a to x = a. Assume , for simplicity, that the shear discontinuity component is zero and that the norma l discontinuity is given by D^(x), x € (—a,a). W h i l e (2.16) could be used to calculate the stress at points along the x-axis result ing from an arb i t rary D D E , C r a w f o r d and C u r r a n [5] have given expressions for the stresses and displacements associated w i th a general two-dimensional D D E (see A p p e n d i x A . ) . U s i n g the expression for ayy, the stress along the x-axis is cryy(x,0) = —^—r- lim (^- - y^-4) f \»-—^DN(e)de. (3.23) 2TT(1 - v) y->o \dy dy2JJ_a(x-e)2 + y2 B y choosing to collocate stresses along the P W L element at noda l posit ions given by the Gauss-Chebyshev integrat ion points, good convergence propert ies are obta ined [5]. Accord ing ly , the l inear D D E is wr i t ten as DN(x) = N!(x)D! + N2(x)D2 (3.24) where Di, D2 are the noda l displacements and Subst i tut ion of (3.24) into (3.23) yields oyy(x, 0) = 2 ? r ( 1 G _ [F(x, -l)Dx + F{x, +l)D2\ (3.25) where F(x, n) = 9 a 9 + nV2 xz — az T h i s result can then be used to construct an approximate solut ion of the pressurized line crack prob lem. T h e construct ion of the P W L C mat r ix equat ion is detai led in A p p e n d i x C . 2a In x — a x + a (3.26) SELF-EFFECT TIP CORRECTIONS FOR PIECEWISE CONSTANT DDE 24 A n easy way to increase the accuracy of the P W L C m e t h o d is to use a crack t ip element, such as (3.21), at each end of the crack. T h e stress influence of the t ip element is given by (3.22). For purposes of numer ica l compar ison, a t ip element of the same length as the linear element wi l l be used. T h e t ip element significantly reduces the error at the crack t ip , as wi l l be shown in C h a p t e r 5. Self-effect Tip Corrections for Piecewise Constant D D E Based on the success of the.crack t ip element, it is clear that numer ica l accuracy at the fracture t ip plays an important role in the overall accuracy of the solut ion. T h e strong dependence of crack w id th accuracy on crack t ip influence coefficients is reasonable since the influence matr ix is an approx imat ion of a F r e d h o l m integral operator which is singular at the crack t ips. In fact, it wi l l be i l lustrated that large changes in the accuracy of the crack w i d t h solution can be achieved by s imply correct ing the t ip influence coefficients. In the most extreme case, on ly a single number in the influence matr ix , the crack t ip self-effect, is changed. T h e correct ion may be interpreted as introduc ing corrective stresses at the t ip (or edge) elements which reduce the crack w i d t h overest imation characteristic of the P W C D D M [27]. M o d e l l i n g a crack wh ich has a length that is not a mult ip le of the constant gr id of size h is also considered. T h i s is often refered to as a fract ional ly m ined crack. The Quarter-grid Correction Scheme U s i n g the "method of rat ional functions," R y d e r and Nap ie r [27] were able to analyt ical ly study the error associated w i th D D M model l ing of cracklike tabular stress problems on a un i fo rm gr id . E r r o r analysis of the piecewise constant D D M solution of a pressurized line crack, for constant pressure, on a constant gr id of size h indicated a quarter gr id discret izat ion error in the crack width . U s i n g this insight, SELF-EFFECT TIP CORRECTIONS FOR PIECEWISE CONSTANT DDE 25 a correct ion strategy was created which mult ip l ied the self-effect influence coeffi-cients corresponding to the crack t ip elements by a factor. T h i s edge/tip correct ion also supported nondiscrete posit ioning of the min ing edges, i.e., fract ional m in ing . Therefore, by s imply correct ing the edge self-effects, the numer ica l m e t h o d cou ld more accurately coincide wi th the locat ion of the excavation edges. T h i s is conve-nient since perfect al ignment of the excavation w i th the gr id is often inconvenient. In addi t ion to 2 D crack problems, the approach was also found to be enhance the numer ica l model l ing of 3 D excavations (see [27], [21]). T h e complete definit ion of this correct ion process was given by R y d e r and Nap ier [27] as follows: P r i o r to prepar ing prob lem input in the fo rm of fixed-size grids, concep-tual ly reduce the area m ined by mov ing al l faces inward by j - g r i d . Le t m be the result ing edge gr id fract ion m ined (0 < m < 1). Numer i ca l l y increase the normal ized edge self-effects by a = (1 — m)/m. For example, for a nonfract ional ly m ined gr id w i th D D E of length h, reduc ing the edge elements by \ implies that m = |. Therefore, and the edge effects must be mult ip l ied by 4 w = l + a = - , i.e., increased by approximate ly 33%. In general, let the gr id be of length h and the crack t ip element of length Xh. a — 1 3 3 4 In this case m = | A and since m G (0,1] therefore A € (0, | ] . T h e self-effect formula gives • ^ q u a r t e r - g r i d l — m m (3.27) SELF-EFFECT TIP CORRECTIONS FOR PIECEWISE CONSTANT DDE 26 so that the self-effect mult ip l ier is ^ q u a r t e r - g r i d — 1 ^ q u a r t e r - g r i d — • Numerically Determined Crack Tip Self-effect Correction Schemes T h e quarter-gr id correct ion scheme was analyt ical ly der ived in order to provide a reduct ion in crack w i d t h error in the case of a constant pressure d is t r ibut ion. In order to find other error reducing t ip element self-effect correct ion schemes a numerica l approach was considered. T h e approach used involves solving for the required self-effect t ip correct ion subject to the min imizat ion of a relevant error. Symmetr i c pressure distr ibutions w i th p(x) > 0 were considered which guarantees that the stress intensity factor is nonzero (see (2.19)) a n d that a squareroot crack t ip of the fo rm (3.14) exists. Far-field Co l locat ion Cor rec t ion Strategy In this case, the error of the D D E closest to x = 0 is min imized . In part icular , discretize the crack using 2N constant D D E each of length h = 2a= JJ. Le t u> = (1+a) be the self-effect mult ip l ier . Use the parametr ized 2N x 2N influence matr ix , A(a), to construct the N x N influence matr ix , B(a), for the posit ive interval x 6 (0,L) (see (B.6)) . T h e n solve (B.5) to find k = 1,...,N. Use - 2Uy(h) as a measure of the error of the solution near the midd le of the crack, where uy is the height of the upper crack face given by (2.8). In this fashion, for a given N, find a = ajv such that . E m i d d i e = 0. Since col location of the solut ion near the midd le of the crack (far f rom the stress singularity) is used to determine the required stress self-effect correct ion at the crack t ip , this approach has been called the "far-field" m i d d l e D: SELF-EFFECT TIP CORRECTIONS FOR PIECEWISE CONSTANT DDE 27 N 5 10 20 40 100 500 Po 0.2635 0.2636 0.2623 0.2621 0.2620 0.2620 P o ( l - f ) 0.2346 0.2476 0.2546 0.2581 0.2603 0.2620 Table 3.1: Alpha values for the far-field correction. correct ion scheme. Numer i ca l l y determined a corrections are shown in Table 3.1 for constant and l inearly decreasing pressure distr ibut ions. Based on these results, the far-field correct ion strategy is to choose ^far-field — 0.26 and therefore increase the t ip element self-effect by 26%. T i p Co l locat ion Cor rec t ion Strategy T h e t ip col locat ion correct ion is constructed like the far-field col locat ion cor-rect ion, w i th the exception that the error min imized is that of the D D E closest to the t ip . Numer i ca l l y determined a corrections for the t ip col locat ion m e t h o d are shown in Table 3.2 for constant and l inearly decreasing pressure distr ibut ions. Based on these results, the t ip col location strategy is to choose ^tip-collocation 0.20 (3.28) and therefore increase the t ip element self-effect by 20%. A Scaled Quarter-grid Correction Unl ike the quarter-gr id correct ion, the numerica l ly determined col locat ion cor-rections do not include the fractional min ing parameter A . However, it was found that the col location corrections prov ided better accuracy than the quarter-gr id cor-A FRACTIONALLY MINED CRACK TIP ELEMENT 28 N 5 10 20 40 100 500 Po 0.2065 0.2027 0.2008 0.1999 0.1993 0.1990 P o ( l - f ) 0.1801 0.1889 0.1938 0.1963 0.1979 0.1987 Table 3.2: Alpha values for the tip collocation correction, rect ion in the case A = 1. In an attempt to construct a col location self-effect correc-t ion that inc luded A dependence^ the following scaled quarter-grid-l ike corrections were constructed: "far-field (A) "tip-collocation (A) ^quarter-grid (A) "quarter-grid (0) ^quarter-grid (A) "far-field "tip-collocation (3.29) ^quarter-grid (0) A Fractionally Mined Crack Tip Element Consider the case of a crack which has a length that is not a mult ip le of a constant gr id of size h. T h i s is refered to as a fractional ly m i n e d crack. H a v i n g a technique, such as the quarter-gr id correct ion scheme, that allows the accurate model l ing of a fractional ly m ined crack on a constant gr id is very useful. In ad-dit ion, the fact that the quarter-gr id scheme only requires the modif icat ion of the single crack t ip self-effect influence coefficient allows the m e t h o d to be generalized to problems of higher dimensions. It would therefore be useful to develop other correct ion strategies that allow fractional min ing yet only require the modif icat ion of local crack t ip influence coefficients. Since crack t ip self-effect correct ion strate-gies were considered in the previous section, now consider construct ing a crack t ip influence correct ion for the neighbour ing element. T h e strategy used wi l l be to construct a t ip element that includes an extension element for the fract ional ly A FRACTIONALLY MINED CRACK TIP ELEMENT 29 mined crack part . T h e extension element wil l produce an extra contr ibut ion to the ord inary constant D D E . B y l imit ing the influence of the extra contr ibut ion to the t ip and neighbouring elements, a local correct ion is achieved. F inal ly , a self-effect correct ion on the t ip element, s imilar to the t ip col location scheme, leads to improved accuracy. Cons ider a piecewise constant t ip element of length Xh = 2Xa for 1 < A < 2. T h e "ma in" part of the element wil l be of w id th — Dy and length h = 2a. T h e "extension" part of the element is of length (A — l)h and lies between the m a i n element and the actual crack t ip. T h e w id th of the extension element is to be determined by the asymptot ic requirement that the D D E width at the midpoints of the m a i n a n d extension elements is consistent w i th the scaling \D\ <x y/x, as x -> 0. (3.30) Specifically, the piecewise constant D D E w i d t h is chosen as Uy(x) = ux(x, 0_) - ux(x, 0+) = < ADy, x 6 (0,(A- l)h), Dy, x e ((A - l)h, Xh) 0, otherwise, where A = V 2 ^ T <3'31> T h e ma in crack opening is therefore located in the x > Xh region. Since the variable length t ip element consists of two piecewise constant D D E , the result ing induced norma l stress m a y be calculated as the s u m of two stress contributions, ind iv idual ly given by (3.10). Denot ing the midpo int of the m a i n element by xm = (A — \)h, the midpoint of the extension element by xe = ^^h and the length of the extension element by he = 2ae, the n o r m a l stress result ing f rom A FRACTIONALLY MINED CRACK TIP ELEMENT 30 Figure 3.1: The fractional mining tip. the crack t ip is ayy(x,0) = (stress from " m a i n " element) + (stress from "extension" element) —G a ^ \ / — G ae , „ 7r(l — v)(x — xm)2 — a2 J V7r(l ~ u) (x ~ xe)2 — a i ^ I \2 f "^V 7r(l - v) \ (x- xm)2 - a2 (x- Xg)2 - a2 } U s i n g this result, the influence coefficients for the crack t ip are readi ly obtainable for a part icular noda l discretization, x = xk, k = 1 , N . A local t ip element is then achieved by only al lowing the extension element to influence the t ip element and the neighbour ing element. T h e self-effect correct ion strategy used only modifies the self-effect of the m a i n element. T h e m a i n element's self-effect is mult ip l ied by co = l + a(A) (3.32) ^tip-collocation = i + A N u m e r i c a l experiments suggested that when the extension element's corrective stresses were l imited to the m a i n element and the neighbour ing element, that m o d -ification of the stress influence of the neighbour ing element on the t ip would also A FRACTIONALLY MINED CRACK TIP ELEMENT 31 aid the correct ion. T h e addit ional corrective stress that that was chosen for the one neighbour ing D D E appl ied to the crack t ip was —G a 7r(l — v) (x — xe)2 — a2 W h i l e this strategy is ad hoc, it is consistent w i th the t ip col location method when A = 1 and good numer ica l solutions are obta ined for A € [1,2]. A l so , since the t ip correct ion scheme affects only the nearest neighbours, the m e t h o d is consistent w i t h the constant D D M in the l imit that N —> oo. Chapter 4 A Galerkin Finite Element Method T h e Ga le rk in approach is to represent the crack w i d t h as the s u m of basis functions wi th compact support and then to use (2.16), or the equivalent equat ion on (—L,L), to construct a discrete relat ion between crack w id th and pressure. A n approximate weak solution of (2.15) may then be constructed by tak ing the inner product of b o t h sides wi th suitable test functions, v(x), and then integrat ing by parts. Specifically, start from (2.15), which is equivalent to (2.16), then take the inner product w i th v(x') (4.1) U s i n g the fact that (4.1) can be integrated by parts: 32 PIECEWISE LINEAR (PWL) BASIS FUNCTIONS 33 Jo \ 2x x-x' ) M pL rL -I dx I dx' i" Jo Jo X — X x + x' X — x' dv(x') j dx' J dv(x') duy(x) dx' dx (4.2) where the condit ion that v(L) = 0 was required. F inal ly , wr i t ten in terms of the ful l crack width , w(x) = 2uy(x), equat ion (4.2) becomes rL/ rL rLi \ p(x')v(x')dx' = — I dx dx' In Jo 4 7 r Jo Jo x + x' x — X' dv(x') dw(x) (4.3) dx' dx Note that while this representation is in terms of crack w id th derivatives, a non-derivative formulat ion can also be constructed using (2.16). However, (4.3) is conve-nient for use w i th the piecewise l inear basis functions given in the following section. Piecewise Linear (PWL) Basis Functions A s s u m e that the interval (0, L) has been part i t ioned by the points {xk}K=Q, where 0 = XQ = Xi < X2 < • • • < Xpf-\ < X J V < XM+\ = L Let Afc = {x\x € (xk,Xk+i)}, k — 0,1,...,AT. N o t e that x0 — 0 a n d Ao = {0} have been introduced for notat ional convenience. T h e set of piecewise l inear basis functions wi l l be {<Pk(x)}k-i for x 6 (0,L). These functions wi l l be consistent w i th the crack w id th boundary condit ion Uy(L) = 0. A suitable choice of basis functions are the hat-like functions mk-i(x - x f c _i) , l e A n , <t>k(x)= { m k ( x k + i - x ) , x e A f c , , fc=l,2,..,iV, (4.4) 0, otherwise, where 0, k = 0, ( x f e + x - X f e ) - 1 , fc = l,...,iV. (4.5) PWL BASIS FUNCTIONS WITH A CRACK-TIP CORRECTION 34 Not ice that <f>(xk) = 1, k = 1, ...,N. Accord ing ly , the crack solution wi l l be approxi-mated at {xk}k=i-PWL Basis Functions with a Crack-tip Correction F r o m the D D M , it is known that incorporat ing the square root dependence of the w id th near the crack t ip (see (2.5)) can lead to a significant improvement in the accuracy of the w id th solution. T h e Ga le rk in approach also benefits from the use of a special square root t ip element near the crack t ip . Le t |0fc(x)| be the new basis that incorporates the special t ip . Specifically, define c ( x ) = <j>k{x), k — 1 , N - 1, (4.6) and mAr_i(x - X J V _ I ) , -. . x e A J V _ I , mN^/(L - x)(L - X J V ) , x e A J V , 0, otherwise (4.7) Matrix Equations U s i n g the basis functions {<f>k(x)}k=i ^or {^ f c ( ; r )} f c ^ *s possible to approxi-mate the crack w id th w(x) by the s u m N n = l T h e n , by lett ing v(x) = </>n(x), n = 1 , 2 , N , (4.3) yields the system of equations I p(x)d>n(x)dx = > — / dx' / dxhx J o y y ' V y ' ^[^Jo Jo B y mak ing the approx imat ion X + x' x — X' d<pn(x') # m ( x ) dx' dx Wm, n = l,2,...,N. (4.8) I p(x)<pn(x)dx « Pn / 4>n(x)dx, Jo Jo STRAIN ENERGY CALCULATION METHOD 35 with PN = p(xn), then (4.8) becomes the linear system of equations MP = SW, (4.9) where M is the N x N d iagonal weight matr ix , Mnm = 6nm / <pn(x)dx, (4.10) Jo S is the N x N stress influence matr ix , E ' fdx' fL Jo Jo Snm = T~ I dx I dxlll 47T X + x' X — X' d(pn(x') d<pm(x) dx' dx ' (4.11) and P = (P1,P2,---,PN)T, w = (W1,W2,---,WN)T. If the square root t ip corrected basis functions are used, all <p(x) functions are replaced by 4>(x). T h e weight and stress matrices are calculated in A p p e n d i x D.. N u m e r i c a l results are presented in the next chapter. Strain Energy Calculation Method W h e n calculat ing the strain energy of a half-crack, a direct approach is to approximate (3.16). For example, this is done in the next chapter when est imating stress intensity factors using the strain energy release rate m e t h o d (see p.42). T h e integral has been est imated as / p(x)w(x)dx « / p(x) I Wk4>k(x) ) dx N . L w y^p(xk)Wk / (j>k(x)dx N = Z~^p{xk)WkMkk fe=i STRAIN ENERGY CALCULATION METHOD 36 where the definit ion of the weighting mat r ix (4.10) has been used. O n a constant gr id, this is consistent w i th the midpoint rule, except for the endpoints which are sl ightly different. W h e n the t ip-corrected basis functions are used, the weighting matr ix changes accordingly. Chapter 5 Comparison of Numerical Methods T h i s chapter presents numerica l results for the various crack w id th methods. Relative Error Notation T h e relative error of a w id th wk at noda l posit ion xk is denned to be p _ Wk - WexactpCfc) _ Wexact(Zfc) where Nmax is the number of nodes for the method . T h e error at the crack t ip wi l l be calculated using Eup = Effmax. T h e error near the middle of the crack wi l l be est imated using Emid = E\. A n average error is calculated using •« m a x Jvmax T h e fact fact that the w id th methods have differing noda l numbers , i V m a x , a n d posit ions should be kept in m i n d when compar ing errors. 37 CRACK WIDTH RESULTS 38 Shortened N a m e Descr ip t ion PWCC piecewise constant collocation D D M (ordinary DDE) PWCC-tip PWCC with squareroot tip element (tip length = ordinary DDE length) PWCC-QGrid PWCC with quarter-grid self-effect correction (a = |) PWCC-far-field PWCC with far-field self-effect correction (a = 0.26) PWCC-tip-coU PWCC with tip collocation self-effect correction (a = 0.20) Galerkin Galerkin with {<pk)k=i basis functions Galerkin-tip Galerkin with {<pk\k-i basis functions P W L C piecewise linear collocation D D M (linear DDE) PWLC-tip P W L C with squareroot tip element (tip length = linear element length) Table 5.1: Naming conventions adopted for crack width methods. Crack Width Results T h i s section presents the results of the crack w i d t h methods for compar ison. For convenience a shortened name has been assigned to each of the methods and is shown in Table 5.1. It is important to compare the methods fairly w i th respect to computat iona l effort. T h i s means that methods wil l be compared based on the number of nodes (degrees of freedom) that are used and not the number of elements. Le t N be the number of nodes on the interval [0, L), w i th the except ion of the P W L C - t i p m e t h o d which has N — 1 nodes. T h e methods are demonstrated for the constant pressure case. T h e accuracy trends that exist for this case tend to also continue for other pressure distr ibutions that have p(x) > 0, and therefore a nonzero stress intensity factor. For the sake of brevity only the constant pressure results are presented. T h e results of the ord inary element based methods ( P W C C , P W C C - t i p , P W C C -QUARTER-GRID METHODS AND FRACTIONAL MINING 39 Emid (% or units of 0.01) N 4 10 20 40 100 PWCC +6.54 +2.54 +1.26 +0.627 +0.250 PWCC-tip +4.55 +1.61 +0.778 +0.382 +0.151 PWCC-QGrid -1.20 -0.476 -0.239 -0.120 -0.0479 PWCC-far-field +0.0686 +0.0201 +0.00835 +0.00371 +0.00137 PWCC-tip-coll +1.25 +0.482 +0.238 +0.118 +0.0471 Galerkin -2.29 -1.05 -0.551 -0.282 -0.114 Galerkin-tip +0.548 +0.0810 +0.0188 +0.00406 +0.000357 P W L C -1.32 -0.837 -0.519 -0.306 -0.146 PWLC-tip +5.00 -0.109 -0.385 -0.291 -0.151 Table 5.2: Mid-crack relative error for width calculation methods. QGrid, PWCC-far-field, PWCC-tip-coll) are plotted in Fig. 5.1 for N = 10 and Fig. 5.2 for N = 20. The low tip error of the PWCC-tip-coll method and the low mid-crack error of the PWCC-far-field method are clearly apparent. The higher-order PWL methods (Galerkin, Galerkin-tip, PWLC, PWLC-tip) are demonstrated in Fig. 5.3 for TV = 10 and Fig. for N = 20. The three relative errors, Emid, Eup and E, are shown in tables 5.2, 5.3 and 5.4, respectively. Note that part of the reason the P W L C method has such a large relative error at its tip node is because the tip node is closer to x = L than the tip nodes of any of the other methods. Quarter-grid Methods and Fractional Mining The quarter-grid method has already been demonstrated, for the case A = 1, in Fig. 5.1 and Fig. 5.2, where it was compared to the other P W C C related QUARTER-GRID METHODS AND FRACTIONAL MINING Etip (in % or units of 0.01) N 4 10 20 40 100 PWCC +27.4 +26.1 +25.7 +25.5 +25.4 PWCC-tip +9.55 +8.59 +8.31 +8.17 +8.08 PWCC-QGrid -11.4 -11.8 -11.9 -12.0 -12.0 PWCC-far-field -5.06 -5.54 -5:69 -5.76 -5.81 PWCC-tip-coll +0.878 +0.732 +0.0813 -0.0130 -0.0989 Galerkin +3.25 +2.25 +2.22 +2.21 +2.20 Galerkin-tip +2.84 +1.08 +0.782 +0.639 +0.555 P W L C +11.7 +11.9 +11.9 +12.0 +12.1 PWLC-tip +6.91 +1.08 +0.782 +0.639 +0.555 Table 5.3: Crack tip relative error for width calculation methods. QUARTER-GRID METHODS AND FRACTIONAL MINING E (in % or units of 0.01) N 4 10 20 40 100 PWCC +13.0 +6.32 +3.59 +2.01 +0.917 PWCC-tip +7.04 +3.47 +1.98 +1.12 +0.515 PWCC-QGrid +4.10 +1.95 +1.08 +0.586 +0.258 PWCC-far-field +1.34 +0.609 +0.318 +0.163 +0.0665 PWCC-tip-coll +1.33 +0.678 +0.407 +0.241 +0.118 Galerkin +2.87 +1.53 +0.962 +0.581 +0.286 Galerkin-tip +1.29 +0.293 +0.122 +0.0541 +0.0196 P W L C +3.91 +2.00 +1.19 +0.702 +0.345 PWLC-tip +6.53 +0.510 +0.387 +0.348 +0.232 Table 5.4: Average relative error for width calculation methods. FRACTIONAL MINING TIP CORRECTION 42 methods. It can be seen that, while the P W C C - Q G r i d result is clearly better t h a n the P W C C result, it sti l l has a round 12% relative error at the crack t ip . However, the scaled quater-gr id methods, which use (3.29), do not have this proper ty since they are identical to the P W C C - f a r - f l e l d and P W C C - t i p - c o l l methods when A = 1. N u m e r i c a l results for the three quarter-gr id related methods are shown in F i g . 5.5 for A =0.8, 1, 1.2. It should be noted that the quarter-gr id correct ion was designed for A G (0, | ] . W h i l e the two scaled quarter-gr id methods are b o t h very accurate for A = 1 and do not produce such large t ip errors as the P W C C - Q G r i d method , it becomes clear for different A that the corrections do not eradicate the error. It should be noted that, a l though the errors are significant for al l three methods, they are sti l l comparable to the error of the P W C C method . Based on the evidence it may be mi ld ly favourable to use the 0.26 scaled quarter-gr id m e t h o d instead of just quarter-gr id. Fractional Mining Tip Correction Some numerica l results for the fractional m in ing t ip correct ion (see p. 28) are presented in F i g . 5.6 for A = 1 . 0 , 1.2, 1.5, 2.0 for N = 10. Results for the quarter-gr id method are also inc luded for A = 1 . 0 and 1.2. Smal l crack t ip error for various A is a feature of the fractional m in ing t ip that is not shared by the quarter-gr id-like corrections. A l t h o u g h the overall error is gett ing worse as A approaches 2, the solution is still reasonably accurate. However, a notable point is that good t ip error results are obtained for al l the A values in the figure. Stress Intensity Factor Estimation T h i s section provides stress intensity factors calculated using various numer ica l crack w i d t h methods. T h e stress intensity factors are est imated using b o t h the STRESS INTENSITY FACTOR ESTIMATION 43 Stress Intensity Factor (units of POA/TTX) N ' 4 10 20 40 100 PWCC 1.2339 1.2455 1.2494 1.2514 1.2525 PWCC-tip 1.0607 1.0722 1.0763 1.0783 1.0795 PWCC-QGrid 0.8576 0.8711 0.8755 0.8777 0.8790 PWCC-far-field 0.9193 0.9327 0.9372 0.9394 0.9407 PWCC-tip-coll 0.9767 0.9901 0.9945 0.9967 0.9981 Galerkin 0.9577 0.9966 1.0093 1.0156 1.0194 Galerkin-tip 0.9549 0.9852 0.9951 1.0001 1.0030 P W L C 1.0967 1.1103 1.1153. 1.1180 1.1198 PWLC-tip 1.0001 0.9921 0.9923 0.9930 0.9937 Exact 0.9682 0.9874 0.9937 0.9969 0.9987 Table 5.5: Stress intensity factors calculated using the asymptotic method. asymptot ic and energy formulations. T h e popular case of a pressurized crack w i t h constant pressure, po, and tota l length 1L is considered. Since the stress intensity factor is given exactly by (2.24), all numer ica l results are stated in units of poy/nL. A lso , in addi t ion to the methods listed in Table 5.1, the exact w id th at the P W C C nodal points is used (called Exac t ) . T h i s " m e t h o d " is used to provide a reasonable indicat ion of the accuracy obtainable w i th the part icular stress intensity factor est imat ion scheme employed. Asymptotic Method For each of the w id th methods, the stress intensity factor is calculated from (3.15) using the noda l w id th magnitude corresponding to the noda l point closest to STRESS INTENSITY FACTOR ESTIMATION 44 the t ip . C lear ly accuracy depends on having smal l Eup and being close to the crack t ip . B u t , as ment ioned elsewhere, not al l methods have t ip nodes the same distance from the crack t ip . T h e results for the various w i d t h methods are presented in Table 5.5. T h e results essentially follow the t ip accuracy trends encountered previously. However, it is clear that the one-term asymptot ic m e t h o d is l imi ted to about a tenth of a percent relative error, for a reasonable number of elements. T h e most accurate methods are P W C C - t i p - c o l l , Ga lerk in- t ip and P W L C - t i p . Strain Energy Method For each of the w id th methods, the stress intensity factor is calculated using (3.18). Since p = po, calculat ion of the stra in energy is simplif ied. St ra in energies are calculated using the midpoint rule for the P W C C and E x a c t methods, G a l e r k i n methods use the weighting matr ix approach, P W L C integrates the l inear elements exactly, and P W L C - t i p is like P W L C except that the t ip element contr ibut ion is in-tegrated exactly. N u m e r i c a l results shown in Table 5.6 are generally more accurate than the asymptot ic results. However, two matr ix inversions are required in the energy m e t h o d compared to only one in the asymptot ic case. T h e most accurate methods are clearly P W C C - f a r - f i e l d and Galerk in-t ip . A part ia l explanat ion for the less than spectacular P W L C - t i p result is that the centred difference derivative approx imat ion uses a step size of 2h, instead of h, and therefore l imits accuracy. STRESS INTENSITY FACTOR ESTIMATION Stress Intensity Factor (units of poy/nL) N 4 10 20 40 100 P W C C 1.030776 1.012423 1.006231 1.003120 1.001249 PWCC-tip 1.019838 1.007790 1.003856 1.001914 1.000760 PWCC-QGrid 0.993992 0.997600 0.998800 0.999400 0.999760 PWCC-far-field 1.000166 1.000065 1.000032 1.000016 1.000006 PWCC-tip-coll 1.005868 1.002349 1.001175 1.000588 1.000235 Galerkin 0.990800 0.995486 0.997555 0.998714 0.999463 Galerkin-tip 1.000063 0.999983 0.999989 0.999994 0.999998 P W L C 0.994618 0.996675 0.997915 0.998748 0.999389 PWLC-tip 0.999698 0.996063 0.997583 0.998580 0.999322 Exact 1.001754 1.000437 1.000154 1.000054 1.000014 Table 5.6: Stress intensity factors calculated using the strain energy method. STRESS INTENSITY FACTOR ESTIMATION Figure 5.1: Plots of constant DD and related methods (N=10). STRESS INTENSITY FACTOR ESTIMATION 47 Figure 5.2: Plots of constant DD and related methods (N=20). STRESS INTENSITY FACTOR ESTIMATION 48 Figure 5.3: Plots of Galerkin and P W C L methods (N=10). STRESS INTENSITY FACTOR ESTIMATION Figure 5.4: Plots of Galerkin and PWCL methods (N=20). STRESS INTENSITY FACTOR ESTIMATION a: o CL a: LU LU > r-0.4 0.5 0.6 X, LAMBDA =1 ill CU •0.02 • •0.04 • •0.06 -0.08 • t o + x O Quarter-grid ' IQu o o Scaled QTuarter-grid 0 + + Scaled Quarter-grid (i .23 0:2 y 0.1 0.2 0.3 0.4 0.5 0.6 X, LAMBDA = 1.2 0.7 0.8 0.9 Figure 5.5: Quarter-grid and scaled variants (N=10). STRESS INTENSITY FACTOR ESTIMATION 0.4 0.5 0.6 X, LAMBDA = 1 0.05 cc O DC CC o UJ tu > <j -0.05 UJ cc -0.1 I •at 1 1 1 4fc i * i * o o o o O o o O O 1 1 i 1 1 i o ) 0.1 0.2 0.3 0.4 X, 0.5 LAMBDA 0.6 = 1.2 0.7 0.8 0.9 1 1 1 1 1 i i i * 1 * 1 * * i 0.05 cc o cc cc 5 UJ cc -0.05 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 X, LAMBDA = 1.5 I 1 1 I I 1 1 1 1 * * * * * * * * CC o CC CC 5 l i l CC -0.05 0.1 0.2 0.3 0.4 0.5 0.6 X, LAMBDA = 2 0.7 0.8 0.9 Figure 5.6: Fractional mining tip results (N=10). Chapter 6 Conclusions T h i s research has addressed the compar ison of D D M and Ga le rk in B E M for the pressurized line crack. T h e construct ion of various t ip correct ion strategies for constant D D E on a gr id of constant spacing was also addressed. For the case of constant D D element discret izat ion of a crack on a constant gr id , numerica l evidence has been found which suggests that various modif icat ions of the crack t ip element's self-effect have the abi l i ty to significantly reduce the crack w i d t h solution's error. T h i s result was known previously, and a similar modi f icat ion cal led the quarter-gr id correct ion had been developed [27]. However, it has been found that other corrections y ie ld comparable or better accuracy. Results indicate that increasing the self-effect by approximate ly 20% reduces the t ip error significantly, while an increase of approximate ly 26% reduces the error near the midd le of the crack. For compar ison, the quarter-gr id correct ion of R y d e r a n d Nap ie r suggests a one th i rd increase in the self-effect. It appears that the correct ion strength is not strongly dependent of the number of D D elements or the pressure d istr ibut ion. For example, in the case of the t ip col location method (20% correct ion), errors are less than hal f of those encountered using a squareroot t ip element. It is hoped that the method may be generalized to more compl icated problems and that benefits of the constant D D element crack p rob lem formulat ion may be exploited w i th greater accuracy. 52 53 W e at tempted to generalize the self-effect correct ion m e t h o d to the fract ional m in ing prob lem. Rescaled quarter-grid-l ike corrections were constructed which cor-responded to the 20% and 26% self-effect increases when there was no fract ional min ing . W h i l e the corrections produced results which had some favourable prop-erties, it was found that relative errors in excess of 5% would have to be to lerated for al l the methods. C e r t a i n scalings produced better errors for different length fractions, but only minor improvement was obtained overall . A n improved crack t ip element correct ion was constructed for the constant D D E , constant gr id spacing fract ional min ing prob lem. T h e m e t h o d corrects not only the crack t ip element's self-effect but also the t ip element's stress influence on the adjacent element and vice versa. T h e correct ion allows the fractional ly m ined crack w i d t h to be calculated wi th errors generally m u c h less t h a n 3%. A notable point is that the t ip error is never m u c h larger than 1%. T h e results of higher-order Ga le rk in and P W L col locat ion ( P W L C ) methods were also compared . It was found that the G a l e r k i n and P W L C methods p roduced solutions wi th comparable accuracy. T h e t ip corrected versions (Galerk in-t ip a n d P W L C - t i p ) were also comparable in accuracy. However, there was a tendency for the t ip-corrected Ga le rk in m e t h o d to be extremely accurate in regions away f rom the crack tips, e.g., w i th 100 nodes, the relative error near the middle of the crack was a few 100 t imes less than comparable Ga lerk in , P W L C a n d t ip-corrected P W L C results. In the absence of t ip corrections, the P W L C m e t h o d is essentially equivalent to the Ga le rk in method , on a constant gr id of equal spacing, w i th the exception that the Ga le rk in m e t h o d has slightly more accuracy at the crack t ip . F inal ly , the methods were also used to calculate stress intesity factors using as-ymptot i c and stra in energy release rate formulas. M e t h o d s that p roduced accurate crack w i d t h est imat ion near the crack t ip, natura l ly p roduced the best stress inten-54 sity factor approximations when using the asymptot ic formula. T h e most accurate methods in that category were the t ip col location m e t h o d (constant D D E w i th 20% self-effect increase) and bo th the G a l e r k i n and P W L C methods w i th squareroot t ip corrections. W h e n the stra in energy release rate formula was used, methods which produced the most accurate overall crack w i d t h were the most accurate. In this category, the Ga le rk in method wi th a squareroot t ip correct ion was the most accurate, followed by the constant D D M wi th the 26% self-effect increase chosen to min imize the error near the middle of the crack. Bibliography [1] S. H . Advani, T. S. Lee and R. H . Dean, "Variational Principles for Hydraulic Fracturing," Journal of Applied Mechanics, 59, pp. 819-826, 1992. [2] G. I. Barenblatt, "Mathematical Theory of Equilibrium Cracks," in Advances in Applied Mechanics, 7, pp. 55-129, 1962. [3] R. S. Carbonell, "Self-Similar Solution of a Fluid-Driven Fracture", Ph.D. The-sis, University of Minnesota, 1996. [4] R. J. Clifton and A. S. Abou-Sayed, On the Computation of the Three-Dimensional Geometry of Hydraulic Fractures, paper SPE 7943, 1979. [5] A. M . Crawford and J. H . Curran, "Higher-order functional variation displace-ment discontinuity elements," Mech. Min. Sci. and Geomech. Abstr., 19, pp. 143-148. [6] S. L. Crouch, Solution of plane elasticity problems by the displacement discon-tinuity method. Int. J. Num. Methods Engng 10, 301-43. [7] S. L. Crouch, Analysis of stresses and displacements around underground ex-cavations: an application of the displacement discontinuity method. Geome-chanics report to the National Science Foundation. Minneapolis: University of Minnesota. [8] S. L. Crouch and A. M . Starfleld, Boundary element methods in solid mechan-ics, Cambridge, Unwin Hyman Inc., 1990. 55 BIBLIOGRAPHY 56 [9] J. Desroches, E. Detournay, B. Lenoach, P. Papanastasiou, J. R. A . Pearson, M . Thiercelin and A. Cheng, "The crack tip region in hydraulic fracturing," Proc. R. Soc. London A, 447, pp. 39-48, 1994. [10] J. Desroches and M . Thiercelin, "Modelling the Propagation and Closure of Micro-Hydraulic Fractures", Int. J. Rock Mech. Min. Sci. &; Geomech. Abstr., 36, No. 7, pp. 1231-1234, 1993. [11] A. H . England and A. E. Green, Some Two-dimensional Punch and Crack Problems in Classical Elasticity, Proc. Cambridge Phil. Soc, London, 59, 489-500, 1963. [12] A. D. Fitt, A. D. Kelly and C. P. Please, "Crack Propagation models for rock fracture in a geothermal energy reservoir," SIAM J. Appl. Math., 55, No. 6, pp. 1592-1608, 1995. [13] L. B. Freund, Dynamic Fracture Mechanics, Cambridge University Press, New York, 1990. [14] J. Geertsma and F. de Klerk, " A Rapid Method of Predicting Width and Extent of Hydraulically Induced Fractures," Journal of Petroleum Technology, pp. 1571-1581, December 1969. [15] A. E. Green and W. Zerna, Theoretical Elasticity, Oxford University Press, London,1968. [16] A. A. Griffith, The Phenomena of Rupture and Flow in Solids, Phil. Trans. Royal Soc. London, Ser. A. 221, 163-198, 1920. [17] A. A. Griffith, The Theory of Rupture, Proc. First international Congress on Applied Mechanics, Delft, 55-63, 1924. [18] M . F. Kanninen and C. H . Popelar, Advanced Fracture Mechanics, Oxford University Press, New York, 1985. [19] S. A. Kristianovitch and Y. P. Zheltov, Formation of Vertical Fractures by Means of Highly Viscous Fluids, Proc. World Pet. Cong., Rome, 2, 579, 1955. BIBLIOGRAPHY 57 [20] N. I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity; Fundamental Equations, Plane Theory of Elasticity, Torsion and Bending, 5th ed., Noordhoff, Groningen, Leyden, The Netherlands, 1953. [21] J. A. L. Napier, Personal communication, 1997. [22] V. Z. Parton and E. M . Morozov, Elastic-Plastic Fracture Mechanics, Mir, Moscow, 1978. [23] R. P. Nordgren, "Propagation of a Vertical Hydraulic Fracture," Society of Petroleum Engineers Journal, pp. 306-314, August 1972. [24] T. K. Perkins and L. R. Kern, "Widths of Hydraulic Fractures," Journal of Petroleum Technology, pp. 937-949, September 1961. [25] S. T . Raveendra and T. A. Cruse, " B E M Analysis of Problems of Fracture Mechanics," in Developments in Boundary Element Methods - 5: Industrial Applications of Boundary Element Methods, P. K. Banerjee and R. B. Wilson (ed.), Elsevier Science, New York, pp. 187-204, 1989. [26] J. R. Rice, "Mathematical Analysis in the Mechanics of Fracture," in Fracture, H . Leibovitz (ed.), Vol. 2, Academic Press, New York, 1968. [27] J. A. Ryder and J. A. L. Napier, "Error analysis and design of a large-scale tab-ular mining stress analyser," Proceedings of the fifth international conference on numerical methods in geomechanics, Nagoya, pp. 1549-1555, 1985. [28] The SCR Geomechanics Group, "On the Modelling of Near Tip Processes in Hydraulic Fractures," Ing. J. Rock Mech. Min. Sci. & Geomech. Abstr., 7, pp. 1127-1134, 1993. [29] R. A. Schmidt, A microcrack model and its significance to hydraulic fracturing and fracture toughness testing, Proc. 21st US Symp. on Rock Mech., pp. 581-590. [30] I. N. Sneddon, Elements of Partial Differential Equations, McGraw-Hill, New York, 1957. BIBLIOGRAPHY 58 [31] I. N. Sneddon, Integral Transform Methods, in Mechanics of Fracture I, Meth-ods of analysis and solutions of crack problems, G. C. Sih (ed.), Nordhoff International, Leyden, 1973. [32] I. S. Sokolnikoff, Mathematical Theory of Elasticity, pp. 179-180, McGraw-Hill, New York, 1956. [33] D. A. Spence and P. Sharp, "Self-similar solutions for elastohydrodynamic cavity flow," Proc. R. Soc. London A, 400, pp. 289-313, 1985. [34] S. P. Timoshenko and J. N. Goodier, Theory of Elasticity, 3rd. ed., McGraw-Hill, New York, 1970. [35] D. J. Unger, Analytical Fracture Mechanics, London, Academic Press, Inc., 1995 [36] P. Valkd and M . J. Economides, Hydraulic Fracture Mechanics, New York, John Wiley and Sons, Inc., 1995. [37] B. N. Whittaker, R. N . Singh and G. Sun, Rock Fracture Mechanics: Princi-ples, Design and Applications, Elsevier Science Publishers, New York, 1992. [38] Y. P. Zheltov and S. A. Kristianovitch, On the Mechanism of Hydraulic Frac-ture of an Oil-Bearing Stratum, Izv. AN SSSR, O T N , (No. 5), 3-41, 1955. Appendices A . Stresses and Displacements for a D D of Finite Length This Appendix restates the results given by Crawford and Curran [5] with due reference to Crouch [6]. Consider a general D D E in the xy-plane of length 2a, located along the x-axis from —a to a. Assume, as usual, that the elastic medium is linear, homogeneous and isotropic. If the normal and shear displacements of the DD are given by DN(X) and Ds(x), respectively, then all stresses and displacements in the medium can be explicitly defined. Defining the normal and shear DD contributions to the stress and displacement fields are given by the following expressions. Normal Displacement Discontinuity, Dj^(e) df(x,y) dx 47r(l — u) x2 + y2 i y 1 X df{x,y) dy 47r(l — u) x2 + y2 u. •x U, •y d 1 ra df(x-e,y) dy DN(e)de 59 B.. INFLUENCE MATRIX SIMPLIFICATION FOR THE SYMMETRIC PROBLEM 60 ®xy — 2G Vdy2 I J —a a df(x-e,y) dx Dpf(e)de Shear Displacement Discontinuity, Ds(e) r° df(x-e,y) <7x 2 ( 1 I * d_ r df(x-e,y) Vdx J_ Ds(e)de dx a df(x-e,y) Ds(e)de +y 2G 2G 2G d_ r df(x-e,y) dy J-a dx " "I r0, J J—a Ds(e)de 2— — dy dy2 Ds(e)de a df(x-e,y) dx Ds(e)de /a -a _ d^_ Vdy2 d_ dy dy2 a df(x--e,y) dx Ds(e)de f J —a df(x-s,y) dy Ds{e)de B. Influence Matrix Simplification for the Symmetric Problem Assume that N is even, p(x) = p(—x), and that the D D E nodal points are sym-metrically located about x — 0, namely xk = —xw-k+i. In this case, the symmetry of the problem implies that the D D E solution with also be symmetric about x = 0, i.e., Dy = Dy~k+1. Let N = 2M. Expressing (3.12) as ' M M -pk = Y,A^Di + lZAk'M^Dy+^ k = l,...,M, (B.l) 3=1 3=1 M M -PM+k = J2AM+^Di + Z-ZAM+k-M+iD™+J' k = l,...,M, (B.2) j=l 3=1 use the fact that pk = p M + ( M + i - f e ) to rewrite (B.2) as M M -Pk = ~P2M+\-k = y^JA2M+i-k,jDj + y^yA2M+i-k,M+jDyA+:', k = l,...,M. (B.3) 3=1 3=1 Adding (B.l) and (B.3) yields M M -2pk = (Akj + A2M+i-k,j) Djy + zZ (Ak,M+3 + A2M+i-k,M+3) D y + i - (B-4) 3=1 3=1 C. PWLC MATRICES 61 Then, express (B.4) in terms of PM+k and DyVI+3-. ^ M ^ M PM+k = g J2 (AM+k,j + AM+1_kJ) Dl + ~Y^ (AM+k,M+j + AM+1_k<M+j) DM+i i=i j=i 1 M = 0 iZ (AM+k,M+l-J + AM+l-k,M+l-j) D¥+1~J + 2J=1 M 2 0 (AM+k,M+3 + AM+l-k,M+j ) D M + j Z 3=1 1 M = ~ ^ (AM+k,M+l-J + AM+l-k,M+l-j) D M + J + 1 J=l 1 M . 9 S~2 {AM+k,M+3 + AM+l-k,M+j) DM+j-1 3=1 where the fact that = D y v + 1 - k was used. Finally, this can be written as M -pk = j2BvDi (B-5) 3=1 where Bkj = g (AM+k,M+j + AM+l-k,M+j + AM+k,M+l-j + ^Uf+l-fc.M+l-j) (B.6) is the influence matrix for the positive interval and the circumflex indicates quan-tities on the positive a>axis, e.g., PM+k = Pk- Note that this approach works for appropriate non-DDM as well. C. P W L C Matrices This section details the construction of the matrix equation for the piecewise linear collocation (PWLC) method. Consider a line crack from x = —L to x — L discretized with 2N PWL D D E of equal length h = 2a = jj. Therefore the line crack has 4iV nodal points along the crack. Using Gauss-Chebyshev collocation point positioning, the nodes are located at xk — —L + a 1 , , * - 2 + ( - l ) fc = l,2,...,4/V, C.. PWLC MATRICES 62 so that -L < X\ < X2 < • • • < X2N < 0 < X2N+1 < • • • < X 4 J V - 1 < X4N < L. Using (3.25), the normal stress at the j t h node is G IN ^ ^ m = l where D™ and D™ are the nodal displacements of the m t h D D E and Xjm = Xj — (centre of m t n DDE). Letting dk be the nodal displacement corresponding to the node located at xk, then the resulting 47V x 47V matrix equation is 4N ( J j — Ajkdk k=l where G F&k, (-1)&) 2n(l-vy K 3 G { a 2w(l-v) \x%-( - l ) f e V2 ^ ^ - + -^ln x2jk - a? 2a Xjk + a and Xjk is the signed distance from Xj to the centre of the element with a node at Xk* Xjk — Xj Xk and xk = -L + ^[2k-l + (-l)h+1] Also, for the symmetric pressurized crack problem, the 4Arx47V matrix equation may be converted to a 2N x 2N equation using the simple procedure given in Appendix D.. GALERKIN MATRICES USING PWL BASIS FUNCTIONS 63 D . Galerkin Matrices using P W L Basis Functions This section concerns the calculation of the stress and weighting matrices, Snm and Mnm respectively. Stress Matrix, Snm In order to calculate the stress matrix, Snm, given by (4.11), the derivatives of the PWL basis functions <fin(x), defined by (4.4) and (4.5), are found to be mn_!, x e A n _ i , 6 A„, 0, otherwise, d<t>n{x) dx where A n = {x\x € (xk, Xk+i)}. Therefore, the stress matrix is Snm — Air E 47T 4TT f dx' r<fain ^ d M x > ) d ^ { x ) Jo J0 x — x' dx' dx - [ dx' [ d x J x + x'\dMx')d4m(x) f dx' J m n _ i m m _ i - m „ m m _ i dx' dx — mn-iirirr [ dx' f A n _i J A / dx' / dx + mnmm I dx' / dx > In JA,i JAM-1 JAN JAM J dx x + x' x — x' 4-K Tfln Tfl m Cn} m } where and Cn,m — Cm,n — / dx / dx In JA„ JAM x + x' x — x' n,m = 1 , 2 , N Co,m = Cm,o = 0, m = 0,l,...,JV. By changing variables with X — Xn ~f" (xn-|-i Xn)fJ D.. GALERKIN MATRICES USING PWL BASIS FUNCTIONS 64 then Cn f dx' f JA„ JAr, dxln Jo Jo dr}]iL + CnrnV - V' where Xn There are three cases to be considered if the symmetry condition, C n > m = C m > n , is used. Case I (m — n)\ When n = m, 2,xn >o, and C n t n simplifies to (ann)2lnan + ( a n n + l ) 2 In ^ ° n n + 2\ + ( 2 f l n n + 3 ) ^ Cn,n = (Xn+1 ~ Xn)2 I dr{ Jo Jo drjlii a n n + V + V' 77-77' / \2 f (ann)2 In arm . / . ^2 i„ ( \/ann + 2^ = ( x n + i - x n ) z <j - h ( a n n + 1) In I ——-j-y- 1 +(2ann + 3) In v^nn + 2} Case II (m = n + 1): In this case, x n x m >o, D.. GALERKIN MATRICES USING PWL BASIS FUNCTIONS 65 >o, Xn+1 and therefore I "nm ~i~ ^ nrnH ~^ V Cn,m = (xn+1 - xn)(xm+x - xm) I drj' j dr/In Jo Jo may be integrated to yield 1 + cnmn - rf Cn,n+i = ^ ( x n + 1 - o ; n ) 2 { ( a + c) 2 lri(l + (a + c)- 1) + ( l + 2(a + c))hi(l + a + c) -(1 +a) 2 ln( l +a) - (1 + c) 2 ln(l + c) + a 2 l n a + c 2lnc} where a = an<n+1 and c = cn>n+1. Case III (m > n + 1): This case requires the most general integration. However, using the fact that bnm > 1, iJjnm + CnmV + V Cn,m = (xn+i ~ xn)(xm+1 - xm) drf dry In m>n+l Jo JO 1 + CnmV - Tj' = ^(xn+i - xnf {{a + c)2 ln(l + (a + c) _ 1 ) + (1 + 2(a + c)) ln(l + a + c) +(b + c - l ) 2 ln(6 + c - 1) - (b + c)2 ln(6 + c) - (6 - l ) 2 ln(6 - 1) -(1 + a) 2 ln(l + a) + b2 In 6 + a 2 lna} where a = a n m , 6 = 6 n m and c = c n O T. Weighting Matrix, Mnm The weighting matrix, Mnm, is given by (4.10) for the case of PWL basis func-tions. This may be simplified using (4.4) and (4.5): Mnm = 6nm / 4>n(x)dx Jo = 5nm\ / [mn-i(x - xn-i) + Sni] dx + / mn(xn+1 - x)dx > [A„_! Jhn ) c 1/1 c \Xn ~~ Xn—1 . c / \ , ^n+l — XT{ = o „ m < ( l - o„i) h o „ i ( x n — X n_lJ -I = -"2— {xn+i ~ xn-i + bnXxn) E.. GALERKIN MATRICES USING TIP-CORRECTED PWL BASIS FUNCTIONS 66 E . Galerkin Matrices using Tip-corrected P W L Basis Functions This section concerns the calculation of the stress and weighting matrices, Snm and Mnm, in the case that the tip-corrected PWL basis functions, {(j>k}k=ii are used. Stress Matrix, Snm Since the 4>k basis functions are identical to the 4>k except for the j = N case, the stiffness matrix, 4TT JO JO da; In x + x' x — x' # n (x ') d<j>m(x) dx' dx will be identical to Snm for n,m ^ N. Also, since Snm is symmetric, it will be sufficient to calculate the m = N column for n = 1 , 2 , i V . From (4.6) and (4.7), the derivatives of the basis functions are d4>n(x) dx m n - l , X £ A n _ i —mn, x e An 0, otherwise, n = l , 2 , . . . , jv - l , and #jy(x) dx mN-i, x e Ajv-i 0, otherwise, Substitution into the Snm expression, when n < N and m = N, yields SnN = — I dx ax In 4TT JO JO x + x' X — X' d<j>n(x') d(f>N(x) dx' dx = E f dx' f ax In ™ 7 A „ _ I + A „ JAN-I+AN X + x' X — X' d(f)n(x') d<f>N(x) dx' dx E' f — I m n - i m j v - i C i - i . j v - i - mn-\mNCn-\,N — mnmN_iCn,N-i +m. E.. GALERKIN MATRICES USING TIP-CORRECTED PWL BASIS FUNCTIONS 67 where C n m is as before and C N , N = CN-, * JAN J A, dx\l E——— In L — x x + x' x — x' In the case n = m = N, E' SNN = {mN-imN-iCN-i,N-i — 2 m j v - i T O j v C W _ i i j v + ^ivWjvC7v ,Ar| where CN, N 4 J AN J A, dx-L — XN :ln x + x' X — X' 4 N AN yJ{L-x')(X -X) The Cn,N integral may be solved. Changing variables and separating the loga-rithm, Cn,N = (xn+1 -xn)(L-xN) 2 (Zn+l - XN)(L - XN) [ dn' f dr](l-T])-i\n Jo Jo {Nn,N — A I . A T } , O-nN + C-nNV + V' where a njv, &njv and c„jv are the same as in the previous section and Nn,N DnM I dr( I dr](l - rj) ^ m (anN + CnNrj + rj') Jo Jo / drf / drj(l-r])~i]n\bnN + cnNri-r]'\ Jo Jo It follows that Nn,N = < 4C l + a + c tanh - l +(3a + 2c)ln( ) +3 l + a + c ln(l + a) - 7 j a + c tanh -l a + c and > i < J V - l 4c 6 + a - 1 tanh - l 6 + c - 1 6 + c tanh -l b + c +(36 +2c) In 6 - 1 31n(6-l) + 7 D N - l l N = - y + 4 v / 3 t a n h - 1 - ^ l n 2 . E.. GALERKIN MATRICES USING TIP-CORRECTED PWL BASIS FUNCTIONS 68 where a = anjv, b = 6njv and c = cn^ are shortened notations. While the CN,N integral has not been solved exactly, it may be simply approx-imated. Introducing a change of variables, CN,N = J [ D X ' [ * JAN JA* dx L — XN y/(L-x')(X-x) x + x' x — X' (L - xN ,2 /•! / dn' f d T / K l - r / X l - f j ) ] * l n Jo Jo aNN+V + v' \v-v'\ where 2x N L — X J V Notice that ajvjv >^ 1 when N 3> 1. Rewriting the CN,N integral as (L-xN)2 /-1 CN,N f drj' ^ [ ( l - ^ l "»/)]"* Jo Jo i l n 1 + ^ - ^ - +luaNN - l n | r / - r / | I L l a N N J J {ri+Ta-ra} with ri = / dr/ / dT7[(l-»/)(l-»7)rihi JO JO T 2 = / dr/ / d»7[(l-»/)(l-»/)]-* lnojvjv, JO JO T 3 = f drf I' dry [(1 - r/)(l - ,7)]-* In |ry - ry'| 7o JO straightforward integration yields T 2 = 41najvjv and Using the expansion r 3 = 81n2-12. 1 1 1 1 ln(l + x) = x - -x2 + -x3 - -xA + -x5 K. GALERKIN MATRICES USING TIP-CORRECTED PWL BASIS FUNCTIONS 69 in the integrand of the T i integral, an approximation can be generated for large 16 _ j _ 1 7 6 _ 2 256 _ 3 _ 7808 _ 4 l l - 3 a N N 4 5 a N N + Q 3 a N N 1 5 7 5 a N N + 346112 _ 5 _ 9029632 _ 6 1507328 _7 n f -s\ 51975 a " N 945945 a " N + 105105 & i V J V + ™ " -In practice, it was found that seven terms of the series were sufficient to approximate T i to six significant figures, assuming a constant grid and N > 5 . Weighting Matrix, M N M The diagonal weighting matrix, MNM, corresponding to the 4>k basis functions, is identical to MNM for n < N. The M J V J V element is given by M N N = j 4>N(x)dx Jo JAN-XXXN-XN-IJ J A N \ L - X N xN - X J V - I 2 2 + 3^-^ )
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A survey of boundary element methods for a two-dimensional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A survey of boundary element methods for a two-dimensional fracture model James, Douglas L. 1997
pdf
Page Metadata
Item Metadata
Title | A survey of boundary element methods for a two-dimensional fracture model |
Creator |
James, Douglas L. |
Date Issued | 1997 |
Description | This thesis surveys boundary element methods for solving a two-dimensional pressurized line crack problem in a homogeneous infinite elastic medium. Starting with an overview of linear elastic fracture mechanics, the exact solution to the pressurized line crack problem is stated and results for several pressure distributions are provided. Various numerical methods for solving the crack problem are then introduced, such as the displacement discontinuity method (DDM) using point collocation and a Galerkin method. Crack tip elements and higher-order DDM are discussed. Self-effect correction methods, stated and developed for the piecewise constant DDM on a uniform grid, are shown to significantly improve the numerical solution. An accurate DDM correction for modelling crack tip element extensions is also presented. Finally, numerical results for the various methods are given and stress intensity factors are also presented for comparison. |
Extent | 2711106 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-25 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080003 |
URI | http://hdl.handle.net/2429/6512 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1997-0552.pdf [ 2.59MB ]
- Metadata
- JSON: 831-1.0080003.json
- JSON-LD: 831-1.0080003-ld.json
- RDF/XML (Pretty): 831-1.0080003-rdf.xml
- RDF/JSON: 831-1.0080003-rdf.json
- Turtle: 831-1.0080003-turtle.txt
- N-Triples: 831-1.0080003-rdf-ntriples.txt
- Original Record: 831-1.0080003-source.json
- Full Text
- 831-1.0080003-fulltext.txt
- Citation
- 831-1.0080003.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080003/manifest