**1. Introduction**

As new energy sources are sought for economic and security reasons, unconventional reservoirs attracted the oil and gas industry's attention. Among the unconventional options, shale gas reservoirs have become conspicuous. It is generally accepted that horizontal drilling and hydraulic fracturing are required to effectively recover hydrocarbons from the shale reservoirs [11]. Creating complex fracture networks by hydraulic fracturing is one of the most efficient ways to produce hydrocarbons from these reservoirs due to very low effective permeability (~500 nano Darcy). However, the numerical modeling of hydraulic fractures in such low permeable reservoirs presents significant challenges in field applications [1].

There remains a need for fast, yet accurate, models that remain consistent to the underlying physics of the problem. For numerical simulations, several researchers have considered a number of issues: the coupling between fracture mechanics and fluid dynamics in the fracture [2], fracture interaction [3-5], proppant transport [6], and others [7-9]. Further, there have been specific codes developed to model complex fracture network development [14-16]. Neverthe‐ less, the available literature within the oil and gas industry often ignores the significance of the crack tip in modeling applications developed for hydraulic fracture design.

The importance of accurate modeling of the stress induced near the crack tip is likely critical in complex geological reservoirs. Multiple propagating crack tips interact with each other along with natural fractures, discontinuities, etc., during stimulation treatments in these reservoirs. Consequently, accurate modeling of the stress ahead of the propagating fracture is required to predict fracture paths in this complex environment. This study investigates the influence of several boundary element numerical techniques, available in the literature [10,12,13], on the stress intensity factor near the crack tip and on the fracture profile, in general. This work is a part of a long-term project in the development of more accurate and efficient numerical simulations for field engineering applications.

shown [1,4,7], the accuracy of tip asymptote is critical in characterizing the stresses in the

The authors examined the numerically derived stress intensity factor for several crack geo‐ metries with and without higher-order elements and with and without special tip elements, to analytical solutions. As expected, they found that the cases with higher-order elements and special tip elements provide more accurate results than the cases with constant displace‐ ment discontinuity and/or no tip elements. However, the numerical technique developed

These results show that numerical simulators can incorporate proper crack-tip treatments ef‐ fectively. In addition, higher-order elements increase computational efficiency by reducing the number of elements instead of simply increasing the discretization of constant displace‐ ment elements. The accurate modeling of stress intensity factors is necessary to better simu‐

**Keywords** displacement discontinuity method, higher order elements, crack tip elements

As new energy sources are sought for economic and security reasons, unconventional reservoirs attracted the oil and gas industry's attention. Among the unconventional options, shale gas reservoirs have become conspicuous. It is generally accepted that horizontal drilling and hydraulic fracturing are required to effectively recover hydrocarbons from the shale reservoirs [11]. Creating complex fracture networks by hydraulic fracturing is one of the most efficient ways to produce hydrocarbons from these reservoirs due to very low effective permeability (~500 nano Darcy). However, the numerical modeling of hydraulic fractures in such low permeable reservoirs presents significant challenges in field applications [1].

There remains a need for fast, yet accurate, models that remain consistent to the underlying physics of the problem. For numerical simulations, several researchers have considered a number of issues: the coupling between fracture mechanics and fluid dynamics in the fracture [2], fracture interaction [3-5], proppant transport [6], and others [7-9]. Further, there have been specific codes developed to model complex fracture network development [14-16]. Neverthe‐ less, the available literature within the oil and gas industry often ignores the significance of

The importance of accurate modeling of the stress induced near the crack tip is likely critical in complex geological reservoirs. Multiple propagating crack tips interact with each other along with natural fractures, discontinuities, etc., during stimulation treatments in these reservoirs. Consequently, accurate modeling of the stress ahead of the propagating fracture is required to predict fracture paths in this complex environment. This study investigates the influence of several boundary element numerical techniques, available in the literature [10,12,13], on the stress intensity factor near the crack tip and on the fracture profile, in general.

the crack tip in modeling applications developed for hydraulic fracture design.

late curved fractures, kinked cracks and interaction between fractures.

near-tip region of a propagating fracture.

832 Effective and Sustainable Hydraulic Fracturing

still proved efficient.

**1. Introduction**

To perform the investigation, we used the displacement discontinuity method (DDM), a version of the boundary element method (BEM). The method was developed for, and has been successfully applied to rock mechanics area such as mining engineering [17,18], fracture analysis [19,20], and wellbore stabilities [12]. We have applied DDM here using both constant and higher-order elements. The higher-order elements use a quadratic variation of displace‐ ment discontinuity, and are based on the use of three collocation points over a three-element patch centered at the source element [10], while the constant elements use a constant variation of displacement discontinuity [12]. Details related to the elements are elaborated on in Shou's work [12]. Further, the authors also applied special crack tip elements [10] to capture the square root displacement variation at the crack tip, expected from Linear Elastic Fracture Mechanics (LEFM). The authors expect that special crack tip elements will provide the necessary flexibility to choose other tip profiles. This flexibility will be instrumental for efficient modeling of the different near-tip displacement profiles exhibited by various regimes in hydraulic fracture propagation (e.g., Viscosity-Dominated or Toughness-Dominated [22,23]). As others have shown, the accuracy of tip asymptote is critical in characterizing the stresses in the near-tip region of a propagating fracture [1].

We examined the numerically derived stress intensity factor for three crack geometries with and without higher-order elements and with and without special tip elements, to analytical solutions. The three crack geometries are a pressurized crack orthogonal to the least principle stress, a slanted straight crack, and a circular arc crack. Several authors selected these specific geometries to justify the use of higher-order or specialized boundary elements [24-27]. However, the quantification of the computational efficiency coupled with the accuracy has been limited. Therefore, we present the following analysis that aids in determining the method that provides the most efficient, yet accurate solutions. Accurate and efficient methods are required for the development of field applications of engineering software packages.

Several other numerical techniques can be implemented within BEM. What we present here is not meant to be a review of possible combinations. We have chosen basic numerical techniques that provide the necessary flexibility to model very complex geometries, yet remain efficient enough for engineering modeling applications. The literature contains numerous examples of refinements to the techniques presented here [24-27]. For example, refinements with respect to the quarter-point method are found in Gray *et al*. [26] and refinements to higherorder elements are suggested by Dong and de Pater [25]. It is expected that implementing more refined methods will increase the efficiency of the numerical calculations. However, this work is primarily concerned with determining the general framework for BEM implementation.

The details related to the crack tip elements are available from a number of sources [10,12]. For brevity, this work will only summarize some basic concepts and mathematical formulas of the higher-order elements and the specialized crack tip elements. The next section of this paper describes the general displacement discontinuity method utilizing constant displacement elements. In section 3, the authors summarize the chosen higher order elements. Section 4 defines the special crack tip elements used in this work. Section 5 compares various combina‐ tions of the presented methods to known solutions of various crack geometries for an estima‐ tion of accuracy of calculations. Section 5 concludes by comparing of the computational efficiencies exhibited by the various methods. Finally, some concluding remarks are provided in Section 6.

#### **2. Displacement discontinuity method**

The displacement discontinuity method (DDM), originally formulated by Crouch [12], is used here. DDM is based on the solution of the stresses and displacements at a point caused by a constant displacement discontinuity (DD) over a line segment in an elastic body under prescribed boundary conditions [12]. Due to the simplicity of mathematical formulas and procedures of DDM (with a constant DD), it has been widely applied to various engineering problems. This paper summarizes some of the important mathematical expressions but limits specificity. The details of DDM are well described in the literature [12].

The 2-D displacements and stresses at a point (*x*, *y*), generated by a displacement discontinuity (*Dx*(*x*), *Dy*(*x*)) on the line segment | *x* | ≤*a*, *y* =0, can be analytically expressed as follows [1-3]:

$$\mu\_x = \left\{ 2(1 - \nu) f\_{,y} - y f\_{,xx} \right\} + \left\{ -(1 - 2\nu) g\_{,x} - y g\_{,xy} \right\} \tag{1}$$

( ) ( ,0 ) ( ,0 ) *Dx ux ux xx x*

( ) ( ,0 ) ( ,0 ) *Dx ux ux yy y*

come out of the integrals, and then Equations (6) and (7) can be simplified as:

2 2

[ ]

( )ln ( ) 2

++ + + -

established in the available literature [12], they are not given herein.

elements increase the number of equations that must be solved.

x

*xa xa y a*

( , ) ln[ ( ) ]


*a <sup>a</sup> I xy x yd*

ò

2 2

 x

arctan arctan ( )ln ( ) () ()

For simplicity, the derivatives of *I*0(*x*, *y*), used to calculate the stresses and displacements (i.e. Equations (1) to (5)), are omitted in this paper. The derivatives are given in Shou *et al.* [10]. Since the numerical procedures of DDM (with constant displacement discontinuities) are well

However, DDM with a constant DD can't accurately calculate the stresses and displacements of the area closer than about one element-length distance from a boundary [12]. To improve the accuracy of calculations in close proximity to the boundaries, Crawford *et al.* developed higher-order displacement elements [24] among others [25]. Although higher-order elements overcame the limitations of constant elements and improved the accuracy of DDM, the method significantly increases the number of degrees of freedom. In other words, the higher-order

To improve the accuracy of DDM without sacrificing the number of degrees of freedom of the overall system, a new higher-order elements method was suggested by Shou *et al.* [10]. The method used collocation points at the centers of the source elements and its two adjacent neighbors, so it could maintain the same degrees of freedom as the constant elements method by sharing the DD of the two adjacent neighbors. Other methods have been suggested by in the literature [15,25] to overcome issues with kinked or intersecting cracks when utilizing

This study uses Shou *et al.*'s method to satisfy one of this research's objectives, which is to develop methods that reduce computation costs while improving accuracy. The next section

<sup>=</sup> - -- - + - +

*y y <sup>y</sup> xa xa y xa xa*

where

0

neighboring elements in calculations.

For constant displacement elements (i.e. constant *Dx*(*x*) and *Dy*(*x*)), the DD components can


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

835

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems


<sup>0</sup> ( , ) ( , ) *<sup>x</sup> f xy I xyD* = (10)

<sup>0</sup> ( , ) ( , ) *<sup>y</sup> gxy I xyD* = (11)

2 2

(12)

$$\mathbf{u}\_y = \left\{ (1 - 2\nu) f\_{,x} - y f\_{,xy} \right\} + \left\{ 2(1 - \nu) \mathbf{g}\_{,y} - y \mathbf{g}\_{,yy} \right\} \tag{2}$$

$$\sigma\_{xx} = \text{2G} \{ \text{2}f\_{,yy} + yf\_{,yyy} \} + \text{2G} \{ \text{g}\_{,yy} + yg\_{,yyy} \} \tag{3}$$

$$\sigma\_{yy} = \text{2G} \{-y f\_{,yy}\} + \text{2G} \{\text{g}\_{,yy} - y \text{g}\_{,yy}\} \tag{4}$$

$$
\sigma\_{xy} = \text{2G} \{ f\_{,yy} + yf\_{,yyy} \} + \text{2G} \{ -yg\_{,xyy} \} \tag{5}
$$

where *f* (*x*, *y*) and *g*(*x*, *y*) are defined as:

$$f(\mathbf{x}, y) = \frac{-1}{4\pi(1 - \nu)} \int\_{-d}^{d} D\_x(\xi) \ln[\sqrt{(\mathbf{x} - \xi)^2 + y^2}] d\xi \tag{6}$$

$$\log(\mathbf{x}, y) = \frac{-1}{4\pi(1 - \nu)} \int\_{-a}^{a} D\_y(\xi) \ln[\sqrt{(\mathbf{x} - \xi)^2 + y^2}] d\xi \tag{7}$$

and where the displacement discontinuity components are defined as:

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems http://dx.doi.org/10.5772/56445 835

$$D\_{\mathbf{x}}(\mathbf{x}) = \mu\_{\mathbf{x}}(\mathbf{x}, 0^{-}) - \mu\_{\mathbf{x}}(\mathbf{x}, 0^{+}) \tag{8}$$

$$D\_y(\mathbf{x}) = \mu\_y(\mathbf{x}, 0^-) - \mu\_y(\mathbf{x}, 0^+) \tag{9}$$

For constant displacement elements (i.e. constant *Dx*(*x*) and *Dy*(*x*)), the DD components can come out of the integrals, and then Equations (6) and (7) can be simplified as:

$$f(\mathbf{x}, y) = I\_0(\mathbf{x}, y)D\_\mathbf{x} \tag{10}$$

$$\log(\mathbf{x}, \mathbf{y}) = I\_0(\mathbf{x}, \mathbf{y}) D\_y \tag{11}$$

where

tions of the presented methods to known solutions of various crack geometries for an estima‐ tion of accuracy of calculations. Section 5 concludes by comparing of the computational efficiencies exhibited by the various methods. Finally, some concluding remarks are provided

The displacement discontinuity method (DDM), originally formulated by Crouch [12], is used here. DDM is based on the solution of the stresses and displacements at a point caused by a constant displacement discontinuity (DD) over a line segment in an elastic body under prescribed boundary conditions [12]. Due to the simplicity of mathematical formulas and procedures of DDM (with a constant DD), it has been widely applied to various engineering problems. This paper summarizes some of the important mathematical expressions but limits

The 2-D displacements and stresses at a point (*x*, *y*), generated by a displacement discontinuity (*Dx*(*x*), *Dy*(*x*)) on the line segment | *x* | ≤*a*, *y* =0, can be analytically expressed as follows [1-3]:

> n

 n

*f yf g yg* (1)

*f yf g yg* (2)

= ++ + (3)

=- + - (4)

= + +- (5)



, , , , [2(1 ) ] [ (1 2 ) ] *<sup>x</sup> y xx x xy u* = - - +- - -

, , , , [(1 2 ) ] [2(1 ) ] *<sup>y</sup> x xy y yy u* =- - + - -

,, , , 2 [2 ] 2 [ ] *xx G f yf G g yg xy xyy yy yyy*

, ,, 2[ ] 2[ ] *yy G yf G g yg xyy yy yyy*

,, , 2[ ] 2[ ] *xy yy yyy G f yf G yg xyy*

<sup>1</sup> 2 2 (,) ( )ln[ ( ) ] 4 (1 ) *a <sup>x</sup> <sup>a</sup> fxy D x yd* x

<sup>1</sup> 2 2 (,) ( )ln[ ( ) ] 4 (1 ) *a <sup>y</sup> <sup>a</sup> gxy D x yd* x

 xx

 xx

specificity. The details of DDM are well described in the literature [12].

n

n

s

s

s

p n-

p n-

and where the displacement discontinuity components are defined as:

where *f* (*x*, *y*) and *g*(*x*, *y*) are defined as:

in Section 6.

834 Effective and Sustainable Hydraulic Fracturing

**2. Displacement discontinuity method**

$$\begin{aligned} I\_0(\mathbf{x}, y) &= \int\_{-a}^{a} \ln[\sqrt{(\mathbf{x} - \boldsymbol{\xi})^2 + y^2}^2] d\boldsymbol{\xi} \\ &= y \Big[ \arctan \frac{y}{\left(\mathbf{x} - a\right)} - \arctan \frac{y}{\left(\mathbf{x} + a\right)} \Big] - (\mathbf{x} - a) \ln \sqrt{\left(\mathbf{x} - a\right)^2 + y^2} \\ &+ (\mathbf{x} + a) \ln \sqrt{\left(\mathbf{x} + a\right)^2 + y^2}^2 - 2a \end{aligned} \tag{12}$$

For simplicity, the derivatives of *I*0(*x*, *y*), used to calculate the stresses and displacements (i.e. Equations (1) to (5)), are omitted in this paper. The derivatives are given in Shou *et al.* [10]. Since the numerical procedures of DDM (with constant displacement discontinuities) are well established in the available literature [12], they are not given herein.

However, DDM with a constant DD can't accurately calculate the stresses and displacements of the area closer than about one element-length distance from a boundary [12]. To improve the accuracy of calculations in close proximity to the boundaries, Crawford *et al.* developed higher-order displacement elements [24] among others [25]. Although higher-order elements overcame the limitations of constant elements and improved the accuracy of DDM, the method significantly increases the number of degrees of freedom. In other words, the higher-order elements increase the number of equations that must be solved.

To improve the accuracy of DDM without sacrificing the number of degrees of freedom of the overall system, a new higher-order elements method was suggested by Shou *et al.* [10]. The method used collocation points at the centers of the source elements and its two adjacent neighbors, so it could maintain the same degrees of freedom as the constant elements method by sharing the DD of the two adjacent neighbors. Other methods have been suggested by in the literature [15,25] to overcome issues with kinked or intersecting cracks when utilizing neighboring elements in calculations.

This study uses Shou *et al.*'s method to satisfy one of this research's objectives, which is to develop methods that reduce computation costs while improving accuracy. The next section will summarize some basic concepts and mathematical formulas for the higher-order elements used in this work. As above, further details of the higher-order elements are available in [10].

Combining Equations (13) through (16) with Equations (6) and (7) gives the following

012

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

012


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

837


ò (19)

3

p n=

p n=

x

definition of these kernels is given by Shou *et al.* [10].

**Figure 2.** Representation of a crack by N elemental displacement discontinuities [12]

1 <sup>1</sup> (,) ( ) (,,) 4 (1 ) *xjj j fxy D FI I I*

3

1 <sup>1</sup> (,) ( ) (,,) 4 (1 ) *yjj j gxy D FI I I*

<sup>012</sup> ( , , ) ( )ln[ ( ) ] , 1 3 *j j F I I I N x y d j to* = -+ =

 xx

which can be expressed in terms of constant, linear, and quadratic kernels (*I*0, *I*1, *I*2). The

Based on these formulas, a crack can be discretized into N elements (see Figure 2) and 2N equations in terms of the DD component unknowns are formed (i.e. 2N unknowns of *Dx* and *Dy*). Under certain boundary conditions, the 2N unknowns can be obtained. Once the 2N unknowns are calculated, the 2-D displacements and stresses at a point (*x*, *y*) can be calculated through Equations (1) to (5). Further, Equations (20) and (21) compute the stress intensity

> \* \* 2 4(1 ) *I n <sup>G</sup> K D <sup>a</sup>*

<sup>=</sup> - (20)

p

n

where the subscript *j* indicates the *j* th collocation node in the three-element patch and

2 2

simplified expressions:

(*I*0, *I*1, *I*2) is defined as

factors at the crack tip.

*Fj*

#### **3. Higher-order elements of displacement discontinuity method**

Higher-order elements, as formulated by Shou *et al.,* use quadratic displacement elements. The calculation of the DD component of a particular element is accomplished by using three collocation points. The center collocation point is within the element of interest, while the bounding collocation points are within the neighboring elements. This configuration forms a three-element "patch" (shown in Figure 1), on which the quadratic formulation is performed. Equation (13) shows how the value of the DD components is formed mathematically.

**Figure 1.** Quadratic collocation for the new higher-order elements [2]

$$D\_i(\xi) = N\_1(\xi)(D\_i)\_1 + N\_2(\xi)(D\_i)\_2 + N\_3(\xi)(D\_i)\_3 \tag{13}$$

where (*Di* ) 1 , (*Di* ) 2 and (*Di* ) 3 are the nodal displacement discontinuities (*i* = *x or y*) and

$$N\_1 = \frac{\xi(\xi - a\_2 - a\_3)}{(a\_1 + a\_2)(a\_1 + 2a\_2 + a\_3)}\tag{14}$$

$$N\_2 = \frac{-(\xi + a\_1 + a\_2)(\xi - a\_2 - a\_3)}{(a\_1 + a\_2)(a\_2 + a\_3)}\tag{15}$$

$$N\_3 = \frac{\xi(\xi + a\_1 + a\_2)}{(a\_2 + a\_3)(a\_1 + 2a\_2 + a\_3)}\tag{16}$$

*N*1, *N*2 and *N*3 are the collocation shape functions whose *a*1, *a*<sup>2</sup> and *a*3 are half length of the three elements of the patch.

Combining Equations (13) through (16) with Equations (6) and (7) gives the following simplified expressions:

will summarize some basic concepts and mathematical formulas for the higher-order elements used in this work. As above, further details of the higher-order elements are available in [10].

Higher-order elements, as formulated by Shou *et al.,* use quadratic displacement elements. The calculation of the DD component of a particular element is accomplished by using three collocation points. The center collocation point is within the element of interest, while the bounding collocation points are within the neighboring elements. This configuration forms a three-element "patch" (shown in Figure 1), on which the quadratic formulation is performed.

Equation (13) shows how the value of the DD components is formed mathematically.

1 12 23 3 ( ) ( )( ) ( )( ) ( )( ) *D NDNDND i iii*

 x

2 3

12 23

1 2

*N*1, *N*2 and *N*3 are the collocation shape functions whose *a*1, *a*<sup>2</sup> and *a*3 are half length of the

 x

1 21 2 3 ( ) ( )( 2 ) *a a <sup>N</sup> a aa a a* x x

1 22 3 ( )( ) ( )( ) *aa aa <sup>N</sup> a aa a* -+ + - -

2 31 2 3 ( ) ( )( 2 ) *a a <sup>N</sup> a aa a a* x x

 x

are the nodal displacement discontinuities (*i* = *x or y*) and

=+ + (13)


<sup>=</sup> + + (15)

+ + <sup>=</sup> + ++ (16)

**Figure 1.** Quadratic collocation for the new higher-order elements [2]

where (*Di*

) 1 , (*Di* ) 2 and (*Di* ) 3

836 Effective and Sustainable Hydraulic Fracturing

three elements of the patch.

xx

1

2

3

x

**3. Higher-order elements of displacement discontinuity method**

$$f(\mathbf{x}, y) = \frac{-1}{4\pi(1 - \nu)} \sum\_{j=1}^{3} (D\_{\mathbf{x}})\_j F\_j(I\_{0}, I\_{1}, I\_{2}) \tag{17}$$

$$\log(\mathbf{x}, \mathbf{y}) = \frac{-1}{4\pi(1-\nu)} \sum\_{j=1}^{3} (D\_y)\_j F\_j(I\_0, I\_1, I\_2) \tag{18}$$

where the subscript *j* indicates the *j* th collocation node in the three-element patch and *Fj* (*I*0, *I*1, *I*2) is defined as

$$F\_j(I\_0, I\_1, I\_2) = \int \mathcal{N}\_j(\xi) \ln[\sqrt{(\mathbf{x} - \xi)^2 + y^2}] d\xi, \quad j = 1 \text{ to } 3 \tag{19}$$

which can be expressed in terms of constant, linear, and quadratic kernels (*I*0, *I*1, *I*2). The definition of these kernels is given by Shou *et al.* [10].

Based on these formulas, a crack can be discretized into N elements (see Figure 2) and 2N equations in terms of the DD component unknowns are formed (i.e. 2N unknowns of *Dx* and *Dy*). Under certain boundary conditions, the 2N unknowns can be obtained. Once the 2N unknowns are calculated, the 2-D displacements and stresses at a point (*x*, *y*) can be calculated through Equations (1) to (5). Further, Equations (20) and (21) compute the stress intensity factors at the crack tip.

**Figure 2.** Representation of a crack by N elemental displacement discontinuities [12]

$$K\_{1} = \frac{G}{4(1-\nu)} \sqrt{\frac{2\pi}{a^{\*}}} D\_{n}^{\*} \tag{20}$$

$$K\_{II} = \frac{G}{4(1-\nu)} \sqrt{\frac{2\pi}{a^\ast}} D\_s^\ast \tag{21}$$

Thus, Shou *et al.* introduced a more sophisticated crack tip elements method to comply with the LEFM physical phenomenon, which is the square root element method summarized below.

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

LEFM predicts that in the vicinity of the crack tip the crack displacement is proportional to the square root of the distance from the crack tip (i.e.*w* ∝ *ξ*). Figure 4 illustrates the basics of the square root crack tip element. Equation (22) shows the representation of *ξ* variation of the

(*ξ*) along the crack tip.

( ) , *D D i xy i ci <sup>a</sup>* x

where *Dci* are the DD values at the center of the crack tip element. Substituting Equation (22) into Equations (1) to (5) the stresses and displacements are resolved in terms of *Dci*

solutions can be expressed in kernel functions, similar to higher-order elements methods. The

To access the accuracy of the selected numerical techniques, this study compares the stress near the crack tips and stress intensity factors for three crack geometries with constant displacement or higher-order elements, both elements will be combined with and without specialized quarter and square root crack tip elements. We chose these particular geometries because analytical solutions are readily available. Further, many previous publications comparing BEMs have chosen these same geometries [25-27]. This research uses following

= = (22)

. The

839

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

x

details of the kernel functions were well documented in previous work [10].

**4.2. Square root element method**

displacement discontinuities *Di*

**Figure 4.** Square root crack tip element [2]

elastic properties: *E* =10<sup>6</sup> *psi* and *ν* =0.2 in the calculations.

**5. Comparison**

where *<sup>a</sup>* \* is the half length of the crack tip element, and *Dn* \* and *Ds* \* are normal and shear DD at the crack tip, respectively.

#### **4. Crack tip elements of displacement discontinuity method**

In addition to advanced elements, Shou *et al.* formulated two special crack tip elements to capture the square root displacement variation at the crack tip, expected from Linear Elastic Fracture Mechanics (LEFM) [10]. One is to use a constrained collocation point one-quarter of an element length away from the end of the crack. This study will refer to this as a quarter element method. The other tip element is simply a prescribed displacement discontinuity proportional to the square root variation at the crack tip. Herein, it is called the square root element method. This study uses their methods to calculate the stresses and displacements at the crack tip.

#### **4.1. Quarter element method**

Shou *et al.* introduced a constrained collocation point one-quarter of the crack tip element length away from the crack tip element. The DDs of the point will be set zero. Figure 3 shows the element.

**Figure 3.** Quarter element method at a crack tip [2]

Numerical implementation of this method is more efficient compared to the square root element method and provides reasonably accurate results. However, the displacement discontinuities near the crack tip do not capture the square root displacement variation at the crack tip, expected from LEFM [10]. Theoretically, this method may give unreliable results. Thus, Shou *et al.* introduced a more sophisticated crack tip elements method to comply with the LEFM physical phenomenon, which is the square root element method summarized below.

#### **4.2. Square root element method**

\* \* 2 4(1 ) *II <sup>s</sup> <sup>G</sup> K D <sup>a</sup>*

<sup>=</sup> - (21)

\* are normal and shear DD

\* and *Ds*

p

In addition to advanced elements, Shou *et al.* formulated two special crack tip elements to capture the square root displacement variation at the crack tip, expected from Linear Elastic Fracture Mechanics (LEFM) [10]. One is to use a constrained collocation point one-quarter of an element length away from the end of the crack. This study will refer to this as a quarter element method. The other tip element is simply a prescribed displacement discontinuity proportional to the square root variation at the crack tip. Herein, it is called the square root element method. This study uses their methods to calculate the stresses and displacements at

Shou *et al.* introduced a constrained collocation point one-quarter of the crack tip element length away from the crack tip element. The DDs of the point will be set zero. Figure 3 shows

Numerical implementation of this method is more efficient compared to the square root element method and provides reasonably accurate results. However, the displacement discontinuities near the crack tip do not capture the square root displacement variation at the crack tip, expected from LEFM [10]. Theoretically, this method may give unreliable results.

n

**4. Crack tip elements of displacement discontinuity method**

where *<sup>a</sup>* \* is the half length of the crack tip element, and *Dn*

at the crack tip, respectively.

838 Effective and Sustainable Hydraulic Fracturing

the crack tip.

the element.

**4.1. Quarter element method**

**Figure 3.** Quarter element method at a crack tip [2]

LEFM predicts that in the vicinity of the crack tip the crack displacement is proportional to the square root of the distance from the crack tip (i.e.*w* ∝ *ξ*). Figure 4 illustrates the basics of the square root crack tip element. Equation (22) shows the representation of *ξ* variation of the displacement discontinuities *Di* (*ξ*) along the crack tip.

$$D\_i(\xi) = D\_{c1} \sqrt{\frac{\xi}{a}} \quad i = \text{x}, y \tag{22}$$

where *Dci* are the DD values at the center of the crack tip element. Substituting Equation (22) into Equations (1) to (5) the stresses and displacements are resolved in terms of *Dci* . The solutions can be expressed in kernel functions, similar to higher-order elements methods. The details of the kernel functions were well documented in previous work [10].

**Figure 4.** Square root crack tip element [2]

#### **5. Comparison**

To access the accuracy of the selected numerical techniques, this study compares the stress near the crack tips and stress intensity factors for three crack geometries with constant displacement or higher-order elements, both elements will be combined with and without specialized quarter and square root crack tip elements. We chose these particular geometries because analytical solutions are readily available. Further, many previous publications comparing BEMs have chosen these same geometries [25-27]. This research uses following elastic properties: *E* =10<sup>6</sup> *psi* and *ν* =0.2 in the calculations.

#### **5.1. Single pressurized crack in an infinite elastic domain**

A single pressurized crack is a basic fracture geometry, and the analytical solutions are well documented [19,28]. Figure 5 shows a schematic diagram for the fracture geometry. The crack is pressurized by a pressure *p* =1000 *psi*. The crack length is 10 inches and is discretized into 10 elements with equal length. According to the specified method, the two crack tips may be replaced by specialized crack tip elements. This research compares half width, stresses at defined locations, and stress intensity factors computed from each method.

analytical solution. The importance of utilizing special crack tip elements is well established

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

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

841

**Figure 6.** Dimensionless crack half width versus Dimensionless distance from the crack center (HDD, CDDCE, CDDCT

Figure 7 is more illustrative for comparing the accuracy of the various methods. It shows the relative error from the computation of the half width compared to the analytical solution (i.e.

*wAM %*). The relative errors of all methods increase as they approach to the crack tip. The majority of the methods demonstrate errors bounded between -5% to 5%, except for CDD, which shows over 20% in close proximity to the crack tip. Computational errors over 20% from

To evaluate the perturbed stress state due to the presence of a pressurized crack we use [30]

CDD methods have been reported in the literature [20].

**Figure 7.** Relative error of the half width

in the literature [10, 26, 27].

and HDDCE overlap)

*w* − *wAM*

Figure 5 illustrates the chosen locations where each of the given methods calculates the stress. The points are arbitrary, but chosen at a location of symmetry with respect to the fracture. In Figure 5, the blue X represents the point orthogonal to the fracture plane at the mid-point of the fracture. The red diamond represents a point ahead of the fracture tip. For convenience, this report uses the following abbreviations to represent each method: AM (analytical method), CDD (only constant displacement discontinuity), HDD (only higher-order elements), CDDCE (constant DD with the quarter element method), CDDCT (constant DD with the square root element method), HDDCE (higher-order elements with the quarter element method), and HDDCT (higher-order elements with the square root element method).

**Figure 5.** A crack under a constant pressure. The blue X represents a point orthogonal to the fracture plane where the induced stress calculated from each method is compared. The red diamond represents the evaluation point ahead of the fracture tip.

Equation (23) is the analytical solution of the dimensionless half width under a constant pressure [29].

$$\frac{w(x)}{c} = \frac{4p}{E'} \sqrt{1 - \left(\frac{x}{c}\right)^2} \tag{23}$$

Figure 6 is a plot of the calculated fracture profile from each method. The analytical solution is the solid blue line. The computed half widths from each of the numerical methods are shown as a series of points. From the results, CDD overestimates the half width particularly near crack tip area as others have found [20]. Conversely, HDD underestimates the fracture width in the proximity of the tip. Methods that use tip elements show fracture width profiles close to the analytical solution. The importance of utilizing special crack tip elements is well established in the literature [10, 26, 27].

**Figure 6.** Dimensionless crack half width versus Dimensionless distance from the crack center (HDD, CDDCE, CDDCT and HDDCE overlap)

Figure 7 is more illustrative for comparing the accuracy of the various methods. It shows the relative error from the computation of the half width compared to the analytical solution (i.e. *w* − *wAM wAM %*). The relative errors of all methods increase as they approach to the crack tip. The majority of the methods demonstrate errors bounded between -5% to 5%, except for CDD, which shows over 20% in close proximity to the crack tip. Computational errors over 20% from CDD methods have been reported in the literature [20].

**Figure 7.** Relative error of the half width

**5.1. Single pressurized crack in an infinite elastic domain**

840 Effective and Sustainable Hydraulic Fracturing

A single pressurized crack is a basic fracture geometry, and the analytical solutions are well documented [19,28]. Figure 5 shows a schematic diagram for the fracture geometry. The crack is pressurized by a pressure *p* =1000 *psi*. The crack length is 10 inches and is discretized into 10 elements with equal length. According to the specified method, the two crack tips may be replaced by specialized crack tip elements. This research compares half width, stresses at

Figure 5 illustrates the chosen locations where each of the given methods calculates the stress. The points are arbitrary, but chosen at a location of symmetry with respect to the fracture. In Figure 5, the blue X represents the point orthogonal to the fracture plane at the mid-point of the fracture. The red diamond represents a point ahead of the fracture tip. For convenience, this report uses the following abbreviations to represent each method: AM (analytical method), CDD (only constant displacement discontinuity), HDD (only higher-order elements), CDDCE (constant DD with the quarter element method), CDDCT (constant DD with the square root element method), HDDCE (higher-order elements with the quarter element method), and

**Figure 5.** A crack under a constant pressure. The blue X represents a point orthogonal to the fracture plane where the induced stress calculated from each method is compared. The red diamond represents the evaluation point ahead of

Equation (23) is the analytical solution of the dimensionless half width under a constant

Figure 6 is a plot of the calculated fracture profile from each method. The analytical solution is the solid blue line. The computed half widths from each of the numerical methods are shown as a series of points. From the results, CDD overestimates the half width particularly near crack tip area as others have found [20]. Conversely, HDD underestimates the fracture width in the proximity of the tip. Methods that use tip elements show fracture width profiles close to the

*cE c* = - ¢ (23)

<sup>4</sup> <sup>2</sup> ( ) <sup>1</sup> ( ) *wx x <sup>p</sup>*

defined locations, and stress intensity factors computed from each method.

HDDCT (higher-order elements with the square root element method).

the fracture tip.

pressure [29].

To evaluate the perturbed stress state due to the presence of a pressurized crack we use [30]

$$\frac{\sigma\_{\text{xx}}}{1-p} = \frac{\left(\frac{\text{x}}{c}\right)}{\sqrt{\binom{\text{x}}{c}^2 - 1}} - 1 \tag{24}$$

These figures also show that CDD overestimates the stress and HDD underestimates it while the other methods give results with less than 1% error. The inaccuracies of the methods without tip elements become significant closer to the crack tip. Similar to the half width results, the results of the methods with tip elements overlap since they are close to the analytical solution.

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

2 2 3 1 1 1

(25)

843

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

*<sup>p</sup>* at a point along the Y-axis (the distance from

+ +

() ()

*L L c c*

() {() }

the crack face is normalized by the crack half-length), which is represented by the blue X in Figure 5. Figure 10 shows the results of the calculation from each numerical method. These results show consistency with the analytical solution, regardless of the numerical method. This is expected, as the location where the stress is calculated is sufficiently far from the crack tip,

**Figure 10.** Dimensionless stress versus Dimensionless distance orthogonal to the fracture plane, i.e. at the blue X. In

Table 1 shows the calculated stress intensity factor, along with the ratio to the analytical solution. For the mode I (or KI), CDD shows the biggest error while CDDCE, HDDCE and HDDCT give around 1% errors. Obviously, the mode II (or KII) stress intensity factor is zero.

Reasonable values from calculations of the stresses and displacements can be achieved at distances greater than the length of one discretized element, regardless of the numerical technique. Close to the crack tip, however, the error of displacements, stresses, and stress intensity factors for CDD elements are significant, whereas CDDCE, HDDCE and HDDCT provide reasonable estimations. This is not surprising; similar results are well documented in

To calculate the stress induced orthogonal to a pressurized crack we use [30]

*p L L c c*

=- + +

*yy*

s

i.e. more than the length of one discretized element [12].

Equation (25) expresses a dimensionless stress *σyy*

this case, all methods overlap

So, Table 1 omits the results.

the literature [12,24,26].

Equation (24) provides a dimensionless stress *σxx <sup>p</sup>* at a point along X-axis (the distance from the crack tip normalized by the crack half-length), which in this case is red diamond in Figure 5.

**Figure 8.** Dimensionless stress versus Dimensionless distance ahead of the crack tip (i.e. at red diamond) (CDDCE, HDDCE and HDDCT overlap)

Figure 8 plots the dimensionless stress ahead of the crack tip. x/c = 1 represents the crack tip. Figure 9 plots the relative stress (the ratio of *σxx p* to the analytical *σxx <sup>p</sup>* ) near the crack tip area for each numerical method.

**Figure 9.** Relative stress versus Dimensionless distance ahead of the crack tip (CDDCE, HDDCE and HDDCT overlap)

These figures also show that CDD overestimates the stress and HDD underestimates it while the other methods give results with less than 1% error. The inaccuracies of the methods without tip elements become significant closer to the crack tip. Similar to the half width results, the results of the methods with tip elements overlap since they are close to the analytical solution.

To calculate the stress induced orthogonal to a pressurized crack we use [30]

2

crack tip normalized by the crack half-length), which in this case is red diamond in Figure 5.

**Figure 8.** Dimensionless stress versus Dimensionless distance ahead of the crack tip (i.e. at red diamond) (CDDCE,

Figure 8 plots the dimensionless stress ahead of the crack tip. x/c = 1 represents the crack tip.

**Figure 9.** Relative stress versus Dimensionless distance ahead of the crack tip (CDDCE, HDDCE and HDDCT overlap)

*p* to the analytical *σxx*

*<sup>p</sup>* ) near the crack tip area

( )

= - - -

*xx*

s

Equation (24) provides a dimensionless stress *σxx*

842 Effective and Sustainable Hydraulic Fracturing

Figure 9 plots the relative stress (the ratio of *σxx*

HDDCE and HDDCT overlap)

for each numerical method.

( )

*x c p x c*

1 1

*<sup>p</sup>* at a point along X-axis (the distance from the

(24)

$$\frac{\sigma\_{yy}}{p} = -\frac{\left(\frac{L}{c}\right)}{\sqrt{\left(\frac{L}{c}\right)^2 + 1}} + 1 + \frac{\left(\frac{L}{c}\right)}{\left\{\sqrt{\left(\frac{L}{c}\right)^2 + 1}\right\}^3} \tag{25}$$

Equation (25) expresses a dimensionless stress *σyy <sup>p</sup>* at a point along the Y-axis (the distance from the crack face is normalized by the crack half-length), which is represented by the blue X in Figure 5. Figure 10 shows the results of the calculation from each numerical method. These results show consistency with the analytical solution, regardless of the numerical method. This is expected, as the location where the stress is calculated is sufficiently far from the crack tip, i.e. more than the length of one discretized element [12].

**Figure 10.** Dimensionless stress versus Dimensionless distance orthogonal to the fracture plane, i.e. at the blue X. In this case, all methods overlap

Table 1 shows the calculated stress intensity factor, along with the ratio to the analytical solution. For the mode I (or KI), CDD shows the biggest error while CDDCE, HDDCE and HDDCT give around 1% errors. Obviously, the mode II (or KII) stress intensity factor is zero. So, Table 1 omits the results.

Reasonable values from calculations of the stresses and displacements can be achieved at distances greater than the length of one discretized element, regardless of the numerical technique. Close to the crack tip, however, the error of displacements, stresses, and stress intensity factors for CDD elements are significant, whereas CDDCE, HDDCE and HDDCT provide reasonable estimations. This is not surprising; similar results are well documented in the literature [12,24,26].


**Table 1.** Stress intensity factors

#### **5.2. A slanted straight crack**

While a straight pressurized crack shows zero KII, a slanted straight crack under a uniform tension can show a variable KII depending on the angle of incidence to the applied tension. Figure 11 illustrates the crack geometry. The stress intensity factors are calculated by [19]

$$K\_{\rm I} = \sigma \sqrt{\pi a} \sin^2(\beta) \qquad K\_{\rm II} = \sigma \sqrt{\pi a} \sin(\beta) \cos(\beta) \tag{26}$$

overlap in Figure 13 and 14. However, the calculation errors of the two stress intensity factors exhibit opposing patterns. For KI, the errors increase as the slanted angle becomes larger.

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

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

845

versus the slanted angle (CDDCE, HDDCE and HDDCT overlap)

Curved cracks may represent more realistic fracture geometry and exhibit complexity in calculation of the stress intensity factors. This study selects a circular arc crack under a far field uniform biaxial tension in order to evaluate the accuracy of the numerical methods. Figure 14 describes the fracture geometry. The uniform biaxial tension is *σ* =1000 *psi*. The analytical

**Figure 13.** Dimensionless KII versus the slanted angle (CDDCE, HDDCE and HDDCT overlap)

solutions of KI and KII for a curved crack are given by [19]

Conversely, the errors of KII are maximal at 45° and at a minimum at 0° and 90°.

**Figure 12.** Dimensionless KI

**5.3. A circular arc crack**

Equation (26) expresses the analytical solution of the stress intensity factors [19]. The uniform tension is *σ* =1000 *psi*. The crack length is 10 inches. It is discretized into 10 elements with equal length. The two crack tips are replaced by crack tip elements according to the applied method.

**Figure 11.** Slanted straight crack under uniform axial tension at infinity [25]

Figure 12 shows the dimensionless KI (i.e. *KI σ πa* ) according to the slanted angle and Figure 13 gives the dimensionless KII (i.e. *KII σ πa* ) according to the slanted angle, respectively. Similar to the previous crack geometry, CDD and CDDCT overestimate KI and KII while HDD underestimates them. CDDCE, HDDCE, and HDDCT show fairly accurate results, so that they overlap in Figure 13 and 14. However, the calculation errors of the two stress intensity factors exhibit opposing patterns. For KI, the errors increase as the slanted angle becomes larger. Conversely, the errors of KII are maximal at 45° and at a minimum at 0° and 90°.

**Figure 12.** Dimensionless KI versus the slanted angle (CDDCE, HDDCE and HDDCT overlap)

**Figure 13.** Dimensionless KII versus the slanted angle (CDDCE, HDDCE and HDDCT overlap)

#### **5.3. A circular arc crack**

**KI [psi √in] KI/ KIAM**

AM 3963.3 1 CDD 4905.6 1.238 HDD 3670.8 0.926 CDDCE 3933.8 0.993 CDDCT 4219 1.065 HDDCE 3933.8 0.993 HDDCT 3918.5 0.989

While a straight pressurized crack shows zero KII, a slanted straight crack under a uniform tension can show a variable KII depending on the angle of incidence to the applied tension. Figure 11 illustrates the crack geometry. The stress intensity factors are calculated by [19]

> sp

Equation (26) expresses the analytical solution of the stress intensity factors [19]. The uniform tension is *σ* =1000 *psi*. The crack length is 10 inches. It is discretized into 10 elements with equal length. The two crack tips are replaced by crack tip elements according to the applied method.

> *KI σ πa*

to the previous crack geometry, CDD and CDDCT overestimate KI and KII while HDD underestimates them. CDDCE, HDDCE, and HDDCT show fairly accurate results, so that they

 b b

(26)

) according to the slanted angle and Figure

) according to the slanted angle, respectively. Similar

<sup>2</sup> sin ( ) sin( )cos( ) *Ka K a <sup>I</sup> II* = =

 b

sp

**Figure 11.** Slanted straight crack under uniform axial tension at infinity [25]

*KII σ πa*

Figure 12 shows the dimensionless KI (i.e.

13 gives the dimensionless KII (i.e.

**Table 1.** Stress intensity factors

**5.2. A slanted straight crack**

844 Effective and Sustainable Hydraulic Fracturing

Curved cracks may represent more realistic fracture geometry and exhibit complexity in calculation of the stress intensity factors. This study selects a circular arc crack under a far field uniform biaxial tension in order to evaluate the accuracy of the numerical methods. Figure 14 describes the fracture geometry. The uniform biaxial tension is *σ* =1000 *psi*. The analytical solutions of KI and KII for a curved crack are given by [19]

$$K\_{\rm I} = \frac{\sigma \sqrt{\pi r \sin(\alpha/2)}}{1 + \sin^2(\alpha/4)} \cos(\alpha/4) \qquad K\_{\rm II} = \frac{\sigma \sqrt{\pi r \sin(\alpha/2)}}{1 + \sin^2(\alpha/4)} \sin(\alpha/4) \tag{27}$$

**Figure 16.** KII versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

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

847

**Figure 18.** Relative KII versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

**Figure 17.** Relative KI

where *r* is the radius of the circular arc [19].

**Figure 14.** Circular arc crack under uniform biaxial tension [8]

Figures (15) and (16) show the values of the calculated KI and KII as a function of the circular crack angle, respectively. The figures depict actual stress intensity factor values under the prescribed biaxial tension, instead of dimensionless values. The unit of KI and KII is *psi in*. The effective half-length of the crack is ambiguous due to variable circular arc length related to the prescribed circular arc angle (α). The circular arc has a 10 inch radius. It is discretized into 20 elements with equal length. The two tip elements are evenly discretized into 10 additional segments to apply the tip elements methods, respectively. Thus, the circular arc has 38 segments (i.e. 18 middle identical elements and 20 identical tip elements).

Figures (17) and (18) show the normalized KI and KII stress intensity factors as a function of the circular crack angle, respectively.

**Figure 15.** KI versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems http://dx.doi.org/10.5772/56445 847

**Figure 16.** KII versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

2 2 sin( / 2) sin( / 2) cos( / 4) sin( / 4) 1 sin ( / 4) 1 sin ( / 4) *<sup>I</sup> II*

sp

Figures (15) and (16) show the values of the calculated KI and KII as a function of the circular crack angle, respectively. The figures depict actual stress intensity factor values under the prescribed biaxial tension, instead of dimensionless values. The unit of KI and KII is *psi in*. The effective half-length of the crack is ambiguous due to variable circular arc length related to the prescribed circular arc angle (α). The circular arc has a 10 inch radius. It is discretized into 20 elements with equal length. The two tip elements are evenly discretized into 10 additional segments to apply the tip elements methods, respectively. Thus, the circular arc has 38

Figures (17) and (18) show the normalized KI and KII stress intensity factors as a function of

segments (i.e. 18 middle identical elements and 20 identical tip elements).

**Figure 15.** KI versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

 a

 a= = + + (27)

 a

*r r K K*

a

sp

846 Effective and Sustainable Hydraulic Fracturing

where *r* is the radius of the circular arc [19].

**Figure 14.** Circular arc crack under uniform biaxial tension [8]

the circular crack angle, respectively.

 a

a

**Figure 17.** Relative KI versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

**Figure 18.** Relative KII versus the circular arc angle (CDDCE, HDDCE and HDDCT overlap)

Based on the results presented in this section, numerical calculations of KI and KII show similar patterns. As the circular crack angle increases, the differences between the analytically derived stress intensity factors increase. The ratio between the analytical and numerical stress intensity factors decreases, however. CDD and CDDCT show the overestimation while HDD method underestimates the stress intensity factors. The other methods (CDDCE, HDDCE, and HDDCT) provide very close results compared to the analytical solution.

#### **5.4. Computational efficiency**

In general, we find that CDDCE, HDDCE, and HDDCT methods significantly increase the accuracy of the computation of stress intensity factors for the geometries presented here. However, the required computational resource varies among the numerical methods. Further, for constant displacement elements, accuracy of the stress calculations ahead of the crack tip can be improved by increasing the number of elements. In order to objectively evaluate the efficiency of the numerical methods, we return to the pressurized crack example from Section 5.1. The following section first compares the accuracy improvement by increasing the number of elements for the constant element method. Then we compare the computation time for calculating the stress intensity factors for this particular problem.

**Figure 20.** Computation time for the solution of the pressurized crack exercise

element choices and increasingly complex crack geometries.

**6. Conclusion**

Finally, the computational efficiency is evaluated by inspecting Figure 20. This figure plots the calculation time of the various numerical methods for the pressurized crack exercise described above. As expected, CDD calculations with fewer elements are completed more quickly. Interestingly the CDDCE, HDDCE, and HDDCT methods show similar computation times while maintaining the highest accuracy of the numerical methods. This result allows for further refinement and evaluation of the CDDCE, HDDCE, and HDDCT methods using more refined

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

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

849

Overall results show that CDD gives prominent errors of calculations of stresses, displace‐ ments, and stress intensity factor compared to the other methods. Particularly, when ap‐ proaching to the crack tips and a fracture is curved, the errors of CDD significantly increase. Replacing tip elements by special crack tip elements can mitigate calculation errors when close to the crack tips. Using higher-order elements helps to reduce errors for the simple straight crack geometries. When a fracture is curved, the efficiency of combining specialized crack tip elements in computational errors in the calculation of KI and KII is more important than for simple fracture geometries. Combination of higher-order elements and crack tip elements give the most accurate calculations, yet retain the necessary efficiency. However, the overall efficiency of CDDCE, HDDCE and HDDCT methods cannot be definitively evaluated using the simple geometries shown here. We reserve that analysis for subsequent publications.

For comparison within this work, the numerical methods maintained a similar number of elements for each fracture geometry. Increasing the number of CDD elements increases the accuracy of the CDD method. However, this requires increased computation time. Thus, the use of higher-order elements and crack tip elements is likely warranted if considering the development of more accurate and efficient numerical simulations for field engineering applications where computation resources are restricted. Evaluating the efficiency of specific

The computer specifications used in this work are as follows: CPU-Intel® Xeon W3670 @ 3.2GHZ, installed memory (RAM)-24 gigabytes, OS- 64-bit Windows 7®, Software- Matlab® R2011b.

This study uses the stress of a horizontal crack near the crack tip to show the computational accuracy of simply increasing the number of CDD elements compared to higher-order elements and/or special tip elements. Figure 19, shows the normalized stress (the ratio of *σxx* / *p* to the analytical *σxx* / *p*) near the crack tip as a function of the number of CDD elements. The number at the legend indicates the number of CDD elements. As expected, increasing the number of elements results in a corresponding improvement in calculation accuracy. Figure 19 illustrates that at least 100 CDD elements are required to show less than 1% error. In other words, an order of magnitude increase in the number of CDD elements is required to match the accuracy derived with CDDCE, HDDCE and HDDCT methods (see Figure 9).

**Figure 19.** Relative stress versus Dimensionless distance ahead of the crack tip (i.e. at red diamond) according to the number of CDD elements used

**Figure 20.** Computation time for the solution of the pressurized crack exercise

Finally, the computational efficiency is evaluated by inspecting Figure 20. This figure plots the calculation time of the various numerical methods for the pressurized crack exercise described above. As expected, CDD calculations with fewer elements are completed more quickly. Interestingly the CDDCE, HDDCE, and HDDCT methods show similar computation times while maintaining the highest accuracy of the numerical methods. This result allows for further refinement and evaluation of the CDDCE, HDDCE, and HDDCT methods using more refined element choices and increasingly complex crack geometries.

#### **6. Conclusion**

Based on the results presented in this section, numerical calculations of KI and KII show similar patterns. As the circular crack angle increases, the differences between the analytically derived stress intensity factors increase. The ratio between the analytical and numerical stress intensity factors decreases, however. CDD and CDDCT show the overestimation while HDD method underestimates the stress intensity factors. The other methods (CDDCE, HDDCE, and

In general, we find that CDDCE, HDDCE, and HDDCT methods significantly increase the accuracy of the computation of stress intensity factors for the geometries presented here. However, the required computational resource varies among the numerical methods. Further, for constant displacement elements, accuracy of the stress calculations ahead of the crack tip can be improved by increasing the number of elements. In order to objectively evaluate the efficiency of the numerical methods, we return to the pressurized crack example from Section 5.1. The following section first compares the accuracy improvement by increasing the number of elements for the constant element method. Then we compare the computation time for

The computer specifications used in this work are as follows: CPU-Intel® Xeon W3670 @ 3.2GHZ, installed memory (RAM)-24 gigabytes, OS- 64-bit Windows 7®, Software- Matlab®

This study uses the stress of a horizontal crack near the crack tip to show the computational accuracy of simply increasing the number of CDD elements compared to higher-order elements and/or special tip elements. Figure 19, shows the normalized stress (the ratio of *σxx* / *p* to the analytical *σxx* / *p*) near the crack tip as a function of the number of CDD elements. The number at the legend indicates the number of CDD elements. As expected, increasing the number of elements results in a corresponding improvement in calculation accuracy. Figure 19 illustrates that at least 100 CDD elements are required to show less than 1% error. In other words, an order of magnitude increase in the number of CDD elements is required to match

**Figure 19.** Relative stress versus Dimensionless distance ahead of the crack tip (i.e. at red diamond) according to the

the accuracy derived with CDDCE, HDDCE and HDDCT methods (see Figure 9).

HDDCT) provide very close results compared to the analytical solution.

calculating the stress intensity factors for this particular problem.

**5.4. Computational efficiency**

848 Effective and Sustainable Hydraulic Fracturing

R2011b.

number of CDD elements used

Overall results show that CDD gives prominent errors of calculations of stresses, displace‐ ments, and stress intensity factor compared to the other methods. Particularly, when ap‐ proaching to the crack tips and a fracture is curved, the errors of CDD significantly increase. Replacing tip elements by special crack tip elements can mitigate calculation errors when close to the crack tips. Using higher-order elements helps to reduce errors for the simple straight crack geometries. When a fracture is curved, the efficiency of combining specialized crack tip elements in computational errors in the calculation of KI and KII is more important than for simple fracture geometries. Combination of higher-order elements and crack tip elements give the most accurate calculations, yet retain the necessary efficiency. However, the overall efficiency of CDDCE, HDDCE and HDDCT methods cannot be definitively evaluated using the simple geometries shown here. We reserve that analysis for subsequent publications.

For comparison within this work, the numerical methods maintained a similar number of elements for each fracture geometry. Increasing the number of CDD elements increases the accuracy of the CDD method. However, this requires increased computation time. Thus, the use of higher-order elements and crack tip elements is likely warranted if considering the development of more accurate and efficient numerical simulations for field engineering applications where computation resources are restricted. Evaluating the efficiency of specific combinations of higher-order elements coupled with specialized crack tip elements requires more complex geometries than presented here.

[2] Naceur, K. B, Thiercelin, M, & Touboul, E. Simulation of Fluid Flow in Hydraulic Fracturing: Implications for 3D Propagation. SPE Production Engineering (1990). SPE

Testing and Review of Various Displacement Discontinuity Elements for LEFM Crack Problems

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

851

[3] Germanovich, L, & Astakhov, D. Multiple fracture model. In: Jeffrey, McLennan, edi‐ tors. Proceedings of the three dimensional and advance hydraulic fracture modeling. Workshop held at 4th North American rock mechanics symposium, Seattle, July

[4] Bunger, A. P, Zhang, X, & Jeffrey, R. G. Parameters Affecting the Interaction Among Closely Spaced Hydraulic Fractures. SPE Journal;17(1):SPE-140426-PA, 292-306. [5] Olson, J. E, & Wu, K. Sequential vs. Simultaneous Multizone Fracturing in Horizon‐ tal Wells: Insights From a Non-Planar, Multifrac Numerical Model. In: SPE Hydraul‐ ic Fracturing Technology Conference. The Woodlands, Texas, USA: Society of

[6] Gu, H, & Siebrits, E. On numerical solutions of hyperbolic proppant transport prob‐ lems. In: Proceedings of the 10th international conference on hyperbolic problems:

[7] Gordeliy, E, & Detournay, E. A fixed grid algorithm for simulating the propagation of a shallow hydraulic fracture with a fluid lag. International Journal for Numerical

[8] Zhang, X, Jeffrey, R. G, & Thiercelin, M. Deflection and propagation of fluid-driven fractures at frictional bedding interfaces: A numerical investigation. Journal of Struc‐

[9] Sellers, E, & Napier, J. A comparative investigation of micro-flaw models fog the sim‐

[10] 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‐

[11] King, G. E. Thirty Years of Gas Shale Fracturing: What Have We Learned? In: SPE Annual Technical Conference and Exhibition. Florence, Italy: Society of Petroleum

[12] Crouch, S. L, & Starfield, A. M. Boundary Element Methods in Solid Mechanics: With Applications in Rock Mechanics and Geological Engineering: George Allen & Unwin;

[13] Crouch, S. L. Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution. International Journal for Numerical Methods in En‐

ulation of brittle fracture in rock. Computational Mechanics (1997).

theory, numerics and applications, Osaka, Japan, September (2004). , 13-17.

16032, 5(2), 133-141.

29-31, (2000). , 45-70.

tural Geology (2007).

Engineers. SPE 133456

1 edition, (1983).

gineering (1976).

Petroleum Engineers. SPE-M, 152602.

ence & Geomechanics Abstracts (1995).

and Analytical Methods in Geomechanics;35(5):602.
