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
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
where
where,
2.2. Stress intensity factor computation
By knowing the crack tip element displacement discontinuities, KI, KII and KIII 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
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, F1, along the crack front y=b for various values of
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 (
The time progression of the
Although
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.
References
- 1.
Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of AmericaOkada Y 1992 82 2 1018 1040 - 2.
Fracture mechanics analysis of joints and veins. PhD Dissertation. Stanford University;Olson J. E 1991 - 3.
A Boundary Element Method for Two Dimensional Linear Elastic Fracture Analysis. International Journal of FractureCahng C. C Mear M. E 1995 74 219 - 4.
Theoretical displacement and stress near fracture in rock: with application to faults, joints, veins, dikes and solution surfaces. In: Atkinson BK. (ed.) Fracture Mechanics of Rock. London: Academic Press;Pollard D. D Segall P 1987 277 350 - 5.
Analysis of stress and strain near the end of a crack traversing a plate. Journal of Applied MechanicsIrwin G. R 1957 24 361 - 6.
On the finite element method in linear fracture mechanics. Engineering Fracture MechanicsChan S. K Tuba I. S Wilson W. K 1970 2 1 - 7.
Crack tip finite elements are unnecessary. International Journal for Numerical Methods in EngineeringHenshell R. D Shaw K. G 1975 9 495 - 8.
On the use of isoparametric finite elements in linear fracture mechanics, International Journal for Numerical Methods in EngineeringBarsoum R. S 1976 10 1 65 88 - 9.
Representation of singularities with isoparametric finite elements., International Journal for Numerical Methods in EngineeringBenzley S. E 1974 8 3 537 545 - 10.
Enhanced blending elements for XFEM applied to linear elastic fracture mechanics. International Journal for Numerical Methods in EngineeringTarancon J. E Vercher A Giner E Fuenmayor F. J 2009 77 1 126 148 - 11.
The optimal XFEM approximation for fracture analysis. IOP Conference Series: Materials Science and EngineeringJiang S Ying Z Du C 2010 - 12.
Re-consideration of Chen’s problem by finite difference method. Engineering Fracture MechanicsLin X Ballmann J 1993 44 5 735 739 - 13.
Numerical computation of dynamic stress intensity factors by a Lagrangian finite-deference method (THE HEMP CODE). Engineering Fracture MechanicsChen Y. M 1975 - 14.
Numerical Methods in Rock Mechanics. West Sussex: John Wiley & Sons Ltd;Pande G. N Beer G Williams J. R 1990 - 15.
An integral equation approach to boundary value problems of classical elastostatics. Quarterly of Applied MathematicsRizzo F. J 1967 25 83 - 16.
Solution of plane elasticity problems by the displacement discontinuity method. International Journal for Numerical Methods in EngineeringCrouch S. L 1976 10 2 301 343 - 17.
Boundary element methods in solid mechanics. London: George Allen & Unwin;Crouch S. L Starfield A. M 1983 - 18.
Stress intensity factor for curved cracks obtained displacement discontinuity method. International Journal of FractureSchultz R. A 1988 R31 R34. - 19.
Effect of mechanical interaction on the development of strike-slip faults with echelon patterns. Journal of Structural GeologyAydin A Schultz R. A 1990 12 1 123 129 - 20.
Higher-order functional variation displacement discontinuity elements. International Journal of Rock Mechanics and Mining Sciences & Geomechanics AbstractCrawford A. M Curran J. H 1982 19 3 143 148 - 21.
The displacement discontinuity method on the analysis of open cracks. MeccanicaScavia C 1991 26 1 27 32 - 22.
Stress intensity factors for cracks emanating from a triangular or square hole in an infinite plate by boundary elements. Engineering Failure AnalysisYan X 2005 12 3 362 375 - 23.
A higher order displacement discontinuity method for analysis of crack problems. International Journal of Rock Mechanics and Mining Sciences & Geomechanics AbstractShou K. J Crouch S. L 1995 32 1 49 55 - 24.
Numerical implementation of displacement discontinuity method and its application in hydraulic fracturing. Computer Methods in Applied Mechanics Engineering de-Dong C. Y Pater C. J 2001 - 25.
A fictitious stress and displacement discontinuity method for dynamic crack problems. In: Brebbia CA. (ed.) Boundary Element Method XVI. Southampton: Computational Mechanics Publications;Wen P. H Aliabadi M. H Rooke D. P 1994 469 476 - 26.
The geometry of echelon fractures in rock-implication from laboratory and numerical experiments. Journal of Structural GeologyThomas A. L Pollard D. D 1993 - 27.
D mechanical analysis of normal fault evolution and joint development in perturbed stress field around normal faults. PhD Dissertation. Stanford University;Kattenhorn S. A. A 1998 - 28.
Willemse EJM Pollard DD. Normal Fault growth: Evolution of tipline shapes and slip distribution. In: Lehner FK, Urai JL. (ed.) Aspects of Tectonic Faulting. Newyork: Springer;2000 193 226 - 29.
Fracture aperture, length and pattern geometry development under biaxial loading: a numerical study with application to natural, cross-jointed systems. Geological Society Special PublicationOlson J. E 2007 289 123 - 30.
Calculation of dyke trajectories from volcanic centers. Journal of Geophysical ResearchMériaux C Lister J. R 2002 B4) 2077-2087. - 31.
Mutlu O Pollard D. D 2006 A complementarity approach for modeling fractures. In: Yale D, Holtz S, Breeds C, and Ozbay U. (eds.) 50 years of Rock Mechanics-Landmarks and Future Challenges: The 41st U.S. Symposium on Rock Mechanics, Paper No: ARMA/USRMS06 1058 June 17-21 2006, Golden, CO, USA. - 32.
Closure of circular arc cracks under general loading: effects on stress intensity factors. International Journal of FractureRitz E Pollard D. D 2011 167 1 3 14 - 33.
Quantitative evaluation of fatigue strength of metals containing various small defects or cracks. Engineering Fracture MechanicsMurakami Y Endo M 1983 17 1 1 15 - 34.
Quantitative evaluation of effects of non-metallic inclusions on fatigue strength of high strength steels. I: Basic fatigue mechanism and evaluation of correlation between the fatigue fracture stress and the size and location of non-metallic inclusions. International Journal of FatigueMurakami Y Kodama S Konuma S 1989 11 5 291 298 - 35.
Irwin G. R 1962 Crack-extension force for a part-through crack in a plate. Journal of Applied Mechanics 1962;29 4 651 654 - 36.
Stress intensity factor of an elliptical crack and semi-elliptical crack in plates subjected to tension. International Journal of fractureNisitani H Murakami Y 1974 10 3 353 368 - 37.
Three dimensional crack analysis. International Journal of Solids and StructuresWeaver J 1977 13 4 321 330 - 38.
rectangular crack subjected to shear loading. International Journal of Solids and Structures 1982;Kassir M. K 1982 A Three-dimensional 18 12 1075 1082 - 39.
A rectangular crack in an infinite solid, a semi-infinite solid and a finite-thickness plate subjected to tension. International Journal of FractureIsida M Yoshida T Noguchi H 1991 52 79 - 40.
Variation of stress intensity factor along the front of a 3D rectangular crack by using a singular integral equation method. International Journal of FractureWang Q Noda N. A Honda M. A Chen M 2001 108 119 - 41.
Stress intensity factor for a three-dimensional rectangular crack. Journal of Applied MechanicsKassir M. K 1981 48 2 309 313 - 42.
Variation of the stress intensity factor along the front of aNoda N. A Kihara T. A 3 D rectangular crack subjected to mixed-mode load. Archive of Applied Mechanics2002 - 43.
Stress intensity factor for a plane crack under normal pressure. International Journal of FractureMastrojannis E. N Keer L. M Mura T. A 1979 15 3 101 114 - 44.
Sublinear scaling of fracture aperture versus length: An exception or the rule. Journal of Geophysical ResearchOlson J. E 2003 B9) 2413-2424. - 45.
An isolated Mode I three-dimensional planar crack: The stress intensity factor is independent to elastic constant. International Journal of FractureMear M. E Rodin G. J 2011 172 2 217 218 - 46.
The distribution of stress in the neighborhood of a crack in an elastic solid. Proceedings of the Royal Society of London. Series A, Mathematical and Physical SciencesSneddon I. N 1946 187 1009 229 260 - 47.
-Erdogan, F, Sih GC. On the crack extension in plates under plane loading and transverse shear. ASME Journal of Basic Engineering 1963;85 519-527 - 48.
Formation and interpretation of dilatant echelon cracks. Geological Society of America BulletinPollard D. D Segall P Delaney P. T 1982 93 1291 - 49.
Subcritical crack growth in geological materials. Journal of geophysical researchAtkinson B. K 1984 B6) 4077-4114 - 50.
Multi-fracture propagation modeling: Application to hydraulic fracturing in shales and tight gas sands. The 42st U.S. Symposium on Rock Mechanics and 2nd U.S.-canada Rock Mechanics Symposium, Paper No: ARMAOlson J. E 08 327 June 29-July 22008 San Francisco, CA, USA. - 51.
Fracturing from highly deviated and horizontal wells: Numerical analysis of non-planar fracture propagation. SPE Rock Mountain Regional/Low Permeability Reservoir Symposium, Paper No. SPE 29573, MarchOlson J. E 20 22 1995 Denver, CO, USA. - 52.
Sequential versus simultaneous multi-zone fracturing in horizontal wells: Insight from a non-planar, multi-frac numerical model. SPE Hydraulic Fracturing Technology Conference, Paper No: SPEOlson J. E Wu K 152602 PP, February 6-82012 The Woodland, TX, USA. - 53.
Analysis of minor fractures associated with joints and faulted joints. Journal of Structural GeologyCruikshank K. M Zhao G Johnson A. M 1991 13 8 865 886