**2. Problem formulation**

#### **2.1. Problem definition**

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

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 particu‐ larly 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

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 independ‐ ent 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

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

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

propagates, and thus is computationally expensive.

726 Effective and Sustainable Hydraulic Fracturing

problems [18-19].

the coupled crack deformation.

**2. Problem formulation**

**2.1. Problem definition**

The stress inside the domain, **σ**, is related to the external loading **F** and the fluid pressure *p* through the equilibrium equations:

$$\begin{aligned} \nabla \cdot \mathbf{o} &= 0, \text{ on } \Omega\\ \mathbf{o} \cdot \mathbf{n} &= \mathbf{F}\_{\prime} \text{ on } \Gamma\_{F} \\ \mathbf{o} \cdot \mathbf{n}^{\prime} &= -\mathbf{o} \cdot \mathbf{n}^{+} = -\mathbf{p} \mathbf{n}^{+} = \mathbf{p} \mathbf{n}^{+}, \text{ on } \Gamma\_{c} \end{aligned} \tag{1}$$

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 condi‐ tions and the crack surfaces separation, read

$$\begin{aligned} \mathbf{e} &= \left(\nabla \, \mathbf{u} + (\nabla \, \mathbf{u})^{\mathrm{T}}\right) / 2 \text{ on } \Omega\\ \mathbf{u} &= 0 \text{ on } \Gamma\_u\\ \mathbf{w} &= \boldsymbol{u}^+ \text{ - } \boldsymbol{u}^+ \text{ on } \Gamma\_c \end{aligned} \tag{2}$$

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

$$
\sigma = \mathbf{D} \colon \varepsilon \tag{3}
$$

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

According to linear elastic fracture mechanics, the criterion that the fracture propagates

is the mode I stress intensity factor and *K*Ic the material fracture toughness.

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

The above equations constitute the complete formulation that can be used to predict the

<sup>∂</sup> *<sup>x</sup>* ) + *g* =0 (7)

http://dx.doi.org/10.5772/56287

729

*w*˙ - ∇<sup>T</sup> (**k**∇ *p*) + *g* =0 (8)

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

*K*I= *K*Ic (9)

*q*|inlet=*Q*<sup>0</sup> (10)

*w*|tip =*q*|tip =0 (11)

dΓ ) =0 (12)

, *δ***u** is an arbitrary virtual

<sup>∂</sup> *<sup>x</sup>* (*<sup>k</sup>* <sup>∂</sup> *<sup>p</sup>*

<sup>12</sup>*μ* is the conductivity. The general form of Eq. (7) reads

∂ *w* <sup>∂</sup> *<sup>t</sup>* - <sup>∂</sup>

continuously in mobile equilibrium (quasi-static) takes the form

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

fracture

where *<sup>k</sup>* <sup>=</sup> *<sup>w</sup>* <sup>3</sup>

where *K*<sup>I</sup>

conditions, i.e.,

*∫*

where **k** is the conductivity tensor.

evolution of the hydraulic fracture.

<sup>Ω</sup>*δε* <sup>T</sup>**σ**d<sup>Ω</sup> - *<sup>∫</sup>*

strain operator **S** as *δ***ε**=**S***δ***u**.

**3. Weak form and FEM discretization**

The weak form of the equilibrium equation is given by

For the fluid pressure on the crack surfaces, we define

<sup>Ω</sup>*δ<sup>u</sup>* <sup>T</sup>**b**d<sup>Ω</sup> - *<sup>∫</sup>*

Γt

where **b** is the body force, t is the applied traction on the boundary Γ<sup>t</sup>

*δu* <sup>T</sup>**t**dΓ - (*∫*

Γ*c* <sup>+</sup>*δ***u***<sup>c</sup>* +T**p***<sup>c</sup>* + dΓ + *∫* Γ*c* -*δ***u***<sup>c</sup>* - <sup>T</sup>**p***<sup>c</sup>* -

displacement and *δ***ε** is the corresponding virtual strain, which is related to *δ***u** through the

where **D** is Hooke's tensor.

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

$$
\Delta q = -\frac{w^3}{12\mu} \frac{\partial \langle p \rangle}{\partial \langle x \rangle} \tag{4}
$$

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

$$
\sigma\_q(\mathbf{x}) = \bar{v}(\mathbf{x})\varpi(\mathbf{x})\tag{5}
$$

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 \tag{6}$$

where the leak-off rate *g*(*x*) accounts for fluid exchange between the fracture and the sur‐ rounding 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) + \mathcal{g} = 0 \tag{7}$$

where *<sup>k</sup>* <sup>=</sup> *<sup>w</sup>* <sup>3</sup> <sup>12</sup>*μ* is the conductivity. The general form of Eq. (7) reads

$$
\dot{w} \cdot \nabla^{\mathrm{T}} (\mathbf{k} \nabla \, p) + \mathcal{g} = 0 \tag{8}
$$

where **k** is the conductivity tensor.

where **u** is the displacement, **w** is the separation between the two faces of the crack, and **ε** is

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

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

where the leak-off rate *g*(*x*) accounts for fluid exchange between the fracture and the sur‐

*<sup>q</sup>* <sup>=</sup> - *<sup>w</sup>* <sup>3</sup> 12*μ* ∂ *p*

∂ *w* <sup>∂</sup> *<sup>t</sup>* + ∂ *q*

**σ**=**D**: **ε** (3)

<sup>∂</sup> *<sup>x</sup>* (4)

*q*(*x*)=*v*¯(*x*)*w*(*x*) (5)

<sup>∂</sup> *<sup>x</sup>* + *g* =0 (6)

the strain.

The isotropic, linear elastic constitutive law is

where **D** is Hooke's tensor.

728 Effective and Sustainable Hydraulic Fracturing

opening *w* (see Figure 2), i.e.,

the fluid may be expressed as

**Figure 2.** Fluid flow within fracture

rounding medium (e.g. porous rock).

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

$$\mathbf{K}\_{\rm I} = \mathbf{K}\_{\rm Ic} \tag{9}$$

where *K*<sup>I</sup> 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.,

$$\left.q\right|\_{\text{inlet}} = Q\_0 \tag{10}$$

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

$$\left.w\right\|\_{\text{tip}} = q \left| \begin{array}{c} \\ \end{array} \right. $$

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

$$\left\{ \mathbf{J}\_{\varepsilon}^{\rm T} \boldsymbol{\delta} \boldsymbol{\varepsilon}\_{\varepsilon}^{\rm T} \mathbf{q} \mathbf{d} \boldsymbol{\Omega} \mathbf{D} \cdot \mathbf{J}\_{\varepsilon}^{\rm T} \mathbf{p} \mathbf{d} \boldsymbol{\Omega} \mathbf{D} \right. \\ \left. \left. \mathbf{J}\_{\varepsilon}^{\rm T} \boldsymbol{\delta} \boldsymbol{u}^{\rm T} \mathbf{t} \mathbf{d} \boldsymbol{\Gamma} \right. \left. \left. \left. \left. \mathbf{J}\_{\varepsilon}^{\rm T} \boldsymbol{\delta} \mathbf{u}^{\rm T} \right\mathbf{p} \right. \right| \boldsymbol{\delta} \boldsymbol{\Gamma} \right. \left. \left. \left. \left. \mathbf{J}\_{\varepsilon}^{\rm T} \boldsymbol{\delta} \mathbf{u}^{\rm T} \right\mathbf{p} \right\mathbf{d} \boldsymbol{\Gamma} \right. \right| \boldsymbol{\delta} \boldsymbol{\Gamma} \right. \right\} = 0 \tag{12}$$

where **b** is the body force, t is the applied traction on the boundary Γ<sup>t</sup> , *δ***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

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

where **N***<sup>i</sup>*

**N***i*

matrix

flow problem

**C**=*Q*<sup>T</sup> =*∫*

in matrix form as:

Γ*c* (*N <sup>p</sup>*)<sup>T</sup>

where

later that the shaped function **N***<sup>i</sup>*

discrete structural problem

where the structural stiffness matrix

and the equivalent nodal force vector

*<sup>u</sup>* according to the relationship Eq. (14).

*<sup>w</sup>* are the approximate crack opening displacement shape function. It will be shown

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

**Ku**˜ - **Qp**˜ - *f*

**K**=*∫*

Ω(*Nu*)<sup>T</sup>

**Q**=*∫* Γ*c* (*Nw*)<sup>T</sup>

**Cu**˜

Γ*c*

(∇*Np*)<sup>T</sup>

**K** -**Q** <sup>0</sup> **<sup>H</sup>** (

The above equations form the basis for the construction of a finite element which couples the

Then, the discrete governing equations for the coupled fluid-fracture problem can be expressed

**k**∇*N <sup>p</sup>*

**u**˜ **p**˜ ) =(*f u f p*

*n* <sup>T</sup>*N* <sup>w</sup>*d*Γ, **H**=*∫*

0 0 **<sup>C</sup>** <sup>0</sup> (

fluid flow within the crack and crack propagation.

**u**˜ ˙ **p**˜ ˙) +

**b***d*Ω + *∫* Γt (*Nu*)<sup>T</sup>

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

**n***N <sup>p</sup>*

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

*f <sup>u</sup>* =*∫*

*<sup>w</sup>* can be expressed in terms of the displacement shape functions

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

http://dx.doi.org/10.5772/56287

731

*<sup>u</sup>* =0 (21)

**t***d*Γ (23)

*d*Γ (24)

˙ + **Hp**˜ - *f <sup>p</sup>* =0 (25)

Γ*c* (*N <sup>p</sup>*)<sup>T</sup>

) (27)

*gd*Γ (26)

*d*Γ, *f <sup>p</sup>* = - *∫*

<sup>Ω</sup>*<sup>B</sup>* <sup>T</sup>**DB***d*<sup>Ω</sup> (22)

The crack opening displacement **w** is given by

$$\mathbf{u} = \mathbf{n}\_c^T \cdot \left(\mathbf{u}\_c^+ \cdot \mathbf{u}\_c^\cdot\right)\_\prime \text{ or } \mathbf{w} = \mathbf{n}\_c \cdot \left(\mathbf{u}\_c^+ \cdot \mathbf{u}\_c^\cdot\right) \cdot \mathbf{n}\_c \tag{14}$$

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

$$\text{d}\_{\text{L}} \text{\u} \xi \,\text{\u} \text{d} \text{\u} \text{\u} \text{\u} \text{\u} \text{\u} \text{\u}^{\text{T}} \text{\u} \text{d} \text{\u} \text{\u} \text{\u} \text{\u} \text{\u}^{\text{T}} \text{\u} \text{\u} \text{\u}^{\text{T}} \text{\u} \text{\u} \text{\u}^{\text{T}} \text{\u} \text{\u} \text{\u}^{\text{T}} \text{\u} \text{\u} \text{\u} \text{\u} \text{\u} = \text{0} \tag{15}$$

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

$$\{\!\!\!\!\!\!\/ \_{\dot{\psi}}\delta p^{\rm T}(\dot{\psi}\text{-}\nabla^{\rm T}(\mathbf{k}\nabla\,p) + \mathbf{g})\mathbf{d}\Gamma = 0\tag{16}$$

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

$$\mathbf{J}\_{\Gamma\_{\varepsilon}} \delta p^{\mathrm{T}} \dot{w} \mathbf{d} \Gamma + \mathbf{J}\_{\Gamma\_{\varepsilon}} \nabla^{\mathrm{T}} (\delta p) \mathbf{k} \nabla\_{p} p \mathbf{d} \Gamma + \mathbf{J}\_{\Gamma\_{\varepsilon}} \delta p^{\mathrm{T}} \mathbf{g} \, \mathrm{d} \Gamma = \mathbf{0} \tag{17}$$

Consider the coupled problem discretized in the standard (displacement) manner with the displacement vector **u** approximated as

$$\mathbf{u} \approx \hat{\mathbf{u}} = \sum\_{i=1}^{n} \mathbf{N}\_{i}^{\mu} \mathbf{u}\_{i} = \mathbf{N}^{\mathbf{u}} \tilde{\mathbf{u}}\_{i} \quad \delta \mathbf{u} \approx \mathbf{N}^{\mathbf{u}} \delta \tilde{\mathbf{u}} \tag{18}$$

and the fluid pressure *p* similarly approximated by

$$p \approx \hat{p} = \sum\_{i=1}^{n} \mathbf{N}\_i \mathbf{^p} p\_i = \mathbf{N}^\mathbf{P} \mathbf{\tilde{p}}\_{\prime} \quad \delta p \approx \mathbf{N}^\mathbf{P} \delta \mathbf{\tilde{p}} \tag{19}$$

where **u***<sup>i</sup>* and *pi* are the nodal displacement and pressure, **N***<sup>i</sup> <sup>u</sup>* and *Ni <sup>p</sup>* are corresponding nodal displacement and fluid pressure shape functions.

The crack opening displacement **w** (or more generally displacement discontinuity) is approxi‐ mated by

$$\mathbf{w} \approx \stackrel{\frown}{\mathbf{w}} = \sum\_{i=1}^{n} \mathbf{N}\_{i}^{w} \mathbf{u}\_{i} = \mathbf{N}^{w} \tilde{\mathbf{u}}\_{\prime} \quad \delta \mathbf{w} \approx \mathbf{N}^{w} \delta \tilde{\mathbf{u}} \tag{20}$$

where **N***<sup>i</sup> <sup>w</sup>* are the approximate crack opening displacement shape function. It will be shown later that the shaped function **N***<sup>i</sup> <sup>w</sup>* can be expressed in terms of the displacement shape functions **N***i <sup>u</sup>* 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}} \cdot f^\prime = 0 \tag{21}$$

where the structural stiffness matrix

$$\mathbf{K} = \mathbb{f}\_{\Omega} \boldsymbol{\mathcal{B}}^{\mathrm{T}} \mathbf{D} \mathbf{B} d \boldsymbol{\Omega} \tag{22}$$

and the equivalent nodal force vector

$$f^{\boldsymbol{u}} = \mathbb{f}\_{\Omega} \Big(\mathcal{N}^{\boldsymbol{u}}\big)^{\mathrm{T}} \mathsf{b}d\,d\Omega + \mathbb{f}\_{\Gamma\_{\boldsymbol{u}}} \Big(\mathcal{N}^{\boldsymbol{u}}\big)^{\mathrm{T}} \mathsf{t}d\,\Gamma \tag{23}$$

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

$$\mathbf{Q} = \mathbf{f}\_{\Gamma\_c} \left( \mathbf{N}^w \right)^{\rm T} \mathbf{n} \mathbf{N}^p d\Gamma \tag{24}$$

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}} \cdot f^p = 0 \tag{25}$$

where

**p**=**p***<sup>c</sup>*

<sup>T</sup> <sup>⋅</sup> (**u***<sup>c</sup>* <sup>+</sup> - **<sup>u</sup>***<sup>c</sup>* -

The crack opening displacement **w** is given by

*∫*

730 Effective and Sustainable Hydraulic Fracturing

*∫* Γc

displacement vector **u** approximated as

yields

where **u***<sup>i</sup>*

mated by

and *pi*

*w* =**n***<sup>c</sup>*

<sup>Ω</sup>*δε* <sup>T</sup>**σ**d<sup>Ω</sup> - *<sup>∫</sup>*

*∫* Γc

*δp* <sup>T</sup>*w*˙ dΓ + *∫*

**u**≈**u** ^ <sup>=</sup> <sup>∑</sup> *i*=1 *n* **N***i*

*p* ≈ *p* ^ <sup>=</sup> <sup>∑</sup> *i*=1 *n Ni*

**w**≈**w** ^ <sup>=</sup> <sup>∑</sup> *i*=1 *n* **N***i*

are the nodal displacement and pressure, **N***<sup>i</sup>*

and the fluid pressure *p* similarly approximated by

displacement and fluid pressure shape functions.

Γc

<sup>+</sup> <sup>=</sup> - **<sup>p</sup>***<sup>c</sup>*

<sup>Ω</sup>*δ<sup>u</sup>* <sup>T</sup>**b**d<sup>Ω</sup> - *<sup>∫</sup>*


), or **w**=**n***<sup>c</sup>* ⋅ (**u***<sup>c</sup>*

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,

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

The crack opening displacement **w** (or more generally displacement discontinuity) is approxi‐

∇<sup>T</sup> (*δp*)**k**∇ *p*dΓ + *∫*

*δu* <sup>T</sup>**t**dΓ - *∫*

Γt


<sup>+</sup> - **<sup>u</sup>***<sup>c</sup>* -

Γc

Γc

*δp* T(*w*˙ - ∇<sup>T</sup> (**k**∇ *p*) + *g*)dΓ =0 (16)

*<sup>u</sup>***u***<sup>i</sup>* <sup>=</sup>*N***uu**˜, *<sup>δ</sup>***u**≈*N***u***δ***u**˜ (18)

*<sup>p</sup> pi* <sup>=</sup>*N*p**p**˜, *<sup>δ</sup><sup>p</sup>* <sup>≈</sup>*N*p*δ***p**˜ (19)

*<sup>u</sup>* and *Ni*

*<sup>w</sup>***u***<sup>i</sup>* <sup>=</sup>*Nw***u**˜, *<sup>δ</sup>***w**≈*Nwδ***u**˜ (20)

<sup>+</sup> (13)

) ⋅**n***<sup>c</sup>* (14)

*δw*T**p**dΓ =0 (15)

*δp* <sup>T</sup>*g*dΓ =0 (17)

*<sup>p</sup>* are corresponding nodal

$$\mathbf{C} = \mathbf{Q}^{\mathrm{T}} = \mathbf{f}\_{\Gamma\_{\varepsilon}} \{ \mathbf{N}^{p} \}^{\mathrm{T}} \mathbf{n}^{\mathrm{T}} \mathbf{N}^{\mathrm{w}} d\Gamma, \ \mathbf{H} = \mathbf{f}\_{\Gamma\_{\varepsilon}} \{ \nabla \mathbf{N}^{p} \}^{\mathrm{T}} \mathbf{k} \nabla \mathbf{N}^{p} d\Gamma,\\ f^{p} = -\mathsf{f}\_{\Gamma\_{\varepsilon}} \{ \mathbf{N}^{p} \}^{\mathrm{T}} \mathbf{g} d\Gamma \tag{26}$$

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

$$
\begin{bmatrix} 0 & 0 \\ \mathbf{C} & 0 \end{bmatrix} \begin{bmatrix} \dot{\tilde{\mathbf{u}}} \\ \dot{\tilde{\mathbf{p}}} \end{bmatrix} + \begin{bmatrix} \mathbf{K} & -\mathbf{Q} \\ 0 & \mathbf{H} \end{bmatrix} \begin{bmatrix} \tilde{\mathbf{u}} \\ \tilde{\mathbf{p}} \end{bmatrix} = \begin{pmatrix} f^{\mu} \\ f^{\rho} \end{pmatrix} \tag{27}
$$

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}(\mathbf{x}) = \sum\_{I \in \mathcal{N}} \mathbf{N}\_I(\mathbf{x}) \mathbf{u}\_I + \sum\_{I \in \mathcal{N}\_{\mathcal{U}}} \widetilde{\mathbf{N}}\_I(\mathbf{x}) [H(\mathbf{x}) \cdot H(\mathbf{x}\_I)] \mathbf{a}\_I + \sum\_{I \in \mathcal{N}\_{\mathcal{U}}} \widetilde{\mathbf{N}}\_I(\mathbf{x}) \sum\_{l=1}^4 \left\{ \mathbf{B}^{\{l\}}(r, \Theta) \cdot \mathbf{B}^{\{l\}}(r\_l, \Theta\_l) \right\} \mathbf{b}\_I^{\{l\}} \tag{28}$$

where *N* is the set of all nodes in the mesh, *Ncr* the set of nodes whose support are bisected by the crack surface Γ*c*, *Ntip* the set of nodes whose support are partially cut by the crack surface, **N***I* (**x**) and **N**˜ *I* (**x**) are the standard finite element shape functions, **u***<sup>I</sup>* are the displacement nodal degrees of freedom, **a***I* and **b***<sup>I</sup>* (*l*) are the additional degrees of freedom for the displacement, and *H* (**x**) and *B*(*l*) (*r*, *θ*) are the appropriate enrichment basis functions which are localized by **N**˜ *I* (**x**). The shape function **N**˜ *I* (**x**) can differ from **N***<sup>I</sup>* (**x**).

The displacement discontinuity given by a crack Γ*c* can be represented by the generalized Heaviside step function

$$H(\mathbf{x}) = H(d(\mathbf{x})) = \operatorname{sign}(d(\mathbf{x})) = \begin{cases} 1 & d(\mathbf{x}) \ge 0 \\ -1 & d(\mathbf{x}) < 0 \end{cases} \tag{29}$$

*p*(**x**)= ∑ *I*∈*Ncr NI*

degrees of freedom **u***<sup>I</sup>* . The additional degrees of freedom **a***I* and **b***<sup>I</sup>*

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

the crack tip field enriched degrees of freedom **b***<sup>I</sup>*

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

where *NI*

regime in a hydraulic fracture [18].

**4.2. Element implementation**

*<sup>p</sup>*(**x**)*pI* **<sup>x</sup>**∈Γ*<sup>c</sup>* (32)

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

(*l*)

are assigned to the four

http://dx.doi.org/10.5772/56287

733

T

18

*<sup>p</sup>*(**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

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

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

1 1 2 2 3 3 4 41 1 2 2 3 3 4 4 3 4

and for the tip element the Heaviside enriched degrees of freedom **a***<sup>I</sup>* need to be replaced by

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

<sup>ˆ</sup> *x x x xx y y y yy y y y x x <sup>e</sup> <sup>x</sup> u u u u u u u ua a a a a a a ap p <sup>u</sup>*ìï üï <sup>=</sup> <sup>í</sup> <sup>ý</sup> ïî ïþ

Standard Heaviside Enriched Coupled

(*l*) .

14444444244444443 14444444244444443 14243 (33)

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:

$$\begin{aligned} \left[ \mathbb{B}^{(0)} \right]\_{l=1}^{4} &= \sqrt{r} \left[ \sin \left( \theta / 2 \right) \quad \cos \left( \theta / 2 \right) \quad \sin \left( \theta / 2 \right) \sin \left( \theta \right) \quad \cos \left( \theta / 2 \right) \sin \left( \theta \right) \right] \end{aligned} \tag{30}$$

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

$$\mathbf{u}(\mathbf{x}) = \mathbf{u}\_c^+(\mathbf{x}) \cdot \mathbf{u}\_c^-(\mathbf{x}) = 2 \sum\_{I \in \mathcal{N}\_{\mathcal{O}}} \widetilde{\mathbf{N}}\_I(\mathbf{x}) \mathbf{a}\_I + 2 \sum\_{I \in \mathcal{N}\_{\mathcal{O}}} \widetilde{\mathbf{N}}\_I(\mathbf{x}) \mathbf{B}^{(1)}(r, \pi) \mathbf{b}\_I^{(1)} \mathbf{x} \in \Gamma\_c \tag{31}$$

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

The fluid pressure field within the crack is approximated by

$$p(\mathbf{x}) = \sum\_{I \in \mathcal{N}\_{\mathcal{O}}} N\_I^{\mathcal{P}}(\mathbf{x}) p\_I \quad \mathbf{x} \in \Gamma\_{\mathcal{E}} \tag{32}$$

where *NI <sup>p</sup>*(**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].

#### **4.2. Element implementation**

**4. The extended finite elemet method and element implementation**

(**x**)(*<sup>H</sup>* (**x**) - *<sup>H</sup>* (**x***<sup>I</sup>* ))**a***<sup>I</sup>* <sup>+</sup> <sup>∑</sup>

(**x**) can differ from **N***<sup>I</sup>*

(*l*)

*I*

where *d*(**x**) is the signed distance of the point **x** to Γ*c*.

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


The fluid pressure field within the crack is approximated by

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

where *N* is the set of all nodes in the mesh, *Ncr* the set of nodes whose support are bisected by the crack surface Γ*c*, *Ntip* 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

crack tip, which are generally chosen as a basis that approximately spans the two-dimensional

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

*I*∈*Ntip* **N**˜ *I* (**x**)*B*(1)

(**x**)**a***<sup>I</sup>* + 2 ∑

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

<sup>4</sup> = *r*{sin (*θ* / 2) cos (*θ* / 2) sin (*θ* / 2)sin (*θ*) cos (*θ* / 2)sin (*θ*)} (30)

*<sup>H</sup>* (**x**)=*<sup>H</sup>* (*d*(**x**))=*sign*(*d*(**x**))={ <sup>1</sup> *<sup>d</sup>*(**x**)≥<sup>0</sup>

plane strain asymptotic crack tip fields in the linear elastic fracture mechanics:

*I*∈*Ntip* **N**˜ *I* (**x**)∑ *l*=1 4 (*B*(*l*)

(**x**) are the standard finite element shape functions, **u***<sup>I</sup>* are the displacement nodal

(*r*, *θ*) are the appropriate enrichment basis functions which are localized by

(**x**).

are the additional degrees of freedom for the displacement, and

(*r*, *θ*) - *B*(*l*)


(*r*, *θ*) are required to model the displacement around the

(*r*, *π*)**b***<sup>I</sup>* (1)

**x**∈Γ*<sup>c</sup>* (31)

(*rI* , *θ<sup>I</sup>* ))**b***<sup>I</sup>*

(*l*) (28)

**4.1. Extended finite element approximation**

(**x**)**u***<sup>I</sup>* + ∑ *I*∈*Ncr* **N**˜ *I*

732 Effective and Sustainable Hydraulic Fracturing

[17]

**N***I*

**N**˜ *I*

**u**(**x**)= ∑ *I*∈*N* **N***I*

(**x**) and **N**˜

*H* (**x**) and *B*(*l*)

*I*

Heaviside step function

{*B*(*l*) }*l*=1

can be obtained as

**w**(**x**)=**u***<sup>c</sup>*

+(**x**) - **<sup>u</sup>***<sup>c</sup>*

degrees of freedom, **a***I* and **b***<sup>I</sup>*

(**x**). The shape function **N**˜

The enrichment basis functions *B*(*l*)

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***<sup>I</sup>* . The additional degrees of freedom **a***I* and **b***<sup>I</sup>* (*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.

**Figure 3.** 2-D 4-node -node plane strain hydraulic fracture elements

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

$$\hat{\boldsymbol{\mu}}^{\varepsilon} = \left\{ \underbrace{\boldsymbol{u}\_{1}^{\boldsymbol{x}} \quad \boldsymbol{u}\_{1}^{\boldsymbol{y}} \quad \boldsymbol{u}\_{2}^{\boldsymbol{x}} \quad \boldsymbol{u}\_{2}^{\boldsymbol{y}} \quad \boldsymbol{u}\_{3}^{\boldsymbol{y}} \quad \boldsymbol{u}\_{4}^{\boldsymbol{y}} \quad \boldsymbol{u}\_{4}^{\boldsymbol{y}} \quad \boldsymbol{u}\_{4}^{\boldsymbol{y}} \quad \boldsymbol{u}\_{4}^{\boldsymbol{y}} \quad \underbrace{\boldsymbol{a}\_{1}^{\boldsymbol{x}} \quad \boldsymbol{a}\_{1}^{\boldsymbol{y}} \quad \boldsymbol{a}\_{2}^{\boldsymbol{x}} \quad \boldsymbol{a}\_{2}^{\boldsymbol{x}} \quad \boldsymbol{a}\_{3}^{\boldsymbol{y}} \quad \boldsymbol{a}\_{4}^{\boldsymbol{x}} \quad \boldsymbol{a}\_{4}^{\boldsymbol{y}} \quad \boldsymbol{a}\_{4}^{\boldsymbol{x}} \quad \boldsymbol{a}\_{4}^{\boldsymbol{y}} \quad \underbrace{\boldsymbol{p}\_{3}}\_{\text{Coulomb}} \right\}\_{18}^{\mathrm{T}} \tag{33}$$

and for the tip element the Heaviside enriched degrees of freedom **a***<sup>I</sup>* need to be replaced by the crack tip field enriched degrees of freedom **b***<sup>I</sup>* (*l*) .

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.

(a) (b)

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

http://dx.doi.org/10.5772/56287

735

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

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

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

tip user-defined element that captures the crack tip singularity correctly.

**Figure 5.** Zero-toughness case: (a) imposed pressure; and (b) imposed opening

along the crack length in the finite element model.

**6. Summary**

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

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems http://dx.doi.org/10.5772/56287 735

**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 zerotoughness 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**

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

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

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

Comparisons of the FEM predictions with the available analytical solutions to the two limiting

(a) (b)

side residual vector, as need in a context of solving the discrete system of equations.

result.

**5. Numerical examples**

734 Effective and Sustainable Hydraulic Fracturing

cases are given in Figures 4 and 5.

the analytical solutions are available in the literature [21].

**Figure 4.** Zero-viscosity case: (a) imposed pressure; and (b) imposed opening

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.

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

scale *L*(*t*) take the explicit forms [21]

Ω¯ *m*0 (*ξ*)= *A*<sup>0</sup>

Π*m*<sup>0</sup> (1) <sup>=</sup> <sup>1</sup> <sup>3</sup>*<sup>π</sup> <sup>B</sup>*( <sup>1</sup> <sup>2</sup> , <sup>2</sup>

where *A*<sup>0</sup> =31/2

, *A*<sup>1</sup> (1)

metric function. Thus, Ω¯

**Acknowledgements**

**Author details**

Zuorong Chen\*

ϵ*m*(*t*)=(*μ*'

(1 - *<sup>ξ</sup>*2)2/3 <sup>+</sup> *<sup>A</sup>*<sup>1</sup>

<sup>3</sup> ) *A*0*F*<sup>1</sup> (- 1 <sup>6</sup> , 1; 1 <sup>2</sup> ;*ξ*2) <sup>+</sup>

*m*0

≅ - 0.156, and *B*(1)

(0)≅1.84.

Ω¯ *k* 0

/ *E*'

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

(1)

For the viscosity scaling, denoted by a subscript *m*, the small parameter

*<sup>t</sup>*)1/3 , *Lm*(*t*)= (*E*'

10 <sup>7</sup> *A*<sup>1</sup> (1) *F*1 (- 7 <sup>6</sup> , 1; 1

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.

CSIRO Earth Science & Resource Engineering, Melbourne, Australia

*Q*0 3 *t* <sup>4</sup> / *μ*'

(1 - *ξ*2)5/3 + *B*(1) 4 1 - *ξ*<sup>2</sup> + 2*ξ*<sup>2</sup>

*<sup>γ</sup><sup>k</sup>* <sup>0</sup> =2 / *<sup>π</sup>*2/3 (38)

An ABAQUS Implementation of the XFEM for Hydraulic Fracture Problems

(*ξ*)= 1 - *ξ*<sup>2</sup> / *π*1/3 (39)

<sup>Π</sup>*<sup>k</sup>* <sup>0</sup> <sup>=</sup>*π*1/3 / <sup>8</sup> (40)

*γm*0≅0.616 (42)

<sup>2</sup> ;*ξ*2) <sup>+</sup> *<sup>B</sup>*(1)

≅0.0663; *B* is Euler beta function, and *F*1 is hypergeo‐

ln| <sup>1</sup> - <sup>1</sup> - *<sup>ξ</sup>* <sup>2</sup> 1 + 1 - *ξ* <sup>2</sup>

ϵ

)1/6 (41)

(*t*) and the length

http://dx.doi.org/10.5772/56287

737


(2 - *π*|*ξ*|) (44)
