The XFEM with an Explicit-Implicit Crack Description for Hydraulic Fracture Problems

Abstract

The Extended Finite Element Method (XFEM) approach is applied to the coupled problem of fluid flow, solid deformation, and fracture propagation. The XFEM model description of hydraulic fracture propagation is part of a joint project in which the developed numerical model will be verified against large-scale laboratory experiments. XFEM forms an important basis towards future combination with heat and mass transport simulators and extension to more complex fracture systems. The crack is described implicitly using three level-sets to evaluate enrichment functions. Additionally, an explicit crack representation is used to update the crack during propagation. The level-set functions are computed exactly from the explicit representation. This explicit/implicit representation is applied to a fluid-filled crack in an impermeable, elastic solid and compared to the early-time solution of a plane-strain hydraulic fracture problem with a fluid lag.

1. Introduction

The large scale conversion of geothermal energy into electrical energy using natural formations as heat exchangers depends on the coincidental occurrence of heat, fluid and permeability. This is valid for only a few locations on earth. Enhanced Geothermal System (EGS) propose to engineer the controlled creation of a heat exchanger between two wells in deep hot rocks, increasing the number of possible locations on earth. Water can be let circulate between the two wells, heat up while passing through the hot rock and be cooled down on the surface for power generation. Yet this engineering of the heat exchanger has to be improved such that the outcome can be predicted within specified uncertainties.

The extension to more complex fracture scenarios as well as the integration with other software for risk assessment simulations requires a computer resource moderate modeling of fracture propagation. The extended finite element method (XFEM) forms a good basis for this. It has been applied to various problems within the area of fracture mechanics. The XFEM allows for the consideration of a priori knowledge about the solution of a hydraulic fracture problem into the approximation space through the addition of enrichment functions [10]. It enables, thereby, the accurate approximation of fields that involve jumps, kinks, singularities, and other non-smooth features within elements [2, 6, 11]. The developed numeric model will be verified against large-scale laboratory experiments. However, the focus of the present paper will lie on the progress in using XFEM for hydraulic fracture modeling. An XFEM approach in combination with an explicit and implicit crack description is applied to a plane-strain hydraulic fracture problem. The implicit description is given by three level-set functions defined in [5] and enables a simple evaluation of the enrichments. In contrast, the explicit crack description is used to perform the crack update. Given the explicit interface, the level-set functions for each propagation step can be calculated straightforward.

The paper is organized as follows: After a short description of the laboratory part of the joint project in Section 2, the governing equations for a hydraulic fracture problem in its basic form are presented in Section 3. Models are discussed for the solid deformation, fluid flow, and fracture propagation. In Section 4, an XFEM formulation with an explicit-implicit interface description is introduced and the discretization of these governing equations is carried out. Numerical results are presented in Section 5. Finally, Section 6 concludes this paper.

2. Laboratory experiments for model verification and optimization

Laboratory experiments for model verification and optimization will be performed on large rock specimen. The large-scale testing facility is currently under construction. Meanwhile preliminary experiments on smaller specimens have been conducted. In our project we focus on fracture creation in basement rock. Therefore mainly rocks like basalt, granite and gneiss are going to be tested for model verification and optimization.

2.1. Large-scale

Blocks of size 30 x 30 x 45 cm³ will be pre-stressed in a massive steel frame to set up realistic primary stress states before hydraulic stimulation (see Figure 1). Stresses in all three directions can be adjusted independently and will be held constant during injection time. After setting up the primary principal stress state with flat-jacks, the injection interval of the borehole will be pressurized by a syringe pump. Injection pressures up to 65 MPa are possible. In order to allow for the verification of the developed numerical model with the experiments, we are going to monitor the borehole-pressure, the deformation of the rock sample and localize acoustic events occurring during crack formation and propagation by means of ultrasonic transducers. Material parameters will be derived from standard rock mechanical tests.

Figure 1.

Large- scale tri-axial testing facility under construction, able to apply up to 30 MPa vertically and 15 MPa horizontally

2.2. Small-scale

Preliminary tests were performed on concrete samples of smaller size (15 x 15 x 15 cm³) and recently been extended to granite and basalt. Acoustic events were recorded during fracture creation and propagation. The experiments indicated the need to lower the compressive energy induced before breakdown. Further, instead of water and light oil, fluids of higher viscosity will be used from now on to enlarge the regime of fracture propagation (see Figure 2c) and to consider the lag of scaling. Optimization of acoustic emission monitoring is continuously ongoing.

Figure 2.

Preliminary small-scale fracturing experiments. a) Fractured concrete sample. b) Located acoustic events, projected onto face D of the sample. Different colors correspond to different time of occurrence. Dashed lines correspond to minimum, mean and maximum (along the direction of projection) height of the two main fracture surfaces, visible on the photo to the left. c) Fluid pressure (black) and flow rate (blue) curve used to fracture the specimen. Colored regions show the time and pressure regimes during which the same-colored acoustic events to the left occurred.

3. Governing equations

Hydraulic fracture propagation is based at least on three physical processes. A fluid flow within the fracture imposes a pressure load on the fracture surfaces. As a result, the rock undergoes a (mechanical) deformation and the fracture starts to propagate when a critical condition is reached [1]. Depending on the different modeling assumptions, this critical condition can be defined by the fracture toughness or another stress-based criterion. The following assumptions are usually made for the hydraulic fracture model [1]: I) the fluid flow is governed by the lubrication theory, II) solid deformation is modeled using the theory of linear elasticity, and III) the propagation criterion is given by the conventional energy- release-rate approach of linear elastic fracture mechanics (LEFM) theory. The crack propagates when the mode I stress intensity factor reaches the fracture toughness. Each physical process is modeled separately and coupled iteratively. The governing equations are given as follows:

3.1. Deformation

A homogenous, isotropic and linear elastic solid is modeled with the concept of equilibrium of forces. Far field stresses and the pressure on the interface are imposed as Neumann boundary conditions, body forces are neglected.

$$\begin{array}{llll}\nabla \cdot \sigma +f\hfill & =0\hfill & \hfill & \text{in }\Omega \hfill \\ u\hfill & =u_0\hfill & \hfill & \text{in }{\Gamma }_d\hfill \\ \sigma \cdot n\hfill & =\widehat{t}\hfill & \hfill & \text{in }{\Gamma }_n\hfill \end{array}$$E1

f denotes the body force vector and u the displacement defined on the region Ω. The traction $\widehat{t}$ is applied at the outer boundary $${\Gamma }_n$$ and Dirichlet boundary conditions are defined on $${\Gamma }_d$$. Hooke’s law of elasticity gives the relation between the stress $$\sigma$$ and the strain $$\epsilon$$

$${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$$E2

where C is the fourth-order stiffness (elasticity) tensor. Since the fracture aperture w is not given directly in this formulation, it has to be determined from the displacement field.

3.2. Fluid flow

The fracturing fluid with the dynamic viscosity $$\mu$$ is modeled as laminar flow between two parallel plates with a constant injection rate Q₀. The fluid flux q then reads

Figure 3.

Sketch of a plane-strain hydraulic fracture with varying aperture w and fluid front position L_f. A fluid lag is shown at the fracture tip. The fluid is injected at the wellbore and flows into the fracture at a constant rate Q₀. It can leak off into the matrix through the fracture surfaces at a rate q_l. Fluid pressure field is denoted by p.

$$q=-\frac{w^3}{12\mu }\nabla p.$$E3

The Reynolds (lubrication) equation is given by

$$\frac{\partial w}{\partial t}-\frac{1}{12\mu }\nabla \cdot \left(w^3\nabla p\right)+q_l=Q_0\delta (x,y)$$E4

and describes the conservation of the fluid mass for a Newtonian fluid. The fluid is injected into the fracture at a constant rate Q₀. For a fracture propagating in an impermeable solid, the leak-off q_l is negligible and, therefore, set to zero. It is assumed that a fluid lag develops between the fluid front L_f and the crack tip. However, for reasons of simplicity the lag size is not part of the solution. Taking into account the symmetry of the problem, the boundary conditions for the fluid flow problem read as follows:

$q\text{=}\frac{{\text{Q}}_{\text{0}}}{\text{2}}at\mathrm{the}fractureinlet$E5

$q=q_{f}atthefluidfront$E6

$p=p_{0}inthefluidlag$E7

The flow condition at the fluid front q_f is determined from the pressure gradient and thus, is part of the solution. The pressure in the lag region is set to a constant value p₀, that is usually chosen to be zero. Finally, the global volume balance condition

$$Q_{0}t={\displaystyle {\int }_{{\Gamma }_f}w}\text{d}\Gamma +{\displaystyle {\int }_0^t{\displaystyle {\int }_{{\Gamma }_f}{q}_l}}\text{d}\Gamma \text{d}\tau =V$$E8

equates the fracture volume V to the volume of injected fluid and the amount lost to the surrounding rock-mass. The integration is performed over the fluid filled part of the crack $${\Gamma }_f$$.

3.3. Propagation condition

Due to the symmetry in loading and geometry, the hydraulic fracture propagates in pure opening mode, i.e. the tensile stress is acting normal to the plane of the crack. The propagation criterion is formulated in the framework of linear elastic fracture mechanics (LEFM) and accounts for the energy required to break the rock. It is characterized by the stress singularity at the tip and a propagation in mobile equilibrium. The LEFM assumption requires the stress intensity factor K_I to be equal to the fracture toughness K_IC of the material [12]

$$K_I=K_{IC}.$$E9

3.4. Asymptotic behavior

The hydraulic fracture problem characterized by a strong fluid-solid coupling that is mainly confined to a small region near the crack tip where rapid variation in the fluid pressure occurs. Analyzing the physical process at the tip by comparing the work done by the fluid in extending a fracture with the energy required to create new crack surfaces leads to understanding of the propagation regime of a fluid-driven fracture. Two limiting regimes can be detected, a toughness- and a viscosity-dominated regime [3]. In the toughness-dominated regime the inverse square root singularity of LEFM captures the effect of the crack tip process on the total fracture. In contrast, the viscosity-dominated regime is characterized by a singularity that is weaker than the singularity predicted by LEFM. Fracture toughness K_IC may become irrelevant [9].

4. XFEM approximation

The extended finite element method (XFEM) allows for the consideration of a priori knowledge about the solution of a hydraulic fracture problem into the approximation space through the addition of enrichment functions [10]. It enables, thereby, the accurate approximation of fields that involve jumps, kinks, singularities, and other non-smooth features within elements [2, 6, 11].

The enrichments, that are realized through the partition of unity (PU) concept, are chosen in such a way that they are able to reproduce the asymptotic behavior near the crack tip. In this work only the toughness-dominated solution is considered. Thus, enrichment functions compatible with the classical square root singularity of LEFM are used to enrich the region near the crack tip.

4.1. Standard formulation

The XFEM formulation with an explicit-implicit crack description used in this work is based on the work done by [5]. The basic idea is recalled in this paper, for further details see the original work. The enriched approximation of the displacements is stated as follows:

$$u^h(x)=\underset{\text{continuous}}{\underbrace{{\displaystyle \sum _{i\in I}{N}_i}(x)u_i}}+\underset{\text{discontinuous}}{\underbrace{{\displaystyle \sum _{j\in I^{step}}{N}_j^{\text{*}}(x)}\cdot {\Psi }_{step}(x)a_j+{\displaystyle \sum _{k\in I^{tip}}{N}_k^{\text{*}}(x)}\cdot ({\displaystyle \sum _{m=1}^4{\Psi }_{tip}^m(r,\theta )}{b}_k^m)}}$$E10

The first term on the right hand side describes the classical FEM-approximation with continuous shape functions N_i(x) and nodal unknowns u_i. The second term accounts for the discontinuity in the displacement field across the crack path by incorporating step-functions Ψ_step with additional nodal unknowns a_i into the enrichment space. The tip region is enriched with a set of enrichment functions $${\Psi }_{tip}^m(r,\theta )$$that consider the singularity according to the dominating regime. They can be defined as [10]

$${\left\{{\Psi }_{tip}^m(r,\theta )\right\}}_{m=1}^4=\left\{r^{\lambda }cos(\lambda \theta ),r^{\lambda }sin(\lambda \theta ),r^{\lambda }sin(\theta )sin(\lambda \theta ),r^{\lambda }sin(\theta )cos(\lambda \theta )\right\}$$E11

where r and θ denote local polar coordinates at the crack tip. When propagating in the toughness-dominated regime the assumption of a square root singularity in LEFM requires to choose λ = 1/2. In the viscosity-dominated regime the weaker singularity is taken into account by λ = 2/3. Additional degrees of freedom $$b_k^m$$ are introduced into the approximation locally within the enriched region.

Figure 4.

The enrichment is acting either along the crack path (dashed field) with the step-enrichment $${\Psi }_{step}$$ or within a specified region near the crack tip by defining the enrichments $${\Psi }_{tip}^m(r,\theta )$$.

The crack opening is obtained through interpolation of the displacement field u(x) at the interface nodes by means of (10). Since the interface represents a discontinuity the interpolation is performed by moving the nodes slightly away in normal direction. The opening is defined as the distance between the positive and negative displacement at the interface (see Figure 5).

$$w(x)=u^{+}(x)-u^{-}(x).$$E12

Figure 5.

Interpolation of the crack opening along the interface.

4.2. Numerical integration

The standard approach in the XFEM for numerical integration is a decomposition of the elements into subelements that align with the discontinuity [6]. A Gauss quadrature is then applied on each of these subelements. For a detailed description of the decomposition method in 2D and 3D the reader is referred to [6].

4.3. Explicit-implicit interface description

The explicit crack description is given by a mesh that is aligned with the interface. For 2D problems the crack is a line and is represented by one dimensional elements in the 2D space. In three dimensions, the crack is a surface and, thus, described through a two dimensional mesh in the 3D space.

Normal and tangential vectors are computed easily on the interface and can be used to define a local coordinate system at the crack tip/front. On the basis of the explicit interface mesh the crack update is applied by simply adding new elements to the crack front according to a given extension vector.

Figure 6.

(a) Arbitrary crack surface with normal vectors on each element. (b) Local coordinate system at the crack front. (c) Crack update according to crack extension vectors at the front.

The implicit interface description is realized by means of the level-set concept. Three level-set functions are defined according to [5]. They are used to define the region to be enriched and to evaluate the enrichment functions.

• $${\Phi }_1(x)$$ is the (un-signed) distance function to the crack path/surface. That is, the level-set value at position x is the shortest distance to the crack path/surface.

• $${\Phi }_2(x)$$ is the (un-signed) distance function to the crack tip(s)/front. That is, the level-set value at position x is the shortest distance to the crack tip(s)/front.

• $${\Phi }_3(x)$$ is a signed distance function to crack path/surface that is extended over the entire domain. The sign is based on the direction of the normal vector of the segment that contains the nearest point.

4.4. Discretization of governing equations

Since the solid deformation and the fluid flow are coupled iteratively, they are solved independently in each iteration step. Solid deformation is discretized with XFEM as follows:

$$\left[{\displaystyle \underset{\Omega }{\int }{B}^{T}CB}d\Omega \right]\cdot u={\displaystyle \underset{{\Gamma }_b}{\int }{N}^T}\widehat{t}d\Gamma +{\displaystyle \underset{{\Gamma }_c}{\int }{N}^T}pd\Gamma$$E13

where the term on the left B^TCB denotes the stiffness matrix with the gradient operator B [4], N the shape and enrichment functions, $\widehat{t}$ the traction on the outer boundary $${\Gamma }_b$$ and p the pressure on the interface $${\Gamma }_c$$. A classical FEM approach is used to solve the fluid flow equation. The pressure is approximated by

$$p^h(x)={\displaystyle \sum _{i\in I}{N}_i}(x)p_i.$$E14

The discretized problem formulation reads

$$\frac{1}{12\mu }{\displaystyle \underset{{\Gamma }_f}{\int }{w}^3\nabla N\nabla N^{T}d\Gamma }={\left[N\frac{w^3}{12\mu }\nabla p\right]}_{{\Gamma }_f}-{\displaystyle \underset{{\Gamma }_f}{\int }N{q}_{l}d\Gamma }-{\displaystyle \underset{{\Gamma }_f}{\int }N\frac{\Delta w}{\Delta t}d\Gamma }$$E15

This formulation is valid for one half of a symmetric crack where $${\Gamma }_f$$ denotes the fluid filled region. The flow boundary conditions at the fluid front and the fracture inlet correspond to the first term on the right-hand side. Fluid leak-off and the change of volume over time are taken into account by the second and third term on the right-hand side.

5. Hydraulic fracture propagation

The problem of a fluid driven fracture in an impermeable elastic solid with a fluid lag is considered here. Simulation results are compared with the asymptotic solutions for zero underpressure/time given in [8]. This solution corresponds to the “beginning” of the fluid-driven fracture evolution and provides initial condition for plane-strain fracture propagations. The propagation regime of a fluid driven fracture is controlled by a parameter representing a dimensionless viscosity M (dimensionless toughness K) defined as

$$M=\frac{{\mu }^{\prime }{Q}_0}{E^{\prime }}{\left(\frac{E^{\prime }}{K^{\prime }}\right)}^4,K=M^{-1/4}.$$E16

This formulation uses effective parameters [6]

$$\begin{array}{l}{E}^{\prime }=\frac{E}{1-{\nu }^2},{\mu }^{\prime }=12\mu ,K^{\prime }=4{\left(\frac{2}{pi}\right)}^{1/2}{K}_{IC},\hfill \end{array}$$E17

where μ′ denotes the fluid viscosity, Q₀ the constant injection rate, E′ the plane-strain elastic modulus with Poisson’s ratio ν and K′ the toughness, respectively. The procedure solving the coupled equations follows that described in [8]. Given a fluid front position L_f, a solution is sought for the pressure distribution and crack opening.

5.1. Numerical algorithm

The simulation process is realized through an iterative coupling of the fluid flow and solid deformation. Starting with an initial solution and a guess for the fluid fraction, the pressure distribution and the crack opening are calculated until convergence is reached. When the propagation condition is met, the crack is updated for the next time step. Otherwise, the fluid front is moved towards the crack tip with a velocity v determined from the fluid flow rate q.

5.2. Numerical results

The numerical results for the crack opening and pressure distribution at the wellbore of a plane-strain hydraulic fracture problem are compared to the similarity solution of a small enough toughness parameter in order to allow a significant fluid lag. The boundary condition of zero displacement at infinity is approximated by a finite body and standard finite elements and a local mesh refinement in the area close to the crack interface. Computational evidence of the validity of approximating the infinite medium with a finite block is provided in [13].

Figure 7.

The numerical results (blue circles) of dimensionless pressure Π (a) and the crack opening Ω (b) at the wellbore as well as the dimensionless crack length γ (c) are compared to the analytical solution (red solid line) for various values of the fluid front position.

The results are scaled to dimensionless quantities in the viscosity scaling. For a detailed description of the scaling for the pressure Π, the opening Ω and the crack length γ see the original publication [13]. The domain and the explicit interface are meshed independently with 5000 and 3000 elements, respectively.

Figures 7(a)-(c) show a good agreement of the similarity solution for various values of the fluid fraction $${\xi }_f=L_f/L$$. However, especially for high fluid fraction values where the fluid front is close to the fracture front the results reveal inaccuracies. Special attention has to be paid to the crack tip behavior in the case for a vanishing fluid lag when the pressure becomes singular. Depending on the propagation regime crack propagation is governed either by the classical singularity of linear elastic fracture mechanics or by viscous fluid effects which would lead to a weaker singularity than given by LEFM.

Figure 8.

The pressure distribution Π(ξ) (a) and the crack opening profile Ω(ξ) (b) for various values of the fluid fraction ξ_f.

The pressure distribution and the crack opening profile along the dimensionless coordinate $$\xi =x/L$$ are shown in Figures 8(a) and (b) in the viscosity scaling for fluid fraction values ξ_f = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.97, 0.99}.

6. Conclusions

The XFEM with an explicit-implicit crack description has been applied to a plane-strain hydraulic fracture problem. The crack is described explicitly by a line (2D)/triangular (3D) mesh that is aligned with the interface and implicitly by three level-set functions. The enrichment functions at the tip can be chosen according to the asymptotic behavior of the hydraulic fracture problem. Depending on the propagation regime the stress singularity can be described either by LEFM or by a singularity, which is weaker than predicted by LEFM. However, in this work a partially filled crack with a significant lag is examined and, therefore, crack propagation is governed by LEFM. The results show a good agreement with the known similarity solutions and can be interpreted as an early-time solution that can be used as a starting point in hydraulic fracture simulations.

Acknowledgments

We acknowledge support for this work from the Federal Ministry for the Environment, Nature Conservation and Nuclear Safety, Germany (FKZ 0325167).