An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

Abstract

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 ${\mathrm{\Gamma }}_c$ propagating in a homogeneous, isotropic, linear elastic, impermeable medium $\mathrm{\Omega }$ under plane strain conditions, see Figure 1. The boundary of the domain consists of ${\mathrm{\Gamma }}_F$ on which prescribed tractions $\mathbf{F}$, are imposed, ${\mathrm{\Gamma }}_u$ on which prescribed displacements (assumed to be zero for simplicity) are imposed, and crack faces ${\mathrm{\Gamma }}_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.

Figure 1.

A two-dimensional domain containing a hydraulic fracture

2.2. Governing equations

The stress inside the domain, $\mathbf{\sigma }$, is related to the external loading $\mathbf{F}$ and the fluid pressure $p$ through the equilibrium equations:

$$\begin{array}{c}\nabla \bullet \mathbf{\sigma }=0,\mathrm{ }\mathrm{o}\mathrm{n}\mathrm{ }\mathrm{\Omega }\\ \mathbf{\sigma }\bullet \mathbf{n}=\mathbf{F},\mathrm{ }\mathrm{o}\mathrm{n}\mathrm{ }{\mathrm{\Gamma }}_F\mathrm{ }\\ \mathbf{\sigma }\bullet {\mathbf{n}}^{-}=-\mathbf{\sigma }\bullet {\mathbf{n}}^{+}=-p{\mathbf{n}}^{+}=p{\mathbf{n}}^{-}, \mathrm{o}\mathrm{n} {\mathrm{\Gamma }}_c \end{array}$$E1

where $\mathbf{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

$$\begin{array}{c}\mathbf{\epsilon }=\left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\mathrm{T}}\right)/2 on \mathrm{\Omega }\\ \mathbf{u}=0\mathrm{ }\mathrm{o}\mathrm{n}\mathrm{ }{\mathrm{\Gamma }}_u\\ \mathbf{w}={\mathbf{u}}^{+}-{\mathbf{u}}^{-} \mathrm{o}\mathrm{n}\mathrm{ }{\mathrm{\Gamma }}_c\end{array}$$E2

where $\mathbf{u}$ is the displacement, $\mathbf{w}$ is the separation between the two faces of the crack, and $\mathbf{\epsilon }$ is the strain.

The isotropic, linear elastic constitutive law is

$$\mathbf{\sigma }=\mathbf{D}:\mathbf{\epsilon }$$E3

where $\mathbf{D}$ is Hooke’s tensor.

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

$$q=-\frac{w^3}{12\mu }\frac{\partial p}{\partial x}$$E4

where $\mu$ 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 $\stackrel{-}{v}$ times the crack opening $w$ (see Figure 2), i.e.,

$$q\left(x\right)=\stackrel{-}{v}\left(x\right)w\left(x\right)$$E5

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

$$\frac{\partial w}{\partial t}+\frac{\partial q}{\partial x}+g=0$$E6

where the leak-off rate $g\left(x\right)$ accounts 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

$$\frac{\partial w}{\partial t}-\frac{\partial }{\partial x}\left(k\frac{\partial p}{\partial x}\right)+g=0$$E7

where $k=\frac{w^3}{12\mu }$ is the conductivity. The general form of Eq. (7) reads

$$\dot{w}-{\nabla }^{\mathrm{T}}\left(\mathbf{k}\nabla p\right)+g=0$$E8

where $\mathbf{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

$$K_{\mathrm{I}}=K_{\mathrm{I}\mathrm{c}}$$E9

where $K_{\mathrm{I}}$ is the mode I stress intensity factor and $K_{\mathrm{I}\mathrm{c}}$ the material fracture toughness.

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

$${\left.q\right|}_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}=Q_0$$E10

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

$${\left.w\right|}_{\mathrm{t}\mathrm{i}\mathrm{p}}={\left.q\right|}_{\mathrm{t}\mathrm{i}\mathrm{p}}=0$$E11

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

$${\int }_{\mathrm{\Omega }}^{}\delta {\mathbf{\epsilon }}^{\mathrm{T}}\mathbf{\sigma }\mathrm{d}\mathrm{\Omega }\mathbf{ }-{\int }_{\mathrm{\Omega }}^{}\delta {\mathbf{u}}^{\mathrm{T}}\mathbf{b}\mathrm{d}\mathrm{\Omega }\mathbf{ }-{\int }_{{\mathrm{\Gamma }}_{\mathrm{t}}}^{}\delta {\mathbf{u}}^{\mathrm{T}}\mathbf{t}\mathrm{d}\mathrm{\Gamma }\mathbf{ }-\left({\int }_{{\mathrm{\Gamma }}_c^{+}}^{}\delta {{\mathbf{u}}_c^{+}}^{\mathrm{T}}{\mathbf{p}}_c^{+}\mathrm{d}\mathrm{\Gamma }\mathbf{ }+{\int }_{{\mathrm{\Gamma }}_c^{-}}^{}\delta {{\mathbf{u}}_c^{-}}^{\mathrm{T}}{\mathbf{p}}_c^{-}\mathrm{d}\mathrm{\Gamma }\mathbf{ }\right)=0$$E12

where $\mathbf{b}$ is the body force, t is the applied traction on the boundary ${\mathrm{\Gamma }}_{\mathrm{t}}$, $\delta \mathbf{u}$ is an arbitrary virtual displacement and $\delta \mathbf{\epsilon }$ is the corresponding virtual strain, which is related to $\delta \mathbf{u}$ through the strain operator $\mathbf{S}$ as $\delta \mathbf{\epsilon }=\mathbf{S}\delta \mathbf{u}$.

For the fluid pressure on the crack surfaces, we define

$$\mathbf{p}={\mathbf{p}}_c^{+}=-{\mathbf{p}}_c^{-}=p{\mathbf{n}}_c=p{\mathbf{n}}_c^{-}=-p{\mathbf{n}}_c^{+}$$E13

The crack opening displacement $\mathbf{w}$ is given by

$$w={\mathbf{n}}_c^{\mathrm{T}}\bullet \left({\mathbf{u}}_c^{+}-{\mathbf{u}}_c^{-}\right), \mathrm{o}\mathrm{r} \mathbf{w}={\mathbf{n}}_c\bullet \left({\mathbf{u}}_c^{+}-{\mathbf{u}}_c^{-}\right)\bullet {\mathbf{n}}_c$$E14

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

$${\int }_{\mathrm{\Omega }}^{}\delta {\mathbf{\epsilon }}^{\mathrm{T}}\mathbf{\sigma }\mathrm{d}\mathrm{\Omega }\mathbf{ }-{\int }_{\mathrm{\Omega }}^{}\delta {\mathbf{u}}^{\mathrm{T}}\mathbf{b}\mathrm{d}\mathrm{\Omega }\mathbf{ }-{\int }_{{\mathrm{\Gamma }}_{\mathrm{t}}}^{}\delta {\mathbf{u}}^{\mathrm{T}}\mathbf{t}\mathrm{d}\mathrm{\Gamma }\mathbf{ }-{\int }_{{\mathrm{\Gamma }}_{\mathrm{c}}}^{}\delta {\mathbf{w}}^{\mathrm{T}}\mathbf{p}\mathrm{d}\mathrm{\Gamma }\mathbf{ }=0$$E15

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

$${\int }_{{\mathrm{\Gamma }}_{\mathrm{c}}}^{}\delta p^{\mathrm{T}}\left(\dot{w}-{\nabla }^{\mathrm{T}}\left(\mathbf{k}\nabla p\right)+g\right)\mathrm{d}\mathrm{\Gamma }\mathbf{ }=0$$E16

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

$${\int }_{{\mathrm{\Gamma }}_{\mathrm{c}}}^{}\delta p^{\mathrm{T}}\dot{w}\mathrm{d}\mathrm{\Gamma }\mathbf{ }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{c}}}^{}{\nabla }^{\mathrm{T}}\left(\delta p\right)\mathbf{k}\nabla p\mathrm{d}\mathrm{\Gamma }\mathbf{ }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{c}}}^{}\delta p^{\mathrm{T}}g\mathrm{d}\mathrm{\Gamma }\mathbf{ }=0$$E17

Consider the coupled problem discretized in the standard (displacement) manner with the displacement vector $\mathbf{u}$ approximated as

$$\mathbf{u}\approx \widehat{\mathbf{u}}=\sum _{i=1}^n{\mathbf{N}}_i^u{\mathbf{u}}_i={\mathbf{N}}^{\mathbf{u}}\tilde{\mathbf{u}}, \mathit{\delta }\mathbf{u}\approx {\mathbf{N}}^{\mathbf{u}}\delta \tilde{\mathbf{u}}$$E18

and the fluid pressure $p$ similarly approximated by

$$p\approx \widehat{p}=\sum _{i=1}^n{N}_i^p{p}_i={\mathbf{N}}^{\mathrm{p}}\tilde{\mathbf{p}},\mathbf{ }\delta p\approx {\mathbf{N}}^{\mathrm{p}}\delta \tilde{\mathbf{p}}$$E19

where ${\mathbf{u}}_i$ and $p_i$ are the nodal displacement and pressure, ${\mathbf{N}}_i^u$ and $N_i^p$ are corresponding nodal displacement and fluid pressure shape functions.

The crack opening displacement $\mathbf{w}$ (or more generally displacement discontinuity) is approximated by

$$\mathbf{w}\approx \widehat{\mathbf{w}}=\sum _{i=1}^n{\mathbf{N}}_i^w{\mathbf{u}}_i={\mathbf{N}}^w\tilde{\mathbf{u}}, \delta \mathbf{w}\approx {\mathbf{N}}^w\delta \tilde{\mathbf{u}}$$E20

where ${\mathbf{N}}_i^w$ are the approximate crack opening displacement shape function. It will be shown later that the shaped function ${\mathbf{N}}_i^w$ can be expressed in terms of the displacement shape functions ${\mathbf{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

$$\mathbf{K}\tilde{\mathbf{u}}-\mathbf{Q}\tilde{\mathbf{p}}-{\mathbf{f}}^u=0$$E21

where the structural stiffness matrix

$$\mathbf{K}={\int }_{\mathrm{\Omega }}^{}{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\mathbf{B}d\mathrm{\Omega }$$E22

and the equivalent nodal force vector

$${\mathbf{f}}^u={\int }_{\mathrm{\Omega }}^{}{\left({\mathbf{N}}^u\right)}^{\mathrm{T}}\mathbf{b}d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{t}}}^{}{\left({\mathbf{N}}^u\right)}^{\mathrm{T}}\mathbf{t}d\mathrm{\Gamma }$$E23

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

$$\mathbf{Q}={\int }_{{\mathrm{\Gamma }}_c}^{}{\left({\mathbf{N}}^w\right)}^{\mathrm{T}}\mathbf{n}{\mathbf{N}}^{p}d\mathrm{\Gamma }$$E24

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

$$\mathbf{C}\dot{\tilde{\mathbf{u}}}+\mathbf{H}\tilde{\mathbf{p}}-{\mathbf{f}}^p=0$$E25

where

$$\mathbf{C}={\mathbf{Q}}^{\mathrm{T}}={\int }_{{\mathrm{\Gamma }}_c}^{}{\left({\mathbf{N}}^p\right)}^{\mathrm{T}}{\mathbf{n}}^{\mathrm{T}}{\mathbf{N}}^{\mathrm{w}}d\mathrm{\Gamma }, \mathbf{H}={\int }_{{\mathrm{\Gamma }}_c}^{}{\left({\nabla \mathbf{N}}^p\right)}^{\mathrm{T}}\mathbf{k}\nabla {\mathbf{N}}^{p}d\mathrm{\Gamma },{\mathbf{f}}^p=-{\int }_{{\mathrm{\Gamma }}_c}^{}{\left({\mathbf{N}}^p\right)}^{\mathrm{T}}gd\mathrm{\Gamma }$$E26

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

$$\left[\begin{array}{cc}0& 0\\ \mathbf{C}& 0\end{array}\right]\left(\begin{array}{c}\dot{\tilde{\mathbf{u}}}\\ \dot{\tilde{\mathbf{p}}}\end{array}\right)+\left[\begin{array}{cc}\mathbf{K}& -\mathbf{Q}\\ 0& \mathbf{H}\end{array}\right]\left(\begin{array}{c}\tilde{\mathbf{u}}\\ \tilde{\mathbf{p}}\end{array}\right)=\left(\begin{array}{c}{\mathbf{f}}^u\\ {\mathbf{f}}^p\end{array}\right)$$E27

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]

$$\mathbf{u}\left(\mathbf{x}\right)=\sum _{I\in N}{\mathbf{N}}_I\left(\mathbf{x}\right){\mathbf{u}}_I+\sum _{I\in N_{cr}}{\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right){\left(H\left(\mathbf{x}\right)-H\left({\mathbf{x}}_I\right)\right)\mathbf{a}}_I+\sum _{I\in N_{tip}}{\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right)\sum _{l=1}^4\left(B^{\left(l\right)}\left(r,\theta \right)-B^{\left(l\right)}\left(r_I,{\theta }_I\right)\right){\mathbf{b}}_I^{\left(l\right)}$$E28

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 ${\mathrm{\Gamma }}_c$, $N_{tip}$ the set of nodes whose support are partially cut by the crack surface, ${\mathbf{N}}_I\left(\mathbf{x}\right)$ and ${\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right)$ are the standard finite element shape functions, ${\mathbf{u}}_I$ are the displacement nodal degrees of freedom, ${\mathbf{a}}_I$ and ${\mathbf{b}}_I^{\left(l\right)}$ are the additional degrees of freedom for the displacement, and $H\left(\mathbf{x}\right)$ and $B^{\left(l\right)}\left(r,\theta \right)$ are the appropriate enrichment basis functions which are localized by ${\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right)$. The shape function ${\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right)$ can differ from ${\mathbf{N}}_I\left(\mathbf{x}\right)$.

The displacement discontinuity given by a crack ${\mathrm{\Gamma }}_c$ can be represented by the generalized Heaviside step function

$$H\left(\mathbf{x}\right)=H\left(d\left(\mathbf{x}\right)\right)=sign\left(d\left(\mathbf{x}\right)\right)=\left\{\begin{array}{c} 1 d\left(\mathbf{x}\right)\ge 0\\ -1 d\left(\mathbf{x}\right)<0\end{array}\right.$$E29

where $d\left(\mathbf{x}\right)$ is the signed distance of the point $\mathbf{x}$ to ${\mathrm{\Gamma }}_c$.

The enrichment basis functions $B^{\left(l\right)}\left(r,\theta \right)$ 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:

$${\left\{B^{\left(l\right)}\right\}}_{l=1}^4=\sqrt{r}\left\{\begin{array}{cccc}\sin \left(\theta /2\right)& \cos \left(\theta /2\right)& \sin \left(\theta /2\right)\sin \left(\theta \right)& \cos \left(\theta /2\right)\sin \left(\theta \right)\end{array}\right\}$$E30

where $\left(r,\theta \right)$ 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

$$\mathbf{w}\left(\mathbf{x}\right)={\mathbf{u}}_c^{+}\left(\mathbf{x}\right)-{\mathbf{u}}_c^{-}\left(\mathbf{x}\right)=2\sum _{I\in N_{cr}}{\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right){\mathbf{a}}_I+2\sum _{I\in N_{tip}}{\tilde{\mathbf{N}}}_I\left(\mathbf{x}\right)B^{\left(1\right)}\left(r,\pi \right){\mathbf{b}}_I^{\left(1\right)}\mathbf{x}\in {\mathrm{\Gamma }}_c\mathrm{ }$$E31

Combination of Eqs. (31) and (20) enables determining the shape function ${\mathbf{N}}^w$.

The fluid pressure field within the crack is approximated by

$$p\left(\mathbf{x}\right)=\sum _{I\in N_{cr}}{N}_I^p\left(\mathbf{x}\right)p_I\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{x}\in {\mathrm{\Gamma }}_c$$E32

where $N_I^p\left(\mathbf{x}\right)$ 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].

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 ${\mathbf{u}}_I$. The additional degrees of freedom ${\mathbf{a}}_I$ and ${\mathbf{b}}_I^{\left(l\right)}$ 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.

Figure 3.

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

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

$${\widehat{u}}^e={\left\{\underset{\mathrm{Standard}}{\underbrace{\begin{array}{cccccccc}{u}_1^x& u_1^y& u_2^x& u_2^y& u_3^x& u_3^y& u_4^x& u_4^y\end{array}}}\underset{\mathrm{Heaviside}Enriched}{\underbrace{\begin{array}{cccccccc}{a}_1^x& a_1^y& a_2^x& a_2^y& a_3^x& a_3^y& a_4^x& a_4^y\end{array}}}\underset{\mathrm{Coupled}}{\underbrace{\begin{array}{cc}{p}_3& p_4\end{array}}}\right\}}_{18}^T$$E33

and for the tip element the Heaviside enriched degrees of freedom ${\mathbf{a}}_I$ need to be replaced by the crack tip field enriched degrees of freedom ${\mathbf{b}}_I^{\left(l\right)}$.

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 $Q_0$ and on the three material parameters $E^{\text{'}}$, $K^{\text{'}}$, and ${\mu }^{\text{'}}$, which are defined as [21]

$E^{\mathrm{\text{'}}}=E/1-{\nu }^2$, $K^{\mathrm{\text{'}}}={\left(32/\pi \right)}^{1/2}{K}_{\mathrm{I}c}$, ${\mu }^{\mathrm{\text{'}}}=12\mu$

For the plane strain KGD hydraulic fracture, the crack opening $w\left(x,t\right)$, crack length (half length) $l\left(t\right)$, and net fluid pressure $p\left(x,t\right)$ can be expressed as [21]

$$\begin{array}{c}w\left(x,t\right)=ϵ\left(t\right)L\left(t\right)\mathrm{\Omega }\left[\xi ,P\left(t\right)\right]=ϵ\left(t\right)L\left(t\right)\gamma \left[P\left(t\right)\right]\stackrel{-}{\mathrm{\Omega }}\left(\xi \right)\\ p\left(x,t\right)=ϵ\left(t\right)E^{\mathrm{\text{'}}}\mathrm{\Pi }\left[\xi ,P\left(t\right)\right]\\ l\left(t\right)=\gamma \left[P\left(t\right)\right]L\left(t\right)\end{array}$$E34

where $\xi =x/l\left(t\right)$ is the scaled coordinate $\left(0\le \xi \le 1\right)$, $ϵ\left(t\right)$ is a small dimensionless parameter, $L\left(t\right)$ denotes a length scale of the same order as the fracture length $l\left(t\right)$, $P\left(t\right)$ is the dimensionless evolution parameter, and $\gamma \left[P\left(t\right)\right]$ is dimensionless fracture length.

The evolution parameter $P\left(t\right)$ can be interpreted as a dimensionless toughness $\kappa$ in the viscosity scaling [21]

$$\kappa =K^{\text{'}}{\left(E^{\text{'}3}{\mu }^{\text{'}}{Q}_0\right)}^{-1/4}$$E35

or as a dimensionless viscosity $\mathcal{M}$ in the toughness scaling [21]

$$\mathcal{M}={\mu }^{\text{'}}{E}^{\text{'}3}{Q}_0/K^{\text{'}}$$E36

For the toughness scaling, denoted by a subscript $k$, the small parameter $ϵ\left(t\right)$ and the length scale $L\left(t\right)$ take the explicit forms [21]

$${ϵ}_k\left(t\right)={\left(K^{\mathrm{\text{'}}4}/E^{\mathrm{\text{'}}4}{Q}_{0}t\right)}^{1/3},L_k\left(t\right)={\left(E^{\mathrm{\text{'}}}{Q}_{0}t/K^{\mathrm{\text{'}}}\right)}^{2/3}$$E37

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

$${\gamma }_{k0}=2/{\pi }^{2/3}$$E38

$${\stackrel{-}{\mathrm{\Omega }}}_{k0}\left(\xi \right)=\sqrt{1-{\mathrm{\xi }}^2}/{\pi }^{1/3}$$E39

$${\mathrm{\Pi }}_{k0}={\pi }^{1/3}/8$$E40

For the viscosity scaling, denoted by a subscript $m$, the small parameter $ϵ\left(t\right)$ and the length scale $L\left(t\right)$ take the explicit forms [21]

$${ϵ}_m\left(t\right)={\left({\mu }^{\mathrm{\text{'}}}/E^{\mathrm{\text{'}}}t\right)}^{1/3}\mathrm{ }\mathrm{ },\mathrm{ }{L}_m\left(t\right)=\mathrm{ }{\left(E^{\mathrm{\text{'}}}{Q}_0^3{t}^4/{\mu }^{\mathrm{\text{'}}}\right)}^{1/6}$$E41

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

$${\gamma }_{m0}\cong 0.616$$E42

$${\stackrel{-}{\mathrm{\Omega }}}_{m0}\left(\xi \right)=A_0{\left(1-{\xi }^2\right)}^{2/3}+A_1^{\left(1\right)}{\left(1-{\xi }^2\right)}^{5/3}+B^{\left(1\right)}\left[4\sqrt{1-{\xi }^2}+2{\xi }^2\mathrm{l}\mathrm{n}\left|\frac{1-\sqrt{1-{\xi }^2}}{1+\sqrt{1-{\xi }^2}}\right|\right]$$E43

$${\mathrm{\Pi }}_{m0}^{\left(1\right)}=\frac{1}{3\pi }B\left(\frac{1}{2},\frac{2}{3}\right)\left[A_0{F}_1\left(-\frac{1}{6},1;\frac{1}{2};{\xi }^2\right)+\frac{10}{7}{A}_1^{\left(1\right)}{F}_1\left(-\frac{7}{6},1;\frac{1}{2};{\xi }^2\right)\right]+B^{\left(1\right)}\left(2-\pi \left|\xi \right|\right)$$E44

where $A_0=3^{1/2}$, $A_1^{\left(1\right)}\cong -0.156$, and $B^{\left(1\right)}\cong 0.0663$; $B$ is Euler beta function, and $F_1$ is hypergeometric function. Thus, ${\stackrel{-}{\mathrm{\Omega }}}_{m0}\left(0\right)\cong 1.84$.

Acknowledgments

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.