## Abstract

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

## 1. Introduction

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

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

There are four general distinctive numerical methods to model fracture propagation problems:

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

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

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

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

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

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

SIF values can be obtained from the displacement discontinuity magnitudes at crack tip elements [17-19]. However, according to [16], DDM consistently overestimates displacement 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 displacement 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 calculation 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

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

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

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

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

## 2. Numerical procedure

### 2.1. Displacement discontinuity method

The general concept of the displacement discontinuity method proposed by [16] is to approximate 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 displacement), a system of *3N* linear algebraic equations can be set up as the following:

where *N* is the total number of elements, *s*, *d*, *n* are the directions of local coordinates depicted in Figure 3, *j*th element, *i*th element, and *A* is the boundary influence coefficient of the stresses tensor. If known values are the displacements of one side of boundary elements, these equations will be modified as:

where, *B* matrix is called the boundary influence coefficient of the displacements tensor.

### 2.2. Stress intensity factor computation

By knowing the crack tip element displacement discontinuities, K_{I}, K_{II} and K_{III} can be directly calculated using Equations 4a, b & c respectively:

where

## 3. Validation of numerical model

### 3.1. Rectangular crack

There is no analytical solution for the stress intensity factor variation along a rectangular crack front. However, rectangular cracks were the subject of several papers where the “Integral Equation” or “Body Force Method” was used to numerically approximate mixed Mode SIF values [37-42]. Results obtained from [39] are in a good agreement with [40] for maximum SIF calculation of rectangular cracks. In addition, [39] investigated how maximum stress intensity factors change in a half-space in terms of crack depth. Because of these reasons, [39] and [40] were selected as reference solutions to which we compare the results from this paper. Studies done by [37], [41] and [43] yield relatively different results for

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.

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 *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

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

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 dimensionless stress intensity factor,

Figure 8 shows the variation of dimensionless stress intensity factor, F_{1}, along the crack front y=b for various values of *x=0,* *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 (*F*_{1}*=1*). As *b/a* increases (equivalent to reducing the crack length *a* relative to *b*), *F*_{1} is reduced. When *b/a*=1.0, the square crack, *F*_{1}=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-uniform mesh.

Considering a rectangular vertical crack in a half-space, and assuming

where

Figure 10-a and b show for greater aspect ratio (

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,

### 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]:

where:

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

### 3.3. Penny-shaped crack

Stress intensity factor at the tip of a circular crack of radius

where:

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

## 4. Fracture propagation

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

where

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.

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]:

## 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

For the propagation cases that follow, we assume

To examine the effect of horizontal well misalignment angle on fracture propagation (Figure 17), first we assume the differential compression in

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 (*K*_{I}) 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 *K*_{I} for the initial fracture segment, where for the 89º misalignment case, the *K*_{I} 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-planar propagation that might also restrict width and hinder infectivity.

The time progression of the *K*_{I} variation along the top fracture front is displayed in Figure 20 for the case *K*_{I} at the initial fracture location (the injection location) grows very slowly in comparison to the curving wings of the fracture.

Although *K*_{I} is restricted in the misaligned portion of the fracture, Mode III or out of plane shear SIF(*K*_{III}) 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 *K*_{III} for varying fracture misalignment based on the simulation of Figure 19.

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

The magnitude of

## 6. Conclusion

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 penny-shaped cracks introduces more error in to the calculation, but the minimum and maximum SIF values can be accurately computed.

## Acknowledgments

Funding for this project is provided by RPSEA through the “Ultra-Deepwater and Unconventional 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 organizations, 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.