**Abstract**

Stress intensity factor determination plays a central role in linearly elastic fracture mechanics (LEFM) problems. Fracture propagation is controlled by the stress field near the crack tip. Because this stress field is asymptotic dominant or singular, it is characterized by the stress intensity factor (SIF). Since many rock types show brittle elastic behaviour under hydrocarbon reservoir conditions, LEFM can be satisfactorily used for studying hydraulic fracture devel‐ opment. The purpose of this paper is to describe a numerical method to evaluate the stress intensity factor in Mode I, II and III at the tip of an arbitrarily-shaped, embedded cracks. The stress intensity factor is evaluated directly based on displacement discontinuities (DD) using a three-dimensional displacement discontinuity, boundary element method based on the equations of proposed in [1]. The boundary element formulation incorporates the fundamental closed-form analytical solution to a rectangular discontinuity in a homogenous, isotropic and linearly elastic half space. The accuracy of the stress intensity factor calculation is satisfactorily examined for rectangular, penny-shaped and elliptical planar cracks. Accurate and fast evaluation of the stress intensity factor for planar cracks shows the proposed procedure is robust for SIF calculation and crack propagation purposes. The empirical constant proposed by [2] relating crack tip element displacement discontinuity and SIF values provides surpris‐ ingly accurate results for planar cracks with limited numbers of constant DD elements. Using the described numerical model, we study how fracturing from misaligned horizontal well‐ bores might results in non-uniform height growth of the hydraulic fracture by evaluating of SIF distribution along the upper front of the fracture.

© 2013 Sheibani and Olson; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### **1. Introduction** %&?>immediately<\$%&?>around<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>in<\$%&?>comparison<\$%&?>with<\$%&?>the<\$%&?> K-

<\$%&?>analysis.

Stress intensity factor determination plays a central role in linearly elastic fracture mechanics problems. Fracture propagation is controlled by the stress field near the crack tip. Because the stress field near the crack tip is asymptotic dominant or singular, it is characterized by the stress intensity factor. The real stress distribution at the vicinity of crack tip and the K-field LEFM approximation can be depicted schematically as in Figure 1. The stress singularity right at the tip of the crack cannot be experienced in real nature because inelastic deformation prevents the crack tip from being perfectly sharp. However, according to small scale yielding of the process zone immediately around the crack tip in comparison with the K-field region (Figure 2), the SIF is the quantity which dictates if/when the crack will propagate. The inaccuracy of the stress field calculation using the SIF based on LEFM is less than 15% of the exact solution over the distance ranging from *r* <*0.01a* to *r* <*0.15a*, where *r* is the radius of Kfield region and *a* is the half length of the crack [4]. field<\$%&?>region<\$%&?>(Figure<\$%&?>2),<\$%&?>the<\$%&?>SIF<\$%&?>is<\$%&?>the<\$%&?>quantity<\$%&?>which<\$%&?>dict ates<\$%&?>if/when<\$%&?>the<\$%&?>crack<\$%&?>will<\$%&?>propagate.<\$%&?>The<\$%&?>inaccuracy<\$%&?>of<\$%&?>the<\$% &?>stress<\$%&?>field<\$%&?>calculation<\$%&?>using<\$%&?>the<\$%&?>SIF<\$%&?>based<\$%&?>on<\$%&?>LEFM<\$%&?>is<\$%& ?>less<\$%&?>than<\$%&?>15%<\$%&?>of<\$%&?>the<\$%&?>exact<\$%&?>solution<\$%&?>over<\$%&?>the<\$%&?>distance<\$%&?>ra nging<\$%&?>from*<\$%&?>r<0.01a<\$%&?>*to<\$%&?>*r<0.15a*,<\$%&?>where<\$%&?>*r*<\$%&?>is<\$%&?>the<\$%&?>radius<\$%&?>of<\$ %&?>Kfield<\$%&?>region<\$%&?>and<\$%&?>*a*<\$%&?>is<\$%&?>the<\$%&?>half<\$%&?>length<\$%&?>of<\$%&?>the<\$%&?>crack<\$%&?>[ 4]. Since<\$%&?>SIF<\$%&?>was<\$%&?>proposed<\$%&?>by<\$%&?>Irwin<\$%&?>[5]<\$%&?>to<\$%&?>express<\$%&?>displacements<\$ %&?>and<\$%&?>stresses<\$%&?>in<\$%&?>the<\$%&?>vicinity<\$%&?>of<\$%&?>crack<\$%&?>tip,<\$%&?>several<\$%&?>analytical< \$%&?>techniques<\$%&?>have<\$%&?>been<\$%&?>developed<\$%&?>for<\$%&?>a<\$%&?>variety<\$%&?>of<\$%&?>common<\$%&? >crack<\$%&?>configurations;<\$%&?>however,<\$%&?>these<\$%&?>analytical<\$%&?>solutions<\$%&?>are<\$%&?>limited<\$%&?>t

eformation<\$%&?>prevents<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>from<\$%&?>being<\$%&?>perfectly<\$%&?>sharp.<\$%&?>H owever,<\$%&?>according<\$%&?>to<\$%&?>small<\$%&?>scale<\$%&?>yielding<\$%&?>of<\$%&?>the<\$%&?>process<\$%&?>zone<\$

?>crack<\$%&?>cannot<\$%&?>be<\$%&?>experienced<\$%&?>in<\$%&?>real<\$%&?>nature<\$%&?>because<\$%&?>inelastic<\$%&?>d eformation<\$%&?>prevents<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>from<\$%&?>being<\$%&?>perfectly<\$%&?>sharp.<\$%&?>H owever,<\$%&?>according<\$%&?>to<\$%&?>small<\$%&?>scale<\$%&?>yielding<\$%&?>of<\$%&?>the<\$%&?>process<\$%&?>zone<\$ %&?>immediately<\$%&?>around<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>in<\$%&?>comparison<\$%&?>with<\$%&?>the<\$%&?>

field<\$%&?>region<\$%&?>(Figure<\$%&?>2),<\$%&?>the<\$%&?>SIF<\$%&?>is<\$%&?>the<\$%&?>quantity<\$%&?>which<\$%&?>dict ates<\$%&?>if/when<\$%&?>the<\$%&?>crack<\$%&?>will<\$%&?>propagate.<\$%&?>The<\$%&?>inaccuracy<\$%&?>of<\$%&?>the<\$% &?>stress<\$%&?>field<\$%&?>calculation<\$%&?>using<\$%&?>the<\$%&?>SIF<\$%&?>based<\$%&?>on<\$%&?>LEFM<\$%&?>is<\$%& ?>less<\$%&?>than<\$%&?>15%<\$%&?>of<\$%&?>the<\$%&?>exact<\$%&?>solution<\$%&?>over<\$%&?>the<\$%&?>distance<\$%&?>ra nging<\$%&?>from*<\$%&?>r<0.01a<\$%&?>*to<\$%&?>*r<0.15a*,<\$%&?>where<\$%&?>*r*<\$%&?>is<\$%&?>the<\$%&?>radius<\$%&?>of<\$

field<\$%&?>region<\$%&?>and<\$%&?>*a*<\$%&?>is<\$%&?>the<\$%&?>half<\$%&?>length<\$%&?>of<\$%&?>the<\$%&?>crack<\$%&?>[

Since<\$%&?>SIF<\$%&?>was<\$%&?>proposed<\$%&?>by<\$%&?>Irwin<\$%&?>[5]<\$%&?>to<\$%&?>express<\$%&?>displacements<\$ %&?>and<\$%&?>stresses<\$%&?>in<\$%&?>the<\$%&?>vicinity<\$%&?>of<\$%&?>crack<\$%&?>tip,<\$%&?>several<\$%&?>analytical< \$%&?>techniques<\$%&?>have<\$%&?>been<\$%&?>developed<\$%&?>for<\$%&?>a<\$%&?>variety<\$%&?>of<\$%&?>common<\$%&? >crack<\$%&?>configurations;<\$%&?>however,<\$%&?>these<\$%&?>analytical<\$%&?>solutions<\$%&?>are<\$%&?>limited<\$%&?>t o<\$%&?>simple<\$%&?>crack<\$%&?>geometries<\$%&?>and<\$%&?>loading<\$%&?>conditions.<\$%&?>For<\$%&?>the<\$%&?>case<

infinite<\$%&?>body,<\$%&?>there<\$%&?>are<\$%&?>less<\$%&?>available<\$%&?>analytical<\$%&?>solutions<\$%&?>for<\$%&?>SIF .<\$%&?>These<\$%&?>exact<\$%&?>analytical<\$%&?>solutions<\$%&?>provide<\$%&?>good<\$%&?>insight<\$%&?>about<\$%&?>fra cture<\$%&?>problems<\$%&?>but<\$%&?>they<\$%&?>are<\$%&?>not<\$%&?>usable<\$%&?>for<\$%&?>general<\$%&?>crack<\$%&?> propagation<\$%&?>modeling<\$%&?>where<\$%&?>the<\$%&?>geometry<\$%&?>of<\$%&?>simultaneously<\$%&?>propagating<\$% &?>cracks<\$%&?>can<\$%&?>be<\$%&?>asymmetrical<\$%&?>and<\$%&?>irregular<\$%&?>and<\$%&?>the<\$%&?>boundary<\$%&?> conditions<\$%&?>can<\$%&?>be<\$%&?>complicated.<\$%&?>Fortunately,<\$%&?>advances<\$%&?>in<\$%&?>numerical<\$%&?>mo deling<\$%&?>procedures<\$%&?>supported<\$%&?>by<\$%&?>the<\$%&?>fast<\$%&?>growing<\$%&?>speed<\$%&?>of<\$%&?>com putational<\$%&?>calculation<\$%&?>have<\$%&?>opened<\$%&?>new<\$%&?>doors<\$%&?>for<\$%&?>fracture<\$%&?>propagation

Figure 1. Schematic<\$%&?>representation<\$%&?>of<\$%&?>stress<\$%&?>distribution<\$%&?>around<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>[3]

http://dx.doi.org/10.5772/56308

743

Inelastic deformation (Process Zone)

K‐field dominates in this "annulas" region

K dominant zone

Actual stress distribution

Figure 2. Process<\$%&?>zone<\$%&?>and<\$%&?>K-filed<\$%&?>representation<\$%&?>[3]

There are four general distinctive numerical methods to model fracture propagation problems: **1.** The boundary element method (BEM) requires discretization and calculation only on boundaries of the domain. The stress resolution is higher in comparison with finite element and finite difference methods because the approximation is imposed only on boundaries of the domain, and there is no further approximation on the solution at interior points. Particularly, for some problems where the ratio of boundary surface to volume is high (for instance for large rock masses), BEM can be advantageous because FEM or other whole-domain-discretizing methods require larger numbers of elements to achieve the

**2.** The Finite Element Method (FEM) has been widely used in fracture mechanics problems since it was implemented by [6] for SIF calculation. Several modifications have proposed to remove its deficiencies in LEFM problem modeling. [7] and [8] devised "quarter point element" or "singularity elements" to improve the accuracy of stress and displacement distributions around the crack and SIF evaluation. To overcome the time consuming process of remeshing in fracture propagation problems, [9] proposed the Extended Finite Element Method (XFEM). XFEM allows fracture propagation without changing the mesh by adding analytical expressions related to the crack tip field to the conventional FE polynomial approximation in what are called "enriched elements". Further work is being done ([10] and [11]) to address the accuracy and stability of XFEM modeling, especially for multiple crack problems and approaching tip elements called "blending elements". **3.** The Finite Difference Method (FDM) requires calculations on a mesh that includes the entire domain. FDM usage in fracture mechanics is mostly limited to dynamic fracture

**4.** The Discrete Element Method (DEM) is mostly applied when continuity cannot be assumed in discontinuous, separated domains. The method apply to describe the behavior of discontinuities between bodies with emphasize on the solution of contact and impact

propagation and dynamic SIF calculation ([12] and [13].)

between multiple bodies [14].

Elastic deformation however, still include perturbation by process zone; perhaps the K‐ field cannot accurately present the stress field

\$%&?>of<\$%&?>3-D<\$%&?>planar<\$%&?>cracks<\$%&?>embedded<\$%&?>in<\$%&?>a<\$%&?>semi-

r

Process Zone

Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method...

Stress field according K‐field

K-

%&?>K-

<\$%&?>analysis.

**Figure 2.** Process zone and K-filed representation [3]

same accuracy.

Yielding Stress

σ

4].

Since SIF was proposed by Irwin [5] to express displacements and stresses in the vicinity of crack tip, several analytical techniques have been developed for a variety of common crack configurations; however, these analytical solutions are limited to simple crack geometries and loading conditions. For the case of 3-D planar cracks embedded in a semi-infinite body, there are less available analytical solutions for SIF. These exact analytical solutions provide good insight about fracture problems but they are not usable for general crack propagation modeling where the geometry of simultaneously propagating cracks can be asymmetrical and irregular and the boundary conditions can be complicated. Fortunately, advances in numerical model‐ ing procedures supported by the fast growing speed of computational calculation have opened new doors for fracture propagation analysis. o<\$%&?>simple<\$%&?>crack<\$%&?>geometries<\$%&?>and<\$%&?>loading<\$%&?>conditions.<\$%&?>For<\$%&?>the<\$%&?>case< \$%&?>of<\$%&?>3-D<\$%&?>planar<\$%&?>cracks<\$%&?>embedded<\$%&?>in<\$%&?>a<\$%&?>semiinfinite<\$%&?>body,<\$%&?>there<\$%&?>are<\$%&?>less<\$%&?>available<\$%&?>analytical<\$%&?>solutions<\$%&?>for<\$%&?>SIF .<\$%&?>These<\$%&?>exact<\$%&?>analytical<\$%&?>solutions<\$%&?>provide<\$%&?>good<\$%&?>insight<\$%&?>about<\$%&?>fra cture<\$%&?>problems<\$%&?>but<\$%&?>they<\$%&?>are<\$%&?>not<\$%&?>usable<\$%&?>for<\$%&?>general<\$%&?>crack<\$%&?> propagation<\$%&?>modeling<\$%&?>where<\$%&?>the<\$%&?>geometry<\$%&?>of<\$%&?>simultaneously<\$%&?>propagating<\$% &?>cracks<\$%&?>can<\$%&?>be<\$%&?>asymmetrical<\$%&?>and<\$%&?>irregular<\$%&?>and<\$%&?>the<\$%&?>boundary<\$%&?> conditions<\$%&?>can<\$%&?>be<\$%&?>complicated.<\$%&?>Fortunately,<\$%&?>advances<\$%&?>in<\$%&?>numerical<\$%&?>mo deling<\$%&?>procedures<\$%&?>supported<\$%&?>by<\$%&?>the<\$%&?>fast<\$%&?>growing<\$%&?>speed<\$%&?>of<\$%&?>com putational<\$%&?>calculation<\$%&?>have<\$%&?>opened<\$%&?>new<\$%&?>doors<\$%&?>for<\$%&?>fracture<\$%&?>propagation

Inelastic deformation (Process Zone)

K‐field dominates in this "annulas" region

Figure 2. Process<\$%&?>zone<\$%&?>and<\$%&?>K-filed<\$%&?>representation<\$%&?>[3]

Figure 1. Schematic<\$%&?>representation<\$%&?>of<\$%&?>stress<\$%&?>distribution<\$%&?>around<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>[3] **Figure 1.** Schematic representation of stress distribution around the crack tip [3]

Elastic deformation however, still include perturbation by process zone; perhaps the K‐ field cannot accurately present the stress field

r

Stress field according K‐field

σ

?>crack<\$%&?>cannot<\$%&?>be<\$%&?>experienced<\$%&?>in<\$%&?>real<\$%&?>nature<\$%&?>because<\$%&?>inelastic<\$%&?>d eformation<\$%&?>prevents<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>from<\$%&?>being<\$%&?>perfectly<\$%&?>sharp.<\$%&?>H owever,<\$%&?>according<\$%&?>to<\$%&?>small<\$%&?>scale<\$%&?>yielding<\$%&?>of<\$%&?>the<\$%&?>process<\$%&?>zone<\$ %&?>immediately<\$%&?>around<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>in<\$%&?>comparison<\$%&?>with<\$%&?>the<\$%&?>

field<\$%&?>region<\$%&?>(Figure<\$%&?>2),<\$%&?>the<\$%&?>SIF<\$%&?>is<\$%&?>the<\$%&?>quantity<\$%&?>which<\$%&?>dict ates<\$%&?>if/when<\$%&?>the<\$%&?>crack<\$%&?>will<\$%&?>propagate.<\$%&?>The<\$%&?>inaccuracy<\$%&?>of<\$%&?>the<\$% &?>stress<\$%&?>field<\$%&?>calculation<\$%&?>using<\$%&?>the<\$%&?>SIF<\$%&?>based<\$%&?>on<\$%&?>LEFM<\$%&?>is<\$%& ?>less<\$%&?>than<\$%&?>15%<\$%&?>of<\$%&?>the<\$%&?>exact<\$%&?>solution<\$%&?>over<\$%&?>the<\$%&?>distance<\$%&?>ra nging<\$%&?>from*<\$%&?>r<0.01a<\$%&?>*to<\$%&?>*r<0.15a*,<\$%&?>where<\$%&?>*r*<\$%&?>is<\$%&?>the<\$%&?>radius<\$%&?>of<\$

field<\$%&?>region<\$%&?>and<\$%&?>*a*<\$%&?>is<\$%&?>the<\$%&?>half<\$%&?>length<\$%&?>of<\$%&?>the<\$%&?>crack<\$%&?>[

Since<\$%&?>SIF<\$%&?>was<\$%&?>proposed<\$%&?>by<\$%&?>Irwin<\$%&?>[5]<\$%&?>to<\$%&?>express<\$%&?>displacements<\$ %&?>and<\$%&?>stresses<\$%&?>in<\$%&?>the<\$%&?>vicinity<\$%&?>of<\$%&?>crack<\$%&?>tip,<\$%&?>several<\$%&?>analytical< \$%&?>techniques<\$%&?>have<\$%&?>been<\$%&?>developed<\$%&?>for<\$%&?>a<\$%&?>variety<\$%&?>of<\$%&?>common<\$%&? >crack<\$%&?>configurations;<\$%&?>however,<\$%&?>these<\$%&?>analytical<\$%&?>solutions<\$%&?>are<\$%&?>limited<\$%&?>t o<\$%&?>simple<\$%&?>crack<\$%&?>geometries<\$%&?>and<\$%&?>loading<\$%&?>conditions.<\$%&?>For<\$%&?>the<\$%&?>case<

infinite<\$%&?>body,<\$%&?>there<\$%&?>are<\$%&?>less<\$%&?>available<\$%&?>analytical<\$%&?>solutions<\$%&?>for<\$%&?>SIF .<\$%&?>These<\$%&?>exact<\$%&?>analytical<\$%&?>solutions<\$%&?>provide<\$%&?>good<\$%&?>insight<\$%&?>about<\$%&?>fra cture<\$%&?>problems<\$%&?>but<\$%&?>they<\$%&?>are<\$%&?>not<\$%&?>usable<\$%&?>for<\$%&?>general<\$%&?>crack<\$%&?> propagation<\$%&?>modeling<\$%&?>where<\$%&?>the<\$%&?>geometry<\$%&?>of<\$%&?>simultaneously<\$%&?>propagating<\$% &?>cracks<\$%&?>can<\$%&?>be<\$%&?>asymmetrical<\$%&?>and<\$%&?>irregular<\$%&?>and<\$%&?>the<\$%&?>boundary<\$%&?> conditions<\$%&?>can<\$%&?>be<\$%&?>complicated.<\$%&?>Fortunately,<\$%&?>advances<\$%&?>in<\$%&?>numerical<\$%&?>mo deling<\$%&?>procedures<\$%&?>supported<\$%&?>by<\$%&?>the<\$%&?>fast<\$%&?>growing<\$%&?>speed<\$%&?>of<\$%&?>com putational<\$%&?>calculation<\$%&?>have<\$%&?>opened<\$%&?>new<\$%&?>doors<\$%&?>for<\$%&?>fracture<\$%&?>propagation

Figure 1. Schematic<\$%&?>representation<\$%&?>of<\$%&?>stress<\$%&?>distribution<\$%&?>around<\$%&?>the<\$%&?>crack<\$%&?>tip<\$%&?>[3]

K dominant zone

Actual stress distribution

\$%&?>of<\$%&?>3-D<\$%&?>planar<\$%&?>cracks<\$%&?>embedded<\$%&?>in<\$%&?>a<\$%&?>semi-

\$%&?>techniques<\$%&?>have<\$%&?>been<\$%&?>developed<\$%&?>for<\$%&?>a<\$%&?>variety<\$%&?>of<\$%&?>common<\$%&? Figure 2. Process<\$%&?>zone<\$%&?>and<\$%&?>K-filed<\$%&?>representation<\$%&?>[3] **Figure 2.** Process zone and K-filed representation [3]

K-

%&?>K-

<\$%&?>analysis.

Yielding Stress

4].

**1. Introduction**

742 Effective and Sustainable Hydraulic Fracturing

K-

%&?>K-

4].

field region and *a* is the half length of the crack [4].

new doors for fracture propagation analysis.

σ

**Figure 1.** Schematic representation of stress distribution around the crack tip [3]

<\$%&?>analysis.

Yielding Stress

Stress intensity factor determination plays a central role in linearly elastic fracture mechanics problems. Fracture propagation is controlled by the stress field near the crack tip. Because the stress field near the crack tip is asymptotic dominant or singular, it is characterized by the stress intensity factor. The real stress distribution at the vicinity of crack tip and the K-field LEFM approximation can be depicted schematically as in Figure 1. The stress singularity right at the tip of the crack cannot be experienced in real nature because inelastic deformation prevents the crack tip from being perfectly sharp. However, according to small scale yielding of the process zone immediately around the crack tip in comparison with the K-field region (Figure 2), the SIF is the quantity which dictates if/when the crack will propagate. The inaccuracy of the stress field calculation using the SIF based on LEFM is less than 15% of the exact solution over the distance ranging from *r* <*0.01a* to *r* <*0.15a*, where *r* is the radius of K-

owever,<\$%&?>according<\$%&?>to<\$%&?>small<\$%&?>scale<\$%&?>yielding<\$%&?>of<\$%&?>the<\$%&?>process<\$%&?>zone<\$

Since SIF was proposed by Irwin [5] to express displacements and stresses in the vicinity of crack tip, several analytical techniques have been developed for a variety of common crack configurations; however, these analytical solutions are limited to simple crack geometries and loading conditions. For the case of 3-D planar cracks embedded in a semi-infinite body, there are less available analytical solutions for SIF. These exact analytical solutions provide good insight about fracture problems but they are not usable for general crack propagation modeling where the geometry of simultaneously propagating cracks can be asymmetrical and irregular and the boundary conditions can be complicated. Fortunately, advances in numerical model‐ ing procedures supported by the fast growing speed of computational calculation have opened

Figure 2. Process<\$%&?>zone<\$%&?>and<\$%&?>K-filed<\$%&?>representation<\$%&?>[3]

Elastic deformation however, still include perturbation by process zone; perhaps the K‐ field cannot accurately present the stress field

\$%&?>of<\$%&?>3-D<\$%&?>planar<\$%&?>cracks<\$%&?>embedded<\$%&?>in<\$%&?>a<\$%&?>semi-

r

Process Zone

Stress field according K‐field

Inelastic deformation (Process Zone)

K‐field dominates in this "annulas" region

K dominant zone

Actual stress distribution

>crack<\$%&?>configurations;<\$%&?>however,<\$%&?>these<\$%&?>analytical<\$%&?>solutions<\$%&?>are<\$%&?>limited<\$%&?>t o<\$%&?>simple<\$%&?>crack<\$%&?>geometries<\$%&?>and<\$%&?>loading<\$%&?>conditions.<\$%&?>For<\$%&?>the<\$%&?>case< There are four general distinctive numerical methods to model fracture propagation problems:

	- **3.** The Finite Difference Method (FDM) requires calculations on a mesh that includes the entire domain. FDM usage in fracture mechanics is mostly limited to dynamic fracture propagation and dynamic SIF calculation ([12] and [13].)

Generally, when the geometry of a problem is changing, whole-domain-discretizing methods like FEM, FDM and DEM are more time-consuming than BEM because of the remeshing process around a propagation fracture. However, BEM loses its advantage when the domain is grossly inhomogeneous.

All of the methods mentioned above including using special crack tip elements or equivalence transformation methods to decrease the error in crack tip element displacement and corre‐ sponding SIF calculation; however, they all need numerical integration and can be more timeconsuming than constant elemental DD approximation. Ref. [2] empirically determined the coincidence between DDM modeling and analytical displacement distribution solution of a straight 2-D crack to remove the error. He showed the margin of error is less than 5% even by using only 2 elements in a 2-D crack. His proposed formula has been widely used in geologic fracture problems [26-29]. This paper extends Olson's method ion [2] to SIF calculation for 3- D homogenous, isotropic and linearly elastic material problems. [30] changed the correction constant. The empirical constant they proposed was used by some researchers afterwards ([31]

Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method...

According to Murakami and [33] and [34] the maximum mode I stress intensity factor appearing at a certain point along the crack front can be estimated by Equation (1) with less

where 'area' is the area of crack projected in the direction of the maximum principal stress. Fortunately, for simple crack geometries like elliptical and circular cracks, there exist analytical formulae for mode I stress intensity factor variation along the crack tip which help us to evaluate the accuracy of the numerical modeling ([35] and [36]). For rectangular defects there are no analytical formulae, but the accuracy of DDM numerical modeling can be examined by

The general concept of the displacement discontinuity method proposed by [16] is to approx‐ imate the distribution of displacement discontinuity of a crack by discretizing it into elements. Knowing the analytical solution for one element, the numerical elastic solution of the whole

The 3-D displacement discontinuity used here is based on the analytical elastic solution of normal and shear displacement of a finite rectangular discontinuity in half-space (Figure 3) proposed by [1]. These equations are closed-form half-space solutions of deformations and deformation derivatives in which most of singularities and mathematical instabilities were

By placing *N* unknown constant displacement elements within the boundaries of the region to be analyzed and knowing the boundary conditions on each element (traction or displace‐

comparing against earlier numerical work using integral equation methods [37-40].

discontinuity can be calculated by adding up the effect of all subdividing elements.

ment), a system of *3N* linear algebraic equations can be set up as the following:

*KI max* =0.50*σ π area* (1)

http://dx.doi.org/10.5772/56308

745

and [32]), but we argue the change does not actually improve SIF accuracy.

than 10% error for an arbitrary-shaped planar crack.

**2. Numerical procedure**

removed.

**2.1. Displacement discontinuity method**

The "Integral equation" approach (also called influence function) and the "displacement discontinuity method" are two types of BEM widely used in LEFM analysis. Both approaches incorporate only boundary data by relating boundary tractions and displacements. In the integral equation technique, superposition of known influence functions (called Green's function) along boundaries generate a system of simultaneous integral equations [15]. In DDM, unknown boundary values are found from a simple system of algebraic equitation [16]. Generally, DDM has the advantage over integral equations in being faster, while integral equations can be more accurate for non-linear problems.

SIF values can be obtained from the displacement discontinuity magnitudes at crack tip elements [17-19]. However, according to [16], DDM consistently overestimates displace‐ ment discontinuities at the tip of the crack (considering element midpoint) by as much as 25%. To improve the accuracy of the solution, some researchers proposed using higher accuracy crack tip element and/or using relatively denser distribution of elements near the crack tip. [20] proposed higher order elements to improve the DDM solution and they used numerical integration to find the fundamental solution of linear and quadratic displace‐ ment discontinuities. [21] proposed another approach called "hybrid displacement discontinuity method" by using parabolic DD for crack tip elements and constant DD for other elements. He concluded increasing the number of elements more than 8-10 times cannot yield more accurate results and the error in mode I stress intensity factor calcula‐ tion for a 2-D straight crack with uniform internal pressure, sporadically changes in a range of 1% to about 10% depending on the ratio of parabolic element length to constant element length. However, [22] used the same combination of DD element and concluded the ratio of crack tip element to constant DD element must be between 1-1.3 to obtain good results with relative error less than 3% in mode II SIF calculation for a straight 2-D crack. [23] presented a new hybrid displacement discontinuity method by using quadratic DD elements and special crack tip elements to show *r* variation of displacement near the crack tip. [24] used the same method with few modifications about the position of collocation points to determine quadratic elemental displacement. They showed the error can be fixed up to 1.5% for Mode I, and about 2% for mode II SIF calculation for a slanted straight crack. [25] took a different approach; instead of direct calculation of stress intensity factors from displacement discontinuities, they proposed a "equivalence transformation method" in which stresses on the crack surface are calculated from displacement discontinuities, and then by using crack line Green's function, the SIF at the crack tip can be obtained from calculated stresses. They implemented the equivalence transformation method to calcu‐ late dynamic stress intensity factors for an isolated 2-D crack in an infinite sheet subject‐ ed to Heaviside loading. By comparison with the exact solution and using 80 DD elements, they inferred the error in mode I SIF is less than 1% and for mode II doesn't exceed 1.5%. All of the methods mentioned above including using special crack tip elements or equivalence transformation methods to decrease the error in crack tip element displacement and corre‐ sponding SIF calculation; however, they all need numerical integration and can be more timeconsuming than constant elemental DD approximation. Ref. [2] empirically determined the coincidence between DDM modeling and analytical displacement distribution solution of a straight 2-D crack to remove the error. He showed the margin of error is less than 5% even by using only 2 elements in a 2-D crack. His proposed formula has been widely used in geologic fracture problems [26-29]. This paper extends Olson's method ion [2] to SIF calculation for 3- D homogenous, isotropic and linearly elastic material problems. [30] changed the correction constant. The empirical constant they proposed was used by some researchers afterwards ([31] and [32]), but we argue the change does not actually improve SIF accuracy.

According to Murakami and [33] and [34] the maximum mode I stress intensity factor appearing at a certain point along the crack front can be estimated by Equation (1) with less than 10% error for an arbitrary-shaped planar crack.

$$K\_{I\\_max} = 0.50\sigma \sqrt{\pi \sqrt{area}}\tag{1}$$

where 'area' is the area of crack projected in the direction of the maximum principal stress.

Fortunately, for simple crack geometries like elliptical and circular cracks, there exist analytical formulae for mode I stress intensity factor variation along the crack tip which help us to evaluate the accuracy of the numerical modeling ([35] and [36]). For rectangular defects there are no analytical formulae, but the accuracy of DDM numerical modeling can be examined by comparing against earlier numerical work using integral equation methods [37-40].
