International Journal of Solids and Structures 43 (2006) 1669–1692 www.elsevier.com/locate/ijsolstr
On the uses of special crack tip elements in numerical rock fracture mechanics
Mohammad Fatehi Marji a,*, Hasan Hosseini_Nasab Amir Hossein Kohsary c,1
b a Faculty of Mining Engineering, Yazd University, Yazd, Iran Department of Industrial Engineering, Yazd University, Yazd, Iran c Department Mining Engineering, Yazd University, Yazd, Iran
b,1
,
Received 27 November 2004; received in revised form 22 April 2005 Available online 1 July 2005
Abstract Numerical methods such as boundary element methods are widely used for the stress analysis in solid mechanics. These methods are also used for crack analysis in rock fracture mechanics. There are singularities for the stresses and displacements at the crack tips in fracture mechanics problem, which decrease the accuracy of the numerical results in areas very close to the crack ends. To overcome this, higher order elements and isoperimetric higher order elements have been used. Recently, special crack tip elements have been proposed and used in most of the numerical fracture mechanics models. These elements can drastically increase the accuracy of the results near the crack tips, but in most of the models only one special crack tip element has been used for each crack end. In this study the uses of higher order crack tip elements are discussed and a higher order displacement discontinuity method is used to investigate the e?ect of these elements on the accuracy of the results in some crack problems. The useful shape functions for two special crack tip elements, are derived and given in the text and appendix for both in?nite and semi-in?nite plane problems. In this analysis both Mode I and Mode II stress intensity factors are computed . Some example problems are solved and the computed results are compared with the results given in the literature. The numerical results obtained here are in good agreement with those cited in the literature. For the curved crack problem, the strain energy release rate, G can be calculated accurately in the vicinity of the crack tips by using the higher order displacement discontinuity method with a quadratic variation of displacement discontinuity elements and with two special crack tip elements at each crack end. ? 2005 Elsevier Ltd. All rights reserved.
Corresponding author. Tel.: +98 351 8225328; fax: +98 351 8216900/0098 351 8210699. E-mail addresses: mfatehi@yazduni.ac.ir (M.F. Marji), hhn@yazduni.ac.ir (H. Hosseini_Nasab), kohsary@yahoo.com (A.H. Kohsary). 1 Fax: +98 351 8210699. 0020-7683/$ - see front matter ? 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2005.04.042
*
1670
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
Keywords: Rock fracture mechanics; Higher order elements; Special crack tip elements; Numerical methods
1. Introduction In the fracture analysis of brittle substances (like most of the rocks), the Mode I and Mode II stress intensity factors can be calculated by numerical methods using ordinary element. The accuracy of the numerical results are increased by using isoparametric elements and crack tip elements (Ingra?ea and Hueze, 1980; Ingra?ea, 1983; Blandford et al., 1982; Ingra?ea, 1987; Guo et al., 1990; Scavia, 1990; Hwang and Ingra?ea, 2004). Boundary element method is one of the powerful numerical methods and has been extensively used in fracture mechanics (Aliabadi and Rooke, 1991; Aliabadi, 1998). Displacement discontinuity method is an indirect boundary element method which has been used for the analysis of crack problems related to rock fracture mechanics and in most cases the problem of crack tip singularities has been improved by the uses of one crack tip element for each crack tip (Guo et al., 1992; Scavia, 1992; Shen and Stephansson, 1994; Scavia, 1995; Tan et al., 1996; Carpinteri and Yang, 1997; Bobet, 2001; Stephansson, 2002). Recently higher order elements have been used to increase the accuracy of the numerical results (Crawford and Curran, 1982; Shou and Crouch, 1995). In this paper the higher order displacement discontinuity method that was originally introduced for ?nite and in?nite crack problems (Shou and Crouch, 1995) is extended to the half plane problems with traction free surfaces (Crouch and Star?eld, 1983). The general series for the crack tip elements is discussed and the required modi?cation for implementation of the higher order special crack tip elements are given in Appendix A. The mixed mode stress intensity factors (i.e. for Mode I and Mode II fractures, which are the most commonly fracture modes occur in rock fracture mechanics) are numerically computed. As most of rocks are brittle and weak under tension the Mode I fracture toughness KIC (under plain strain condition) together with the maximum tangential stress fracture criterion (r-criterion) introduced by Erdogan and Sih are used to predict the crack propagation direction (Erdogan and Sih, 1963). The three fundamental fracture criteria, the maximum tangential stress criterion (or r-criterion), the maximum strain energy release rate criterion (or G-criterion) and the minimum strain energy density criterion (or S-criterion) or any modi?ed form of these three criteria (e.g. F-criterion which is a modi?ed form of G-criterion) have been mostly used to study the fracture behaviour of brittle materials (Ingra?ea, 1983; Broek, 1989; Whittaker et al., 1992; Shen and Stephansson, 1994). All of these criteria have demonstrated that a crack in a plate under a general in-plane load does not initiate and propagate in its original plane, but rather crack initiation take place at an angle with respect to it. The Mode II fracture toughness, KIIC predicted by r-criterion and G-criterion are smaller than the Mode I fracture toughness, KIC, but it is generally larger than the Mode I fracture toughness, KIC predicted by S-criterion depending on the material parameter of the Poisson?s ratio m (because S-criterion depends on m). The mixed modes I–II fracture problems in compression have been shown to be more complicated and also quite di?erent from those under tension. Various existing fracture criteria have been applied to study the fracture problems in compression but the results are poorly correlated to the existing experimental data (Whittaker et al., 1992). Although in brittle substances like rocks Mode II fracture initiation and propagation plays an important role under certain loading conditions and Mode I fracture toughness KIC is less than Mode II fracture toughness, KIIC, but due to the weakness (low strength) of rock under tension, the rock breaks due to tensile and in most cases the condition of KIC will prevail to that of KIIC under pure tensile, pure shear, tension-shear and compression-shear loading conditions. Recently, a lot of work has been done on the application of Mode I, Mode II and mixed Mode fracture theories for rock type materials (Ingra?ea, 1981; Atkinson et al., 1982; Huang and Wang, 1985; Zipf and Bieniawski, 1987; Ouchterlony, 1988; Swartz et al., 1988; Sun et al., 1990; Fowell, 1995; Pang,
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692 Table 1 Mode I fracture toughness KIC and Mode II fracture toughness KIIC for some typical rocks (after Whittaker et al. (1992)) Rock type Basalt Newhurst Granite Welsh Limestone Coarse-grained Sandstone Fine-grained Sandstone Dark gray Syenite Grayish white Syenite KIC (MPa m1/2) 2.27 1.72 0.85 0.28 0.38 1.75 1.36
1671
KIIC (MPa m1/2) 1.878 1.750 0.960 0.360 0.420 1.180 0.830
1995; Stephansson et al., 2001; Backers et al., 2002, 2003, 2004; Rao et al., 2003; Shen et al., 2004). The Mode I fracture toughness KIC and Mode II fracture toughness KIIC for some typical rocks are given in Table 1. Most of the existing numerical tools are based on the continuum assumption and rock failure is predicted by means of plastic deformation. Some recent numerical codes simulate the e?ect of existing fractures explicitly, and some of them are designed to model the fracture initiation and propagation of individual cracks. The recent fracture codes like FRACOD have been used (based on some useful mixed fracture criteria like F-criterion) to model the fracture propagation mechanism in brittle materials like rock (Shen and Stephansson, 1994). In the present work a general higher order displacement discontinuity method implementing two crack tip elements for each crack end is used and based on the linear elastic fracture mechanics principles the mixed mode r-criterion (Erdogan and Sih, 1963) is implemented in this numerical model to handle the fracture initiation and propagation mechanism in rock type material considering the ?nite, in?nite, and semi-in?nite bodies. The emphasises is made on the uses of two special crack tip elements which increase the accuracy of the computed Mode I and Mode II stress intensity factors near the crack ends. Two equal crack tip elements have been used for each crack tip and the results seem to be more accurate than using only one crack tip element (which has been used in most of the existing numerical codes). The formulation given in Appendix A are in a concised form for both in?nite and semi-in?nite plane problems. The displacement discontinuity solution based on one and two special crack tip elements are thoroughly explained and the required formulations are derived and given in the text and the crack propagation analysis is accomplished by using the r-criterion. This criterion compares the computed Mode I and Mode II stress intensity factors KI and KII with their corresponding material properties i.e. Mode I and Mode II fracture toughnesses KIC and KIIC (Guo et al., 1992; Scavia, 1992). It should be noted that any fracture criterion can be implemented in the proposed method to study the fracture behaviour of brittle materials. A 45° circular arc crack under biaxial tension has been solved to show the validity of the results obtained by using the higher order displacement discontinuity program TDQCR2 in which uses two special crack tip elements at each crack end. Comparing the numerical and analytical values (i.e. the values of strain energy release rate, G) calculated for this problem proves the validity and accuracy of the numerical results for curved crack problems too.
2. Higher order displacement discontinuity method A displacement discontinuity element of length 2a along the x-axis is shown in Fig. 1(a), which is characterized by a general displacement discontinuity distribution u(e). By taking the ux and uy components of the general displacement discontinuity u(e) to be constant and equal to Dx and Dy respectively, in the
1672
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
y
? u (ε )
x
y
Dy Dx 2a b
x
-a
a < ε < +a
a
+a
Fig. 1. (a) Displacement discontinuity element and the distribution of u(e); (b) constant element displacement discontinuity.
interval (?a + a) as shown in Fig. 1(b), two displacement discontinuity element surfaces can be distinguished, one on the positive side of y(y = 0+) and another one on the negative side (y = 0?). The displacement undergoes a constant change in value when passing from one side of the displacement discontinuity element to the other side. Therefore the constant element displacement discontinuities Dx and Dy can be written as Dx ? ux ?x; 0? ? ? ux ?x; 0? ?; Dy ? uy ?x; 0? ? ? uy ?x; 0? ? ?1? The positive sign convention of Dx and Dy is shown in Fig. 1(b) and demonstrates that when the two surfaces of the displacement discontinuity overlap Dy is positive, which leads to a physically impossible situation. This conceptual di?culty is overcome by considering that the element has a ?nite thickness, in its undeformed state which is small compared to its length, but bigger than Dy (Crouch, 1976; Crouch and Star?eld, 1983). 2.1. Quadratic element formulation The quadratic element displacement discontinuity is based on analytical integration of quadratic collocation shape functions over collinear, straight-line displacement discontinuity elements (Shou and Crouch, 1995). Fig. 2 shows the quadratic displacement discontinuity distribution, which can be written in a general form as Di ?e? ? N 1 ?e?D1 ? N 2 ?e?D2 ? N 3 ?e?D3 ; i i i where, D1 ; i D2 , i and D3 i N 1 ?e? ? e?e ? 2a1 ?=8a2 ; 1 i ? x; y N 3 ?e? ? e?e ? 2a1 ?=8a2 1 ?2?
are the quadratic nodal displacement discontinuities, and N 2 ?e? ? ??e2 ? 4a2 ?=4a2 ; 1 1 ?3?
are the quadratic collocation shape functions using a1 = a2 = a3. A quadratic element has three nodes, which are at the centers of its three sub-elements (see Fig. 2). The displacements and stresses for a line crack in an in?nite body along the x-axis, in terms of single harmonic functions g(x, y) and f(x, y), are given by Crouch and Star?eld (1983) as
y
Element
Di1 1 2a1
Di2 2 2a2
Di3 3 2a3
ε
Fig. 2. Quadratic collocations for the higher order displacement discontinuity elements.
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1673
ux ? ?2?1 ? m?f;y ? yf ;xx ? ? ???1 ? 2m?g;x ? yg.xy ? uy ? ??1 ? 2m?f;x ? yf ;xy ? ? ?2?1 ? m?g;y ? yg.yy ? and the stresses are rxx ? 2l?2f ;xy ? yf ;xyy ? ? 2l?g;yy ? yg;yyy ? ryy ? 2l??yf ;xyy ? ? 2l?g;yy ? yg;yyy ? rxy ? 2l?2f ;yy ? yf ;yyy ? ? 2l??yg;xyy ?
?4?
?5?
l is shear modulus and f,x, g,x, f,y, g,y, etc. are the partial derivatives of the single harmonic functions f(x, y) and g(x, y) with respect to x and y. These potential functions (for a quadratic variation of displacement discontinuity along the element) can be ?nd from f ?x; y? ?
3 X ?1 Dj F j ?I 0 ; I 1 ; I 2 ?; 4p?1 ? m? j?1 x
g?x; y? ?
3 X ?1 Dj F j ?I 0 ; I 1 ; I 2 ? 4p?1 ? m? j?1 y
?6?
the common function Fj, is de?ned as Z 1 F j ?I 0 ; I 1 ; I 2 ? ? N j ?e? ln??x ? e? ? y 2 ?2 de; the integrals I0, and I1 and I2 are expressed as I 0 ?x; y? ? I 1 ?x; y? ? I 2 ?x; y? ? Z
a 2
j ? 1–3
?7?
ln??x ? e? ? y 2 ?2 de ? y?h1 ? h2 ? ? ?x ? a? ln?r1 ? ? ?x ? a? ln?r2 ? ? 2a e ln??x ? e? ? y 2 ?2 de ? xy?h1 ? h2 ? ? 0.5?y 2 ? x2 ? a2 ? ln
2
1
1
?8:a? ?8:b?
Z?a a Z?a a
r1 ? ax r2
1 y 1 e2 ln??x ? e?2 ? y 2 ?2 de ? ?3x2 ? y 2 ??h1 ? h2 ? ? ?3xy 2 ? x3 ? a3 ? ln?r1 ? 3 3 ?a 2 1 2a 2 a x ? y2 ? ? ?3xy 2 ? x3 ? a3 ? ln?r2 ? ? 3 3 3
?8:c?
where, the terms h1, h2, r1 and r2 are de?ned as y y h1 ? arctan ; h2 ? arctan ; x?a x?a r1 ? ??x ? a? ? y 2 ?2
2
1
?9?
1
and r2 ? ??x ? a? ? y 2 ?2
2
2.2. Higher order displacement discontinuity in a half-plane Half-plane problems in solid mechanics can also be solved by in?nite boundary element methods explained before, however, a more accurate and economic way for solving semi-in?nite problems with a traction free surface, using the method of images, is originally introduced by Crouch and Star?eld (1983) for the constant element displacement discontinuity method. They used the analytical solution to a constant element displacement discontinuity, over the line segment jxj 6 a, y = 0 in the semi-in?nite region y 6 0 as shown in Fig. 3. The displacements and stresses due to the actual displacement discontinuity are denoted
1674
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
y
cx
D y = Dy
y
cy
Dx = ? Dx ?β x
t x = σ yx = 0
x
x
cy
y
Dy
cx
Dx
β
Fig. 3. Actual and image displacement discontinuities in half-plane y 6 0, with a traction-free surface (Crouch and Star?eld, 1983).
by uA and rA , those due to its image by uI and rI , and those resulting from the supplementary solution by i ij i ij uS and rS . The complete solution for the semi-in?nite plane y 6 0 can be written as i ij u i ? uA ? uI ? uS i i i and rij ? rA ? rI ? rS ij ij ij ?10?
Based on these formulas and using quadratic element formulations explained in the previous section, a two dimensional semi-in?nite displacement discontinuity computer program can be developed for the analysis of rock fracture mechanics problems, in which uses two special crack tip elements explained in the next section of this paper.
3. Higher order crack tip element formulation and stress intensity factor computation The displacement discontinuity method permits the crack surfaces to be discretized and computes the crack opening displacement (normal displacement discontinuity), and crack sliding displacement (shear displacement discontinuity) directly as a part of the solution for each element. Due to the singularity variap p tions 1/ r, and r for the stresses and displacements near the crack ends, the accuracy of the displacement discontinuity method at the vicinity of the crack tip decreases, and usually a special treatment of the crack at the tip is necessary to increase the accuracy and make the method more e?cient. In this study the hybrid elements are implemented in a general higher order displacement discontinuity method (i.e. the quadratic displacement discontinuity elements and two special crack tip elements for each crack end). Using a special crack tip element of length 2a, as shown in Fig. 4, the displacement discontinuity variations along this element are given as e 1 e 1 2 2 Dy ?e? ? Dy ?a? and Dx ?e? ? Dx ?a? ?11? a a where e is the distance from the crack tip and Dy(a) and Dx(a) are the opening and sliding displacement discontinuities at the center of the special crack tip element.
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1675
y
v A Crack B
r θ u x
D y (ε )
Dx (ε ) a
ε
Fig. 4. Displacement correlation technique for the special crack tip element.
The potential functions fC(x, y) and gC(x, y) for the crack tip element can be expressed as Z a 1 ?1 Dx ?a? 1 e2 ln??x ? e?2 ? y 2 ?2 de 4p?1 ? m? ?a a1 2 Z a 1 ?1 Dy ?a? 1 2 gC ?x; y? ? e2 ln??x ? e? ? y 2 ?2 de 4p?1 ? m? ?a a1 2 fC ?x; y? ? These equations have a common integral of the following form: Z a 1 1 IC ? e2 ln??x ? e?2 ? y 2 ?2 de
?a
?12?
?13?
For the higher order crack tip elements the following series can be used (Crouch and Star?eld, 1983): Di ?e? ? C 1 e2 ? C 2 e2 ? C 3 e2 ? ? ? ?
1 3 5
?14?
In order to use two crack tip elements the ?rst two terms of Eq. (14) should be considered, which can be arranged in the following form: Di ?e? ? ?N C1 ?e??D1 ?a? ? ?N C2 ?e??D2 ?a? i i ?15?
The crack tip element has a length a = a2 + a1, and the shape functions NC1(e) and NC2(e) can be obtained as 1 3 1 3 a2 e 2 ? e 2 a1 e 2 ? e 2 N C1 ?e? ? 1 and N C2 ?e? ? ? 1 ?16? a2 ?a2 ? a1 ? a2 ?a2 ? a1 ? 1 2 The potential functions fC(x, y) and gC(x, y) for the crack tip can be expressed as Z a 1 ?1 fC ?x; y? ? Dx ?e? ln??x ? e?2 ? y 2 ?2 de 4p?1 ? m? ?a Z a 1 ?1 1 2 Dy ?e?e2 ln??x ? e? ? y 2 ?2 de gC ?x; y? ? 4p?1 ? m? ?a
?17:a? ?17:b?
1676
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
inserting Eq. (16) in Eqs. (15) and (17), gives Z a Z a 1 1 ?1 2 2 fC ?x; y? ? N 1 ?e? ln??x ? e? ? y 2 ?2 de D1 ? N 2 ?e? ln??x ? e? ? y 2 ?2 de D2 ?18:a? x x 4p?1 ? m? ?a ?a 1 ?18:b? fC ?x; y? ? ?I C1 ?x; y?D1 ? I C2 ?x; y?D2 ? ? x x 4p?1 ? v? 2 X 1 fC ?x; y? ? ? Dj F Cj ?I C1 ; I C2 ? ?18:c? 4p?1 ? v? j?1 x and similarly gC ?x; y? ? ? in which F Cj ?I C1 ; I C2 ? ? Z
a ?a 2 X 1 Dj F Cj ?I C1 ; I C2 ? 4p?1 ? v? j?1 y
?19?
N Cj ?e? ln??x ? e?2 ? y 2 ?2 de;
1
j ? 1; 2
?20?
From Eq. (20) the following integrals are obtained: I C1 ?x; y? ? Z
a
a2 e 2 ? e 2 a1 ?a2 ? a1 ?
a
1 2
1
3
ln??x ? e? ? y 2 ?2 de;
3
2
1
I C2 ?x; y? ?21?
?a
??
Z
a1 e 2 ? e 2 a2 ?a2 ? a1 ?
1 2
1
ln??x ? e? ? y 2 ?2 de
2
1
?a
These integrals can be arranged as I C1 ?x; y? ?
1
a2 1 1 2 IC ? IC d C1 d C1
1
and
I C2 ?x; y? ?
1
?a1 1 1 2 IC ? IC d C2 d C2
?a
?22?
1
where d C1 ? a2 ?a2 ? a1 ?; d C2 ? a2 ?a2 ? a1 ?; I C ? 1 2 IC ?
2
Ra
e2 ln??x ? e?2 ? y 2 ?2 de, and
1
Z
a
e2 ln??x ? e?2 ? y 2 ?2 de
1 2
3
1
?23?
?a
The derivations of integrals I C and I C which are used in the solution of the displacement discontinuities near the crack tip for both in?nite and semi-in?nite plane problems are given in Appendix A. Several mixed mode fracture criteria has been used in literature to investigate the crack initiation direction and its path (Whittaker et al., 1992; Shen and Stephansson, 1994; Rao et al., 2003). The maximum tangential stress criterion (or r-criterion), is the used here to investigate the crack initiation and propagation direction. This is a widely used mixed mode fracture criterion which is well ?tted with the experimental results (Ingra?ea, 1983; Guo et al., 1990; Scavia, 1992; Whittaker et al., 1992). Based on LEFM theory, the Mode I and Mode II stress intensity factors KI and KII can be written in terms of the normal and shear displacement discontinuities as (Shou and Crouch, 1995): 1 1 l 2p 2 l 2p 2 Dy ?a? and K II ? Dx ?a? ?24? KI ? 4?1 ? m? a 4?1 ? m? a
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1677
4. Veri?cation of higher order displacement discontinuity Veri?cation of this method is made through the solution of several example problems i.e. a pressurized crack in an in?nite body, a circular arc crack under biaxial tension and circular holes with emanating cracks in in?nite and semi-in?nite bodies. These simple examples are used here because they have simple analytical solutions and have also been solved numerically by other researches, so that the computed numerical results in this paper can be compared and the validity of the proposed computer programs can be proved. 4.1. A pressurized crack in an in?nite plane Because of its simple solution, the problem of a pressurized crack have been used for the veri?cation of the numerical methods developed here. The analytical solution of this problem has been derived and explained by Sneddon (1951). Based on Fig. 5, the analytical solution for the normal displacement discontinuity Dy along the crack boundary, and the normal stress ry near the crack tip (jxj > b), can be written as Dy ? ?
1 2?1 ? m?P 2 ?b ? x2 ?2 ; jxj > b l
and
ry ?
Px ?x2 ? b2 ?2
1
? P ; jxj < b
?25?
m is the Poisson?s ratio of the body. Consider a pressurized crack of a half length b = 1 meter (m), under a normal pressure P = ?10 MPa, with a modulus of elasticity E = 2.2 GPa, and a Poisson?s ratio m = 0.1 (Fig. 5). The normalized displacement discontinuity distribution Dy/b · 103 along the surface of the pressurized crack are given in Table 2 using the constant (ordinary) displacement discontinuity program TWODD with di?erent crack tip elements (Crouch and Star?eld, 1983). As shown in this table by using only one special crack tip element, the percent error of displacement discontinuity Dy at a distance x = 0.05b from the crack tip, reduces from 26.12 to 8.59, and by using two special crack tip elements it reduces to 5.07. The same problem is solved by the higher order displacement discontinuity method (i.e. TWODQ program) using only 20 quadratic elements (60 nodes). The percentage error of displacement discontinuity Dy at a distance x = 0.05b from the crack tip will be about 5.82% (with out using any special crack tip element) in which the results obtained by TWODD program with the same number of nodes (60 constant elements) gives an error of about 24.17%. Therefore using this higher order program for calculating the normalized
y
P x P 2b
Fig. 5. A pressurized crack in an in?nite body.
1678
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
Table 2 Displacement discontinuity Dy(a) at the center x = a of an element at the crack tip using ordinary elements, one special crack tip element and two special crack tip elements Number of elements Distance from crack tip (Dy(a)/b) · 103 Analytical solution ?1.1906 ?0.7846 ?0.5621 ?0.4000 Ordinary elements ?1.5463 ?0.9964 ?0.7089 ?0.5029 One special crack tip element ?1.3275 ?0.8569 ?0.6104 ?0.4332 Two special crack tip elements ?1.2150 ?0.8202 ?0.5906 ?0.4212
4 10 20 40
0.25 0.1 0.05 0.025
displacement discontinuity distribution Dy/b · 103 along the surface of the pressurized crack, gives more accurate results (see Table 3). The normalized normal stress (ry/P) near the crack tip and along the x-axis of the pressurized crack is presented in Table 4. The over all results show that the program using quadratic elements (i.e. TWODQ) gives more accurate results compared to the program using constant elements (i.e. TWODD). 4.2. Curved cracks Curved and kink cracks may occur in cracked bodies (Shou and Crouch, 1995). The proposed method is applied to the problem of a 45° circular arc crack under far ?eld biaxial tension (Fig. 6). The program TDQCR2 (using quadratic displacement discontinuity elements with two special crack tip elements at each crack end) and the program TDQCR1 (using quadratic displacement discontinuity elements with only one special crack tip element at each crack end) have been developed for the analysis of the crack problems. Analytical values of the Mode I and Mode II stress intensity factors, KI and KII, and strain energy release rate, G for a general circular arc crack under biaxial tension given by Cotterell and Rice (1980) as " #1 " #1 2 2 a pr sin a a pr sin a 1 ? m2 2 2 2 ?K I ? K 2 ? ; K II ? r sin and G ? ?26? K I ? r cos II 4 1 ? sin2 a 4 1 ? sin2 a E 4 4
Table 3 Comparison of variation of the displacement discontinuity Dy/b · 103 along the surface of a pressurized crack using ordinary (constant) elements and higher order elements displacement discontinuity methods without using any special crack tip element x/b Distance from crack tip 0.0 0.086 0.171 0.257 0.343 0.429 0.514 0.600 0.686 0.771 0.857 0.950 Dy/b · 1000 Const. Elems. (TWODD) 4.046 4.032 3.989 3.916 3.812 3.673 3.497 3.276 3.000 2.654 2.203 1.536 Analytical results 3.980 3.946 3.901 3.827 3.720 3.578 3.396 3.168 2.882 2.520 2.040 1.237 Quad. Elems. (TWODQ) 3.980 3.966 3.922 3.848 3.741 3.599 3.419 3.193 2.909 2.556 2.094 1.309
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692 Table 4 Variation of the normalized normal stress, ry/P near the crack tip and along the x-axis of the pressurized crack problem (x ? b)/b ry/P TWODD 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 0.36 0.40 3.022 1.783 1.294 1.021 0.844 0.718 0.623 0.549 0.490 0.441 ANALYTIC 2.641 1.648 1.221 0.973 0.809 0.691 0.602 0.532 0.475 0.429
1679
Higher order elements 2.845 1.719 1.259 0.998 0.826 0.704 0.612 0.541 0.483 0.435
σ yy = σ xx
E = 10.0GPa
ν = 0.2
r = 1.0m
y
α
2
σ = 10MPa
α
2
r
x
o
σ xx = σ
Fig. 6. Circular arc crack under uniform biaxial tension (Shou and Crouch, 1995).
Considering a circular arc crack (a = 45°) under biaxial tension, r = 10 MPa with a radius of r = 1 m. The modulus of elasticity and Poisson?s ratio of the material are taken as E = 10 GPa and m = 0.2. The analytical values of the problem are obtained from Eq. (26) as G = 11.47 · 10?3 and based on the r-criterion, the crack propagation angle, h0 = 20.90° (degrees) (Whittaker et al., 1992). The numerical solution of this problem have been obtained by using two higher order displacement discontinuity programs TDQCR1 (using one special crack tip element at each crack end), and TDQCR2 (using two special crack tip elements at each crack end). The numerical values for G and h0 are given in Tables 5 and 6; and Figs. 7 and 8, considering the two cases: (i) using a small crack tip element length L equal to 0.5° (i.e L/b = 0.011) and different number of nodes along the crack and (ii) using 60 nodes along the crack and di?erent crack tip element length ratios (i.e. di?erent L/b ratios). Comparing the numerical results of these two cases with the analytical results show that the results obtained by the program TDQCR2 are somewhat superior to those obtained by the program TDQCR1 specially when using relatively smaller number of nodes along the crack. The e?ect of L/b ratio on the results given by the program TDQCR2 is also negligible. Although the results obtained by both programs are in good agreement with the analytical results but in most cases the results obtained by the program TDQCR2 are preferred and can be used for the analysis of crack problems in rock fracture mechanics.
1680
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
Table 5 The numerical values of the strain release rate, G for a circular arc crack (a = 45°), using di?erent number of nodes along the crack with a 0.5° crack tip Number of nodes The strain energy release rate G · 10?3 TDQCR2 6 12 18 24 30 36 42 48 54 60 18.92 13.99 12.55 11.94 11.61 11.43 11.32 11.27 11.21 11.20 TDQCR1 23.82 16.44 14.25 13.24 12.68 12.34 12.11 11.94 11.82 11.73 Crack initiation angle h0 (degrees) TDQCR2 23.13 21.57 21.17 21.01 20.93 20.90 20.88 20.86 20.85 20.85 TDQCR1 22.58 21.57 21.26 21.12 21.05 20.99 20.96 20.93 20.92 20.91
Table 6 The numerical values of the strain release rate, G for a circular arc crack (a = 45°), using di?erent L/b ratios and 60 nodes along the crack L/b ratio The strain energy release rate G · 10?3 TDQCR2 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 13.47 11.75 11.34 11.20 11.15 11.13 11.13 11.14 TDQCR1 15.53 12.92 12.18 11.85 11.68 11.58 11.51 11.47
The strain energy release rate G for a 45 degree circular arc crack 25 Strain energy release rate (G/1000)
TDQCR2
20
TDQCR1 ANALYTIC
15
10
5
0 6 12 18 24 30 36 42 48 54 60 Number of nodes along the crack
Fig. 7. The strain energy release rate, G for a circular arc crack (a = 45°).
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1681
20
The strain energy release rate G for a 45 degrees circular arc crack using different L/b ratios
TDQCR2 TDQCR1 ANALYTIC
strain energy release rate (G/1000)
18 16 14 12 10 8 6 4 0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
L/b ratio
Fig. 8. The strain energy release rate, G for a circular arc crack (a = 45°).
4.3. Circular holes with two emanating cracks In order to show the bene?t of both higher order elements and special crack tip elements explained above the example problem shown in Fig. 9 is solved numerically by the higher order displacement discontinuity method using quadratic displacement discontinuity elements for crack analysis (i.e. the programs TDQCR1 and TDQCR2). The following assumptions are made to solve this problem numerically: the far ?eld stress r = 10 MPa, the hole radius R = 1 m, modulus of elasticity E = 10 GPa, Poisson?s ratio t = 0.2, and Mode I fracture toughness KIC = 2 MPa m1/2 (for a typical hard rock under plane strain condition). p?????? The analytical value of the normalized stress intensity factor K I =?r pb? obtained from the solutions given by Sihp?????? (1973) is about 1.96, for b/R = 0.4. The numerical result of the normalized stress intensity factor K I =?r pb?, for b/R = 0.4 using 90 nodes along the hole and 49 nodes along each crack, for di?erent L/b ratios (i.e. the ratio of crack tip element length L to the crack length b) are shown graphically in Fig. 10.
y b x R b
σx =σ
Fig. 9. A circular hole with two emanating cracks of length b under far ?eld uniform tension.
1682
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
The stress intensity factor for the problem of the circular hole with two cracks, for b/R= 0.4
TDQCR2
2.1 Normalized stress intensity factors
2.05
TDQCR1 ANALYTIC
2
1.95
1.9
1.85 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 l/b ratio
p????? ? Fig. 10. The normalized stress intensity factor K I =?r pb?, for b/R = 0.4 and di?erent l/b ratios.
The numerical results show that with L/b ratios between 0.1 and 0.2, the numerical results are accurate (in most cases the error is less than about 0.5%). 4.4. Crack problems in semi-in?nite plane The problem shown in Fig. 11 has been solved with the semi-in?nite displacement discontinuity method using quadratic elements and special crack tip elements at each crack tip , considering the two limiting cases of pressurized circular holes having four cracks with (i) no pressure penetration (empty cracks), and (ii) full pressure penetration (fully pressurized cracks). The results of the normalized Mode I and Mode II stress
y
x
b
C
p R b
Fig. 11. A pressurized circular hole with four radial cracks at depth C in a semi-in?nite body.
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1683
p????????? p????????? intensity factors K I =?p pqR?, and K II =?p pqR?, and the crack propagation angle h0, for these two extreme cases (i.e. empty cracks and fully pressurized cracks) against di?erent normalized depths from the free surface of the half plane (C/R ratio), have been calculated numerically and given in Tables 7 and 8, respectively. All the results are numerically calculated by using 20 quadratic elements along the hole and 10 quadratic elements along each crack, with a constant value of q ? b?R ? 2.5, R = 1 m, and p = rp????????? MPa. = 10 R The analytical values of the normalized Mode I and Mode II stress intensity factors K I =?p pqR? and p????????? K II =?p pqR? and the crack initiation angle h0, for the problem of a pressurized circular hole????????? four symp with metricp????????? cracks in a an in?nite plane are given by Ouchterlony (1983) as: K I =?p pqR? ?? 0.1966, radial p????????? p???????? K II =?p pqR? ? 0.0, and h0 = 0.0 for the empty cracks; and K I =?p pqR? ? 0.9085, K II =?p pqR? ? 0.0, and h0 = 0.0 for the fully pressurized cracks, respectively.
Table 7 p????????? p????????? The normalized stress intensity factors K I =?p pqR?, and K II =?p pqR?, and the crack propagation angle h0, for a pressurized hole with four empty radial cracks at di?erent depths (C/R ratios) p????????? p????????? C/R ratios K I =?p pqR? K II =?p pqR? h0 (degrees) Up. crack 2.0 2.25 2.5 2.75 3.0 3.25 3.5 3.75 4.0 4.25 4.5 4.75 5.0 0.8473 0.7068 0.6053 0.5196 0.4527 0.3959 0.3592 0.3245 0.3086 0.2808 0.2686 0.2565 0.2477 Lo. crack 0.4247 0.3668 0.3354 0.3141 0.2977 0.2792 0.2703 0.2628 0.2570 0.2510 0.2478 0.2441 0.2410 Up. crack 0.1056 0.0644 0.0633 0.0593 0.0533 0.0462 0.0418 0.0326 0.0267 0.0227 0.0176 0.0128 0.0106 Lo. crack 0.0008 ?0.0196 ?0.0250 ?0.0270 ?0.0290 ?0.0256 ?0.0254 ?0.0250 ?0.0240 ?0.0230 ?0.0208 ?0.0199 ?0.0187 Up. crack ?13.8 ?10.2 ?11.7 ?12.7 ?13.1 ?13.0 ?12.9 ?11.3 ?10.4 ?9.1 ?7.4 ?5.7 ?4.9 Lo. crack ?0.2 6.1 8.4 10.0 10.1 10.3 10.6 10.7 10.5 10.3 9.5 9.2 8.8
Table 8 p????????? p????????? The normalized stress intensity factors K I =?p pqR?, and K II =?p pqR? and the crack propagation angle h0, for a pressurized hole under the uniform inside pressure p, with four fully pressurized radial cracks at di?erent depths (C/R ratios) p????????? p????????? C/R ratios K I =?p pqR? K II =?p pqR? h0 (degrees) Up. crack 2.0 2.25 2.5 2.75 3.0 3.25 3.5 3.75 4.0 4.25 4.5 4.75 5.0 2.2036 1.8940 1.6993 1.5298 1.3937 1.2690 1.1927 1.1116 1.0344 1.0114 0.9863 0.9582 0.9390 Lo. crack 1.3534 1.2366 1.1766 1.1352 1.1022 1.0578 1.0392 1.0230 1.0085 0.9960 0.9903 0.9816 0.9745 Up. crack 0.1954 0.1230 0.1285 0.1229 0.1109 0.0949 0.0862 0.0609 0.0470 0.0362 0.0223 0.0091 0.0044 Lo. crack 0.0194 ?0.0307 ?0.0465 ?0.0565 ?0.0626 ?0.0552 ?0.0573 ?0.0587 ?0.0580 ?0.0573 ?0.0523 ?0.0510 ?0.0487 Up. crack ?10.0 ?7.4 ?8.6 ?9.1 ?9.0 ?8.5 ?8.2 ?6.2 ?5.5 ?4.1 ?2.5 ?1.1 ?0.5 Lo. crack ?1.6 2.8 4.5 5.7 6.5 5.9 6.3 6.5 6.5 6.5 6.0 5.9 5.7
1684
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
Figs. 12 and 13 are based on the normalized stress intensity factors given in the Tables 7 and 8, which compare the di?erent results obtained for the upper cracks (the cracks near to the free surface of the half plane), and the lower cracks, with the available analytical results of the circular pressurized hole in an in?nite plane (Ouchterlony, 1983). The numerical results show that as the depth of the circular pressurized hole (C/R ratio) increases the mixed mode stress intensity factors KI and KII, and the crack propagation angle h0 tend to their corresponding analytical values of the circular pressurized hole in an in?nite plane.
Effect of depth on the mixed mode stress intensity factors, for the 4 empty radial crack problem 0.9 Normalized Stress intensity factors 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 2 2.25 2.5 2.75 3 3.25 3.5 3.75 C/R ratio 4 4.25 4.5 4.75 5 Mode II Mode I
Up. Crack Lo. Crack Up. Crack Lo. Crack Analytic
p????????? p????????? Fig. 12. The normalized stress intensity factors K I =?p pqR? and K II =?p pqR?, for di?erent C/R ratios, for a pressurized hole with four empty radial cracks in a semi-in?nite rock mass.
Effect of depth on the mixed mode stress intensity factors, for the fully pressurized 4 radial cracks problem 2.5 Normalized Stress intensity factors 2 1.5 1 0.5 Mode II 0
2 2.25 2.5 2.75 3 3.25 3.5 3.75 4 4.25 4.5 4.75 5 Up. Crack Lo. Crack Up. Crack Lo. Crack Analytic
Mode I
-0.5 C/R ratio
p????????? p????????? Fig. 13. The normalized stress intensity factors K I =?p pqR? and K II =?p pqR?, for fully pressurized radial cracks emanating from a pressurized hole at di?erent depths, in a semi-in?nite rock mass.
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1685
5. Conclusion The e?ect of using one and two special crack tip elements and also the e?ect of using higher order displacement discontinuity in two dimensional in?nite and semi-in?nite crack problems have been investigated. The complete and proper solutions and formulations are explained and given in the text and Appendix A. This method is modi?ed to include, the ?nite, in?nite, and semi-in?nite problems using higher order (quadratic) displacement discontinuity elements. Several simple and mostly used example problems are selected and explained to verify the proposed modi?cations and improvements in the higher order displacement discontinuity method. At the end a somewhat more complicated and useful semi-in?nite crack problem is solved and the numerical results are compared with the analytical values of the same problem in an in?nite plane. It has been shown that as the depth of pressurized hole increases the problem changes to that of the in?nite plane case. Therefore, it is concluded that the numerical results obtained by using the modi?ed semi-in?nite program gives very accurate results. This method is also modi?ed for crack initiation and propagation analysis. Although any mixed mode fracture criterion can be implemented to this numerical code, but in this work, based on the linear elastic fracture mechanics (LEFM) principles, the maximum tangential stress criterion or r-criterion is employed to investigate the crack initiation and propagation direction.
Appendix A. The integrals and their derivatives used for one and two special crack tip elements for both in?nite and semi-in?nite plane problems A.1. The integral and its derivatives for one special crack tip element Z Z Z
a
1
Ic ?
e2 ln??x ? e? ? y 2 ?2 de
a
1
2
?a
I c;x ?
e2 ?x ? e? ??x ? e? ? y 2 ? e2 y ??x ? e?2 ? y 2 ?
1
1
?a a
2
de ? xA1 ? A2
I c;y ?
de ? yA1
?a
I c;xy ? yA1;x I c;xx ? ?A1 ? yA1;y ? ?I c;yy I c;xyy ? A1;x ? yA1;xy I c;yyy ? 2A1;y ? yA1;yy The following formulae (derivatives) are also required for treatment of the traction free half-plane problems: I c;xyyy ? 2A1;xy ? yA1;xyy I c;yyyy ? 3A1;yy ? yA1;yyy
1686
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
where, A1, A2, and the derivatives of A1, are de?ned as Z a 1 e2 A1 ? de 2 2 ?a ??x ? e? ? y ? p????? 3 2 x 2a ? 2 2aq cos u ? q2 x 6 0.5 cos u ? y sin u ln 2a ? 2p?????q cos u ? q2 ? sin u ? y cos u 7 2a 6 7 ! p????? ? q?1 6 7 4 5 2 2aq sin u ? arctan 2 ? 2a q Z a 3 2 e A2 ? de 2 ??x ? e? ? y 2 ? ?a p????? 3 2 x 2a ? 2 2aq cos u ? q2 x p????? ? sin u ? sin u ln cos u 7 0.5 cos u ? 6 y y 2a ? 2 ! q cos u ? q2 2a 6 7 p????? ? q6 7 4 5 2 2aq sin u ? arctan q2 ? 2a where q ? ?x2 ? y 2 ?4
1
and u ? 0.5 arctan
y x
and the derivatives of A1 are x A1.x ? q?1 A1x1 ? 4 A1 2q y A1.y ? q?1 A1y1 ? 4 A1 2q y xy x A1.xy ? q?1 A1x2 ? 5 A1x1 ? 8 A1 ? 4 A1;y 2q q 2q y x2 ? y 2 y A1.yy ? q?1 A1y2 ? 5 A1y1 ? A1 ? 4 A1;y 2q 2q 2q8 y 2x2 ? 3y 2 x2 ? 3y 2 2xy x A1.xyy ? q?1 A1x3 ? 5 A1x2 ? A1x1 ? A1 ? 8 A1;y ? 4 A1;yy q q 2q 2q9 q12 y 2x2 ? 3y 2 y?3x2 ? y 2 ? x2 ? y 2 y A1.yyy ? q?1 A1y3 ? 5 A1y2 ? A1y1 ? A1 ? A1;y ? 4 A1;yy 9 12 8 q q 2q 2q q where A1x1 ? CN 1 ? FLX ? FL ? CN 1X ? CN 2 ? TX ? TC ? CN 2X A1y1 ? CN 1 ? FLY ? FL ? CN 1Y ? CN 2 ? TY ? TC ? CN 2Y A1x2 ? CN 1 ? FLXY ? FLX ? CN 1Y ? FL ? CN 1XY ? FLY ? CN 1X ? CN 2 ? TXY ? TX ? CN 2Y ? TY ? CN 2X ? TC ? CN 2XY A1y2 ? CN 1 ? FLYY ? 2 ? FLY ? CN 1Y ? FL ? CN 1YY ? 2 ? TY ? CN 2Y ? CN 2 ? TYY ? TC ? CN 2YY A1x3 ? CN 1 ? FLXYY ? 2 ? FLXY ? CN 1Y ? FLX ? CN 1YY ? 2 ? FLY ? CN 1XY ? FL ? CN 1XYY ? FLYY ? CN 1X ? 2 ? CN 2Y ? TXY ? CN 2 ? TXYY ? TX ? CN 2YY ? 2 ? CN 2XY ? TY ? TYY ? CN 2X ? TC ? CN 2XYY A1y3 ? 3 ? CN 1Y ? FLYY ? CN 1 ? FLYYY ? 3 ? FLY ? CN 1YY ? FL ? CN 1YYY ? 3 ? CN 2Y ? TYY ? 3 ? TY ? CN 2YY ? CN 2 ? TYYY ? TC ? CN 2YYY
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1687
The three form of terms CNs, FL?s, and T ?s are necessary to be de?ned, for the above equations, and are given, respectively, in the following formulae. (i) The terms CN1, CN2, CN1X, CN2X, etc. are de?ned as x CN 1 ? 0.5 cos u ? sin u y x CN 2 ? sin u ? cos u y CN 1X ? ?CN 1?;X ? 0.5?CP 2 ? ?SNP ? x ? SP 2?=y? CN 1Y ? ?CN 1?;y ? 0.5?CP 3 ? x ? SP 3=y ? x ? SNP =y 2 ? CN 1XY ? ?CN 1?;xy ? 0.5?CP 4 ? ?y ? SP 3 ? SNP ? x ? SP 2?=y 2 ? x ? SP 4=y?; etc:; and CN 2X ? ?CN 2?;X ? SP 2 ? ?CSP ? x ? CP 2?=y CN 2Y ? ?CN 1?;y ? SP 3 ? x ? ?CSP ? y ? CP 3?=y 2 CN 2XY ? ?CN 2?;xy ? SP 4 ? ?CSP ? CP 2?=y 2 ? ?CP 3 ? x ? CP 4?=y; etc.; where h y i SNP ? sin u ? sin 0.5 arctan x h y i CSP ? cos u ? cos 0.5 arctan x CP 2 ? ?cos u?;x ? y ? SNP =?2q4 ? CP 3 ? ?cos u?;y ? ?x ? SNP =?2q4 ? CP 4 ? ?cos u?;xy ? ?2?x2 ? y 2 ? ? SNP ? xy ? CSP ?=?4q8 ?; etc.; and SP 2 ? ?sin u?;x ? ?y ? CSP =?2q4 ? SP 3 ? ?sin u?;y ? x ? CSP =?2q4 ? SP 4 ? ?sin u?;xy ? ?2?y 2 ? x2 ? ? CSP ? xy ? SNP ?=?4q8 ?; etc. (ii) The terms FL, FLX, FLY, etc. are de?ned as FL ? ln?DL1=DL2? FLX ? FL1X ? FL2X FLY ? FL1Y ? FL2Y FLXY ? FL1XY ? FL2XY ; etc.; where FL1X ? DL1X =DL1 FL2X ? DL2X =DL2 FL1XY ? ?FL1X ?;y ? ?DL1Y ? DL1X =?DL1?2 ? DL1XY =DL1 FL2XY ? ?FL2X ?;y ? ?DL2Y ? DL2X =?DL2?2 ? DL2XY =DL2; etc.;
1688
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
with p????? DL1 ? 2a ? 2 2aq cos u ? q2 p????? DL1X ? ?DL1?;x ? ? 2a?x ? CSP =q3 ? y ? SNP =q3 ? ? x=q2 p????? DL1Y ? ?DL1?;y ? ? 2a?y ? CSP =q3 ? x ? SNP =q3 ? ? y=q2 p????? DL1XY ? ?DL1?;xy ? ? 2a??2xy ? CSP ? ?x2 ? y 2 ? ? SNP ?=?2q7 ? ? xy=q6 ; etc:; and p????? DL2 ? 2a ? 2 2aq cos u ? q2 p????? DL2X ? ?DL2?;x ? 2a?x ? CSP =q3 ? y ? SNP =q3 ? ? x=q2 p????? DL2Y ? ?DL2?;y ? 2a?y ? CSP =q3 ? x ? SNP =q3 ? ? y=q2 p????? DL2XY ? ?DL2?;xy ? 2a??2xy ? CSP ? ?x2 ? y 2 ? ? SNP ?=?2q7 ? ? xy=q6 ; etc: (iii) The terms TC, TX, TY, TXY, TYY, etc. are de?ned as p????? 2 2aq sin u TC ? arctan ; and q2 ? 2a TX ? DT 1=DT TY ? DT 2=DT TXY ? ?TX ?;y ? ?DTY ? DT 1=?DT ? ? DT 1Y =DT TYY ? ?TY ?;y ? ?DTY ? DT 2=?DT ?2 ? DT 2Y =DT ; etc:; where DT ? q4 ? 4aq2 cos 2u ? 4a2 DTY ? 2y ? 4a?y cos 2u ? x sin 2u?=q2 DTYY ? 2;
p????? DT 1 ? 2 2a? RY 1 ? PY 1 ? RY 2 ? PY 2? p????? DT 1Y ? ?DT 1?;y ? 2 2a??RY 11 ? PY 1 ? RY 1 ? PY 11 ? RY 21 ? PY 2 ? RY 2 ? PY 21? and DT 1YY ? ?DT 1?;yy p????? ? 2 2a??RY 12 ? PY 1 ? 2RY 11 ? PY 11 ? RY 1 ? PY 12 ? RY 22 ? PY 2 ? 2RY 21 ? PY 21 ? RY 2 ? PY 22? p????? DT 2 ? 2 2a? RY 1 ? PY 4 ? RY 2 ? PY 3? p????? DT 2Y ? ?DT 2?;y ? 2 2a?RY 11 ? PY 4 ? RY 1 ? PY 41 ? RY 21 ? PY 3 ? RY 2 ? PY 31? with DT 2YY ? ?DT 2?;yy p????? ? 2 2a?RY 12 ? PY 4 ? 2RY 11 ? PY 41 ? RY 1 ? PY 42 ? RY 22 ? PY 3 ? 2RY 21 ? PY 31 ? RY 2 ? PY 32?
2
PY 1 ? y ? CSP ? x ? SNP PY 2 ? y ? CSP ? x ? SNP PY 3 ? x ? CSP ? y ? SNP PY 4 ? x ? CSP ? y ? SNP
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1689
PY 12 ? CSP ? x ? PY 4=?2q4 ? 4xy ? PY 4 ? x2 ? PY 1 4q8 4 PY 21 ? CSP ? x ? PY 3=?2q ? PY 12 ? ?x ? SNP =q4 ? 4xy ? PY 3 ? x2 ? PY 2 4q8 4 PY 31 ? SNP ? x ? PY 2=?2q ? xy ? PY 2 x ? PY 21 ? PY 32 ? x ? CSP =?2q4 ? ? q8 2q4 4 PY 41 ? ?SNP ? x ? PY 1=?2q ? xy ? PY 1 and PY 42 ? ?x ? ?CSP ? PY 11?=?2q4 ? ? q8 1 ?y ?2q4 ? 5y 2 ; RY 11 ? 5 ; RY 12 ? RY 1 ? 2q 4q 8q9 a ?3ay ?3a?2q4 ? 7y 2 ? RY 2 ? 3 ; RY 21 ? ; RY 22 ? 7 q 2q 4q11 PY 22 ? ?x ? SNP =q4 ? A.2. The integrals and their derivatives for two special crack tip elements Starting with the integrals I C and I C given in Section 3, the ?rst integral I C ? I c and its derivatives are the 2 same as those given in A.1 of this appendix. The integral I C and its derivatives are as follows: Z a 2 1 3 2 e2 ln??x ? e? ? y 2 ?2 de IC ? I C;x ? I C;y ?
2 2 2 2 2 2 1 2 1
Z
?a a
e2 ?x ? e? ??x ? e?2 ? y 2 ? e2 y ??x ? e? ? y 2 ?
2
3
3
de ? x de ? y
Z Z
a
e2 ??x ? e?2 ? y 2 ? e2 ??x ? e? ? y 2 ?
2
3
3
de ?
Z
a ?a
e2 ??x ? e?2 ? y 2 ?
5
de ? xA2 ? A3
Z
?a a
?a a
de ? yA2
?a
?a
I C;xy ? yA2;x I C;yy ? yA2;y ? A2 I C;xyy ? yA2;xy ? A2;x I C;yyy ? yA2;yy ? 2A2;y and for the semi-in?nite case the following derivatives are also needed: I C;xyyy ? yA2;xyy ? 2A2;xy I C;yyyy ? yA2;yyy ? 3A2;yy where, A1 and its derivatives are given in section A of this appendix, A3 and, A2 and its derivatives which needed for in?nite and semi-in?nite plane problems, are given as
2 2
1690
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
?a ??x ? e? ? and A2 and its derivatives can be written as Z a 3 1 1 e2 x2 ? y 2 x2 ? y 2 ?h1 ? h2 ? ? 2xA1 ? 2?2a?2 ? I 0;y ? 2xA1 de ? 2?2a?2 ? A2 ? 2 2 2y 2y ?a ??x ? e? ? y ? x2 ? y 2 x I 0;xy ? I 0;y ? 2xA1;x ? 2A1 A2;x ? y 2y x2 ? y 2 y 2 ? x2 I 0;yy ? A2;y ? I 0;y ? 2xA1;y 2y 2y 2 x2 ? y 2 y 2 ? x2 x x I 0;xyy ? A2;xy ? I 0;xy ? I 0;yy ? 2 I 0;y ? 2xA1;xy ? 2A1;y y y 2y 2y 2 x2 ? y 2 y 2 ? x2 x2 I 0;yyy ? 2 A2;yy ? I 0;yy ? 3 I 0;y ? 2xA1;yy 2y 2y 2 y 2 2 x ?y y 2 ? x2 x2 x x I 0;xyyy ? 2 A2;xyy ? I 0;xyy ? 3 I 0;xy ? I 0;yyy ? 2 2 I 0;yy ? 2xA1;xyy ? 2A1;yy y y 2y 2y 2 y x2 ? y 2 y 2 ? x2 x2 3x2 I 0;yyyy ? 3 A2;yyy ? I 0;yyy ? 3 3 I 0;yy ? 4 I 0;y ? 2xA1;yyy 2y 2y 2 y y and for the semi-in?nite case the following derivatives are also needed: x2 ? y 2 y 2 ? x2 3x2 3x2 x 3x I 0;xyyyy ? 3 I 0;xyyy ? 3 I 0;xyy ? 4 I 0;xy ? I 0;yyyy ? 2 I 0;yyy A2;xyyy ? 2 y y 2y 2y y y 6x 6x ? 3 I 0;yy ? 4 I 0;y ? 2xA1;xyyy ? 2A1;yyy y y x2 ? y 2 y 2 ? x2 6x2 12x2 12x2 I 0;yyyyy ? 4 A2;yyyy ? I 0;yyyy ? 3 I 0;yyy ? 4 I 0;yy ? 5 I 0;y ? 2xA1;yyyy 2y 2y 2 y y y
A3 ?
Z
a
e2
2
5
y2?
3 2 de ? ? ?2a?2 ? 2xA2 ? ?x2 ? y 2 ?A1 3
where I0 ?
Z
a 2 ?a
ln??x ? e? ? y 2 ?2 de ? y?h1 ? h2 ? ? ?x ? a? ln?r1 ? ? ?x ? a? ln?r2 ? ? 2a
1
@I 0 ? ln r1 ? ln r2 I 0;x ? @x I 0;y ? h1 ? h2 y y I 0;xy ? ? 2 ? 2 r r2 1 ?x ? a? ?x ? a? ? I 0;xx ? ? ? ?I 0;yy r2 r2 1 2
! 2 2 ?x ? a? ? y 2 ?x ? a? ? y 2 I 0;xyy ? ? ? ? ?I 0;xxx r4 r4 1 2 ?x ? a? ?x ? a? I 0;yyy ? ?2y ? ? ?I 0;xxy r4 r4 1 2 ?x ? a??r2 ? 4y 2 ? ?x ? a??r2 ? 4y 2 ? 1 1 I 0;xyyy ? ?2 ? r6 r6 1 2 ! 2 2 3?x ? a? ? y 2 3?x ? a? ? y 2 I 0;yyyy ? 2y ? r6 r6 1 2
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
1691
and for the semi-in?nite case the following derivatives are also needed: ! ?x ? a???x ? a?2 ? y 2 ? ?x ? a???x ? a?2 ? y 2 ? I 0;xyyyy ? 24y ? r8 r8 1 2 ! 2 2 2 I 0;yyyy 5?x ? a? ? y 5?x ? a? ? y 2 I 0;yyyyy ? ? 8y 2 ? y r8 r8 1 2 The terms h1, h2, r1 and r2 are de?ned in Eq. (9) in Section 2.1.
References
Aliabadi, M.H., 1998. Fracture of Rocks. Computational Mechanics Publications, Southampton, UK. Aliabadi, M.H., Rooke, D.P., 1991. Numerical Fracture Mechanics. Computational Mechanics Publications, Southampton, UK. Atkinson, C., Smelser, R.E., Sanchez, J., 1982. Combined mode fracture via the cracked Brazilian disc test. Int. J. Fract. 18, 279–291. Backers, T., Stephansson, O., Rybacki, E., 2002. Rock fracture toughness testing in mode II—punch-through shear test. Int. J. Rock Mech. Miner. Sci. 41 (3). Backers, T., Fardin, N., Dresen, G., Stephansson, O., 2003. E?ect of loading rate on mode I fracture toughness, roughness and micromechanics of sandstone. Int. J. Rock Mech. Miner. Sci. 41 (3). Backers, T., Dresen, G., Rybacki, E., Stephansson, O., 2004. New data on Mode II fracture toughness of rock from the punch-through shear test. In: SINOROCK 2004 Symposium, Paper 1A 01, Int. J. Rock Mech. Miner. Sci. 41(3). Blandford, G.E., Ingra?ea, A.R., Ligget, J.A., 1982. Two dimensional stress intensity factor computations using the boundary element method. Int. J. Numer. Methods Eng. 17, 387–404. Bobet, A., 2001. A hybridized displacement discontinuity method for mixed mode I–II–III loading, Int. J. Rock Mech. Miner. Sci. 38, 1121–1134. Broek, D., 1989. The Practical Use of Fracture Mechanics, fourth ed. Kluwer Academic Publishers, The Netherlands. Carpinteri, A., Yang, G., 1997. Size e?ects in brittle specimen with microcrack interaction. Comput. Struct. 63, 429–437. Cotterell, B., Rice, J.R., 1980. Slightly curved or kinked cracks. Int. J. Fract. Mech. 16, 155–169. Crawford, A.M., Curran, J.H., 1982. Higher order functional variation displacement discontinuity elements. Int. J. Rock Mech. Miner. Sci. Geomech. Abstr. 19, 143–148. Crouch, S.L., 1976. Solution of plane elasticity problems by the displacement discontinuity method. Int. J. Numer. Methods Eng. 10, 301–343. Crouch, S.L., Star?eld, A.M., 1983. Boundary Element Methods in Solid Mechanics. Allen and Unwin, London. Erdogan, F., Sih, G.C., 1963. On the crack extension in plates under plane loading and transverse shear. J. Basic Eng. 85, 519–527. Fowell, R.J., 1995. Suggested methods for determining mode I fracture toughness using Chevren notched Brazilian disc specimen. Int. J. Rock Mech. Miner. Sci. 32 (1), 57–64. Guo, H., Aziz, N.I., Schmidt, R.A., 1990. Linear elastic crack tip modeling by displacement discontinuity method. Eng. Fract. Mech. 36, 933–943. Guo, H., Aziz, N.I., Schmidt, R.A., 1992. Rock cutting study using linear elastic fracture mechanics. Eng. Fract. Mech. 41, 771–778. Hwang, C.G., Ingra?ea, A.R., 2004. Shape prediction and stability analysis of Mode-I planar cracks. Eng. Fract. Mech. 71, 1751–1777. Huang, J., Wang, S., 1985. An experimental investigation concerning the comprehensive fracture toughness of some brittle rocks. Int. J. Rock Mech. Miner. Sci. Geomech. Abstr. 22, 99–104. Ingra?ea, A.R., 1981. Mixed-mode fracture initiation in Indiana Sandstone and Westerly Granite. In: Einstein, H.H. (Ed.), Rock Mechanics from Research to Application. Proceedings of the 22nd US Symposium in Rock Mechanics. MIT, Cambridge, MA, pp. 199–204. Ingra?ea, A.R., 1983. Numerical modeling of fracture propagation. In: Rossmanith, H.P. (Ed.), Rock Fracture Mechanics. SpringerVerlag, Wien, New York, pp. 151–208. Ingra?ea, A.R., 1987. Theory of crack initiation and propagation in rock. In: Atkinston, B.K. (Ed.), Fracture Mechanics of Rock, Geology Series. Academy Press, New York, pp. 71–110, pp. 151–208. Ingra?ea, A.R., Hueze, F.E., 1980. Finite element models for rock fracture mechanics. Int. J. Numer. Anal. Methods Geomech. 4, 25. Ouchterlony, F., 1983. Analysis of cracks. In: Rossmanith, H.P. (Ed.), Rock Fracture Mechanics. Springer-Verlag, Wien, New York, pp. 31–67. Ouchterlony, F., 1988. Suggested methods for determining the fracture toughness of rock, ISRM commission on testing methods. Int. J. Rock Mech. Miner. Sci. Geomech. Abstr. 25, 71–96.
1692
M.F. Marji et al. / International Journal of Solids and Structures 43 (2006) 1669–1692
Pang, H.L.J., 1995. Mixed mode fracture analysis and toughness of adhesive joints. Eng. Fract. Mech. 51, 575–583. Rao, Q., Sun, Z., Stephansson, O., Li, C., Stillborg, B., 2003. Shear fracture (Mode II) of brittle rock. Int. J. Rock Mech. Miner. Sci. 40, 355–375. Scavia, C., 1990. Fracture mechanics approach to stability analysis of crack slopes. Eng. Fract. Mech. 35, 889–910. Scavia, C., 1992. A numerical technique for the analysis of cracks subjected to normal compressive stresses. Int. J. Numer. Methods Eng. 33, 929–942. Scavia, C., 1995. A method for the study of crack propagation in rock structures. Geotechnics 45, 447–463. Shen, B., Stephansson, O., 1994. Modi?cation of the G-criterion for crack propagation subjected to compression. Eng. Fract. Mech. 47 (2), 177–189. Shen, B., Stephansson, O., Rinne, M., Lee, H.-S., Jing, L., Rosho?, K., 2004. A fracture propagation code and its application to nuclear waste disposal. Int. J. Rock Mech. Miner. Sci. 41 (3), 448–453. Shou, K.J., Crouch, S.L., 1995. A higher order displacement discontinuity method for analysis of crack problems. Int. J. Rock Mech. Miner. Sci. Geomech. Abstr. 32, 49–55. Sih, G.C., 1973. Methods of analysis and solutions of crack ProblemsMech. Fract., vol. 1. Noordho? International Publishing, Leyden. Sneddon, I.N., 1951. Fourier Transforms. McGraw-Hill, New York. Stephansson, O., 2002. Recent rock fracture mechanics developments. In: 1st Iranian Rock Mechanics Conference, pp. 675–698. Stephansson, O., Backers, T., Dresen, G., Rybacki, E., 2001. Shear fracture mechanics of rocks and a new testing method for KIIC. In: Sarkka P. & Eloranta P., (Eds.), Rocks Mechanics a Challenge for Society, Proceedings of the ISRM Regional Symposium, EUROCK 2001, Finland: EPSOO, pp. 163–168. Sun, G.X., Whittaker, B.N., Sing, R.N., 1990. Use of Brazilian disc test for determination of the mixed Mode I–II Fracture toughness envelope of rock. In: International Conference on Mechanical of Jointed and Fractured Rocks, Vienna, Austria, pp. 18–20. Swartz, S.E., Lu, L.W., Tang, L.D., Refai, T.M.E., 1988. Mode II fracture parameters estimates for concrete from beam specimens. Exp. Mech. 28, 46–53. Tan, X.C., Kou, S.Q., Lindqvist, P.A., 1996. Simulation of rock fragmentation by indenters using DDM and fracture mechanics. In: Aubertin, M., Hassani, F., Mitri, H. (Eds.), Rock Mechanics, Tools and Techniques. Balkema, Roterdam. Whittaker, B.N., Singh, R.N., Sun, G., 1992. Rock Fracture Mechanics, Principles, Design and Applications. Elsevier, The Netherlands. Zipf, R.K., Bieniawski, Z.T., 1987. Development of the mixed mode testing system for geological materials. In: Shah, S.P., Swartz, S.E., (Eds.), SEM/RILEM International Conference on Fracture of Concrete and Rock, Houston, TX, pp. 338–352.