Extended Finite Element Method for Predicting Productivity of Multifractured Horizontal Wells

Based on the theory of the extended finite elementmethod (XFEM), which was first proposed byMoë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.


Introduction
Prediction of productivity of multifractured horizontal wells has always been a research hotspot.Various methods including the analytical method [1,2], the hybrid numericalanalytical 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 numericalanalytical 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][8][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][14][15][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 Horizontal wellbore 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-fromwellbore region of the porous matrix, radial flow in the nearwellbore region of the matrix, linear flow in the away-fromwellbore region of fracture, and radial flow in the nearwellbore 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.

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.

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 Γ  on which the point  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 (Γ +  and Γ −  ).The fluid is single-phase.Every flow regime is the isothermal event and obeys Darcy's law.The internal boundary Γ  relates to the wellbore.n Γ  is the unit vector that points from the wellbore to the fracture tip.n Γ  is the unit normal vector to the discontinuity pointing to . 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 (Figures 2-3).
The essential boundary conditions are imposed on the external and the internal boundary of the model as 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 ⟦ * ⟧ = * + − * − represents the difference between the corresponding values at the two fracture faces.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].

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 strongform equations.
The continuity equation for the porous flow in the awayfrom-wellbore region of the reservoir matrix is written as The continuity equation for the linear flow in the awayfrom-wellbore region of the fracture is given as follows The continuity equation for the radial flow in the nearwellbore region of the fracture can be written as follows: Similar to (7), the continuity equation for the porous flow in the near-wellbore region of the reservoir matrix can be derived as The symbol ∇ denotes the vector gradient operator in the total Cartesian coordinate system (, ), while ∇  denotes that in the local Cartesian coordinate system (  ,   ).To account for the feature of radial flow, we use another local Cartesian coordinate (,   ) and ∇  , respectively, and replace (  ,   ) and ∇  .

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 The derivation is as follows: Divergence theorem is applied in the derivation of (12).  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 The derivation is as follows: 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: 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 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 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 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.

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 where   (x) are the standard finite element shape functions;  is the set of all nodes in the mesh.  denotes the classical degree of freedom at node ; nodes in  enr are such that their support bisected by the fracture; (x) is the shape enrichment function and   is the corresponding additional degree of freedom at node .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: (x) represents the level set function.Due to the general Heaviside function (x) in ( 21), the gradient of the enrichment function is discontinuous across the fracture surfaces, while the enrichment function itself is continuous.

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.

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]: 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: Substitute the new operator  into (25), and we obtain the following:

𝛿( 𝜕𝑝 𝜕𝑟 )
T   Γ. (28) Go on to substitute the original functional Π into (28) and we can finally get the governing equation of coupled flow within the well element: 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: Skin factor  can be easily taken into account by substituting the effective wellbore radius   for the actual wellbore radius   .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)

Numerical Examples
In this section, Case 1 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 m 3 /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.
The details of the different meshes and the corresponding solutions are shown in Table 2.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.The pressure distribution in the steady state (under Plan D) is shown in Figure 4.The discontinuity of the 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 .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 m 3 /d when the pressure drawdown is 9.1 MPa.The reservoir properties and well data for case 2 are summarized in Table 3.
The calculated production rate of Maoping well 1 based on unstructured mesh is 20.90 m 3 /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 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 Table 4), the error of the calculated bottomhole flowing pressure tends to be zero and the production rate converges well.
The pressure distribution (under Plan D) of Case 2 is shown in Figure 6.The discontinuity of the x-component of the pressure gradient across the fracture lines (under Plan D) is illustrated in Figure 7.
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.

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.(3) Cases 1 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

Figure 2 :
Figure 2: The local Cartesian coordinate systems in the fracture.

Figure 3 :
Figure 3: Flow regimes in the hydraulically fractured reservoir.
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 nearwellbore 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.

Table 1 :
Reported reservoir and well data (Case 1).Name of parameters (unit) ValueReservoir thickness (m) 14Length in the  direction (m) 600Length in the  direction (m

Table 2 :
Comparison between the numerical results for different meshes.

Table 3 :
Reported reservoir and well data (Case 2).

Table 4 :
Comparison between the numerical results for different meshes.