**2. Numerical procedure**

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

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

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%.

is grossly inhomogeneous.

744 Effective and Sustainable Hydraulic Fracturing

equations can be more accurate for non-linear problems.

#### **2.1. Displacement discontinuity method**

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 discontinuity can be calculated by adding up the effect of all subdividing elements.

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 removed.

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‐ ment), a system of *3N* linear algebraic equations can be set up as the following:

$$\begin{aligned} \sigma\_s^{\ j} &= \sum\_{j=1}^N A\_{\text{ss}}^{\ ij} \mathbf{D}\_s^{\ j} + \sum\_{j=1}^N A\_{\text{sd}}^{\ ij} \mathbf{D}\_d^{\ j} + \sum\_{j=1}^N A\_{\text{sn}}^{\ ij} \mathbf{D}\_n^{\ j} \ \ \ a \\ \sigma\_d^{\ j} &= \sum\_{j=1}^N A\_{\text{ss}}^{\ ij} \mathbf{D}\_s^{\ j} + \sum\_{j=1}^N A\_{\text{sd}}^{\ ij} \mathbf{D}\_d^{\ j} + \sum\_{j=1}^N A\_{\text{sn}}^{\ ij} \mathbf{D}\_n^{\ j} \ \ \ \ \ \ \ & \tag{2} \\ \sigma\_n^{\ j} &= \sum\_{j=1}^N A\_{\text{ns}}^{\ i j} \mathbf{D}\_s^{\ j} + \sum\_{j=1}^N A\_{\text{nd}}^{\ i j} \mathbf{D}\_d^{\ j} + \sum\_{j=1}^N A\_{\text{sn}}^{\ i j} \mathbf{D}\_n^{\ j} \ \ \ \ \ \ \end{aligned} \tag{2}$$

$$\begin{aligned} \boldsymbol{u}\_{s}^{\,j} &= \sum\_{j=1}^{N} \boldsymbol{B}\_{ss}^{\,ij} \boldsymbol{\mathcal{D}}\_{s}^{\,j} + \sum\_{j=1}^{N} \boldsymbol{B}\_{sd}^{\,ij} \boldsymbol{\mathcal{D}}\_{d}^{\,j} + \sum\_{j=1}^{N} \boldsymbol{B}\_{sn}^{\,ij} \boldsymbol{\mathcal{D}}\_{n}^{\,j} \quad \boldsymbol{a} \\ \boldsymbol{u}\_{d}^{\,j} &= \sum\_{j=1}^{N} \boldsymbol{B}\_{sd}^{\,ij} \boldsymbol{\mathcal{D}}\_{s}^{\,j} + \sum\_{j=1}^{N} \boldsymbol{B}\_{dd}^{\,ij} \boldsymbol{\mathcal{D}}\_{d}^{\,j} + \sum\_{j=1}^{N} \boldsymbol{B}\_{dn}^{\,ij} \boldsymbol{\mathcal{D}}\_{n}^{\,j} \quad \boldsymbol{b} \\ \boldsymbol{u}\_{n}^{\,j} &= \sum\_{j=1}^{N} \boldsymbol{B}\_{sn}^{\,ij} \boldsymbol{\mathcal{D}}\_{s}^{\,j} + \sum\_{j=1}^{N} \boldsymbol{B}\_{nd}^{\,ij} \boldsymbol{\mathcal{D}}\_{d}^{\,j} + \sum\_{j=1}^{N} \boldsymbol{B}\_{nn}^{\,ij} \boldsymbol{\mathcal{D}}\_{n}^{\,j} \quad \boldsymbol{c} \end{aligned} \tag{3}$$

$$\begin{aligned} K\_I &= \mathbb{C} \frac{D\_\nu E \sqrt{\pi}}{4(1 \cdot \nu^2) \sqrt{P}} \quad a\\ K\_{II} &= \mathbb{C} \frac{D\_\nu E \sqrt{\pi}}{4(1 \cdot \nu^2) \sqrt{P}} \quad b\\ K\_{III} &= \mathbb{C} \frac{D\_\nu E \sqrt{\pi}}{4(1 + \nu) \sqrt{P}} \quad c \end{aligned} \tag{4}$$

stress intensity factor at the edge-midpoints of a square crack as *FI* =0.768 which the error is about 1%.

Considering a rectangular crack as shown in Figure 4, the following dimensionless parameter is proposed to demonstrate the result of stress intensity factor of a rectangular crack. *FI* is called the dimensionless stress intensity factor along the crack front *y* =*b*:

$$\,\_{1}F\_{1} = \frac{K\_{1}\begin{pmatrix} x & y \end{pmatrix} \, \big|\, x = x, \quad y = x \, b}{\sigma\_{n}\sqrt{\pi ab}} \tag{5}$$

**Figure 5.** Error in strain energy calculation as a function of number of elements on each side of a pressurized square

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

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

749

The error in strain energy calculation is mainly related to the largest error occurring at the corners of the square crack where the displacement gradient is highest. Figure 6 shows the stress intensity factor variation along the half-length of the crack tip using DDM compared with the integral equation solution suggested by [40]. The total number of elements used in the simulation was 22×22 to be consistent with the number of colocation points used in [40]. The difference between these two solutions is negligible for all elements but the corners (element No. 11). However, the corner elements of rectangular cracks don't play an important role in fracture propagation problems because the level of SIF is the lowest there and unlikely

crack (*a*=*b*)

to control the initiation of crack propagation.

The stability of the solution can be examined by investigation of the strain energy variation through increasing the number of elements. Figure 5-b shows that strain energy (*U* ) linearly varies with <sup>1</sup> *<sup>n</sup>* and has an asymptotic behavior with respect to *n*, where *n* is the number of element on each side of a square crack shown in Figure 5-b. The area of the square crack is *A* under constant pressure *p*.Assuming the error in strain energy calculation approaches zero if *<sup>n</sup>* <sup>→</sup>*<sup>∞</sup>* ( <sup>1</sup> *<sup>n</sup>* →0), the correct answer for error estimation in the strain energy calculation can be obtained from Figure 5-b. Figure 5a shows the error calculation in strain energy. The displace‐ ment discontinuity method always overestimates the strain energy (or displacement across the crack surface) but it yields more accurate results closer to the exact solution when the number of elements increases. The error changes from 48.8% using a 3×3 mesh to about 1.99% for a mesh including 71×71 elements. In comparison with the two dimensional analysis of a straight crack [16], the rate of convergence is faster, but the error in strain energy calculation is higher using the same number of elements to divide one side of a crack.

**Figure 4.** Approximation of the exact solution of strain energy for a square pressurized crack

stress intensity factor at the edge-midpoints of a square crack as *FI* =0.768 which the error is

Considering a rectangular crack as shown in Figure 4, the following dimensionless parameter is proposed to demonstrate the result of stress intensity factor of a rectangular crack. *FI* is called

The stability of the solution can be examined by investigation of the strain energy variation through increasing the number of elements. Figure 5-b shows that strain energy (*U* ) linearly

element on each side of a square crack shown in Figure 5-b. The area of the square crack is *A* under constant pressure *p*.Assuming the error in strain energy calculation approaches zero if

obtained from Figure 5-b. Figure 5a shows the error calculation in strain energy. The displace‐ ment discontinuity method always overestimates the strain energy (or displacement across the crack surface) but it yields more accurate results closer to the exact solution when the number of elements increases. The error changes from 48.8% using a 3×3 mesh to about 1.99% for a mesh including 71×71 elements. In comparison with the two dimensional analysis of a straight crack [16], the rate of convergence is faster, but the error in strain energy calculation

is higher using the same number of elements to divide one side of a crack.

**Figure 4.** Approximation of the exact solution of strain energy for a square pressurized crack

*<sup>n</sup>* and has an asymptotic behavior with respect to *n*, where *n* is the number of

*<sup>n</sup>* →0), the correct answer for error estimation in the strain energy calculation can be

*<sup>σ</sup><sup>n</sup> <sup>π</sup><sup>b</sup>* (5)

*FI* <sup>=</sup> *KI* (*x*, *<sup>y</sup>*)|*<sup>x</sup>* <sup>=</sup> *<sup>x</sup>*, *<sup>y</sup>* <sup>=</sup> <sup>±</sup> *<sup>b</sup>*

the dimensionless stress intensity factor along the crack front *y* =*b*:

about 1%.

748 Effective and Sustainable Hydraulic Fracturing

varies with <sup>1</sup>

*<sup>n</sup>* <sup>→</sup>*<sup>∞</sup>* ( <sup>1</sup>

**Figure 5.** Error in strain energy calculation as a function of number of elements on each side of a pressurized square crack (*a*=*b*)

The error in strain energy calculation is mainly related to the largest error occurring at the corners of the square crack where the displacement gradient is highest. Figure 6 shows the stress intensity factor variation along the half-length of the crack tip using DDM compared with the integral equation solution suggested by [40]. The total number of elements used in the simulation was 22×22 to be consistent with the number of colocation points used in [40]. The difference between these two solutions is negligible for all elements but the corners (element No. 11). However, the corner elements of rectangular cracks don't play an important role in fracture propagation problems because the level of SIF is the lowest there and unlikely to control the initiation of crack propagation.

subdivision number. Figure 7 shows that the most reliable value of *FI max* for a square crack is 0.7607, which is slightly different (0.6%) than the value reported by [39] using body force

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

Figure 8 shows the variation of dimensionless stress intensity factor, F1, along the crack front

Figure 9 shows the maximum dimensionless stress intensity factor (*FI max*) at the location (*x=0, y* =*b*). When *b/a* <1, the crack tip at *y=b* is represents the longer edge of a rectangular crack, whereas when *b/a*>1 the crack tip at *y=b* represents the shorter tip. The dimensionless SIF is referenced to the plane strain SIF for a crack with half-length *b* for all *b/a*. The results show that at *b/a*=0.125, the maximum SIF (at location *x=0, y=b*) has reached the plane strain value (*F1=1*). As *b/a* increases (equivalent to reducing the crack length *a* relative to *b*), *F1* is reduced. When *b/a*=1.0, the square crack, *F1*=0.75. A penny-shaped crack has more restricted opening, and has the ratio of 0.64 to the plane strain SIF. Reducing *a* further such that *b/a*>1 makes *a* the short dimension of the crack and thus the limiting dimension for crack opening and SIF value. The SIF at y=b will then go to 0 as *a*→0. In comparing to the solution of [Wang et al 2001], it is evident that the distribution of SIF near the *x=a* crack tip is more accurate when *b/a* <1, but the maximum value of SIF is a good match for all cases. Using higher element density around the rectangular crack front and a coarser mesh at the center was investigated, but we found a uniform mesh yielded more accurate results using fewer elements in comparison with non-

**Figure 8.** Dimensionless stress intensity factor variation along the crack front *y* =*b*.

*<sup>a</sup>* , using 22×22 elements, a mesh refinement consistent with [40].

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

751

method.

uniform mesh.

y=b for various values of *<sup>b</sup>*

**Figure 6.** Dimensionless stress intensity factor variation along the half length of a square crack front

**Figure 7.** Extrapolation of FI max for a square crack in an infinite body

It is always desirable to use a coarser mesh to save computation time, but the accuracy of DDM depends strongly on mesh refinement. Figure 7 shows the extrapolation of maximum dimen‐ sionless stress intensity factor, *FI max*(which occurs at side-midpoint of a square crack) as a function of <sup>1</sup> *<sup>n</sup>* . It shows the numerical result of *FI max* is parabolic with the reciprocal of the subdivision number. Figure 7 shows that the most reliable value of *FI max* for a square crack is 0.7607, which is slightly different (0.6%) than the value reported by [39] using body force method.

Figure 8 shows the variation of dimensionless stress intensity factor, F1, along the crack front y=b for various values of *<sup>b</sup> <sup>a</sup>* , using 22×22 elements, a mesh refinement consistent with [40]. Figure 9 shows the maximum dimensionless stress intensity factor (*FI max*) at the location (*x=0, y* =*b*). When *b/a* <1, the crack tip at *y=b* is represents the longer edge of a rectangular crack, whereas when *b/a*>1 the crack tip at *y=b* represents the shorter tip. The dimensionless SIF is referenced to the plane strain SIF for a crack with half-length *b* for all *b/a*. The results show that at *b/a*=0.125, the maximum SIF (at location *x=0, y=b*) has reached the plane strain value (*F1=1*). As *b/a* increases (equivalent to reducing the crack length *a* relative to *b*), *F1* is reduced. When *b/a*=1.0, the square crack, *F1*=0.75. A penny-shaped crack has more restricted opening, and has the ratio of 0.64 to the plane strain SIF. Reducing *a* further such that *b/a*>1 makes *a* the short dimension of the crack and thus the limiting dimension for crack opening and SIF value. The SIF at y=b will then go to 0 as *a*→0. In comparing to the solution of [Wang et al 2001], it is evident that the distribution of SIF near the *x=a* crack tip is more accurate when *b/a* <1, but the maximum value of SIF is a good match for all cases. Using higher element density around the rectangular crack front and a coarser mesh at the center was investigated, but we found a uniform mesh yielded more accurate results using fewer elements in comparison with nonuniform mesh.

**Figure 8.** Dimensionless stress intensity factor variation along the crack front *y* =*b*.

**Figure 6.** Dimensionless stress intensity factor variation along the half length of a square crack front

It is always desirable to use a coarser mesh to save computation time, but the accuracy of DDM depends strongly on mesh refinement. Figure 7 shows the extrapolation of maximum dimen‐ sionless stress intensity factor, *FI max*(which occurs at side-midpoint of a square crack) as a

*<sup>n</sup>* . It shows the numerical result of *FI max* is parabolic with the reciprocal of the

**Figure 7.** Extrapolation of FI max for a square crack in an infinite body

750 Effective and Sustainable Hydraulic Fracturing

function of <sup>1</sup>

**Figure 9.** Maximum dimensionless stress intensity factor along the crack front *y* =*b*.

Considering a rectangular vertical crack in a half-space, and assuming *ν* =0.3, the dimension‐ less stress intensity factor at midpoints of crack fronts nearest (*A*1) and farthest (*A*2) from the free surface are presented in Figure 10a and b respectively, as a function of *b* / *a*and *b* / *d*. *F*<sup>1</sup> *max* and *F*<sup>2</sup> *max* are the dimensionless stress intensity factors corresponding to points *A*1and *A*<sup>2</sup> respectively and can be defined as the following:

$$\begin{aligned} F\_{1\,\max} &= \frac{K\_I\left(\mathbf{x},\ \mathbf{y}\right)\,\vert\,A\_1}{\sigma\_\pi\sqrt{\pi b}} \,a\\ F\_{2\,\max} &= \frac{K\_I\left(\mathbf{x},\ \mathbf{y}\right)\,\vert\,A\_2\\ F\_{3\,\max} &= \frac{K\_I\left(\mathbf{x},\ \mathbf{y}\right)\,\vert\,A\_2\\ &\quad \sigma\_\pi\sqrt{\pi b}}{\sigma\_\pi\sqrt{\pi b}} \,b \end{aligned} \tag{6}$$

Figure 10-a and b show for greater aspect ratio (*b* / *a*grater or taller crack) SIF is less affected

**[43]** Mastrojannis EN, Keer LM, Mura TA. Stress intensity factor for a plane crack under normal

**[44]** Olson JE. Sublinear scaling of fracture aperture versus length: An exception or the rule. Journal

**[45]** Mear ME, Rodin GJ. An isolated Mode three-dimensional planar crack: The stress intensity factor is independent to elastic constant. International Journal of Fracture 2011;172(2) 217-218. **[46]** Sneddon, IN. The distribution of stress in the neighborhood of a crack in an elastic solid. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences

**[47]** -Erdogan, F, Sih GC. On the crack extension in plates under plane loading and transverse shear.

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

753

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

**[48]** Pollard DD, Segall P, Delaney PT. Formation and interpretation of dilatant echelon cracks.

**[49]** Atkinson BK. Subcritical crack growth in geological materials. Journal of geophysical research

**[50]** Olson JE. Multi-fracture propagation modeling: Application to hydraulic fracturing in shales and tight gas sands. The 42st U.S. Symposium on Rock Mechanics and 2nd U.S.-canada Rock Mechanics

**[51]** Olson JE. Fracturing from highly deviated and horizontal wells: Numerical analysis of nonplanar fracture propagation. SPE Rock Mountain Regional/Low Permeability Reservoir Symposium,

**[52]** Olson JE, Wu K. Sequential versus simultaneous multi-zone fracturing in horizontal wells: Insight from a non-planar, multi-frac numerical model. SPE Hydraulic Fracturing Technology

**[53]** Cruikshank KM, Zhao G, Johnson AM. Analysis of minor fractures associated with joints and

pressure. International Journal of Fracture 1979;15(3) 101-114.

of Geophysical Research 2003;108(B9) 2413-2424.

ASME Journal of Basic Engineering 1963;85 519-527

Geological Society of America Bulletin 1982;93 1291-1303

Paper No. SPE 29573, March 20-22 1995, Denver, CO, USA.

faulted joints. Journal of Structural Geology 1991;13(8) 865-886.

1946;187(1009) 229-260.

1984;89(B6) 4077-4114

by the depth. Both *F*<sup>1</sup> *max* and *F*<sup>2</sup> *max* increase as the crack approaches the surface of solid. The

mode I stress intensity factor along the crack fronts of a rectangular discontinuity in an infinite

Symposium, Paper No: ARMA 08-327, June 29-July 2 2008, San Francisco, CA, USA.

Conference, Paper No: SPE 152602-PP, February 6-8 2012, The Woodland, TX, USA.

body is independent of Young's modulus [45]. Figure 11a and 11b show that Poisson's ratio *ν*

variation has a slight effect on *F*<sup>1</sup> *max* and *F*<sup>2</sup> *max*, but only for cracks close to the free surface.

**Figure 10.** a Dimensionless stress intensity factor, *F*<sup>1</sup> *max*as a function of b/a and b/d for a rectangular crack in halfspace (ν =0.3); b Dimensionless stress intensity factor *F*<sup>2</sup> *max* as a function of b/a and b/d for a rectangular crack in half-

27

space (ν =0.3)

where *σn* is the normal pressure at the surface of crack. For every combination of *b* / *a*and *b* / *d*, the stress intensity factor along the side nearest to the free surface is greater than the side farthest away.

**[47]** -Erdogan, F, Sih GC. On the crack extension in plates under plane loading and transverse shear.

**[49]** Atkinson BK. Subcritical crack growth in geological materials. Journal of geophysical research

tight gas sands. The 42st U.S. Symposium on Rock Mechanics and 2nd U.S.-canada Rock Mechanics

planar fracture propagation. SPE Rock Mountain Regional/Low Permeability Reservoir Symposium,

**[53]** Cruikshank KM, Zhao G, Johnson AM. Analysis of minor fractures associated with joints and

**[43]** Mastrojannis EN, Keer LM, Mura TA. Stress intensity factor for a plane crack under normal

**[44]** Olson JE. Sublinear scaling of fracture aperture versus length: An exception or the rule. Journal

**[45]** Mear ME, Rodin GJ. An isolated Mode three-dimensional planar crack: The stress intensity factor is independent to elastic constant. International Journal of Fracture 2011;172(2) 217-218. **[46]** Sneddon, IN. The distribution of stress in the neighborhood of a crack in an elastic solid.

pressure. International Journal of Fracture 1979;15(3) 101-114.

Geological Society of America Bulletin 1982;93 1291-1303

Paper No. SPE 29573, March 20-22 1995, Denver, CO, USA.

of Geophysical Research 2003;108(B9) 2413-2424.

Figure 10-a and b show for greater aspect ratio (*b* / *a*grater or taller crack) SIF is less affected ASME Journal of Basic Engineering 1963;85 519-527 **[48]** Pollard DD, Segall P, Delaney PT. Formation and interpretation of dilatant echelon cracks.

by the depth. Both *F*<sup>1</sup> *max* and *F*<sup>2</sup> *max* increase as the crack approaches the surface of solid. The 1984;89(B6) 4077-4114 **[50]** Olson JE. Multi-fracture propagation modeling: Application to hydraulic fracturing in shales and

mode I stress intensity factor along the crack fronts of a rectangular discontinuity in an infinite **[51]** Olson JE. Fracturing from highly deviated and horizontal wells: Numerical analysis of non-

Symposium, Paper No: ARMA 08-327, June 29-July 2 2008, San Francisco, CA, USA.

Conference, Paper No: SPE 152602-PP, February 6-8 2012, The Woodland, TX, USA.

body is independent of Young's modulus [45]. Figure 11a and 11b show that Poisson's ratio *ν* **[52]** Olson JE, Wu K. Sequential versus simultaneous multi-zone fracturing in horizontal wells: Insight from a non-planar, multi-frac numerical model. SPE Hydraulic Fracturing Technology

**Figure 9.** Maximum dimensionless stress intensity factor along the crack front *y* =*b*.

respectively and can be defined as the following:

752 Effective and Sustainable Hydraulic Fracturing

farthest away.

Considering a rectangular vertical crack in a half-space, and assuming *ν* =0.3, the dimension‐ less stress intensity factor at midpoints of crack fronts nearest (*A*1) and farthest (*A*2) from the free surface are presented in Figure 10a and b respectively, as a function of *b* / *a*and *b* / *d*. *F*<sup>1</sup> *max* and *F*<sup>2</sup> *max* are the dimensionless stress intensity factors corresponding to points *A*1and *A*<sup>2</sup>

*a*

*b*

(6)

27

*<sup>F</sup>*<sup>1</sup> *max* <sup>=</sup> *KI* (*x*, *<sup>y</sup>*)|*A*<sup>1</sup> *σ<sup>n</sup> πb*

*<sup>F</sup>*<sup>2</sup> *max* <sup>=</sup> *KI* (*x*, *<sup>y</sup>*)|*A*<sup>2</sup> *σ<sup>n</sup> πb*

where *σn* is the normal pressure at the surface of crack. For every combination of *b* / *a*and *b* / *d*, the stress intensity factor along the side nearest to the free surface is greater than the side

variation has a slight effect on *F*<sup>1</sup> *max* and *F*<sup>2</sup> *max*, but only for cracks close to the free surface. faulted joints. Journal of Structural Geology 1991;13(8) 865-886.

**Figure 10.** a Dimensionless stress intensity factor, *F*<sup>1</sup> *max*as a function of b/a and b/d for a rectangular crack in halfspace (ν =0.3); b Dimensionless stress intensity factor *F*<sup>2</sup> *max* as a function of b/a and b/d for a rectangular crack in halfspace (ν =0.3)

**Figure 12.** Effect of Poisson's ratio on Mode II dimensionless stress intensity factor for a rectangular crack in an infinite

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

For an elliptical crack embedded in an infinite body, the stress intensity factor variation along

sin<sup>2</sup> *<sup>θ</sup>* <sup>+</sup> *<sup>b</sup>* <sup>2</sup>

*E*(*k*)is the complete elliptical integral of the second kind while *a* is the major axis and *b*is the minor axis of ellipse. The maximum and minimum stress intensity factor at the end of minor

> <sup>2</sup> ) <sup>=</sup> *<sup>σ</sup><sup>n</sup> <sup>π</sup><sup>b</sup> <sup>E</sup>*(*<sup>k</sup>* ) *a*

> > *b <sup>a</sup> b*

(*<sup>θ</sup>* <sup>=</sup> *<sup>π</sup>*

(*<sup>θ</sup>* =0)= *<sup>σ</sup><sup>n</sup> <sup>π</sup><sup>b</sup> E*(*k* )

*<sup>a</sup>* <sup>4</sup> cos<sup>2</sup> *<sup>θ</sup>*

) 1 4

(7)

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

755

(9)

*<sup>a</sup>* <sup>2</sup> cos<sup>2</sup> *<sup>θ</sup>*

the crack edge can be obtained from the following analytical solution [36]:

1 2 *<sup>E</sup>*(*<sup>k</sup>* ) ( sin<sup>2</sup> *<sup>θ</sup>* <sup>+</sup> *<sup>b</sup>* <sup>4</sup>

(*θ*)= *<sup>σ</sup>n*(*πb*)

*a* 2

and major axes, respectively, can be calculated using Equations 8a and 8b:

(*<sup>K</sup> <sup>I</sup>* )*max* <sup>=</sup> *KI*

(*<sup>K</sup> <sup>I</sup>* )*min* <sup>=</sup> *KI*

*KI*

*<sup>b</sup>* <sup>2</sup> =1 and,

sin<sup>2</sup> *<sup>θ</sup>*)*d<sup>θ</sup> and <sup>k</sup>* =1 - *<sup>b</sup>* <sup>2</sup>

space.

where:

*<sup>θ</sup>* =tan-1 *<sup>y</sup>*

*E*(*k*)=*∫* 0 *π* <sup>2</sup> (1 - *k* <sup>2</sup>

*<sup>x</sup>* , *<sup>x</sup>* <sup>2</sup> *<sup>a</sup>* <sup>2</sup> <sup>+</sup> *<sup>y</sup>* <sup>2</sup>

**3.2. Elliptical crack**

**Figure 11.** a. Effect of Poisson's ratio on dimensionless stress intensity factor, *F*<sup>1</sup> *max*for a rectangular crack in halfspace; b. Effect of Poisson's ratio on dimensionless stress intensity factor *F*<sup>2</sup> *max* for a rectangular crack in half-space

In contrast to Mode I, for mode II and III stress intensity factor of a crack in an infinite body is dependent on elastic constants. By defining the dimensionless stress intensity factor for mode II, *FII* <sup>=</sup> *KII* (*x*, *<sup>y</sup>*)|*<sup>x</sup>* <sup>=</sup> *<sup>x</sup>*, *<sup>y</sup>* <sup>=</sup> <sup>±</sup> *<sup>b</sup> <sup>τ</sup>zx <sup>π</sup>b* and assuming a frictionless surface crack, Figure 13 shows the maximum dimensionless stress intensity factor along the rectangular crack front *y* =*b*subject to front-perpendicular shear stress *τzx*. The figure shows increasing Poisson's ratio will increase mode II stress intensity factor at the tip of a rectangular crack embedded in an infinite space. Results were satisfactory compared with [38].

**Figure 12.** Effect of Poisson's ratio on Mode II dimensionless stress intensity factor for a rectangular crack in an infinite space.

#### **3.2. Elliptical crack**

For an elliptical crack embedded in an infinite body, the stress intensity factor variation along the crack edge can be obtained from the following analytical solution [36]:

$$\mathcal{K}\_I\{\Theta\} = \frac{\sigma\_n(\pi b)^{\frac{1}{2}}}{E(k)} \left( \frac{\sin^2 \theta + \frac{b^4}{a^4 \cos^2 \theta}}{\sin^2 \theta + \frac{b^2}{a^2 \cos^2 \theta}} \right)^{\frac{1}{4}} \tag{7}$$

where:

(a)

754 Effective and Sustainable Hydraulic Fracturing

(b)

II, *FII* <sup>=</sup> *KII* (*x*, *<sup>y</sup>*)|*<sup>x</sup>* <sup>=</sup> *<sup>x</sup>*, *<sup>y</sup>* <sup>=</sup> <sup>±</sup> *<sup>b</sup>*

space. Results were satisfactory compared with [38].

**Figure 11.** a. Effect of Poisson's ratio on dimensionless stress intensity factor, *F*<sup>1</sup> *max*for a rectangular crack in halfspace; b. Effect of Poisson's ratio on dimensionless stress intensity factor *F*<sup>2</sup> *max* for a rectangular crack in half-space

In contrast to Mode I, for mode II and III stress intensity factor of a crack in an infinite body is dependent on elastic constants. By defining the dimensionless stress intensity factor for mode

maximum dimensionless stress intensity factor along the rectangular crack front *y* =*b*subject to front-perpendicular shear stress *τzx*. The figure shows increasing Poisson's ratio will increase mode II stress intensity factor at the tip of a rectangular crack embedded in an infinite

*<sup>τ</sup>zx <sup>π</sup>b* and assuming a frictionless surface crack, Figure 13 shows the

$$\begin{array}{l} \Theta = \tan^{-1} \frac{y}{x} \; , \; \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \text{ and} \\\\ E\ (k) = \int\_0^{\frac{\pi}{2}} (1 - k^2 \sin^2 \theta) d\theta \; \text{and} \; k = 1 - \frac{b^{\frac{3}{2}}}{a^2} \end{array}$$

*E*(*k*)is the complete elliptical integral of the second kind while *a* is the major axis and *b*is the minor axis of ellipse. The maximum and minimum stress intensity factor at the end of minor and major axes, respectively, can be calculated using Equations 8a and 8b:

$$\begin{aligned} \left(K\_{\,\,I}\right)\_{\max} &= K\_{\,\,I} \left(\Theta = \frac{\pi}{2}\right) = \frac{\sigma\_{\pi} \sqrt{\pi b}}{E(k)} \quad a\\ \left(K\_{\,\,I}\right)\_{\min} &= K\_{\,\,I} \left(\Theta = 0\right) = \frac{\sigma\_{\pi} \sqrt{\pi b}}{E(k)} \sqrt{\frac{b}{a}} \quad b \end{aligned} \tag{9}$$

Figure 13a and b show dimensionless stress intensity factor variation along the elliptical crack front using analytical solutions and DDM numerical modeling. Totally 154 DD elements were used in the model depicted in Figure 13a, and 628 elements in Figure 13b. Whereas SIF is proportional to the area of planar crack, the area of boundary element mesh in both cases is almost equal to the area of the modeled ellipse. For both models, the aspect ratio of the ellipse is *<sup>b</sup> <sup>a</sup>* =2 and *F*<sup>1</sup> <sup>=</sup> *<sup>K</sup>*1(*θ*) *σn*(*πb*) 1 2 . Both figures show that the trend of stress

intensity factor variation can be appropriately modeled by DDM. Oscillation in SIF is because of stepwise mesh boundary used to define the geometry of the ellipse using rectangular elements. However, by using the average of SIF of the neighboring circumfer‐ ential elements, the accuracy improves for both models and the maximum error decreas‐ es from about 24% to 9% for the first model and from 28% to 10% for the second model, as compared to the analytical solution in [36]. Using 20 elements along the major axis and <sup>10</sup> along the minor axis of the ellipse results in good agreement for F1 at *<sup>θ</sup>* =0 and *<sup>π</sup>* <sup>2</sup> (Figure 13-a). For *θ* ≥60° , the rectangular mesh deviates less from the ellipse, and the error in dimensionless stress intensity factor is non-oscillatory and small. Increasing the number of elements doesn't improve the accuracy (Figure 13-b).

#### **3.3. Penny-shaped crack**

Stress intensity factor at the tip of a circular crack of radius *a*in an infinite solid under uniaxial tension *σn* is [46]:

$$\frac{2}{\pi}\sigma\_n\sqrt{\pi a} \tag{10}$$

the crack front is stepwise. Similar to elliptical cracks, using the average SIF of neighbor circumferential elements considerably increases the accuracy of SIF distribution along the

model No. 1 including 154 elements; b. Dimensionless SIF variation along an elliptical crack front using analytical solu‐

*<sup>a</sup>* =2.0),

**Figure 13.** a. Dimensionless SIF variation along an elliptical crack front using analytical solution and DDM ( *<sup>b</sup>*

(b)

(a)

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

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

757

crack front of the penny-shaped discontinuity.

*<sup>a</sup>* =2.0), model No. 2 including 628 elements

17

tion and DDM ( *<sup>b</sup>*

where:

$$\theta = \tan^{-1} \frac{y}{x} \text{ , } \propto x^2 + y^2 = a^2$$

Two different size meshes were considered to calculate dimensionless stress intensity factor variation along the tip of a circular crack as depicted in Figures 14a and 14b. The first model includes 76 elements and the second one has 308 elements. According to Figure 7, for a rectangular crack using 9×9 elements, the error in stress intensity factor is about 3 percent. For the penny-shaped crack, as with the elliptical crack, the error is a strong function of location. Because of the symmetry, error calculations are shown only for one eighth of the circle. The main reason of error in stress intensity factor along the crack front is jagged geometrical definition of the circle by using rectangular displacement discontinu‐ ity elements. The error in SIF can reach up to 20% along the crack front; however, the results are better for *<sup>θ</sup>* =0 *or <sup>π</sup>* <sup>2</sup> - about 2.5% for the coarser model and almost zero for the finer model. Figure 15 compares *F*1variation along the quarter front of the penny-shaped crack for two DD models as well as analytical solution. The figure shows the finer mesh helps to increase the accuracy where the crack front is straight, but is not helpful where Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method... http://dx.doi.org/10.5772/56308 757

Figure 13a and b show dimensionless stress intensity factor variation along the elliptical crack front using analytical solutions and DDM numerical modeling. Totally 154 DD elements were used in the model depicted in Figure 13a, and 628 elements in Figure 13b. Whereas SIF is proportional to the area of planar crack, the area of boundary element mesh in both cases is almost equal to the area of the modeled ellipse. For both models, the aspect

intensity factor variation can be appropriately modeled by DDM. Oscillation in SIF is because of stepwise mesh boundary used to define the geometry of the ellipse using rectangular elements. However, by using the average of SIF of the neighboring circumfer‐ ential elements, the accuracy improves for both models and the maximum error decreas‐ es from about 24% to 9% for the first model and from 28% to 10% for the second model, as compared to the analytical solution in [36]. Using 20 elements along the major axis and

dimensionless stress intensity factor is non-oscillatory and small. Increasing the number of

Stress intensity factor at the tip of a circular crack of radius *a*in an infinite solid under uniaxial

Two different size meshes were considered to calculate dimensionless stress intensity factor variation along the tip of a circular crack as depicted in Figures 14a and 14b. The first model includes 76 elements and the second one has 308 elements. According to Figure 7, for a rectangular crack using 9×9 elements, the error in stress intensity factor is about 3 percent. For the penny-shaped crack, as with the elliptical crack, the error is a strong function of location. Because of the symmetry, error calculations are shown only for one eighth of the circle. The main reason of error in stress intensity factor along the crack front is jagged geometrical definition of the circle by using rectangular displacement discontinu‐ ity elements. The error in SIF can reach up to 20% along the crack front; however, the

finer model. Figure 15 compares *F*1variation along the quarter front of the penny-shaped crack for two DD models as well as analytical solution. The figure shows the finer mesh helps to increase the accuracy where the crack front is straight, but is not helpful where

2

, the rectangular mesh deviates less from the ellipse, and the error in

. Both figures show that the trend of stress

*<sup>π</sup> σ<sup>n</sup> πa* (10)

<sup>2</sup> - about 2.5% for the coarser model and almost zero for the

<sup>2</sup> (Figure

*<sup>a</sup>* =2 and *F*<sup>1</sup> <sup>=</sup> *<sup>K</sup>*1(*θ*)

elements doesn't improve the accuracy (Figure 13-b).

*σn*(*πb*) 1 2

<sup>10</sup> along the minor axis of the ellipse results in good agreement for F1 at *<sup>θ</sup>* =0 and *<sup>π</sup>*

ratio of the ellipse is *<sup>b</sup>*

756 Effective and Sustainable Hydraulic Fracturing

13-a). For *θ* ≥60°

tension *σn* is [46]:

where:

*<sup>θ</sup>* =tan-1 *<sup>y</sup>*

**3.3. Penny-shaped crack**

*<sup>x</sup>* , *<sup>x</sup>* <sup>2</sup> <sup>+</sup> *<sup>y</sup>* <sup>2</sup> <sup>=</sup>*<sup>a</sup>* <sup>2</sup>

results are better for *<sup>θ</sup>* =0 *or <sup>π</sup>*

**Figure 13.** a. Dimensionless SIF variation along an elliptical crack front using analytical solution and DDM ( *<sup>b</sup> <sup>a</sup>* =2.0), model No. 1 including 154 elements; b. Dimensionless SIF variation along an elliptical crack front using analytical solu‐ tion and DDM ( *<sup>b</sup> <sup>a</sup>* =2.0), model No. 2 including 628 elements

the crack front is stepwise. Similar to elliptical cracks, using the average SIF of neighbor circumferential elements considerably increases the accuracy of SIF distribution along the crack front of the penny-shaped discontinuity.

17

0.5

0 20 40 60 80

**q(degree)**

**Figure 15.** Comparison between dimensionless SIF for two DDM models with analytical solution of a penny-shaped

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

For vertical fractures, lateral kinking propagation is modeled based on maximum circumfer‐ ential stress criteria [47], which states growth should occur at the crack tip along a radial path

Model#1

Model#2

Analytical

Average Model # 1

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

759

Average, Model #2

0.55

0.6

**F1**

crack stress intensity factor

**4. Fracture propagation**

perpendicular to the direction of greatest tension:

0.65

0.7

**Figure 14.** a Error in dimensionless calculation along a penny-shaped crack front, Model 1 including 76 elements; b. Error in dimensionless calculation along a penny-shaped crack front, Model 2 containing 308 elements

Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method... http://dx.doi.org/10.5772/56308 759

**Figure 15.** Comparison between dimensionless SIF for two DDM models with analytical solution of a penny-shaped crack stress intensity factor

#### **4. Fracture propagation**

(a)

758 Effective and Sustainable Hydraulic Fracturing

(b)

**Figure 14.** a Error in dimensionless calculation along a penny-shaped crack front, Model 1 including 76 elements; b.

Error in dimensionless calculation along a penny-shaped crack front, Model 2 containing 308 elements

For vertical fractures, lateral kinking propagation is modeled based on maximum circumfer‐ ential stress criteria [47], which states growth should occur at the crack tip along a radial path perpendicular to the direction of greatest tension:

$$\begin{aligned} \tan\frac{\theta\_0}{2} &= \frac{1}{4} \left[ \frac{K\_I}{K\_{II}} - \text{Sgn}\left(K\_{II}\right) \sqrt{\left(\frac{K\_I}{K\_{II}}\right)^2 + 8} \right] \text{ a} \\ K\_{eq} &= K\_I \cos^3\frac{\theta\_0}{2} - \frac{3}{2} K\_{II} \cos\frac{\theta\_0}{2} \sin\frac{\theta\_0}{2} \quad \text{b} \end{aligned} \tag{11}$$

where *θ*<sup>0</sup> is the angle of kinking and *Sgn*(*KII* )denotes the sign of *K II* . Equation 10-a is used to calculate the equivalent opening mode stress intensity factor in the direction of crack extension (*Keq*) as in formula (10b).

Our model takes into account the height growth as pure Mode I propagation. Any contribution of Mode III or out of plane shear on vertical propagation is neglected; however, the possibility of fringe crack generation based on Mode I+III combination will be studied by Mode III SIF evaluation along the upper front of the fracture. The angle of twisting is dependent on the magnitude of Mode III and Mode I SIF as well as mechanical properties [48] and can be calculated using Equation 11. Higher values of Mode III SIF (or lower opening mode) result in bigger twisting angle.

$$\alpha = \frac{1}{2} \tan^{-1} \left[ \frac{\mathbb{K}\_{\text{II}}}{\mathbb{K}\_{\text{I}} \left( \big| \begin{array}{c} \cdot \\ \cdot \\ \end{array} \right)} \right] \tag{12}$$

**Figure 16.** Ideal longitudinal fractured horizontal well with hydraulic fracture perpendicular to σ*hmin*.

stresses), *P frac*is constant and equal to 20.0 *MPa*, the remote compression differential stress is

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

initial fracture length and height are assumed to be 3 meters (a square crack), subdivided by 9 DDM elements. The fracture is assumed to remain rectangular during the propagation (i.e., the height is uniform along the entire length, but the crack path in plan-view can be non-

To examine the effect of horizontal well misalignment angle on fracture propagation (Figure

misalignment angle, *β*, especially for extreme cases. The starter fracture is centered at (0,0) and is rotated counterclockwise by *β*. The smallest misalignment *β* =10° is the closest to planar

Nonplanar propagation has an impact on height growth (Figures 18 and 19). For the smaller misalignment cases (*β* <=45º), crack height keeps pace with crack length growth for our imposed rectangular shape (Figure 18). For our stronger misalignment cases of *β* >45º, the crack height growth is somewhat hindered to only ~80% of the length. Looking at the opening mode

) distribution along the top edge of the fracture is more interesting, however, since our propagation algorithm responds only to the average crack tip SIF. The more severe the fracture reorientation, the lower the *KI* for the initial fracture segment, where for the 89º misalignment case, the *KI* at the center of the crack is 50% lower than it would be for a planar fracture. This implies that at the wellbore, there could be a restriction in fracture height because of the non-

*<sup>r</sup>* ) =2.0 *MPa*, the propagation velocity exponent is *n* =2, *ν* =0.25and *E* =30.0 *GPa*. The

*<sup>r</sup>* = 15 MPa (where *r*denotes remote

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

761

*<sup>r</sup>* - *Sxx*

*<sup>r</sup>* ). Fracture path non-planarity is strongly affected by the initial

*<sup>r</sup>* ) is 40*%* of the net

For the propagation cases that follow, we assume *Shmin* =*Sxx*

17), first we assume the differential compression in *y*direction (*S yy*

planar propagation that might also restrict width and hinder infectivity.

(*σ yy <sup>r</sup>* - *<sup>σ</sup>xx*

planar).

SIF (*KI*

injection pressure (*P frac* - *Sxx*

fracture and *β* =89° is the most curved path.

Fracture front propagation velocity defines which edge extends first. Charles power law [49] was used to relate the equivalent opening Mode stress intensity factor at the tip of the crack to the propagation velocity as the following [49]:

$$V = AK\_{eq}^{n} \tag{13}$$

#### **5. Application: Fracture misalignment and height growth**

Figure 16 shows the ideal alignment of horizontal well and longitudinal hydraulic fracture system where the horizontal well is perpendicular to the minimum remote horizontal stress *Shmin* =*S*3 and the wellbore lies in the principal remote stress plane, parallel to *SHmax* =*S*2. However, hydraulic fractures may not necessarily start perpendicular to the minimum horizontal remote stress because of the lack of alignment between the wellbore and the principal stresses, local stress perturbation, or natural fracture adjacent to a horizontal well [50]. The geometry of a hydraulic fracture could be further complicated by lateral propagation which is non-planar and height growth that is non-uniform. The non-planarity of the fracture path and its resultant near-wellbore width restriction and excessive treating pressure were considered by [51] and [52] using 2-D and pseudo 3-D displacement discontinuity modeling, respectively. In this paper, we study the effect of misalignment angle on the possibility of irregular height growth as well as fringe fracture generation by contemplating the stress intensity factor distribution around the periphery of misaligned hydraulic fracture. Wellbore stress effects are not considered in this study.

**Figure 16.** Ideal longitudinal fractured horizontal well with hydraulic fracture perpendicular to σ*hmin*.

tan *θ*0 <sup>2</sup> <sup>=</sup> <sup>1</sup> 4 *KI*

(*Keq*) as in formula (10b).

760 Effective and Sustainable Hydraulic Fracturing

in bigger twisting angle.

*Keq* <sup>=</sup> *KI* cos<sup>3</sup> *<sup>θ</sup>*<sup>0</sup>

*<sup>α</sup>* <sup>=</sup> <sup>1</sup>

**5. Application: Fracture misalignment and height growth**

to the propagation velocity as the following [49]:

stress effects are not considered in this study.

<sup>2</sup> tan-1 *KIII KI* ( 1 2

*V* = *AKeq*

Fracture front propagation velocity defines which edge extends first. Charles power law [49] was used to relate the equivalent opening Mode stress intensity factor at the tip of the crack

Figure 16 shows the ideal alignment of horizontal well and longitudinal hydraulic fracture system where the horizontal well is perpendicular to the minimum remote horizontal stress *Shmin* =*S*3 and the wellbore lies in the principal remote stress plane, parallel to *SHmax* =*S*2. However, hydraulic fractures may not necessarily start perpendicular to the minimum horizontal remote stress because of the lack of alignment between the wellbore and the principal stresses, local stress perturbation, or natural fracture adjacent to a horizontal well [50]. The geometry of a hydraulic fracture could be further complicated by lateral propagation which is non-planar and height growth that is non-uniform. The non-planarity of the fracture path and its resultant near-wellbore width restriction and excessive treating pressure were considered by [51] and [52] using 2-D and pseudo 3-D displacement discontinuity modeling, respectively. In this paper, we study the effect of misalignment angle on the possibility of irregular height growth as well as fringe fracture generation by contemplating the stress intensity factor distribution around the periphery of misaligned hydraulic fracture. Wellbore


*<sup>n</sup>* (13)

*KII* - *Sgn*(*KII*

<sup>2</sup> - <sup>3</sup>

<sup>2</sup> *KII* cos

where *θ*<sup>0</sup> is the angle of kinking and *Sgn*(*KII* )denotes the sign of *K II* . Equation 10-a is used to calculate the equivalent opening mode stress intensity factor in the direction of crack extension

Our model takes into account the height growth as pure Mode I propagation. Any contribution of Mode III or out of plane shear on vertical propagation is neglected; however, the possibility of fringe crack generation based on Mode I+III combination will be studied by Mode III SIF evaluation along the upper front of the fracture. The angle of twisting is dependent on the magnitude of Mode III and Mode I SIF as well as mechanical properties [48] and can be calculated using Equation 11. Higher values of Mode III SIF (or lower opening mode) result

) ( *KI KII* ) 2 + 8 *a*

> *θ*0 <sup>2</sup> sin *θ*0 <sup>2</sup> *b*

(11)

For the propagation cases that follow, we assume *Shmin* =*Sxx <sup>r</sup>* = 15 MPa (where *r*denotes remote stresses), *P frac*is constant and equal to 20.0 *MPa*, the remote compression differential stress is (*σ yy <sup>r</sup>* - *<sup>σ</sup>xx <sup>r</sup>* ) =2.0 *MPa*, the propagation velocity exponent is *n* =2, *ν* =0.25and *E* =30.0 *GPa*. The initial fracture length and height are assumed to be 3 meters (a square crack), subdivided by 9 DDM elements. The fracture is assumed to remain rectangular during the propagation (i.e., the height is uniform along the entire length, but the crack path in plan-view can be nonplanar).

To examine the effect of horizontal well misalignment angle on fracture propagation (Figure 17), first we assume the differential compression in *y*direction (*S yy <sup>r</sup>* - *Sxx <sup>r</sup>* ) is 40*%* of the net injection pressure (*P frac* - *Sxx <sup>r</sup>* ). Fracture path non-planarity is strongly affected by the initial misalignment angle, *β*, especially for extreme cases. The starter fracture is centered at (0,0) and is rotated counterclockwise by *β*. The smallest misalignment *β* =10° is the closest to planar fracture and *β* =89° is the most curved path.

Nonplanar propagation has an impact on height growth (Figures 18 and 19). For the smaller misalignment cases (*β* <=45º), crack height keeps pace with crack length growth for our imposed rectangular shape (Figure 18). For our stronger misalignment cases of *β* >45º, the crack height growth is somewhat hindered to only ~80% of the length. Looking at the opening mode SIF (*KI* ) distribution along the top edge of the fracture is more interesting, however, since our propagation algorithm responds only to the average crack tip SIF. The more severe the fracture reorientation, the lower the *KI* for the initial fracture segment, where for the 89º misalignment case, the *KI* at the center of the crack is 50% lower than it would be for a planar fracture. This implies that at the wellbore, there could be a restriction in fracture height because of the nonplanar propagation that might also restrict width and hinder infectivity.

**Figure 19.** *KI*

**Figure 20.** *KI*

for the case *β* =80° . The *KI*

along the upper front of hydraulic fracture implying height growth restriction around the wellbore due

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

distribution variation normalized by the absolute maximum opening mode SIF during propagation along

the upper front of a hydraulic fracture perforated from a misaligned horizontal wellbore. Misalignment angle, β =80° .

at the initial fracture location (the injection location) grows very

The time progression of the *KI* variation along the top fracture front is displayed in Figure 20

<sup>2</sup> =25.0 *m*)

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

763

to misalignment normalized to SIF of planar fracture at *x* =0(upper front propagation, <sup>∆</sup> *<sup>H</sup>*

slowly in comparison to the curving wings of the fracture.

**Figure 17.** Map view of non-planar fracture paths (upper front propagation, <sup>∆</sup> *<sup>H</sup>* <sup>2</sup> =25.0 *m*)

**Figure 18.** Vertical versus lateral growth of the hydraulic fracture

Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method... http://dx.doi.org/10.5772/56308 763

**Figure 19.** *KI* along the upper front of hydraulic fracture implying height growth restriction around the wellbore due to misalignment normalized to SIF of planar fracture at *x* =0(upper front propagation, <sup>∆</sup> *<sup>H</sup>* <sup>2</sup> =25.0 *m*)

The time progression of the *KI* variation along the top fracture front is displayed in Figure 20 for the case *β* =80° . The *KI* at the initial fracture location (the injection location) grows very slowly in comparison to the curving wings of the fracture.

**Figure 17.** Map view of non-planar fracture paths (upper front propagation, <sup>∆</sup> *<sup>H</sup>*

762 Effective and Sustainable Hydraulic Fracturing

**Figure 18.** Vertical versus lateral growth of the hydraulic fracture

<sup>2</sup> =25.0 *m*)

**Figure 20.** *KI* distribution variation normalized by the absolute maximum opening mode SIF during propagation along the upper front of a hydraulic fracture perforated from a misaligned horizontal wellbore. Misalignment angle, β =80° .

Although *KI* is restricted in the misaligned portion of the fracture, Mode III or out of plane shear SIF(*KIII*) is accentuated. This twisting SIF could cause the fracture to break down into several en echelon cracks, causing further propagation hindrance in the vertical direction. Figure 21 depicts the distribution of *KIII* for varying fracture misalignment based on the simulation of Figure 19.

**Figure 21.** Mode III SIF along the upper front of hydraulic fracture normalized to SIF of planar fracture at *x* =0, imply‐ ing height growth restriction around the wellbore due to misalignment (upper front propagation, <sup>∆</sup> *<sup>H</sup>* <sup>2</sup> =25.0 *m*)

Fracture path is affected by remote stresses as well as near-tip stress distribution and is quantifies by ratio *R*[53] assuming compression is positive:

$$\mathbf{R} = \frac{\begin{pmatrix} \boldsymbol{\sigma}\_{Hmax} \ \cdot \ \boldsymbol{\sigma}\_{hmin} \end{pmatrix}}{\begin{Bmatrix} \boldsymbol{P}\_{fuc} \ \cdot \ \boldsymbol{\sigma}\_{hmin} \end{Bmatrix}} = \frac{\begin{pmatrix} \boldsymbol{\sigma}\_{yy}^{r} \ \cdot \ \boldsymbol{\sigma}\_{xx}^{r} \end{pmatrix}}{\begin{pmatrix} \boldsymbol{P}\_{fuc} \ \cdot \ \boldsymbol{\sigma}\_{xx}^{r} \end{pmatrix}} \tag{14}$$

**Figure 22.** R ratio effect on fracture path. Upper front propagation, <sup>∆</sup> *<sup>H</sup>*

SIF values can be accurately computed.

**Acknowledgements**

Numerical methods are necessary for the SIF evaluation of 3-D planar cracks because analytical solutions are limited to simple geometries with special boundary conditions. In this paper, the capability of DDM using constant rectangular discontinuity elements and considering the empirical constant proposed by Olson (1991) was satisfactory examined for cracks with simple geometry. The accuracy of the model is excellent especially for rectangular and square shaped cracks. The stepwise shape of the mesh boundary when representing elliptical or pennyshaped cracks introduces more error in to the calculation, but the minimum and maximum

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

Funding for this project is provided by RPSEA through the "Ultra-Deepwater and Unconven‐ tional Natural Gas and Other Petroleum Resources" program authorized by the U.S. Energy Policy Act of 2005. RPSEA (www.rpsea.org) is a nonprofit corporation whose mission is to provide a stewardship role in ensuring the focused research, development and deployment of safe and environmentally responsible technology that can effectively deliver hydrocarbons from domestic resources to the citizens of the United States. RPSEA, operating as a consortium of premier U.S. energy research universities, industry, and independent research organiza‐

**6. Conclusion**

<sup>2</sup> =25.0 *m*, β =80°and (σ *yy*

*<sup>r</sup>* - <sup>σ</sup>*xx*

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

765

*<sup>r</sup>* )=2.0 *MPa*.

The magnitude of *R*shows how fast the misaligned fracture will be aligned with maximum horizontal stress. Figure 20 present the bigger the magnitude of *R*ratio, the faster the fracture will be rotated to be aligned perpendicular to minimum horizontal stress. Because the differential remote stress is kept constant for these 3 cases, smaller magnitude of ratio *R*means the dominance of fracture driving stresses results in a straighter fracture path.

Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method... http://dx.doi.org/10.5772/56308 765

**Figure 22.** R ratio effect on fracture path. Upper front propagation, <sup>∆</sup> *<sup>H</sup>* <sup>2</sup> =25.0 *m*, β =80°and (σ *yy <sup>r</sup>* - <sup>σ</sup>*xx <sup>r</sup>* )=2.0 *MPa*.

#### **6. Conclusion**

Although *KI*

simulation of Figure 19.

764 Effective and Sustainable Hydraulic Fracturing

is restricted in the misaligned portion of the fracture, Mode III or out of plane

shear SIF(*KIII*) is accentuated. This twisting SIF could cause the fracture to break down into several en echelon cracks, causing further propagation hindrance in the vertical direction. Figure 21 depicts the distribution of *KIII* for varying fracture misalignment based on the

**Figure 21.** Mode III SIF along the upper front of hydraulic fracture normalized to SIF of planar fracture at *x* =0, imply‐

Fracture path is affected by remote stresses as well as near-tip stress distribution and is

The magnitude of *R*shows how fast the misaligned fracture will be aligned with maximum horizontal stress. Figure 20 present the bigger the magnitude of *R*ratio, the faster the fracture will be rotated to be aligned perpendicular to minimum horizontal stress. Because the differential remote stress is kept constant for these 3 cases, smaller magnitude of ratio *R*means

*<sup>r</sup>* - *<sup>σ</sup>xx r* ) (*P frac* - *σxx*

<sup>2</sup> =25.0 *m*)

*<sup>r</sup>* ) (14)

ing height growth restriction around the wellbore due to misalignment (upper front propagation, <sup>∆</sup> *<sup>H</sup>*

*<sup>R</sup>* <sup>=</sup> (*σHmax* - *<sup>σ</sup>hmin*)

the dominance of fracture driving stresses results in a straighter fracture path.

(*<sup>P</sup> frac* - *<sup>σ</sup>hmin*) <sup>=</sup> (*<sup>σ</sup> yy*

quantifies by ratio *R*[53] assuming compression is positive:

Numerical methods are necessary for the SIF evaluation of 3-D planar cracks because analytical solutions are limited to simple geometries with special boundary conditions. In this paper, the capability of DDM using constant rectangular discontinuity elements and considering the empirical constant proposed by Olson (1991) was satisfactory examined for cracks with simple geometry. The accuracy of the model is excellent especially for rectangular and square shaped cracks. The stepwise shape of the mesh boundary when representing elliptical or pennyshaped cracks introduces more error in to the calculation, but the minimum and maximum SIF values can be accurately computed.

#### **Acknowledgements**

Funding for this project is provided by RPSEA through the "Ultra-Deepwater and Unconven‐ tional Natural Gas and Other Petroleum Resources" program authorized by the U.S. Energy Policy Act of 2005. RPSEA (www.rpsea.org) is a nonprofit corporation whose mission is to provide a stewardship role in ensuring the focused research, development and deployment of safe and environmentally responsible technology that can effectively deliver hydrocarbons from domestic resources to the citizens of the United States. RPSEA, operating as a consortium of premier U.S. energy research universities, industry, and independent research organiza‐ tions, manages the program under a contract with the U.S. Department of Energy's National Energy Technology Laboratory. Authors gratefully appreciate RPSEA for providing the funding for this work.

[10] Tarancon, J. E, Vercher, A, Giner, E, & Fuenmayor, F. J. Enhanced blending elements for XFEM applied to linear elastic fracture mechanics. International Journal for Nu‐

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

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

767

[11] Jiang, S, Ying, Z, & Du, C. The optimal XFEM approximation for fracture analysis.

[12] Lin, X, & Ballmann, J. Re-consideration of Chen's problem by finite difference meth‐

[13] Chen, Y. M. Numerical computation of dynamic stress intensity factors by a Lagran‐ gian finite-deference method (THE HEMP CODE). Engineering Fracture Mechanics

[14] Pande, G. N, Beer, G, & Williams, J. R. Numerical Methods in Rock Mechanics. West

[15] Rizzo, F. J. An integral equation approach to boundary value problems of classical

[16] Crouch, S. L. Solution of plane elasticity problems by the displacement discontinuity method. International Journal for Numerical Methods in Engineering (1976). , 10(2),

[17] Crouch, S. L, & Starfield, A. M. Boundary element methods in solid mechanics. Lon‐

[18] Schultz, R. A. Stress intensity factor for curved cracks obtained displacement discon‐

[19] Aydin, A, & Schultz, R. A. Effect of mechanical interaction on the development of strike-slip faults with echelon patterns. Journal of Structural Geology (1990). , 12(1),

[20] Crawford, A. M, & Curran, J. H. Higher-order functional variation displacement dis‐ continuity elements. International Journal of Rock Mechanics and Mining Sciences &

[21] Scavia, C. The displacement discontinuity method on the analysis of open cracks.

[22] Yan, X. Stress intensity factors for cracks emanating from a triangular or square hole in an infinite plate by boundary elements. Engineering Failure Analysis (2005). ,

[23] Shou, K. J, & Crouch, S. L. A higher order displacement discontinuity method for analysis of crack problems. International Journal of Rock Mechanics and Mining Sci‐

merical Methods in Engineering (2009). , 77(1), 126-148.

IOP Conference Series: Materials Science and Engineering (2010).

od. Engineering Fracture Mechanics (1993). , 44(5), 735-739.

elastostatics. Quarterly of Applied Mathematics (1967). , 25-83.

tinuity method. International Journal of Fracture (1988). RR34., 31.

Sussex: John Wiley & Sons Ltd; (1990).

don: George Allen & Unwin; (1983).

Geomechanics Abstract (1982). , 19(3), 143-148.

ences & Geomechanics Abstract (1995). , 32(1), 49-55.

Meccanica (1991). , 26(1), 27-32.

(1975).

301-343.

123-129.

12(3), 362-375.
