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

### 2.2. Governing equations

The stress inside the domain,

where

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

where

The isotropic, linear elastic constitutive law is

where

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

where

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

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

where

where

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

where

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

where

For the fluid pressure on the crack surfaces, we define

The crack opening displacement

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

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

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

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

and the fluid pressure

where

The crack opening displacement

where

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

where

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

The displacement discontinuity given by a crack

where

The enrichment basis functions

where

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

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

The fluid pressure field within the crack is approximated by

where

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

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

(33) |

and for the tip element the Heaviside enriched degrees of freedom

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.

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

For the plane strain KGD hydraulic fracture, the crack opening

where

The evolution parameter

or as a dimensionless viscosity

For the toughness scaling, denoted by a subscript

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

For the viscosity scaling, denoted by a subscript

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

where