Extended Finite Element Method for Predicting Productivity of Multifractured Horizontal Wells

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.

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 Γ_𝑑 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.

Figure 1 

Fluid flow in hydraulically multi-fractured reservoir.

𝐧_Γ𝑐 is the unit vector that points from the wellbore to the fracture tip. 𝐧_Γ𝑐 is the unit normal vector to the discontinuity pointing to Ω⁺(𝐧_Γ𝑑=𝐧_Γ⁻_𝑑=−𝐧_Γ⁺_𝑑). 𝐧_Γ 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).

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 𝑝=𝑝_𝑒onΓ,𝑝=𝑝_𝑤𝑓onΓ_𝑤. (1)

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 𝑊̇_𝑤𝐧_Γ𝑑=𝑞_𝑤𝑑onΓ_𝑑. (2)

⟦∗⟧=∗⁺−∗⁻ 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 ∇𝑊̇_𝑤−𝑘=∇^𝑚𝜇𝑘∇𝑝=−^𝑚𝜇𝜕²𝑝𝜕𝑥²+𝜕²𝑝𝜕𝑦²=0. (3)

The continuity equation for the linear flow in the away-from-wellbore region of the fracture is given as follows ∇𝑊̇_𝑤−𝑘=∇^𝑓𝜇𝑘∇𝑝=−^𝑓𝜇𝜕²𝑝𝜕𝑥^′2+𝜕²𝑝𝜕𝑦^′2=0. (4)

The continuity equation for the radial flow in the near-wellbore region of the fracture can be written as follows: 𝑟𝜕𝜈_𝑟𝜕𝑟+𝜈_𝑟+𝑟𝜕𝜈_𝑦′𝜕𝑦^′=0,𝜕𝑟𝜈_𝑟+𝜕𝑟𝜕𝑟𝜈_𝑦′𝜕𝑦^′=0,∇^′′𝑊̇_𝑤=𝜕𝜕𝜕𝑟𝜕𝑦^′⋅𝑘−𝑟^𝑓𝜇𝜕𝑝𝑘𝜕𝑟−𝑟^𝑓𝜇𝜕𝑝𝜕𝑦^′=0. (5)(6)(7)

Similar to (7), the continuity equation for the porous flow in the near-wellbore region of the reservoir matrix can be derived as ∇^′′𝑊̇_𝑤=𝜕𝜕𝜕𝑟⋅𝑘𝜕𝑦−𝑟^𝑚𝜇𝜕𝑝𝑘𝜕𝑟−𝑟^𝑚𝜇𝜕𝑝𝜕𝑦=0. (8)

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 ∇^′.

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 𝑘_𝑚𝜇_Ω∇𝜁_𝑤∇𝑝𝑑Ω−_Γ𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ𝑑𝑑Γ_𝑑=0. (9) The derivation is as follows: _Ω𝜁_𝑤∇𝑊̇_𝑤𝑑Ω=0,_Ω∇𝜁_𝑤𝑊̇_𝑤−∇𝜁_𝑤⋅𝑊̇_𝑤𝑑Ω=0,_Ω∇𝜁_𝑤𝑊̇_𝑤𝑑Ω=_Γ𝜁_𝑤𝑊̇_𝑤𝐧_Γ=𝑑Γ_Γ+𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ⁺_𝑑𝑑Γ+_Γ−𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ⁻_𝑑=𝑑Γ_Γ+𝑑𝜁_𝑤𝑊̇_𝑤−𝐧_Γ⁻_𝑑+𝑑Γ_Γ−𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ⁻_𝑑𝑑Γ=−_Γ𝑑𝜁_𝑤𝑊̇⁺_𝑤−𝑊̇⁻_𝑤𝐧_Γ𝑑𝑑Γ_𝑑=−_Γ𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ𝑑𝑑Γ_𝑑. (10)(11)(12)

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 𝑘_𝑓𝜇_Γ𝑑𝑤_(Γ𝑑)𝜕𝜁_𝑤𝜕𝑥^′𝜕𝑝𝜕𝑥^′𝑑Γ_𝑑+_Γ𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ𝑑𝑑Γ_𝑑=0. (13) The derivation is as follows: _Ω′𝜁_𝑤∇^′𝑊̇_𝑤𝑑Ω^′=0,_Γ𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ𝑑𝑑Γ_𝑑+𝑘_𝑓𝜇_Ω^′∇^′𝜁_𝑤∇^′𝑝𝑑Ω^′=0,_Ω′∇^′𝜁_𝑤∇^′𝑝𝑑Ω^′=_Γ𝑑^ℎ_−ℎ𝜕𝜁_𝑤𝜕𝑥^′𝜕𝑝𝜕𝑥^′+𝜕𝜁_𝑤𝜕𝑦^′𝜕𝑝𝜕𝑦^′𝑑𝑦^′𝑑𝑥^′=_Γ𝑑^ℎ_−ℎ𝜕𝜁_𝑤𝜕𝑥^′𝜕𝑝𝜕𝑥^′=𝑑𝑡𝑑𝑟_Γ𝑑𝑤_(Γ𝑑)𝜕𝜁_𝑤𝜕𝑥^′𝜕𝑝𝜕𝑥^′𝑑Γ_𝑑. (14)

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: 𝑘_𝑚Ω∇𝜁_𝑤∇𝑝𝑑Ω+𝑘_𝑓Γ𝑑𝑤_(Γ𝑑)𝜕𝜁_𝑤𝜕𝑥^′𝜕𝑝𝜕𝑥^′𝑑Γ_𝑑=0. (15)

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 𝑘_𝑓𝜇_Γ𝑑𝑤_(Γ𝑑)𝜕𝜁_𝑤𝑟𝜕𝑟𝜕𝑝𝜕𝑟𝑑Γ_𝑑+_Γ𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ𝑑𝑑Γ=0. (16)

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 𝑘_𝑚𝜇_Ω∇^′′𝜁_𝑤⋅𝑟⋅∇^′′𝑝𝑑Ω−_Γ𝑑𝜁_𝑤𝑊̇_𝑤𝐧_Γ𝑑𝑑Γ=0. (17)

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 𝑘_𝑚Ω∇^′′𝜁_𝑤⋅𝑟⋅∇^′′𝑝𝑑Ω+𝑘_𝑓Γ𝑑𝑤_(Γ𝑑)𝜕𝜁_𝑤𝑟𝜕𝑟𝜕𝑝𝜕𝑟𝑑Γ_𝑑=0. (18)

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 𝑝^ℎ(𝐱)=_𝑖∈𝑁𝑁_𝑖(𝐱)𝑃_𝑖+_{𝑗∈𝑁enr𝑠}𝑁_𝑗(𝐱)𝐷(𝐱)−𝐷𝐱_𝑗𝑎_𝑗, (19) where 𝑁_𝑖(𝐱) 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; 𝐷(𝐱) 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: 𝐷(𝐱)=𝜙(𝐱),∇𝐷(𝐱)⋅𝐧_Γ𝑑=∇𝜙(𝐱)⋅𝐻𝜙(𝐱),𝐻(𝐱)=1if𝐱−𝐱^∗⋅𝐧_Γ𝑑>0−1else. (20)(21)(22)

𝜙(𝐱) represents the level set function. Due to the general Heaviside function 𝐻(𝐱) 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]: 𝐄(𝐮)=𝑝−𝑝_𝑤𝑓=0,Π^∗=Π+_Γ𝜆^T𝐄(𝐮)𝑑Γ,𝛿Π^∗=𝛿Π+_Γ𝛿𝜆^T𝐄(𝐮)𝑑Γ+_Γ𝜆^T𝛿𝐄(𝐮)𝑑Γ=0. (23)(24)(25)

𝐄(𝐮) 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)

Substitute the new operator 𝜆 into (25), and we obtain the following: 𝛿Π^∗=𝛿Π+_Γ𝛿𝝀^T𝐄(𝐮)𝑑Γ+_Γ𝝀^T𝛿𝐄(𝐮)𝑑Γ=0=𝛿Π+_Γ𝛿−𝑟_𝑤𝑘_𝑓𝜇𝜕𝑝𝜕𝑟^T𝑝−𝑝_𝑤𝑓+𝑑Γ_Γ−𝑟_𝑤𝑘_𝑓𝜇𝜕𝑝𝜕𝑟^T𝛿𝑝−𝑝_𝑤𝑓𝑑Γ,𝑟𝛿Π−^𝑤𝑘_𝑓𝜇_Γ𝑤𝛿𝜕𝑝𝜕𝑟^T𝑝+𝜕𝑝𝜕𝑟^T𝑟𝛿𝑝𝑑Γ=−^𝑤𝑘_𝑓𝜇_Γ𝑤𝛿𝜕𝑝𝜕𝑟^T𝑝_𝑤𝑓𝑑Γ. (27)(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: 𝑘_𝑚Ω𝑟∇^′′𝜁_𝑤∇^′′𝑝𝑑Ω+𝑘_𝑓Γ𝑑𝑤_(Γ𝑑)𝜕𝜁_𝑤𝑟𝜕𝑟𝜕𝑝𝜕𝑟𝑑Γ_𝑑+𝑘_𝑓Γ𝑤𝑟_𝑤𝛿𝜕𝑝𝜕𝑟^T𝑝+𝜕𝑝𝜕𝑟^T𝛿𝑝𝑑Ω=𝑘_𝑓Γ𝑤𝛿𝜕𝑝𝜕𝑟^T𝑟_𝑤𝑝_𝑤𝑓𝑑Ω. (29)

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: 𝑄=^2𝑛𝑓_𝑖𝜋𝑟_𝑤𝑘_𝑓𝑤_𝑓𝜇𝐵𝜕𝑝𝜕𝑟_{(𝑟=𝑟𝑤)𝑖}. (30)

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)

5. 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³/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 (Case 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.

Table 2 

Comparison between the numerical results for different meshes.

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.

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 m³/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 (Case 2).

The calculated production rate of Maoping well 1 based on unstructured mesh is 20.90 m³/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 Table 4), 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.

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.

Figure 6 

Matrix pore pressure distribution.

Figure 7 

The distribution of the 𝑥-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 (Case 3).

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

Nomenclature

𝜈𝑟:	Radial velocity in the near-wellbore region of fracture, m/s
𝜈𝑦′:	Fluid velocity in the direction normal to the fracture line, m/s
𝐵:	Oil volume factor, rb/stb
𝐻𝑓:	Reservoir thickness, m
𝑊̇𝑤:	Velocity vector of the pore fluid, m/s
𝑘𝑓:	Fracture permeability, m−2
𝑘𝑚:	Reservoir permeability, m−2
𝑝:	Reservoir pressure, Pa
𝜇:	Oil viscosity, Pa⋅s
𝜁𝑤:	Test function of fluid pressure, Pa
𝑤(Γ𝑑):	Fracture width along the direction parallel to the fracture line, m
ℎ:	Half of fracture width, m
𝑝𝑤𝑓:	Flowing bottomhole pressure, Pa
𝑄:	Production rate, m3/s
𝑟𝑤:	Radius of wellbore, m
𝑛𝑓:	Number of fractures.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The research is supported by the National Natural Sciences Foundation of China (51344005), National Major Special Science and Technology Project of China (2011ZX05014-006), and Program for Innovative Research Team of Southwest Petroleum University (2013XJT001).