Stress intensity factors
Abstract
The numerical modeling of hydraulic fractures in unconventional reservoirs presents significant challenges for field applications. There remains a need for accurate models that field personnel can use, yet remains consistent to the underlying physics of the problem [1]. For numerical simulations, several authors have considered a number of issues: the coupling between fracture mechanics and fluid dynamics in the fracture [2], fracture interaction [35], proppant transport [6], and others [79]. However, the available literature within the oil and gas industry often ignores the importance of the crack tip in modeling applications developed for engineering design. The importance of accurate modeling of the stress induced near the crack tip is likely critical in complex geological reservoirs where multiple propagating crack tips are interacting with natural fractures. This study investigates the influence of various boundary element numerical techniques on the accuracy of the calculated stress intensity factor near the crack tip and on the fracture profile, in general. The work described here is a part of a longterm project in the development of more accurate and efficient numerical simulations for field engineering applications.
Keywords
 displacement discontinuity method
 higher order elements
 crack tip elements
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 [35], proppant transport [6], and others [79]. Further, there have been specific codes developed to model complex fracture network development [1416]. Nevertheless, 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 longterm project in the development of more accurate and efficient numerical simulations for field engineering applications.
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 higherorder elements. The higherorder elements use a quadratic variation of displacement discontinuity, and are based on the use of three collocation points over a threeelement 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 neartip displacement profiles exhibited by various regimes in hydraulic fracture propagation (e.g., ViscosityDominated or ToughnessDominated [22,23]). As others have shown, the accuracy of tip asymptote is critical in characterizing the stresses in the neartip region of a propagating fracture [1].
We examined the numerically derived stress intensity factor for three crack geometries with and without higherorder 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 higherorder or specialized boundary elements [2427]. 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 [2427]. For example, refinements with respect to the quarterpoint method are found in Gray
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 higherorder 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 combinations of the presented methods to known solutions of various crack geometries for an estimation 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 2D displacements and stresses at a point
where
and where the displacement discontinuity components are defined as:
For constant displacement elements (i.e. constant
where
For simplicity, the derivatives of
However, DDM with a constant DD can’t accurately calculate the stresses and displacements of the area closer than about one elementlength distance from a boundary [12]. To improve the accuracy of calculations in close proximity to the boundaries, Crawford
To improve the accuracy of DDM without sacrificing the number of degrees of freedom of the overall system, a new higherorder elements method was suggested by Shou
This study uses Shou
3. Higherorder elements of displacement discontinuity method
Higherorder elements, as formulated by Shou
where
Combining Equations (13) through (16) with Equations (6) and (7) gives the following simplified expressions:
where the subscript
which can be expressed in terms of constant, linear, and quadratic kernels
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
where
4. Crack tip elements of displacement discontinuity method
In addition to advanced elements, Shou
4.1. Quarter element method
Shou
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
4.2. Square root element method
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.
where
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 higherorder 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 [2527]. This research uses following elastic properties:
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
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 midpoint 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 higherorder elements), CDDCE (constant DD with the quarter element method), CDDCT (constant DD with the square root element method), HDDCE (higherorder elements with the quarter element method), and HDDCT (higherorder elements with the square root element method).
Equation (23) is the analytical solution of the dimensionless half width under a constant pressure [29].
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 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.
To evaluate the perturbed stress state due to the presence of a pressurized crack we use [30]
Equation (24) provides a dimensionless stress
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
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]
Equation (25) expresses a dimensionless stress
Table 1 shows the calculated stress intensity factor, along with the ratio to the analytical solution. For the mode I (or K_{I}), CDD shows the biggest error while CDDCE, HDDCE and HDDCT give around 1% errors. Obviously, the mode II (or K_{II}) 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].
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 
5.2. A slanted straight crack
While a straight pressurized crack shows zero K_{II}, a slanted straight crack under a uniform tension can show a variable K_{II} depending on the angle of incidence to the applied tension. Figure 11 illustrates the crack geometry. The stress intensity factors are calculated by [19]
Equation (26) expresses the analytical solution of the stress intensity factors [19]. The uniform tension is
Figure 12 shows the dimensionless K_{I} (i.e.
5.3. A circular arc crack
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
where
Figures (15) and (16) show the values of the calculated K_{I} and K_{II} 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 K_{I} and K_{II} is
Figures (17) and (18) show the normalized K_{I} and K_{II} stress intensity factors as a function of the circular crack angle, respectively.
Based on the results presented in this section, numerical calculations of K_{I} and K_{II} 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.
The computer specifications used in this work are as follows: CPUIntel^{®} Xeon W3670 @ 3.2GHZ, installed memory (RAM)24 gigabytes, OS 64bit 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 higherorder elements and/or special tip elements. Figure 19, shows the normalized stress (the ratio of
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
Overall results show that CDD gives prominent errors of calculations of stresses, displacements, and stress intensity factor compared to the other methods. Particularly, when approaching 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 higherorder 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 K_{I} and K_{II} is more important than for simple fracture geometries. Combination of higherorder 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 higherorder 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 higherorder elements coupled with specialized crack tip elements requires more complex geometries than presented here.
Nomenclature
Displacement  
Stress  
Young’s modulus  
Poisson’s ratio  
Shear modulus  
Plane strain, 

Fracture halflength or half height  
Fracture length or height 
Acknowledgments
The authors thank Baker Hughes for supporting this research and for permission to publish this paper, Randy LaFollette, for sincere advice and encouragement, and Russell Maharidge, PhD for considerate reviews.References
 1.
Computer simulation of hydraulic fractures. International Journal Of Rock Mechanics And Mining SciencesAdachi A Siebrits E Peirce A Desroches J 2007 44 5 739 757  2.
Simulation of Fluid Flow in Hydraulic Fracturing: Implications for 3D Propagation. SPE Production EngineeringNaceur K. B Thiercelin M Touboul E 1990 5 2 133 141 SPE 16032  3.
Multiple fracture model. In: Jeffrey, McLennan, editors. Proceedings of the three dimensional and advance hydraulic fracture modeling. Workshop held at 4th North American rock mechanics symposium, Seattle, July 2931,Germanovich L Astakhov D 2000 45 70  4.
Parameters Affecting the Interaction Among Closely Spaced Hydraulic Fractures. SPE Journal;17(1):Bunger A. P Zhang X Jeffrey R. G 292 306 SPE140426PA  5.
Sequential vs. Simultaneous Multizone Fracturing in Horizontal Wells: Insights From a NonPlanar, Multifrac Numerical Model. In: SPE Hydraulic Fracturing Technology Conference. The Woodlands, Texas, USA: Society of Petroleum Engineers. SPEOlson J. E Wu K 152602 M  6.
On numerical solutions of hyperbolic proppant transport problems. In: Proceedings of the 10th international conference on hyperbolic problems: theory, numerics and applications, Osaka, Japan, SeptemberGu H Siebrits E 13 17 2004  7.
A fixed grid algorithm for simulating the propagation of a shallow hydraulic fracture with a fluid lag. International Journal for Numerical and Analytical Methods in Geomechanics;35(5):602.Gordeliy E Detournay E  8.
Deflection and propagation of fluiddriven fractures at frictional bedding interfaces: A numerical investigation. Journal of Structural GeologyZhang X Jeffrey R. G Thiercelin M 2007  9.
A comparative investigation of microflaw models fog the simulation of brittle fracture in rock. Computational MechanicsSellers E Napier J 1997  10.
A higher order displacement discontinuity method for analysis of crack problems. International Journal of Rock Mechanics and Mining Science & Geomechanics AbstractsShou K. J Crouch S. L 1995  11.
Thirty Years of Gas Shale Fracturing: What Have We Learned? In: SPE Annual Technical Conference and Exhibition. Florence, Italy: Society of Petroleum Engineers. SPE 133456King G. E  12.
Boundary Element Methods in Solid Mechanics: With Applications in Rock Mechanics and Geological Engineering: George Allen & Unwin; 1 edition,Crouch S. L Starfield A. M 1983  13.
Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution. International Journal for Numerical Methods in EngineeringCrouch S. L 1976  14.
Predicting fracture swarms the influence of subcritical crack growth and the cracktip process zone on joint spacing in rock. Geological Society, London, Special PublicationsOlson J. E 2004 231 1 73 88  15.
Modeling of Interaction of Hydraulic Fractures in Complex Fracture Networks. In: SPE Hydraulic Fracturing Technology Conference. The Woodlands, Texas, USA: Society of Petroleum Engineers. SPE152052Wu R Kresse O Weng X Cohen Ce Gu H  16.
Modeling of HydraulicFractureNetwork Propagation in a Naturally Fractured Formation. SPE Production & Operations;26(4):Weng X Kresse O Cohen CE Wu R Gu H 368 380 SPE140253  17.
Napier JAL. Modelling of fracturing near deep level gold mine excavations using a displacement discontinuity approach. Mechanics of jointed and faulted rock. Rossmanith (ed), Balkema: Rotterdam:709 716 1990  18.
Napier JAL. A Spectral Multipole Method For Efficient Solution Of LargeScale BoundaryElement Models In Elastostatics. International Journal For Numerical Methods In EngineeringPeirce A. P 1995 38 23 4009 4034  19.
The Stress Analysis of Cracks Handbook. 3 ed: American Society of Mechanical Engineering;Tada H Paris P. C Irwin G. R 2000  20.
A G2 constant displacement discontinuity element for analysis of crack problems. Computational Mechanics;Exadaktylos G Xiroudakis G 45 4 245 261  21.
Napier JAL. A twodimensional linear variation displacement discontinuity method for threelayered elastic media. International Journal of Rock Mechanics and Mining SciencesShou KJ 1999  22.
Near tip processes of a fluiddriven fracture. ASME J Appl MechGaragash D. I Detournay E 2000 67 183 92  23.
Detournay E 2004 Propagation regimes of fluid driven fractures in impermeable rocks. Int. J. Geomechanics4 1 1 11  24.
Higherorder functional variation displacement discontinuity elements. International Journal of Rock Mechanics and Mining Sciences & Geomechanics AbstractsCrawford A. M Curran J. H 1982  25.
Dong C. Y De Pater C. J 2001 Numerical implementation of displacement discontinuity method and its application in hydraulic fracturing. Computer Methods in Applied Mechanics and Engineering;191:745.  26.
Improved quarterpoint crack tip element. Engineering Fracture MechanicsGray L. J Phan A. V Paulino G. H Kaplan T 2003 70 2 269 283  27.
A special crack tip displacement discontinuity element. Mechanics Research CommunicationsYan X. Q 2004 31 6 651 659  28.
Theoretical displacements and stresses near fractures in rock: with applications to faults, joints, veins, dikes and solution surfaces. In: ATKINSON, B K. Fracture Mechanics of Rock. Academic Press, London,Pollard D Segall P 277 350 1987  29.
Hydraulic Fracture Mechanics. 1 ed: Wiley;Valkó P Economides M. J 1995  30.
The Distribution Of Stress In The Neighbourhood Of A Crack In An Elastic Solid. Proceedings Of The Royal Society Of London Series AMathematical And Physical SciencesSneddon I. N 1946 187 1009 229 260