# Extended Finite Element Method for Predicting Productivity of Multifractured Horizontal  Wells

Authors: Youshi Jiang; Jinzhou Zhao; Yongming Li; Hu Jia; Liehui Zhang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Based on the theory of the extended finite element method (XFEM), which was first proposed by Moës for dealing with the problem characterized by discontinuities, an extended finite element model for predicting productivity of multifractured horizontal well has been established. The model couples four main porous flow regimes, including fluid flow in the away-from-wellbore region of reservoir matrix, radial flow in the near-wellbore region of reservoir matrix, linear flow in the away-from-wellbore region of fracture, and radial flow in the near-wellbore region of fracture by considering mass transfer between fracture and matrix. The method to introduce the interior well boundary condition into the XFEM is proposed, and therefore the model can be highly adaptable to the complex and asymmetrical physical conditions. Case studies indicate that this kind of multiflow problems can be solved with high accuracy by the use of the XFEM.


---

## Body

## 1. Introduction

Prediction of productivity of multifractured horizontal wells has always been a research hotspot. Various methods including the analytical method [1, 2], the hybrid numerical-analytical method [3, 4], and the numerical method [5, 6] have been utilized to tackle that problem, but all of them have shortages. The analytical method is simple and convenient for engineering applications, but it is usually based on the oversimplified model and many questionable assumptions which may cause erroneous results. Both the hybrid numerical-analytical method and the numerical method (including the classical FEM) have advantages for problems with multiscale flow in the fractured reservoir. However, the grid in the local domain surrounding the fracture must be defined to get the satisfactory result and then results in the unexpected computational consumption.In recent decades, the XFEM has been continuously developing and already has become an important and effective technique for the analysis of the problems characterized by discontinuities, singularities, and complex geometries, such as modeling crack growth and complex flow in fractured reservoir [7–9]. It has many advantages over other numerical methods, and two of the advantages based on the theory of partition of unity [10] are the most important. The first is that the mesh is simplified and independent of the morphology of the fractures. The second is that the local field can be accurately reproduced. Belytschko and Black first present the essential ideas of modeling crack and crack growth by finite elements with no remeshing and introducing a discontinuous enrichment function into the approximation [11]. Moës et al. further improved Belytschko and Black’s work by bringing in the step function and the near-tip asymptotic functions and named it the extended finite element method [12]. Since then, many a research has done a lot to develop the XFEM [13–16].Naturally, the XFEM was developed for simulation of strong discontinuities in fracture mechanics. It was later extended to solve problems for weak discontinuity such as interface problems and fluid flow in fractured reservoir. Sukumar et al. presented a method that coupled the level set method to the XFEM to predict the weak discontinuity due to material interfaces in composite [17]. Lamb et al. combined the dual permeability method (DPM) and the XFEM to analyze coupled deformation and fluid flow in fractured porous media [18]. To model different subsurface features in porous media flow problem, Huang et al. proposed several corresponding enrichment shape functions [19]. Mohammadnejad and Khoei developed a fully coupled numerical model for modeling hydraulic fracture propagation using the XFEM [20]. Sheng et al. built an XFEM model for multiscale flow in fractured shale gas reservoir which accounted for mass transfer between matrix and fractures, four kinds of flow regimes, and pressure dependence [21]. However, they totally neglected the radial flow regime in the near-wellbore region of fracture to the horizontal wellbore. In addition, the formers usually located the wellbore on the outer boundary of the whole domain rather than its interior according to the assumption of symmetry which did not always make sense, because both the complex fracture network and the reservoir were more likely to be asymmetrical in relation to the wellbore.This paper aims to establish an extended finite element model for predicting of productivity of multifractured horizontal well, which couples fluid flow in the away-from-wellbore region of the porous matrix, radial flow in the near-wellbore region of the matrix, linear flow in the away-from-wellbore region of fracture, and radial flow in the near-wellbore region of fracture by considering mass transfer between fracture and matrix. The method to load the interior well boundary condition is also proposed, and therefore the model can be completely adaptable in the complex and asymmetrical physical conditions.


## 2. Governing Equations

In this section we demonstrate the derivation of the governing equations for fluid flow in multifractured horizontal well. Firstly, we sketch the physical model and the main assumptions.


### 2.1. Physical Model and Assumptions

#### 2.1.1. Physical Model

Consider a two-dimensional rectangular domainΩ bounded by the outer boundary Γ (Figure 1). A fractured horizontal well is located in the domain. The whole fracture reflects the weak discontinuity Γd on which the point O represents the wellbore. Oil in reservoir initially flows within the porous matrix and then from the matrix into the two-wing fracture and finally from the fracture into the wellbore. Because the width of the fracture is far less than its length, the fracture seems more like a discontinuous line in the view of the whole domain. But on a local scale the fracture is actually two-dimensional and is enveloped by two boundary lines (Γd+ and Γd-). The fluid is single-phase. Every flow regime is the isothermal event and obeys Darcy’s law. The internal boundary Γw relates to the wellbore.Figure 1
Fluid flow in hydraulically multi-fractured reservoir.n
Γ
c is the unit vector that points from the wellbore to the fracture tip. nΓc is the unit normal vector to the discontinuity pointing to Ω+(nΓd=nΓd-=-nΓd+). nΓ is the unit outward normal vector to the external boundary.Fluid flow in fractured horizontal well involves four main kinds of flow mechanisms, including fluid porous flow in the away-from-wellbore region of reservoir matrix (main reservoir flow), radial flow in the near-wellbore region of reservoir matrix (reservoir radial flow), linear flow in the away-from-wellbore region of fracture (fracture linear flow), and radial flow in the near-wellbore region of fracture (fracture radial flow). Different flow regimes are coupled by considering mass transfer between fracture and matrix (Figures2-3).Figure 2
The local Cartesian coordinate systems in the fracture.Figure 3
Flow regimes in the hydraulically fractured reservoir.The essential boundary conditions are imposed on the external and the internal boundary of the model as(1)p=peon⁡Γ,p=pwfon⁡Γw.Mass transfer is considered to couple the fluid flow in the fracture and the porous matrix. Regarding this, the complementary boundary condition is written as(2)⟦W˙w⟧nΓd=qwdon⁡Γd.⟦
*
⟧
=
*
+
-
*
- represents the difference between the corresponding values at the two fracture faces.


#### 2.1.2. The Main Assumptions of the Model Include the Following

(a) Fluid flow in the physical model is in steady state; (b) fracture radial flow exists within a circle in the vertical plane with the center at the wellbore and diameter equal to reservoir thickness; (c) reservoir radial flow region exists in the near-wellbore region surrounding the near-wellbore fracture, and its range is similar to fracture radial flow; (d) the fluid pressure within the horizontal wellbore keeps constant; (e) the fluid pressure within the fracture keeps constant along the inward direction that is perpendicular to the fracture line, while the pressure gradient and the fluid flow are discontinuous.Based on the forth main assumption, this paper uses the same pressure field and the same mesh to approximate the fluid pressure in both the fracture and the matrix. This treatment is different from that in [18, 21].


### 2.2. Strong Form

In order to derive the weak form of the governing equations, we must first constitute the respective continuity equations of the four flow regimes as the strong-form equations.The continuity equation for the porous flow in the away-from-wellbore region of the reservoir matrix is written as(3)∇W˙w=∇(-kmμ∇p)=-kmμ(∂2p∂x2+∂2p∂y2)=0.The continuity equation for the linear flow in the away-from-wellbore region of the fracture is given as follows(4)∇W˙w=∇(-kfμ∇p)=-kfμ(∂2p∂x′2+∂2p∂y′2)=0.The continuity equation for the radial flow in the near-wellbore region of the fracture can be written as follows:(5)r∂νr∂r+νr+r∂νy′∂y′=0,(6)∂rνr∂r+∂rνy′∂y′=0,(7)∇′′W˙w=[∂∂r∂∂y′]·[-rkfμ∂p∂r-rkfμ∂p∂y′]=0.Similar to (7), the continuity equation for the porous flow in the near-wellbore region of the reservoir matrix can be derived as
(8)∇′′W˙w=[∂∂r∂∂y]·[-rkmμ∂p∂r-rkmμ∂p∂y]=0.The symbol∇ denotes the vector gradient operator in the total Cartesian coordinate system (x,y), while ∇′ denotes that in the local Cartesian coordinate system (x′,y′). To account for the feature of radial flow, we use another local Cartesian coordinate (r,y′) and ∇′′, respectively, and replace (x′,y′) and ∇′.


### 2.3. Weak Form

Based on the strong-form equations of the four flow regimes, we can derive the weak form equations and then couple them to get the final two kinds of governing equations. There are four kinds of weak-form equations and the derivation processes of the first two are similar to that in [20].The weak form of the continuity equation of fluid flow in the away-from-wellbore region of the porous medium is given by(9)kmμ∫Ω∇ζw∇pdΩ-∫Γdζw⟦W˙w⟧nΓddΓd=0.
The derivation is as follows:
(10)∫Ωζw∇W˙wdΩ=0,(11)∫Ω[∇(ζwW˙w)-∇ζw·W˙w]dΩ=0,(12)∫Ω∇(ζwW˙w)dΩ=∫ΓζwW˙wnΓdΓ=∫Γd+ζwW˙wnΓd+dΓ+∫Γd-ζwW˙wnΓd-dΓ=∫Γd+ζwW˙w(-nΓd-)dΓ+∫Γd-ζwW˙wnΓd-dΓ=-∫Γdζw[W˙w+-W˙w-]nΓddΓd=-∫Γdζw⟦W˙w⟧nΓddΓd.Divergence theorem is applied in the derivation of (12). ζw is the test function for the fluid pressure and is equal to zero on the exterior boundary. Substitute (4) and (12) into (11), and then we get (9).The weak form of the continuity equation of fluid flow in the away-from-wellbore region of the fracture is given by(13)kfμ∫Γdw(Γd)∂ζw∂x′∂p∂x′dΓd+∫Γdζw⟦W˙w⟧nΓddΓd=0.
The derivation is as follows:
(14)∫Ω′ζw∇′W˙wdΩ′=0,∫Γdζw⟦W˙w⟧nΓddΓd+kfμ∫Ω′∇′ζw∇′pdΩ′=0,∫Ω′∇′ζw∇′pdΩ′=∫Γd∫-hh(∂ζw∂x′∂p∂x′+∂ζw∂y′∂p∂y′)dy′dx′=∫Γd∫-hh∂ζw∂x′∂p∂x′dtdr=∫Γdw(Γd)∂ζw∂x′∂p∂x′dΓd.Equation (13) is derived by the application of divergence theorem and the assumption (e). Couple the main reservoir flow and the fracture linear flow by substituting (13) into (9), and the first kind of the final governing equation (15) can be derived:
(15)km∫Ω∇ζw∇pdΩ+kf∫Γdw(Γd)∂ζw∂x′∂p∂x′dΓd=0.Similar to the derivation of (13), the weak form equation of fluid flow in the near-wellbore region of the fracture can be derived by introducing the feature of fracture radial flow [22] as
(16)kfμ∫Γdw(Γd)∂ζw∂rr∂p∂rdΓd+∫Γdζw⟦W˙w⟧nΓddΓ=0.In the near-wellbore region of the reservoir matrix surrounding the fracture radial flow region, the weak form equation of fluid flow can be derived by reference to the derivation of (16) as
(17)kmμ∫Ω∇′′ζw·r·∇′′pdΩ-∫Γdζw⟦W˙w⟧nΓddΓ=0.Couple the reservoir radial flow and the fracture radial flow by substituting (17) into (16), and the second kind of the final governing equation (18) is
(18)km∫Ω∇′′ζw·r·∇′′pdΩ+kf∫Γdw(Γd)∂ζw∂rr∂p∂rdΓd=0.Here we refer to the element with the wellbore located on its edge as “well element.” Equation (18) needs to be further developed by taking into account the internal boundary condition based on the method of Lagrangian multiplier. The details are specified in Section 4.


## 2.1. Physical Model and Assumptions

### 2.1.1. Physical Model

Consider a two-dimensional rectangular domainΩ bounded by the outer boundary Γ (Figure 1). A fractured horizontal well is located in the domain. The whole fracture reflects the weak discontinuity Γd on which the point O represents the wellbore. Oil in reservoir initially flows within the porous matrix and then from the matrix into the two-wing fracture and finally from the fracture into the wellbore. Because the width of the fracture is far less than its length, the fracture seems more like a discontinuous line in the view of the whole domain. But on a local scale the fracture is actually two-dimensional and is enveloped by two boundary lines (Γd+ and Γd-). The fluid is single-phase. Every flow regime is the isothermal event and obeys Darcy’s law. The internal boundary Γw relates to the wellbore.Figure 1
Fluid flow in hydraulically multi-fractured reservoir.n
Γ
c is the unit vector that points from the wellbore to the fracture tip. nΓc is the unit normal vector to the discontinuity pointing to Ω+(nΓd=nΓd-=-nΓd+). nΓ is the unit outward normal vector to the external boundary.Fluid flow in fractured horizontal well involves four main kinds of flow mechanisms, including fluid porous flow in the away-from-wellbore region of reservoir matrix (main reservoir flow), radial flow in the near-wellbore region of reservoir matrix (reservoir radial flow), linear flow in the away-from-wellbore region of fracture (fracture linear flow), and radial flow in the near-wellbore region of fracture (fracture radial flow). Different flow regimes are coupled by considering mass transfer between fracture and matrix (Figures2-3).Figure 2
The local Cartesian coordinate systems in the fracture.Figure 3
Flow regimes in the hydraulically fractured reservoir.The essential boundary conditions are imposed on the external and the internal boundary of the model as(1)p=peon⁡Γ,p=pwfon⁡Γw.Mass transfer is considered to couple the fluid flow in the fracture and the porous matrix. Regarding this, the complementary boundary condition is written as(2)⟦W˙w⟧nΓd=qwdon⁡Γd.⟦
*
⟧
=
*
+
-
*
- represents the difference between the corresponding values at the two fracture faces.


### 2.1.2. The Main Assumptions of the Model Include the Following

(a) Fluid flow in the physical model is in steady state; (b) fracture radial flow exists within a circle in the vertical plane with the center at the wellbore and diameter equal to reservoir thickness; (c) reservoir radial flow region exists in the near-wellbore region surrounding the near-wellbore fracture, and its range is similar to fracture radial flow; (d) the fluid pressure within the horizontal wellbore keeps constant; (e) the fluid pressure within the fracture keeps constant along the inward direction that is perpendicular to the fracture line, while the pressure gradient and the fluid flow are discontinuous.Based on the forth main assumption, this paper uses the same pressure field and the same mesh to approximate the fluid pressure in both the fracture and the matrix. This treatment is different from that in [18, 21].


## 2.1.1. Physical Model

Consider a two-dimensional rectangular domainΩ bounded by the outer boundary Γ (Figure 1). A fractured horizontal well is located in the domain. The whole fracture reflects the weak discontinuity Γd on which the point O represents the wellbore. Oil in reservoir initially flows within the porous matrix and then from the matrix into the two-wing fracture and finally from the fracture into the wellbore. Because the width of the fracture is far less than its length, the fracture seems more like a discontinuous line in the view of the whole domain. But on a local scale the fracture is actually two-dimensional and is enveloped by two boundary lines (Γd+ and Γd-). The fluid is single-phase. Every flow regime is the isothermal event and obeys Darcy’s law. The internal boundary Γw relates to the wellbore.Figure 1
Fluid flow in hydraulically multi-fractured reservoir.n
Γ
c is the unit vector that points from the wellbore to the fracture tip. nΓc is the unit normal vector to the discontinuity pointing to Ω+(nΓd=nΓd-=-nΓd+). nΓ is the unit outward normal vector to the external boundary.Fluid flow in fractured horizontal well involves four main kinds of flow mechanisms, including fluid porous flow in the away-from-wellbore region of reservoir matrix (main reservoir flow), radial flow in the near-wellbore region of reservoir matrix (reservoir radial flow), linear flow in the away-from-wellbore region of fracture (fracture linear flow), and radial flow in the near-wellbore region of fracture (fracture radial flow). Different flow regimes are coupled by considering mass transfer between fracture and matrix (Figures2-3).Figure 2
The local Cartesian coordinate systems in the fracture.Figure 3
Flow regimes in the hydraulically fractured reservoir.The essential boundary conditions are imposed on the external and the internal boundary of the model as(1)p=peon⁡Γ,p=pwfon⁡Γw.Mass transfer is considered to couple the fluid flow in the fracture and the porous matrix. Regarding this, the complementary boundary condition is written as(2)⟦W˙w⟧nΓd=qwdon⁡Γd.⟦
*
⟧
=
*
+
-
*
- represents the difference between the corresponding values at the two fracture faces.


## 2.1.2. The Main Assumptions of the Model Include the Following

(a) Fluid flow in the physical model is in steady state; (b) fracture radial flow exists within a circle in the vertical plane with the center at the wellbore and diameter equal to reservoir thickness; (c) reservoir radial flow region exists in the near-wellbore region surrounding the near-wellbore fracture, and its range is similar to fracture radial flow; (d) the fluid pressure within the horizontal wellbore keeps constant; (e) the fluid pressure within the fracture keeps constant along the inward direction that is perpendicular to the fracture line, while the pressure gradient and the fluid flow are discontinuous.Based on the forth main assumption, this paper uses the same pressure field and the same mesh to approximate the fluid pressure in both the fracture and the matrix. This treatment is different from that in [18, 21].


## 2.2. Strong Form

In order to derive the weak form of the governing equations, we must first constitute the respective continuity equations of the four flow regimes as the strong-form equations.The continuity equation for the porous flow in the away-from-wellbore region of the reservoir matrix is written as(3)∇W˙w=∇(-kmμ∇p)=-kmμ(∂2p∂x2+∂2p∂y2)=0.The continuity equation for the linear flow in the away-from-wellbore region of the fracture is given as follows(4)∇W˙w=∇(-kfμ∇p)=-kfμ(∂2p∂x′2+∂2p∂y′2)=0.The continuity equation for the radial flow in the near-wellbore region of the fracture can be written as follows:(5)r∂νr∂r+νr+r∂νy′∂y′=0,(6)∂rνr∂r+∂rνy′∂y′=0,(7)∇′′W˙w=[∂∂r∂∂y′]·[-rkfμ∂p∂r-rkfμ∂p∂y′]=0.Similar to (7), the continuity equation for the porous flow in the near-wellbore region of the reservoir matrix can be derived as
(8)∇′′W˙w=[∂∂r∂∂y]·[-rkmμ∂p∂r-rkmμ∂p∂y]=0.The symbol∇ denotes the vector gradient operator in the total Cartesian coordinate system (x,y), while ∇′ denotes that in the local Cartesian coordinate system (x′,y′). To account for the feature of radial flow, we use another local Cartesian coordinate (r,y′) and ∇′′, respectively, and replace (x′,y′) and ∇′.


## 2.3. Weak Form

Based on the strong-form equations of the four flow regimes, we can derive the weak form equations and then couple them to get the final two kinds of governing equations. There are four kinds of weak-form equations and the derivation processes of the first two are similar to that in [20].The weak form of the continuity equation of fluid flow in the away-from-wellbore region of the porous medium is given by(9)kmμ∫Ω∇ζw∇pdΩ-∫Γdζw⟦W˙w⟧nΓddΓd=0.
The derivation is as follows:
(10)∫Ωζw∇W˙wdΩ=0,(11)∫Ω[∇(ζwW˙w)-∇ζw·W˙w]dΩ=0,(12)∫Ω∇(ζwW˙w)dΩ=∫ΓζwW˙wnΓdΓ=∫Γd+ζwW˙wnΓd+dΓ+∫Γd-ζwW˙wnΓd-dΓ=∫Γd+ζwW˙w(-nΓd-)dΓ+∫Γd-ζwW˙wnΓd-dΓ=-∫Γdζw[W˙w+-W˙w-]nΓddΓd=-∫Γdζw⟦W˙w⟧nΓddΓd.Divergence theorem is applied in the derivation of (12). ζw is the test function for the fluid pressure and is equal to zero on the exterior boundary. Substitute (4) and (12) into (11), and then we get (9).The weak form of the continuity equation of fluid flow in the away-from-wellbore region of the fracture is given by(13)kfμ∫Γdw(Γd)∂ζw∂x′∂p∂x′dΓd+∫Γdζw⟦W˙w⟧nΓddΓd=0.
The derivation is as follows:
(14)∫Ω′ζw∇′W˙wdΩ′=0,∫Γdζw⟦W˙w⟧nΓddΓd+kfμ∫Ω′∇′ζw∇′pdΩ′=0,∫Ω′∇′ζw∇′pdΩ′=∫Γd∫-hh(∂ζw∂x′∂p∂x′+∂ζw∂y′∂p∂y′)dy′dx′=∫Γd∫-hh∂ζw∂x′∂p∂x′dtdr=∫Γdw(Γd)∂ζw∂x′∂p∂x′dΓd.Equation (13) is derived by the application of divergence theorem and the assumption (e). Couple the main reservoir flow and the fracture linear flow by substituting (13) into (9), and the first kind of the final governing equation (15) can be derived:
(15)km∫Ω∇ζw∇pdΩ+kf∫Γdw(Γd)∂ζw∂x′∂p∂x′dΓd=0.Similar to the derivation of (13), the weak form equation of fluid flow in the near-wellbore region of the fracture can be derived by introducing the feature of fracture radial flow [22] as
(16)kfμ∫Γdw(Γd)∂ζw∂rr∂p∂rdΓd+∫Γdζw⟦W˙w⟧nΓddΓ=0.In the near-wellbore region of the reservoir matrix surrounding the fracture radial flow region, the weak form equation of fluid flow can be derived by reference to the derivation of (16) as
(17)kmμ∫Ω∇′′ζw·r·∇′′pdΩ-∫Γdζw⟦W˙w⟧nΓddΓ=0.Couple the reservoir radial flow and the fracture radial flow by substituting (17) into (16), and the second kind of the final governing equation (18) is
(18)km∫Ω∇′′ζw·r·∇′′pdΩ+kf∫Γdw(Γd)∂ζw∂rr∂p∂rdΓd=0.Here we refer to the element with the wellbore located on its edge as “well element.” Equation (18) needs to be further developed by taking into account the internal boundary condition based on the method of Lagrangian multiplier. The details are specified in Section 4.


## 3. The Extended Finite Element Method

### 3.1. Approximation for Weak Discontinuity

The fluid flow jump across the fracture line reflects the feature of weak discontinuity which means that the pressure field is continuous, whereas its gradient is discontinuous. We can use the absolute value of the signed distance function to enrich the classical finite element approximation [17, 21].The fluid pressure is approximated as the linear combination of the standard and enriched shape function as(19)ph(x)=∑i∈NNi(x)Pi+∑j∈NsenrNj(x)[D(x)-D(xj)]aj,
where Ni(x) are the standard finite element shape functions; N is the set of all nodes in the mesh. Pi denotes the classical degree of freedom at node i; nodes in Nenr are such that their support bisected by the fracture; D(x) is the shape enrichment function and aj is the corresponding additional degree of freedom at node j. The first term on the right-hand side of (19) represents the classical finite element approximation, while the second one represents the weak discontinuous enrichment for the fracture.The shape enrichment function and its gradient are as follows:(20)D(x)=|ϕ(x)|,(21)∇D(x)·nΓd=∇ϕ(x)·H(ϕ(x)),(22)H(x)={1if(x-x*)·nΓd>0-1else.ϕ
(
x
) represents the level set function. Due to the general Heaviside function H(x) in (21), the gradient of the enrichment function is discontinuous across the fracture surfaces, while the enrichment function itself is continuous.


### 3.2. Quadrature Rules for Different Types of Elements

In the XFEM theory, split nodes are defined as nodes with their support intersecting the fracture. According to the number of the split nodes that an element contains, all the elements can be subsumed under three types, including the standard finite element, the split element, and the blending element. The quadrature rule differs by the element type.To cope with the standard finite element, we directly use the standard second-order Gauss quadrature rule, and there are 4 Gauss points per element. With regard to the split element, we first divide the element into triangular subdomains, then apply Hammer quadrature rule to every subdomain, and finally add together the integration of all the Gauss points in the element. Each split element generally contains 104 Gauss points.In order to avoid the decrease in the convergence rate of the blending element, we adopt the method of constructing the shifted enrichment function [23] and increasing the order of Gauss integration. Each spending element generally contains 36 Gauss points.


## 3.1. Approximation for Weak Discontinuity

The fluid flow jump across the fracture line reflects the feature of weak discontinuity which means that the pressure field is continuous, whereas its gradient is discontinuous. We can use the absolute value of the signed distance function to enrich the classical finite element approximation [17, 21].The fluid pressure is approximated as the linear combination of the standard and enriched shape function as(19)ph(x)=∑i∈NNi(x)Pi+∑j∈NsenrNj(x)[D(x)-D(xj)]aj,
where Ni(x) are the standard finite element shape functions; N is the set of all nodes in the mesh. Pi denotes the classical degree of freedom at node i; nodes in Nenr are such that their support bisected by the fracture; D(x) is the shape enrichment function and aj is the corresponding additional degree of freedom at node j. The first term on the right-hand side of (19) represents the classical finite element approximation, while the second one represents the weak discontinuous enrichment for the fracture.The shape enrichment function and its gradient are as follows:(20)D(x)=|ϕ(x)|,(21)∇D(x)·nΓd=∇ϕ(x)·H(ϕ(x)),(22)H(x)={1if(x-x*)·nΓd>0-1else.ϕ
(
x
) represents the level set function. Due to the general Heaviside function H(x) in (21), the gradient of the enrichment function is discontinuous across the fracture surfaces, while the enrichment function itself is continuous.


## 3.2. Quadrature Rules for Different Types of Elements

In the XFEM theory, split nodes are defined as nodes with their support intersecting the fracture. According to the number of the split nodes that an element contains, all the elements can be subsumed under three types, including the standard finite element, the split element, and the blending element. The quadrature rule differs by the element type.To cope with the standard finite element, we directly use the standard second-order Gauss quadrature rule, and there are 4 Gauss points per element. With regard to the split element, we first divide the element into triangular subdomains, then apply Hammer quadrature rule to every subdomain, and finally add together the integration of all the Gauss points in the element. Each split element generally contains 104 Gauss points.In order to avoid the decrease in the convergence rate of the blending element, we adopt the method of constructing the shifted enrichment function [23] and increasing the order of Gauss integration. Each spending element generally contains 36 Gauss points.


## 4. Introduction of the Internal Boundary Condition

In the context of the XFEM, the mesh no longer needs to match the geometry of the fracture. So if the whole fracture is located in the interior of the domain, we cannot load the internal boundary condition as simply as the outer boundary condition. Here, a solution based on the method of Lagrangian multiplier is put forward [24]:
(23)E(u)=p-pwf=0,(24)Π*=Π+∫ΓλTE(u)dΓ,(25)δΠ*=δΠ+∫ΓδλTE(u)dΓ+∫ΓλTδE(u)dΓ=0.E
(
u
) represents the essential boundary condition on the internal boundary. Π denotes the original functional before the internal boundary condition is brought in and Π* denotes the new functional. It can be confirmed that the physical meaning of the operator λ is the negative fluid flux on the internal boundary. By further taking into account the feature of the fracture radial flow, the new form of λ is as follows:
(26)λ=-rwkfμ∂p∂r.Substitute the new operatorλ into (25), and we obtain the following:
(27)δΠ*=δΠ+∫ΓδλTE(u)dΓ+∫ΓλTδE(u)dΓ=0=δΠ+∫Γδ(-rwkfμ∂p∂r)T(p-pwf)dΓ+∫Γ(-rwkfμ∂p∂r)Tδ(p-pwf)dΓ,(28)δΠ-rwkfμ∫Γw[δ(∂p∂r)Tp+(∂p∂r)Tδp]dΓ=-rwkfμ∫Γwδ(∂p∂r)TpwfdΓ.Go on to substitute the original functionalΠ into (28) and we can finally get the governing equation of coupled flow within the well element:
(29)km∫Ωr∇′′ζw∇′′pdΩ+kf∫Γdw(Γd)∂ζw∂rr∂p∂rdΓd+kf∫Γwrw[δ(∂p∂r)Tp+(∂p∂r)Tδp]dΩ=kf∫Γwδ(∂p∂r)TrwpwfdΩ.In the specific implementation process, in order to reduce the difficult of numerical simulation, the wellbore is located on the element edge.The total productivity of multifractured horizontal well can be predicted by adding together the flow of the different fractures:(30)Q=∑i2nfπrwkfwfμB∂p∂r|(r=rw)i.Skin factorS can be easily taken into account by substituting the effective wellbore radius rwe for the actual wellbore radius rw. However, it may be unreasonable to assume a skin factor existing within the fracture, because the main aim of hydraulic fracturing technology is to eliminate pollution surrounding the wellbore:
(31)rwe=rwe-S.


## 5. Numerical Examples

In this section, Case1 and Case 2 are utilized to determine the accuracy and reliability of the extended finite element model, while Case 3 is utilized to show the capability of solving problems for the complicated asymmetrical underground condition. The basic data of the three cases are from [2, 25]. The skin factor is set equal to zero as we think it is unreasonable to assume a skin factor existing within the fracture.Case 1.
A horizontal well with three vertical hydraulic fractures (Shuiping Well 1) is producing oil on the condition of constant wellbore flowing pressure. The actual flow rate of the well is 6.52 m3/d with a pressure drawdown of 10.14 MPa after 12 months of production. The reservoir properties and well data for case 1 are summarized in Table 1.Table 1
Reported reservoir and well data (Case1).
Name of parameters (unit)
Value
Reservoir thickness (m)
14
Length in thex direction (m)
600
Length in they direction (m)
400
Matrix permeability (mD)
1.3
Fracture permeability (D)
30
Oil viscosity (mPa·s)
3.5
Oil volume factor
1.13
Wellbore radius (m)
0.068
Skin factor
0
Half fracture lengths
Fracture 1 (m)
127.5
Fracture 2 (m)
105.0
Fracture 3 (m)
107.5
Fracture width (m)
0.005
Average distance between fractures (m)
95.4
Pressure on the external boundary (MPa)
16.41
Bottomhole flowing pressure (MPa)
6.2742The details of the different meshes and the corresponding solutions are shown in Table2. With the mesh from coarse to fine (from left to right in Table 2), the error of the calculated bottomhole flowing pressure tends to be zero. As the average side length of well element decreases, the production rate converges reasonably well. The convergence result seems questionable because it is only about 77.1% of the reported production rate, but we tend to regard the deviation as the normal difference caused by the different flow states. As a result of the very low matrix permeability and porosity, it may take a long time for flow in a reservoir to change from unsteady state to steady state. The flow in the reservoir in this case after 12 months of production may be still in unsteady state, so the actual production is normally higher than the calculated result in our steady-state model.Table 2
Comparison between the numerical results for different meshes.
Key parameters and results
Plan A
Plan B
Plan C
Plan D
Mesh strategy
Structured, uniform
Structured, uniform
Structured, uniform
Unstructured
Number of elements
10000
40000
90000
16107
Average side length of well element (m)
5.00
2.50
1.67
0.10
Calculated bottomhole flowing pressure (MPa)
2.7567
5.1523
5.7034
6.2715
Error of bottomhole flowing pressure (%)
-
56.06
-
17.88
-
9.10
-
0.04
Calculated production rate (m3/d)
6.52
5.52
5.28
5.03The pressure distribution in the steady state (under Plan D) is shown in Figure4. The discontinuity of the x-component of the pressure gradient across the fracture lines (under Plan D) is illustrated in Figure 5 and it should be noted that the range of the magnitude of the pressure gradient reflected by the Colorbar is only part of the real complete range which we have narrowed down to make the gradient jump more obvious in the figure.Figure 4
Matrix pore pressure distribution.Figure 5
The distribution of the x-component of pressure gradient.Case 2.
A horizontal well with four vertical hydraulic fractures (Maoping Well 1) is producing oil on the condition of constant wellbore flowing pressure. The actual flow rate of the well is 20.35 m3/d when the pressure drawdown is 9.1 MPa. The reservoir properties and well data for case 2 are summarized in Table 3.Table 3
Reported reservoir and well data (Case2).
Name of parameters (unit)
Value
Reservoir thickness (m)
11.9
Length in thex direction (m)
700
Length in they direction (m)
260
Matrix permeability (mD)
7.5
Fracture permeability (D)
30
Oil viscosity (mPa·s)
4.8
Oil volume factor
1.084
Wellbore radius (m)
0.058
Skin factor
0
Half fracture length (m)
75
Fracture width (m)
0.005
Average distance between fractures (m)
110.95
Pressure on the external boundary (MPa)
17.94
Bottomhole flowing pressure (MPa)
8.82The calculated production rate of Maoping well 1 based on unstructured mesh is 20.90 m3/d at the same pressure drawdown, and the relative error is 2.7%. In Case 2, the convergence result and the real production rate are nearly equal, which did not happen in Case 1. For that problem, we infer that the relatively high permeability and porosity accelerate the shift of the flow from unsteady state to steady state. So the flow of the reservoir in Case 2 may have been in steady state after 12 months of production.In this coupled model, the accuracy is closely related to the size of the elements surrounding wellbore. Due to the relatively large size of the well element, the numerical results under Plan A are completely wrong. With the mesh from coarse to fine (from left to right in Table4), the error of the calculated bottomhole flowing pressure tends to be zero and the production rate converges well.Table 4
Comparison between the numerical results for different meshes.
Key parameters and results
Plan A
Plan B
Plan C
Plan D
Mesh strategy
Structured, uniform
Structured, uniform
Structured, locally refined
Unstructured
Number of elements
90000
160000
151096
19404
Average side length of well element (m)
1.60
1.20
0.45
0.10
Calculated bottomhole flowing pressure (MPa)
-
7.7470
3.0947
7.9298
8.7675
Error of bottomhole flowing pressure (%)
-
187.83
-
64.91
-
10.09
-
0.60
Calculated production rate (m3/d)
-
46.93
32.75
22.49
20.90The pressure distribution (under Plan D) of Case2 is shown in Figure 6. The discontinuity of thex-component of the pressure gradient across the fracture lines (under Plan D) is illustrated in Figure 7.Figure 6
Matrix pore pressure distribution.Figure 7
The distribution of thex-component of pressure gradient.Case 3.
An asymmetrical model is established to analyze the flow regime and the pressure distribution in the hydraulically fractured reservoir. The computation parameters are basically the same as those reported by Li et al. [25] and Guo et al. [2]. One difference is that the formation is 440 m × 400 m. In addition, the location and the lengths of the two fracture wings are adjusted to make the model asymmetrical. The specific parameters of the fractures are summarized in Table 5. The mesh of the extended finite element model is shown in Figure 8. The result reflects an asymmetrical distribution of fluid pressure (see Figure 9), so the extended finite element model for multiscale fluid flow can no longer be based on the assumption of symmetry.Table 5
Reservoir and well data (Case3).
Fracture parameter
Fracture 1
Fracture 2
Location of wellbore
(−100 m, 52 m)
(50 m, 0 m)
Length of the upper wing
80 m
80 m
Length of the lower wing
100 m
80 m
Conductivity
infinite
infiniteFigure 8
The mesh of the extended finite element model.Figure 9
The asymmetrical pressure distribution.


## 6. Concluding Remarks

Based on the above, the following conclusions can be drawn.(1)
An effective extended finite element model has been formulated for prediction of productivity of multifractured horizontal well, and it shows favorable prospect in future engineering application.(2)
Within the frame of the XFEM, the model couples four kinds of the flow regimes, including fluid flow in the away-from-wellbore region of porous matrix (main reservoir flow), radial flow in the near-wellbore region of porous matrix surrounding the fracture (reservoir radial flow), linear flow in the away-from-wellbore region of fracture (fracture linear flow), and radial flow in the near-wellbore region of fracture (fracture radial flow).(3)
Cases1 and 2 illustrate the accuracy of the extended finite element model. Case 3 shows that, in solving the multiscale flow problem in multifractured horizontal well, the XFEM allows the mesh to be relatively coarse with respect to the classical FEM and independent of the fracture.(4)
Furthermore the simulation results show that the extended finite element model constructed in this paper is more appropriate for the complicated asymmetrical physical condition.



---

*Source: 810493-2014-07-23.xml*