An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

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 ele‐ ment 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 de‐ formation. 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 compar‐ ing the finite element results with the analytical solutions available in the literature.


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][4][5][6][7][8][9][10][11][12][13][14][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.

Problem definition
Consider a two-dimensional hydraulically driven fracture Γ c propagating in a homogeneous, isotropic, linear elastic, impermeable medium Ω under plane strain conditions, see Figure 1. The boundary of the domain consists of Γ F on which prescribed tractions F, are imposed, Γ u on which prescribed displacements (assumed to be zero for simplicity) are imposed, and crack faces Γ c subject to fluid pressure. The fracture propagation is driven by injection of an incompressible Newtonian fluid at constant volumetric rate Q 0 at 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.

Governing equations
The stress inside the domain, σ, is related to the external loading F and the fluid pressure p through the equilibrium equations: where n is 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 where u is the displacement, w is the separation between the two faces of the crack, and ε is the strain.
The isotropic, linear elastic constitutive law is where D is 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 g(x) accounts for fluid exchange between the fracture and the surrounding medium (e.g. porous rock).
12μ is the conductivity. The general form of Eq. (7) reads where k is 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 K I is the mode I stress intensity factor and K Ic the 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., w| tip = q| tip = 0 (11) The above equations constitute the complete formulation that can be used to predict the evolution of the hydraulic fracture.

Weak form and FEM discretization
The weak form of the equilibrium equation is given by where b is the body force, t is the applied traction on the boundary Γ t , δu is an arbitrary virtual displacement and δε is the corresponding virtual strain, which is related to δu through the strain operator S as δε = Sδu.
For the fluid pressure on the crack surfaces, we define The crack opening displacement w is given by 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 u approximated as and the fluid pressure p similarly approximated by where u i and p i are the nodal displacement and pressure, N i u and N i p are corresponding nodal displacement and fluid pressure shape functions.
The crack opening displacement w (or more generally displacement discontinuity) is approximated by where N i w are the approximate crack opening displacement shape function. It will be shown later that the shaped function N i w can be expressed in terms of the displacement shape functions N i u according 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 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.

Extended finite element approximation
The XFEM approximation of the displacement field for the crack problem can be expressed as [17] u(x) = ∑ where N is the set of all nodes in the mesh, N cr the set of nodes whose support are bisected by the crack surface Γ c , N tip the set of nodes whose support are partially cut by the crack surface, The displacement discontinuity given by a crack Γ c can be represented by the generalized

Heaviside step function
where d (x) is the signed distance of the point x to Γ c .
The enrichment basis functions B (l ) (r, θ) 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: { B (l ) } l =1 4 = r {sin (θ / 2) cos (θ / 2) sin (θ / 2)sin (θ) cos (θ / 2)sin (θ)} 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 Combination of Eqs. (31) and (20) enables determining the shape function N w .
The fluid pressure field within the crack is approximated by where N I p (x) are 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].

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 u I . The additional degrees of freedom a I and b I (l ) are 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.  u u u u u a a a a a a a a p  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.

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

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.