Open access peer-reviewed chapter

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

By Zuorong Chen

Submitted: October 16th 2012Reviewed: March 1st 2013Published: May 17th 2013

DOI: 10.5772/56287

Downloaded: 3990


A new finite element has been implemented in ABAQUS to incorporate the extended finite element method (XFEM) for the solution of hydraulic fracture problems. The proposed element includes the desired aspects of the XFEM so as to model crack propagation without explicit remeshing. In addition, the fluid pressure degrees of freedom have been defined on the element to describe the fluid flow within the crack and its contribution to the crack deformation. Thus the fluid flow and resulting crack propagation are fully coupled in a natural way and are solved simultaneously. Verification of the element has been made by comparing the finite element results with the analytical solutions available in the literature.

1. Introduction

Hydraulic fracturing is a powerful technology for enhancing conventional petroleum production. It is playing a central role in fast growing development of unconventional gas and geothermal energy. The fully 3-D numerical simulation of the hydraulic fracturing process is of great importance to understand the complex, multiscale mechanics of hydraulic fracturing, to the efficient application of this technology, and to develop innovative, advanced hydraulic fracture technologies for unconventional gas production. The accurate numerical simulation of hydraulic fracture growth remains a significant challenge because of the strong nonlinear coupling between the viscous flow of fluid inside the fracture and fracture propagation (a moving boundary), complicated by the need to consider interactions with existing natural fractures and with rock layers with different properties.

Great effort has been devoted to the numerical simulation of hydraulic fractures with the first 3D modelling efforts starting in the late 1970s [1-2]. Significant progress has been made in developing 2-D and 3-D numerical hydraulic fracture models [3-15]. Boundary integral equation methods or displacement discontinuity techniques have generally been employed to investigate the propagation of simple hydraulic fractures such as radial or plane-strain fractures in a homogeneous, infinite or semi-infinite elastic medium where the appropriate fundamental solutions are available. The finite element method has been used and is particularly useful in modelling the hydraulic fracture propagation in inhomogeneous rocks which may include nonlinear mechanical properties and may be subjected to complex boundary conditions. However, the standard finite element model requires remeshing after every crack propagation step and the mesh has to conform exactly to the fracture geometry as the fracture propagates, and thus is computationally expensive.

By adding special enriched shape functions in conjunction with additional degrees of freedom to the standard finite element approximation within the framework of partition of unity, the extended finite element method [16-17] (XFEM) overcomes the inherent drawbacks associated with use of the conventional finite element methods and enables the crack to be represented without explicitly meshing crack surfaces, and so the crack geometry is completely independent of the mesh and remeshing is not required, allowing for the convenient simulation of the fracture propagation. The XFEM has been employed to investigate the hydraulic fracture problems [18-19].

In this paper, we explore the application of the extended finite element method to hydraulic fracture problems. By taking good advantage of the XFEM and the flexible functionality of user subroutines provided in ABAQUS [20], a user-defined 2-D quadrilateral plane strain element has been coded in Fortran to incorporate the extended finite element capabilities in 2-D hydraulic fracture problems. The user-defined element includes the desired aspects of the XFEM so as to model crack propagation without explicit remeshing. In addition, the extended fluid pressure degrees of freedom are assigned to the appropriate nodes of the proposed elements in order to describe the viscous flow of fluid inside the crack and its contribution to the coupled crack deformation.

2. Problem formulation

2.1. Problem definition

Consider a two-dimensional hydraulically driven fracture Γcpropagating in a homogeneous, isotropic, linear elastic, impermeable medium Ωunder plane strain conditions, see Figure 1. The boundary of the domain consists of ΓFon which prescribed tractions F, are imposed, Γuon which prescribed displacements (assumed to be zero for simplicity) are imposed, and crack faces Γcsubject to fluid pressure. The fracture propagation is driven by injection of an incompressible Newtonian fluid at constant volumetric rate Q0at a fixed injection point. It is assumed that the fracture propagation is quasi-static, and that the fracture is completely filled with the injected fluid, i.e., there is no lag between the fluid front and the fracture tip. The solution of the problem consists of determining the evolution of the fracture length, as well as the fracture opening, the fluid pressure, and the deformations and stresses inside the domain as functions of both position and time.

Figure 1.

A two-dimensional domain containing a hydraulic fracture

2.2. Governing equations

The stress inside the domain, σ, is related to the external loading Fand the fluid pressure pthrough the equilibrium equations:

σ=0, on Ωσn=F, on ΓF σn-=-σn+=-pn+=pn-, on Γc E1

where nis the unit normal vector.

Under the assumptions of small strains and displacements, the kinematic equations, which include the strain-displacement relationship, the prescribed displacement boundary conditions and the crack surfaces separation, read

ε=u+uT/2 on Ωu=0 on Γuw=u+-u- on ΓcE2

where uis the displacement, wis the separation between the two faces of the crack, and εis the strain.

The isotropic, linear elastic constitutive law is


where Dis Hooke’s tensor.

The fluid flow in the crack is modelled using lubrication theory, given by Poiseuille’s law


where μis the dynamic viscosity of the fracturing fluid, q, the flow rate inside the crack per unit extend of the crack in the direction of x, is equal to the average velocity v-times the crack opening w(see Figure 2), i.e.,


The fracturing fluid is considered to be incompressible, so the mass conservation equation for the fluid may be expressed as


where the leak-off rate gxaccounts for fluid exchange between the fracture and the surrounding medium (e.g. porous rock).

Figure 2.

Fluid flow within fracture

Substituting of Eq. (4) into Eq. (6) leads to the governing equation for the fluid flow within the fracture


where k=w312μis the conductivity. The general form of Eq. (7) reads


where kis the conductivity tensor.

According to linear elastic fracture mechanics, the criterion that the fracture propagates continuously in mobile equilibrium (quasi-static) takes the form


where KIis the mode I stress intensity factor and KIcthe material fracture toughness.

At the inlet, the fluid flux is equal to the injection rate, i.e.,


At the crack tip, the boundary conditions are given by the zero fracture opening and zero flow conditions, i.e.,


The above equations constitute the complete formulation that can be used to predict the evolution of the hydraulic fracture.

3. Weak form and FEM discretization

The weak form of the equilibrium equation is given by

ΩδεTσdΩ -ΩδuTbdΩ -ΓtδuTtdΓ -Γc+δuc+Tpc+dΓ +Γc-δuc-Tpc-dΓ =0E12

where bis the body force, t is the applied traction on the boundary Γt, δuis an arbitrary virtual displacement and δεis the corresponding virtual strain, which is related to δuthrough the strain operator Sas δε=Sδu.

For the fluid pressure on the crack surfaces, we define


The crack opening displacement wis given by

w=ncTuc+-uc-, or w=ncuc+-uc-ncE14

Then the weak form of the equilibrium equation can be expressed in a more compact form as

ΩδεTσdΩ -ΩδuTbdΩ -ΓtδuTtdΓ -ΓcδwTpdΓ =0E15

The weak form of the governing equation for the fluid flow within the fracture can be written as

ΓcδpTw˙-Tkp+gdΓ =0E16

which, after integration by parts and substitution of the boundary conditions described above, yields

ΓcδpTw˙dΓ +ΓcTδpkpdΓ +ΓcδpTgdΓ =0E17

Consider the coupled problem discretized in the standard (displacement) manner with the displacement vector uapproximated as

uu^=i=1nNiuui=Nuu~, δuNuδu~E18

and the fluid pressure psimilarly approximated by

pp^=i=1nNippi=Npp~, δpNpδp~E19

where uiand piare the nodal displacement and pressure, Niuand Nipare corresponding nodal displacement and fluid pressure shape functions.

The crack opening displacement w(or more generally displacement discontinuity) is approximated by

ww^=i=1nNiwui=Nwu~, δwNwδu~E20

where Niware the approximate crack opening displacement shape function. It will be shown later that the shaped function Niwcan be expressed in terms of the displacement shape functions Niuaccording to the relationship Eq. (14).

Substitution of the displacement and pressure approximations (Eqs. (18)-(20)) and the constitutive equation (Eq. (3)) into Eq. (15) yields a system of algebraic equations for the discrete structural problem


where the structural stiffness matrix


and the equivalent nodal force vector


and the coupling term arise due to the pressure (tractions) on the crack surface through the matrix


By substituting Eqs. (19) and (20) into Eq. (17), the standard discretization applied to the weak form of the fluid flow equation leads to a system of algebraic equations for the discrete fluid flow problem



C=QT=ΓcNpTnTNwdΓ, H=ΓcNpTkNpdΓ,fp=-ΓcNpTgdΓ E26

Then, the discrete governing equations for the coupled fluid-fracture problem can be expressed in matrix form as:


The above equations form the basis for the construction of a finite element which couples the fluid flow within the crack and crack propagation.

4. The extended finite elemet method and element implementation

4.1. Extended finite element approximation

The XFEM approximation of the displacement field for the crack problem can be expressed as [17]


where Nis the set of all nodes in the mesh, Ncrthe set of nodes whose support are bisected by the crack surface Γc, Ntipthe set of nodes whose support are partially cut by the crack surface, NIxand N~Ixare the standard finite element shape functions, uIare the displacement nodal degrees of freedom, aIand bIlare the additional degrees of freedom for the displacement, and Hxand Blr,θare the appropriate enrichment basis functions which are localized by N~Ix. The shape function N~Ixcan differ from NIx.

The displacement discontinuity given by a crack Γccan be represented by the generalized Heaviside step function

Hx=Hdx=signdx=  1       dx0-1       dx<0E29

where dxis the signed distance of the point xto Γc.

The enrichment basis functions Blr,θare required to model the displacement around the crack tip, which are generally chosen as a basis that approximately spans the two-dimensional plane strain asymptotic crack tip fields in the linear elastic fracture mechanics:


where r,θare the local polar coordinates at the crack tip.

According to Eq. (28), the displacement discontinuity between the two surfaces of the crack can be obtained as

wx=uc+x-uc-x=2INcrN~IxaI+2INtipN~IxB1r,πbI1xΓc  E31

Combination of Eqs. (31) and (20) enables determining the shape function Nw.

The fluid pressure field within the crack is approximated by

px=INcrNIpxpI    xΓcE32

where NIpxare the standard finite element shape functions. In some cases, it can also be chosen as a special function so as to allow for the pressure singularity at the crack tip and the associated near-tip asymptotic fracture opening associated with a zero-lag viscosity dominated regime in a hydraulic fracture [18].

4.2. Element implementation

As shown in Figure 3, the two-dimensional 4-node plane strain channel and tip elements have been constructed for the hydraulic fracture problem. Each node has the standard displacement degrees of freedom uI. The additional degrees of freedom aIand bIlare assigned to the four nodes of channel and tip elements, respectively. In addition, the virtual degree of freedom of fluid pressure has been assigned to nodes 3 and 4 so as to represent the internal fluid pressure within the crack. It should be pointed out that the nodes 3 and 4 physically do not have fluid pressure degrees of freedom because here the fluid flow is confined within the crack, and the integral calculation of the related element matrixes and equivalent nodal forces (e.g. Eq. (26)) must be correctly carried out along the true crack path within the element.

Figure 3.

2-D 4-node -node plane strain hydraulic fracture elements

So, the active degrees of freedom for the channel element are


and for the tip element the Heaviside enriched degrees of freedom aIneed to be replaced by the crack tip field enriched degrees of freedom bIl.

Gauss quadrature is used to calculate the system matrix and equivalent nodal force. Since the discontinuous enrichment functions are introduced in approximating the displacement field, integration of discontinuous functions is needed when computing the element stiffness matrix and equivalent nodal force. In order to ensure the integral accuracy, it is necessary to modify the quadrature routine. Both the channel and tip elements are partitioned by the crack surface into two quadrature sub-cells where the integrands are continuous and differentiable. Then Gauss integration is carried out by a loop over the sub-cells to obtain an accurate integration result.

Due to the flexibility, the user subroutine of UEL provided in the finite element package ABAQUS [20] has been employed in implementing the proposed elements in Fortran code. The main purpose of UEL is to provide the element stiffness matrix as well as the right hand side residual vector, as need in a context of solving the discrete system of equations.

5. Numerical examples

The proposed user element together with the structural elements provided in the ABAQUS element library are used to establish a finite element model to investigate a plane strain hydraulic fracture problem in an infinite impermeable elastic medium. The far-field boundary conditions are modelled by using infinite elements. The initial testing of this new element formulation involves using boundary value problems of an imposed fluid pressure and an imposed fracture opening. These problems are used to test for both of the two limiting cases of a toughness-dominated and viscosity-dominated plane-strain hydraulic fracture for which the analytical solutions are available in the literature [21].

Comparisons of the FEM predictions with the available analytical solutions to the two limiting cases are given in Figures 4 and 5.

Figure 4.

Zero-viscosity case: (a) imposed pressure; and (b) imposed opening

Figure 5.

Zero-toughness case: (a) imposed pressure; and (b) imposed opening

The simulation results for a plane strain toughness-dominated KGD hydraulic fracture are shown in Figure 4. The corresponding analytical solutions for the zero-viscosity case are also shown for comparison. The crack opening is obtained by imposing a given pressure calculated according to the analytical solution (Eq. (41)) on the crack surface of the finite element model. While the fluid pressure is obtained by applying an opening profile calculated from the analytical solution (Eq. (40)) to the crack surface of the finite element model. For the zero-toughness case (Figure 5), the crack opening and the fluid pressure are obtained by imposing the analytical solution of pressure (Eq. (45)) and crack opening (Eq. (44)) to the crack surface of the finite element model, respectively. Only twenty channel elements in total are meshed along the crack length in the finite element model.

It can be seen that the XFEM predictions generally compare well with the analytical solutions for crack openings, while for the fluid pressure the XFEM predictions differ from the analytical solutions at the region close to the crack tip. One main reason for the deviation of the predicted fluid pressure from the analytical solutions near the tip is likely to be because the user-defined element is assumed to be cut through by the crack and no tip element is included in the finite element model. Another reason could be that a static fracture rather than a propagating fracture is simulated here. Improved prediction can be expected with the implementation of a crack tip user-defined element that captures the crack tip singularity correctly.

6. Summary

The application of the extended finite element method to the hydraulic fracture problems has been presented. The discrete governing equations for the coupled fluid-fracture problem have been derived. A user element based on the XFEM has been implemented in ABAQUS, which includes the desired aspects of the XFEM so as to model crack propagation without explicit remeshing. In addition, the fluid pressure degrees of freedom have been introduced and assigned to the appropriate nodes of the proposed element to describe the fluid flow within the crack and its contribution to the crack deformation. Verification of the user-defined element has been made by comparing the FEM predictions with the analytical solutions available in the literature. The preliminary result presented here is a first attempt to the promising application of the XFEM to the hydraulic fracture simulation.

Apendix: Analytical solutions for plane strain Kristianovic-Geertsma-de Klerk (KGD) hydraulic fractures

The solution of a plane strain KGD hydraulic fracture in an infinite elastic body depends on the injection rate Q0and on the three material parameters E', K', and μ', which are defined as [21]

E'=E/1-ν2, K'=32/π1/2KIc, μ'=12μ

For the plane strain KGD hydraulic fracture, the crack opening wx,t, crack length (half length) lt, and net fluid pressure px,tcan be expressed as [21]


where ξ=x/ltis the scaled coordinate 0ξ1, ϵtis a small dimensionless parameter, Ltdenotes a length scale of the same order as the fracture length lt, Ptis the dimensionless evolution parameter, and γPtis dimensionless fracture length.

The evolution parameter Ptcan be interpreted as a dimensionless toughness κin the viscosity scaling [21]


or as a dimensionless viscosity Min the toughness scaling [21]


For the toughness scaling, denoted by a subscript k, the small parameter ϵtand the length scale Lttake the explicit forms [21]


The solution for the zero viscosity case is given by [21]


For the viscosity scaling, denoted by a subscript m, the small parameter ϵtand the length scale Lttake the explicit forms [21]

ϵmt=μ'/E't1/3  , Lmt= E'Q03t4/μ'1/6E41

The first order approximation of the zero toughness solution is [21]


where A0=31/2, A11-0.156, and B10.0663; Bis Euler beta function, and F1is hypergeometric function. Thus, Ω-m001.84.


The author would like to thank Dr. Rob Jeffrey for the support of this work. Furthermore, the author thanks CSIRO CESRE for support and for granting permission to publish.

© 2013 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Zuorong Chen (May 17th 2013). An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems, Effective and Sustainable Hydraulic Fracturing, Andrew P. Bunger, John McLennan and Rob Jeffrey, IntechOpen, DOI: 10.5772/56287. Available from:

chapter statistics

3990total chapter downloads

16Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Stress Intensity Factor Determination for Three-Dimensional Crack Using the Displacement Discontinuity Method with Applications to Hydraulic Fracture Height Growth and Non- Planar Propagation Paths

By Farrokh Sheibani and Jon Olson

Related Book

First chapter

Introduction to Infrared Spectroscopy

By Theophile Theophanides

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us