In this paper we present a numerical model for hydraulic fracturing purposes. The rock formation is modelled as a poroelastic material based on Biot’s Theory. A fracture is represented in a discrete manner using the eXtended Finite Element Method (X-FEM). The fluid flow is governed by a local mass balance. This means that there is an equilibrium between the opening of the fracture, the tangential fluid flow, and the fluid leakage. The mass balance in the fracture is solved with a separate equation by including an additional degree of freedom for the pressure in the fracture. The fracture can grow in arbitrary directions by using an average stress criterion. We show a result of hydraulic fracture propagation for a 2D circular borehole. The fracture direction is consistent with the expected direction

Huyghe, J.M.R.J.

Remij, E.W.

Remmers, J.J.C.

Smeulders, D.M.J.

English

Pure OAI Repository

 Numerical modelling of hydraulic fracturingCitation for published version (APA):Remij, E. W., Remmers, J. J. C., Huyghe, J. M. R. J., & Smeulders, D. M. J. (2015). Numerical modelling ofhydraulic fracturing. In F. Oka, A. Murakami, R. Uzuoka, & S. Kimoto (Eds.), Computer Methods and RecentAdvances in Geomechanics - Proceedings of the 14th Int. Conference of International Association for ComputerMethods and Recent Advances in Geomechanics, IACMAG 2014 (pp. 1719-1723). Taylor & Francis/ Balkema.Document status and date:Published: 01/01/2015Document Version:Accepted manuscript including changes made at the peer-review stagePlease check the document version of this publication:• A submitted manuscript is the version of the article upon submission and before peer-review. There can beimportant differences between the submitted version and the official published version of record. Peopleinterested in the research are advised to contact the author for the final version of the publication, or visit theDOI to the publisher's website.• The final author version and the galley proof are versions of the publication after peer review.• The final published version features the final layout of the paper including the volume, issue and pagenumbers.Link to publicationGeneral rightsCopyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright ownersand it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.            • Users may download and print one copy of any publication from the public portal for the purpose of private study or research.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying the publication in the public portal.If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, pleasefollow below link for the End User Agreement:www.tue.nl/taverneTake down policyIf you believe that this document breaches copyright please contact us at:openaccess@tue.nlproviding details and we will investigate your claim.Download date: 04. Oct. 2023Numerical modelling of hydraulic fracturingE.W. Remij & J.J.C. Remmers & J.M. Huyghe & D.M.J. SmeuldersDepartment of Mechanical EngineeringEindhoven University of Technology, PO BOX 513, 5600 MB EindhovenABSTRACT: In this paper we present a numerical model for hydraulic fracturing purposes. The rockformation is modelled as a poroelastic material based on Biot’s Theory. A fracture is represented in adiscrete manner using the eXtended Finite Element Method (X-FEM). The fluid flow is governed by alocal mass balance. This means that there is an equilibrium between the opening of the fracture, thetangential fluid flow, and the fluid leakage. The mass balance in the fracture is solved with a separateequation by including an additional degree of freedom for the pressure in the fracture. The fracture cangrow in arbitrary directions by using an average stress criterion. We show a result of hydraulic fracturepropagation for a 2D circular borehole. The fracture direction is consistent with the expected direction.1 INTRODUCTIONHydraulic fracturing is the process in which a frac-ture propagates by applying a high pressure insidethe fracture. In geo-mechanics this process is usedto stimulate oil and gas reservoirs by fracturingthe underground rock formation. The induced frac-tures remain open, due to the addition of a prop-pant to the fracturing fluid, and therefore greatlyenhance the permeability in the formation. Numer-ical models can be used to obtain more insight inthe fracture process and may eventually be usedto optimize the fracture process.Boone & Ingraffea (1990) developed a numericalmodel based on the finite element method (FEM)for the poroelastic material where a cohesive zonedescription was used for the fracture. Using a meshadaptation scheme in a poroelastic FEM, Schref-fler and co-workers (Schrefler et al. 2006, Secchiet al. 2007, Secchi and Schrefler 2012) simulatedpropagating hydraulic fractures in a arbitrary di-rections. Carrier & Granet (2012) used a similarapproach with a priori placed interface elementsthat contained an additional degree of freedom forthe pressure in the fracture.The eXtended Finite Element Method (X-FEM)was used by Mohammadnejad & Khoei (2012) forcohesive crack growth in multiphase porous ma-terials. They also successfully applied their modelfor hydraulic fracturing simulations (Mohammad-nejad & Khoei 2013).X-FEM is a common technique in solid mechan-ics and has as an important advantage comparedto the previously mentioned fracture models; afracture can grow in arbitrary directions withoutthe need to remesh (Belytschko and Black 1999,Moës et al. 1999). By exploiting the partition ofunity property of finite element shape functions(Melenk & Babuška 1996), a fracture can be in-serted in a finite element mesh by adding addi-tional degrees of freedom to the nodes surround-ing the fracture (Figure 1). In this contributionwe present a similar X-FEM model as was usedby Mohammadnejad & Khoei (2013) but we alsoinclude the pressure in the fracture as a separatedegree of freedom. Such a formulation ensures thatall the pumping fluid goes into the fracture. Thenew model is refereed to as the Enhanced LocalPressure (ELP) model.Enhanced nodeRegular nodeΩ+Ω−Figure 1: A two dimensional finite element meshcrossed by a fracture represented by the black line.The black nodes surrounding the fracture are en-hanced with additional degrees of freedom. Thegrey elements therefore contain additional termsin the stiffness matrix and the force vector.2 KINEMATIC RELATIONSThe total displacement field of the solid skeletoncan, at any time t, be described by a continuousdisplacement field û(x, t) and a discontinuous dis-placement field ũ(x, t) (Belytschko & Black 1999,Moës, Dolbow, & Belytschko 1999)u(x, t) = û(x, t) +HΓd(x)ũ(x, t), (1)We can also write the pressure field as a sum of thecontinuous field p̂(x, t) and a discontinuous fieldp̃(x, t).p(x, t) = p̂(x, t) +HsΓd(x)p̃(x, t). (2)The functions HΓd and HsΓd are Heaviside stepfunctions defined over the fracture asHΓd ={1 if x ∈ Ω+0 if x ∈ Ω−, (3)HsΓd ={12if x ∈ Ω+−12if x ∈ Ω−. (4)With these formulations the pressure inside thefracture is undetermined. We therefore introduce athird pressure pd, representing the pressure in thefracture. The continuous fields are present in allthe finite element nodes. The discontinuous fieldsand the pressure in the fracture are only includedin the enhanced nodes (Figure 1).The pressure in the fracture is a one dimensionalvariable that can be used to describe the fluid massbalance in the fracture. Integrating the local massbalance givesq+Γd · nd − q−Γd· nd + [u̇]n + un〈δu̇sδs〉− (5)unδδs(kdδpdδs) +unMfṗd = 0,with q+Γd and q−Γdbeing the fluid flow from thefracture into formation for the fracture side of theΩ+ and the Ω− domain, respectively, [u̇]n denotingthe time derivative of the normal opening of thefracture, us being the shear opening of the fracture,and kd being the permeability in the fracture.Assuming a Couette flow, the latter is given bykd =u2n12µ, (6)where µ is the viscosity of the fluid, and un theopening of the fracture. The derivation of thisequation can be found in Irzal et al. (2013).Using the cohesive zone approach, the damagingof the material is described by a traction acting onthe fracture. This traction is coupled to the hydro-static pressure in the crack. Assuming continuityof stress from the formation to the fracture, we canwrite the local momentum balance asσ · nd = td − pdnd, (7)3 NUMERICAL FORMULATIONThe time discretization is performed with an im-plicit Euler scheme. The non-linear system of equa-tions is finally solved by using a Newton-Raphsoniterative procedure. We consider the porous mate-rial to be fully saturated with a fluid and subjectedto small variations in the displacement gradient.The bulk poroelasticity is based on Biot theory.The propagation time and the propagation di-rection of the fracture are based on an averagestress criterion. The stress values surrounding thetip of the cohesive zone are weighted with a Gaus-sian function (Jirásek 1998). The average stress isthen the weighted sum of stresses in these integra-tion points.xy~n ~sθηFigure 2: Schematic representation of a materialwith at the crack tip the global x-y coordinatesystem and the local coordinate system, describedwith a normal unit vector n and a tangential unitvector s.σav =nint∑i=1wiwtotσe,i with wtot =nint∑j=1wj. (8)Here nint is the number of integration points in thedomain, σe,i is the current effective stress state inintegration point i which has a weight factor widefined aswi =(2π)23l3ae−r2i2l2a , (9)with ri being the distance between the tip of thecohesive zone and the integration point ni, and labeing a length scale parameter defining how fastthe weight factor decays as a function of the dis-tance between the integration points and the co-hesive zone tip. Based on the average stress, anequivalent traction is calculated (Camacho & Or-tiz 1996)teq(θ) =√< tn >2 +1βt2s , (10)with< tn >={0 if tn ≤ 0tn if tn > 0, (11)where tn and ts respectively are the normal andshear tractiontn = nTσavn ts = sTσavn. (12)Here n is the normal vector and s is the tangentvector to an axis η which is rotated by an angle θwith respect to the x-axis (Figure 2). If the maxi-mum equivalent traction exceeds the yield strengthτult of the material the fracture is extended in thedirection of angle θ through one element. The frac-ture is assumed to open as a cleavage crack. Thetraction separation relation is given by (Camacho& Ortiz 1996)tn = tn0(1− vnvncr), (13)ts = ts0(1− vnvncr)sgn(vs).Here vn and vs are respectively the normal andsliding the displacement, sgn(·) is the Signum func-tion. The parameter vncr is the length of the fullydeveloped traction-free crack. This parameter de-pends on the fracture toughness Gc and the initialnormal traction tn0-101t n/tn00 1vn/vn,crFigure 3: The normalized tractions across the frac-ture as a function of the displacement jump in.vncr =2Gctn0. (14)The initial tractions, tn0 and ts0 , are taken equal tothe normal and shear tractions at the time of prop-agation, respectively. Self-contact of the fracture issimulated by using a penalty stiffness method. Adetailed description of the numerical implementa-tion can be found in the work (Remmers, de Borst,& Needleman 2008) and the usage of X-FEM inporous materials is given in (Kraaijeveld & Huyghe2011, Réthoré, Borst, & Abellan 2007).Numerical solutionAnalytical solution-3-2.5-2-1.5-1-0.500 2 4 6 8 10 12 14 16σrr[MPa]x [m](a) Radial stressNumerical solutionAnalytical solution-1.75-1.7-1.65-1.6-1.55-1.5-1.45-1.4-1.35-1.3-1.250 2 4 6 8 10 12 14 16σθθ[MPa]x [m](b) Circumferential stressFigure 4: The results of the numerical solution andanalytical solution for the Kirsch problem plottedover the x-axis.4 EXAMPLESThe performance of the numerical model is ad-dressed by considering a 2D borehole problem(Figure 5). The rock formation is assumed to be asquared specimen with a dimension of 10.0m. Thecircular hole has a radius of 1.0m. The confiningstress is taken as σ0=5.0[Mpa] in y-direction, andis taken twice as large in the x-direction. The ma-terial parameters are given in Table 1.Applying the confining stress leads to consolida-tion effects in the bulk material due to a poroelas-tic response. In order to obtain an initial equili-brated state we perform one time step with a pre-scribed zero pressure in all the nodes. After thisfirst step we release the pressure, place the initialfractures and start the fluid inflow q0.The initial stress field is compared with Kirschanalytical solution for a circular hole in an infi-nite elastic solid (Kirsch 1898). The radial stressand the circular stress over the x-axis are shownin Figure 4.Inintialfracturesσ0σ02σ02σ0p = 0p = 0yxFigure 5: Scheme of the borehole fracture problem.Table 1: Model Parameters used in the boreholefracture problem.Name Symbol ValueYoung’s modulus E 17.0 GPaPoisson’s ratio ν 0.2Ultimate strength τult 1.25 MPaFracture toughness Gc 120.0 Nm−1Fluid viscosity µ 1.0e−4 Pa · sFluid injection rate q0 0.05 ms−1Bulk modulus solid Ks 36.0 GPaBulk modulus fluid Kf 3.0 GPaBiot coefficient α 0.75The preferred propagation direction of a hy-draulic fracture is perpendicular to the minimumconfining stress (Hubbert & David 1972). In or-der to investigate if the numerical model can pre-dict this behaviour we also perform the simulationwith only the initial fracture in the y-direction.The time evolution of this fracture is shown inFigure 6a. If we consider both initial fractures, in-deed only the fracture that is perpendicular to theminimum confining stress propagates, as shown inFigure 6b.5 CONCLUSIONSThe performance of the ELP model is addressedin a 2D circular borehole problem. In the first caseonly the initial fracture perpendicular to the high-est confining stress is included in the simulation.This is the unfavourable direction of fracture prop-agation.a bFigure 6: Contour plots of the hydrostatic pressurefor the borehole problem at 3 different times. Thedeformed mesh is magnified 100 times. (a) A prop-agating fracture turns in the direction of the min-imum confining stress. (b) Only the fracture per-pendicular to the minimum confining stress propa-gates. The second fracture in the y-direction is notvisible because it is closed.The propagated fracture indeed turns in the direc-tion of the highest confining stress. In the secondcase also the fracture that is already in the di-rection of the highest confining stress is included.Now only the second fracture grows while the frac-ture that is initially perpendicular to the highestconfining stress remains closed. This behaviour isconsistent with findings in the literature (Hubbert& David 1972).ACKNOWLEDGEMENTThis research project is supported by the DutchTKI Gas Foundation.REFERENCESBelytschko, T. & T. Black (1999). Elastic crack growth infinite elements with minimal remeshing. Internationaljournal for numerical methods in engineering 45 (5),601–620.Boone, T. & A. Ingraffea (1990). A numerical procedurefor simulation of hydraulically-driven fracture propaga-tion in poroelastic media. International Journal for Nu-merical and Analytical Methods in Geomechanics 14 (1),27–47.Camacho, G. & M. Ortiz (1996). Computational modellingof impact damage in brittle materials. InternationalJournal of solids and structures 33 (20), 2899–2938.Carrier, B. & S. Granet (2012). Numerical modeling of hy-draulic fracture problem in permeable medium using co-hesive zone model. Engineering Fracture Mechanics 79,312–328.Hubbert, M. & G. W. David (1972). Mechanics of hydraulicfracturing. US Geological Survey 210, 153–168.Irzal, F., J. Remmers, J. Huyghe, & R. de Borst (2013). Alarge deformation formulation for fluid flow in a progres-sively fracturing porous material. Computer Methods inApplied Mechanics and Engineering , 29–37.Jirásek, M. (1998). Embedded crack models for concretefracture. In Computational Modelling of Concrete Struc-tures, EURO C-98, Volume 1, pp. 291–300.Kirsch, G. (1898). Die Theorie der Elastizität und dieBedürfnisse der Festigkeitslehre. Zeitschrift des Vereinesdeutscher Ingenieure 42, 797–807.Kraaijeveld, F. & J. Huyghe (2011). Propagating cracks insaturated ionized porous media. Multiscale Methods inComputational Mechanics, 425–442.Melenk, J. & I. Babuška (1996). The partition of unityfinite element method: basic theory and applications.Computer methods in applied mechanics and engineer-ing 139 (1), 289–314.Moës, N., J. Dolbow, & T. Belytschko (1999). A finite el-ement method for crack growth without remeshing. In-ternational journal for numerical methods in engineer-ing 46 (1), 131–150.Mohammadnejad, T. & A. Khoei (2012). Hydro-mechanicalmodeling of cohesive crack propagation in multiphaseporous media using the extended finite element method.International Journal for Numerical and AnalyticalMethods in Geomechanics 37 (10), 1247–1279.Mohammadnejad, T. & A. Khoei (2013). An extended finiteelement method for hydraulic fracture propagation indeformable porous media with the cohesive crack model.Finite Elements in Analysis and Design 73, 77–95.Remmers, J., R. de Borst, & A. Needleman (2008). Thesimulation of dynamic crack propagation using the co-hesive segments method. Journal of the Mechanics andPhysics of Solids 56 (1), 70–92.Réthoré, J., R. Borst, & M. Abellan (2007). A two-scaleapproach for fluid flow in fractured porous media. In-ternational Journal for Numerical Methods in Engineer-ing 71 (7), 780–800.Schrefler, B., S. Secchi, & L. Simoni (2006). On adaptiverefinement techniques in multi-field problems includingcohesive fracture. Computer methods in applied mechan-ics and engineering 195 (4), 444–461.Secchi, S., L. Simoni, & B. Schrefler (2007). Mesh adapta-tion and transfer schemes for discrete fracture propaga-tion in porous materials. International Journal for Nu-merical and Analytical Methods in Geomechanics 31 (2),331–345.Secchi, S. & B. Schrefler (2012). A method for 3D hy-draulic fracturing simulation. International journal offracture 178 (1-2), 245–258.

Numerical modelling of hydraulic fracturing

Remij, EW Ernst

Remmers, JJC Joris

Huyghe, JMRJ Jacques

Smeulders, DMJ David

Repository TU/e

 Numerical modelling of hydraulic fracturingCitation for published version (APA):Remij, E. W., Remmers, J. J. C., Huyghe, J. M. R. J., & Smeulders, D. M. J. (2015). Numerical modelling ofhydraulic fracturing. In F. Oka, A. Murakami, R. Uzuoka, & S. Kimoto (Eds.), Computer Methods and RecentAdvances in Geomechanics - Proceedings of the 14th Int. Conference of International Association for ComputerMethods and Recent Advances in Geomechanics, IACMAG 2014 (pp. 1719-1723). Taylor & Francis/ Balkema.Document status and date:Published: 01/01/2015Document Version:Accepted manuscript including changes made at the peer-review stagePlease check the document version of this publication:• A submitted manuscript is the version of the article upon submission and before peer-review. There can beimportant differences between the submitted version and the official published version of record. Peopleinterested in the research are advised to contact the author for the final version of the publication, or visit theDOI to the publisher's website.• The final author version and the galley proof are versions of the publication after peer review.• The final published version features the final layout of the paper including the volume, issue and pagenumbers.Link to publicationGeneral rightsCopyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright ownersand it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.            • Users may download and print one copy of any publication from the public portal for the purpose of private study or research.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying the publication in the public portal.If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, pleasefollow below link for the End User Agreement:www.tue.nl/taverneTake down policyIf you believe that this document breaches copyright please contact us at:openaccess@tue.nlproviding details and we will investigate your claim.Download date: 17. Nov. 2023Numerical modelling of hydraulic fracturingE.W. Remij & J.J.C. Remmers & J.M. Huyghe & D.M.J. SmeuldersDepartment of Mechanical EngineeringEindhoven University of Technology, PO BOX 513, 5600 MB EindhovenABSTRACT: In this paper we present a numerical model for hydraulic fracturing purposes. The rockformation is modelled as a poroelastic material based on Biot’s Theory. A fracture is represented in adiscrete manner using the eXtended Finite Element Method (X-FEM). The fluid flow is governed by alocal mass balance. This means that there is an equilibrium between the opening of the fracture, thetangential fluid flow, and the fluid leakage. The mass balance in the fracture is solved with a separateequation by including an additional degree of freedom for the pressure in the fracture. The fracture cangrow in arbitrary directions by using an average stress criterion. We show a result of hydraulic fracturepropagation for a 2D circular borehole. The fracture direction is consistent with the expected direction.1 INTRODUCTIONHydraulic fracturing is the process in which a frac-ture propagates by applying a high pressure insidethe fracture. In geo-mechanics this process is usedto stimulate oil and gas reservoirs by fracturingthe underground rock formation. The induced frac-tures remain open, due to the addition of a prop-pant to the fracturing fluid, and therefore greatlyenhance the permeability in the formation. Numer-ical models can be used to obtain more insight inthe fracture process and may eventually be usedto optimize the fracture process.Boone & Ingraffea (1990) developed a numericalmodel based on the finite element method (FEM)for the poroelastic material where a cohesive zonedescription was used for the fracture. Using a meshadaptation scheme in a poroelastic FEM, Schref-fler and co-workers (Schrefler et al. 2006, Secchiet al. 2007, Secchi and Schrefler 2012) simulatedpropagating hydraulic fractures in a arbitrary di-rections. Carrier & Granet (2012) used a similarapproach with a priori placed interface elementsthat contained an additional degree of freedom forthe pressure in the fracture.The eXtended Finite Element Method (X-FEM)was used by Mohammadnejad & Khoei (2012) forcohesive crack growth in multiphase porous ma-terials. They also successfully applied their modelfor hydraulic fracturing simulations (Mohammad-nejad & Khoei 2013).X-FEM is a common technique in solid mechan-ics and has as an important advantage comparedto the previously mentioned fracture models; afracture can grow in arbitrary directions withoutthe need to remesh (Belytschko and Black 1999,Moës et al. 1999). By exploiting the partition ofunity property of finite element shape functions(Melenk & Babuška 1996), a fracture can be in-serted in a finite element mesh by adding addi-tional degrees of freedom to the nodes surround-ing the fracture (Figure 1). In this contributionwe present a similar X-FEM model as was usedby Mohammadnejad & Khoei (2013) but we alsoinclude the pressure in the fracture as a separatedegree of freedom. Such a formulation ensures thatall the pumping fluid goes into the fracture. Thenew model is refereed to as the Enhanced LocalPressure (ELP) model.Enhanced nodeRegular nodeΩ+Ω−Figure 1: A two dimensional finite element meshcrossed by a fracture represented by the black line.The black nodes surrounding the fracture are en-hanced with additional degrees of freedom. Thegrey elements therefore contain additional termsin the stiffness matrix and the force vector.2 KINEMATIC RELATIONSThe total displacement field of the solid skeletoncan, at any time t, be described by a continuousdisplacement field û(x, t) and a discontinuous dis-placement field ũ(x, t) (Belytschko & Black 1999,Moës, Dolbow, & Belytschko 1999)u(x, t) = û(x, t) +HΓd(x)ũ(x, t), (1)We can also write the pressure field as a sum of thecontinuous field p̂(x, t) and a discontinuous fieldp̃(x, t).p(x, t) = p̂(x, t) +HsΓd(x)p̃(x, t). (2)The functions HΓd and HsΓd are Heaviside stepfunctions defined over the fracture asHΓd ={1 if x ∈ Ω+0 if x ∈ Ω−, (3)HsΓd ={12if x ∈ Ω+−12if x ∈ Ω−. (4)With these formulations the pressure inside thefracture is undetermined. We therefore introduce athird pressure pd, representing the pressure in thefracture. The continuous fields are present in allthe finite element nodes. The discontinuous fieldsand the pressure in the fracture are only includedin the enhanced nodes (Figure 1).The pressure in the fracture is a one dimensionalvariable that can be used to describe the fluid massbalance in the fracture. Integrating the local massbalance givesq+Γd · nd − q−Γd· nd + [u̇]n + un〈δu̇sδs〉− (5)unδδs(kdδpdδs) +unMfṗd = 0,with q+Γd and q−Γdbeing the fluid flow from thefracture into formation for the fracture side of theΩ+ and the Ω− domain, respectively, [u̇]n denotingthe time derivative of the normal opening of thefracture, us being the shear opening of the fracture,and kd being the permeability in the fracture.Assuming a Couette flow, the latter is given bykd =u2n12µ, (6)where µ is the viscosity of the fluid, and un theopening of the fracture. The derivation of thisequation can be found in Irzal et al. (2013).Using the cohesive zone approach, the damagingof the material is described by a traction acting onthe fracture. This traction is coupled to the hydro-static pressure in the crack. Assuming continuityof stress from the formation to the fracture, we canwrite the local momentum balance asσ · nd = td − pdnd, (7)3 NUMERICAL FORMULATIONThe time discretization is performed with an im-plicit Euler scheme. The non-linear system of equa-tions is finally solved by using a Newton-Raphsoniterative procedure. We consider the porous mate-rial to be fully saturated with a fluid and subjectedto small variations in the displacement gradient.The bulk poroelasticity is based on Biot theory.The propagation time and the propagation di-rection of the fracture are based on an averagestress criterion. The stress values surrounding thetip of the cohesive zone are weighted with a Gaus-sian function (Jirásek 1998). The average stress isthen the weighted sum of stresses in these integra-tion points.xy~n ~sθηFigure 2: Schematic representation of a materialwith at the crack tip the global x-y coordinatesystem and the local coordinate system, describedwith a normal unit vector n and a tangential unitvector s.σav =nint∑i=1wiwtotσe,i with wtot =nint∑j=1wj. (8)Here nint is the number of integration points in thedomain, σe,i is the current effective stress state inintegration point i which has a weight factor widefined aswi =(2π)23l3ae−r2i2l2a , (9)with ri being the distance between the tip of thecohesive zone and the integration point ni, and labeing a length scale parameter defining how fastthe weight factor decays as a function of the dis-tance between the integration points and the co-hesive zone tip. Based on the average stress, anequivalent traction is calculated (Camacho & Or-tiz 1996)teq(θ) =√< tn >2 +1βt2s , (10)with< tn >={0 if tn ≤ 0tn if tn > 0, (11)where tn and ts respectively are the normal andshear tractiontn = nTσavn ts = sTσavn. (12)Here n is the normal vector and s is the tangentvector to an axis η which is rotated by an angle θwith respect to the x-axis (Figure 2). If the maxi-mum equivalent traction exceeds the yield strengthτult of the material the fracture is extended in thedirection of angle θ through one element. The frac-ture is assumed to open as a cleavage crack. Thetraction separation relation is given by (Camacho& Ortiz 1996)tn = tn0(1− vnvncr), (13)ts = ts0(1− vnvncr)sgn(vs).Here vn and vs are respectively the normal andsliding the displacement, sgn(·) is the Signum func-tion. The parameter vncr is the length of the fullydeveloped traction-free crack. This parameter de-pends on the fracture toughness Gc and the initialnormal traction tn0-101t n/tn00 1vn/vn,crFigure 3: The normalized tractions across the frac-ture as a function of the displacement jump in.vncr =2Gctn0. (14)The initial tractions, tn0 and ts0 , are taken equal tothe normal and shear tractions at the time of prop-agation, respectively. Self-contact of the fracture issimulated by using a penalty stiffness method. Adetailed description of the numerical implementa-tion can be found in the work (Remmers, de Borst,& Needleman 2008) and the usage of X-FEM inporous materials is given in (Kraaijeveld & Huyghe2011, Réthoré, Borst, & Abellan 2007).Numerical solutionAnalytical solution-3-2.5-2-1.5-1-0.500 2 4 6 8 10 12 14 16σrr[MPa]x [m](a) Radial stressNumerical solutionAnalytical solution-1.75-1.7-1.65-1.6-1.55-1.5-1.45-1.4-1.35-1.3-1.250 2 4 6 8 10 12 14 16σθθ[MPa]x [m](b) Circumferential stressFigure 4: The results of the numerical solution andanalytical solution for the Kirsch problem plottedover the x-axis.4 EXAMPLESThe performance of the numerical model is ad-dressed by considering a 2D borehole problem(Figure 5). The rock formation is assumed to be asquared specimen with a dimension of 10.0m. Thecircular hole has a radius of 1.0m. The confiningstress is taken as σ0=5.0[Mpa] in y-direction, andis taken twice as large in the x-direction. The ma-terial parameters are given in Table 1.Applying the confining stress leads to consolida-tion effects in the bulk material due to a poroelas-tic response. In order to obtain an initial equili-brated state we perform one time step with a pre-scribed zero pressure in all the nodes. After thisfirst step we release the pressure, place the initialfractures and start the fluid inflow q0.The initial stress field is compared with Kirschanalytical solution for a circular hole in an infi-nite elastic solid (Kirsch 1898). The radial stressand the circular stress over the x-axis are shownin Figure 4.Inintialfracturesσ0σ02σ02σ0p = 0p = 0yxFigure 5: Scheme of the borehole fracture problem.Table 1: Model Parameters used in the boreholefracture problem.Name Symbol ValueYoung’s modulus E 17.0 GPaPoisson’s ratio ν 0.2Ultimate strength τult 1.25 MPaFracture toughness Gc 120.0 Nm−1Fluid viscosity µ 1.0e−4 Pa · sFluid injection rate q0 0.05 ms−1Bulk modulus solid Ks 36.0 GPaBulk modulus fluid Kf 3.0 GPaBiot coefficient α 0.75The preferred propagation direction of a hy-draulic fracture is perpendicular to the minimumconfining stress (Hubbert & David 1972). In or-der to investigate if the numerical model can pre-dict this behaviour we also perform the simulationwith only the initial fracture in the y-direction.The time evolution of this fracture is shown inFigure 6a. If we consider both initial fractures, in-deed only the fracture that is perpendicular to theminimum confining stress propagates, as shown inFigure 6b.5 CONCLUSIONSThe performance of the ELP model is addressedin a 2D circular borehole problem. In the first caseonly the initial fracture perpendicular to the high-est confining stress is included in the simulation.This is the unfavourable direction of fracture prop-agation.a bFigure 6: Contour plots of the hydrostatic pressurefor the borehole problem at 3 different times. Thedeformed mesh is magnified 100 times. (a) A prop-agating fracture turns in the direction of the min-imum confining stress. (b) Only the fracture per-pendicular to the minimum confining stress propa-gates. The second fracture in the y-direction is notvisible because it is closed.The propagated fracture indeed turns in the direc-tion of the highest confining stress. In the secondcase also the fracture that is already in the di-rection of the highest confining stress is included.Now only the second fracture grows while the frac-ture that is initially perpendicular to the highestconfining stress remains closed. This behaviour isconsistent with findings in the literature (Hubbert& David 1972).ACKNOWLEDGEMENTThis research project is supported by the DutchTKI Gas Foundation.REFERENCESBelytschko, T. & T. Black (1999). Elastic crack growth infinite elements with minimal remeshing. Internationaljournal for numerical methods in engineering 45 (5),601–620.Boone, T. & A. Ingraffea (1990). A numerical procedurefor simulation of hydraulically-driven fracture propaga-tion in poroelastic media. International Journal for Nu-merical and Analytical Methods in Geomechanics 14 (1),27–47.Camacho, G. & M. Ortiz (1996). Computational modellingof impact damage in brittle materials. InternationalJournal of solids and structures 33 (20), 2899–2938.Carrier, B. & S. Granet (2012). Numerical modeling of hy-draulic fracture problem in permeable medium using co-hesive zone model. Engineering Fracture Mechanics 79,312–328.Hubbert, M. & G. W. David (1972). Mechanics of hydraulicfracturing. US Geological Survey 210, 153–168.Irzal, F., J. Remmers, J. Huyghe, & R. de Borst (2013). Alarge deformation formulation for fluid flow in a progres-sively fracturing porous material. Computer Methods inApplied Mechanics and Engineering , 29–37.Jirásek, M. (1998). Embedded crack models for concretefracture. In Computational Modelling of Concrete Struc-tures, EURO C-98, Volume 1, pp. 291–300.Kirsch, G. (1898). Die Theorie der Elastizität und dieBedürfnisse der Festigkeitslehre. Zeitschrift des Vereinesdeutscher Ingenieure 42, 797–807.Kraaijeveld, F. & J. Huyghe (2011). Propagating cracks insaturated ionized porous media. Multiscale Methods inComputational Mechanics, 425–442.Melenk, J. & I. Babuška (1996). The partition of unityfinite element method: basic theory and applications.Computer methods in applied mechanics and engineer-ing 139 (1), 289–314.Moës, N., J. Dolbow, & T. Belytschko (1999). A finite el-ement method for crack growth without remeshing. In-ternational journal for numerical methods in engineer-ing 46 (1), 131–150.Mohammadnejad, T. & A. Khoei (2012). Hydro-mechanicalmodeling of cohesive crack propagation in multiphaseporous media using the extended finite element method.International Journal for Numerical and AnalyticalMethods in Geomechanics 37 (10), 1247–1279.Mohammadnejad, T. & A. Khoei (2013). An extended finiteelement method for hydraulic fracture propagation indeformable porous media with the cohesive crack model.Finite Elements in Analysis and Design 73, 77–95.Remmers, J., R. de Borst, & A. Needleman (2008). Thesimulation of dynamic crack propagation using the co-hesive segments method. Journal of the Mechanics andPhysics of Solids 56 (1), 70–92.Réthoré, J., R. Borst, & M. Abellan (2007). A two-scaleapproach for fluid flow in fractured porous media. In-ternational Journal for Numerical Methods in Engineer-ing 71 (7), 780–800.Schrefler, B., S. Secchi, & L. Simoni (2006). On adaptiverefinement techniques in multi-field problems includingcohesive fracture. Computer methods in applied mechan-ics and engineering 195 (4), 444–461.Secchi, S., L. Simoni, & B. Schrefler (2007). Mesh adapta-tion and transfer schemes for discrete fracture propaga-tion in porous materials. International Journal for Nu-merical and Analytical Methods in Geomechanics 31 (2),331–345.Secchi, S. & B. Schrefler (2012). A method for 3D hy-draulic fracturing simulation. International journal offracture 178 (1-2), 245–258.

https://pure.tue.nl/ws/files/4003655/59396300220197.pdf

Numerical modelling of hydraulic fracturing

Abstract

Similar works

Full text

Available Versions

Pure OAI Repository

Repository TU/e

Pure OAI Repository