Next Article in Journal
Battery Energy Storage System for Aggregated Inertia-Droop Control and a Novel Frequency Dependent State-of-Charge Recovery
Next Article in Special Issue
Reservoir Properties of Low-Permeable Carbonate Rocks: Experimental Features
Previous Article in Journal
Effect of Existing Façade’s Construction and Orientation on the Performance of Low-E-Based Retrofit Double Glazing in Tropical Climate
Previous Article in Special Issue
An Analytical Solution for Transient Productivity Prediction of Multi-Fractured Horizontal Wells in Tight Gas Reservoirs Considering Nonlinear Porous Flow Mechanisms
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Productivity-Index Behavior for a Horizontal Well Intercepted by Multiple Finite-Conductivity Fractures Considering Nonlinear Flow Mechanisms under Steady-State Condition

1
School of Computer & Information Technology, Northeast Petroleum University, Daqing 163318, China
2
Department of Well Logging & Remote Sensing Technology, Research Institute of Petroleum Exploration and Development, PetroChina, Beijing 100083, China
*
Author to whom correspondence should be addressed.
Energies 2020, 13(8), 2015; https://doi.org/10.3390/en13082015
Submission received: 23 March 2020 / Revised: 10 April 2020 / Accepted: 13 April 2020 / Published: 17 April 2020
(This article belongs to the Special Issue Development of Unconventional Reservoirs 2020)

Abstract

:
In this paper, a mathematical model is proposed to investigate the effect of nonlinear flow mechanisms on productivity-index (PI) behavior in hydraulically fractured reservoirs during steady-state condition. This approach focuses on the fact that PI approaches a constant value at a certain time, indicating the beginning of steady state. In this model, the reservoirs are considered as an elliptical-shaped drainage with constant-pressure boundary, which is depleted by a multiple-fractured horizontal well (MFHW), and various nonlinear flow mechanisms, such as the non-Darcy flow effect and pressure-dependency effect, control flow patterns in the hydraulic fractures. Then, an exact algorithm of solving the resulting nonlinear equations is developed to obtain the PI of MFHW using a semi-analytical approach. Next, type curves are generated to investigate the influences of flow mechanisms and fracture properties on PI. The most interesting points in this study are the following: (1) PI is determined by the properties of MHFW (i.e., dimensions and configuration), the reservoir geometry, and flow mechanism; (2) PI is deteriorated by non-Darcy flow caused by inertial forces; and (3) PI is reduced under the influence of pressure sensitivity caused by the degradation of dynamic conductivity. Generally, this study provides a significant insight into understanding the factors affecting the productivity of a MFHW with nonlinear flow mechanisms.

1. Introduction

With the success of Bakken tight oil resource in the United States, tight oil and gas reservoirs have become a significant source of hydrocarbon supply globally [1]. Subsequently, the technology of horizontal drilling combined with hydraulically fracturing has effectively enabled commercial production from the reservoir with extremely low permeability [2,3]. In practice, the characteristics for a stimulated well are identified by flow regimes; for example, the approach of pressure-derivative is used to identify flow-regimes for pressure transient analysis. For a multi-fractured horizontal well, the sequence of flow-regimes may be more complex (Figure 1) [4]. At some point after compound linear flow caused by the fracture interference, a compound elliptical flow regime (regime 5 in Figure 1a) would appear as the transient pressure waves move past the ends of the fractures [5,6]. Put another way, the compound elliptical response is the pseudo-steady state in an infinite drainage area. It is noted that the isopotential line approximates to an elongated elliptical shape for a long horizontal wellbore, as shown in Figure 1b.
In physics, the performance of a fractured reservoir is determined by numerous factors such as formation petrophysical properties, fluid properties, and reservoir and wellbore configuration [7]. Numerous methods are provided to capture the physics of fluid flowing in the fractured reservoirs. Originating from the pioneer work suggested by Warren and Root [8], the fractures are represented by a 3D uniformly spaced fracture network, and the matrix blocks are contributing sources. It is a fictitious homogeneous porous medium. As a result, multiple-continuum models are used to compute the transient response of uniform naturally (continuously) fractured reservoirs. In analytical modeling, these models consist of a set of one-dimensional linear flow regions, and the main difference is the way in which these regions are coupled [9,10,11]. Therefore, it is difficult for the multiple-continuum model to capture the intrinsic behavior of complexity of nonuniform naturally fractured reservoirs, such as hydraulically fractured reservoirs. To capture the strong heterogeneity between fracture and matrix, various improved techniques are presented to subdivise the grids, such as the embedded discrete fracture model (EDFM), Galerkin finite-element method (FEM), and discrete-element method (DEM) [12,13,14]. However, it is still a huge challenge to obtain accurate transmissibility by computing the irregular connections between fracture segments. Considering the high computational burden and the great difficulties in gridding, the mesh-free semi-analytical method is necessarily developed as an alternative to the numerical simulation method. On the basis of the works suggested by Cinco-Ley et al. [15], various semi-analytical models, which have the flexibility of accounting for the individual fracture properties, are proposed to investigate the influence factors on the performance of the fractured reservoirs [16,17,18,19]. It is proven that the semi-analytical methods have huge potential in analyzing transient performance response and evaluating productivity of wells.
From the viewpoint of petroleum engineering, the effects of these parameters on reservoir performance can be represented by the productivity index (PI). It is defined as the ratio of production rate to a certain pressure drop. PI is determined by the flow regime, with different behavior characteristics with time. It starts out with a high value at an early time when the transient response is dominant. PI declines with time to a small value, and approximates to a constant value, known as stabilized PI, at a late time period when the compound elliptical response is reached [20]. Currently, most studies have put the focus on the transient productivity index to identify the flow regimes [21,22,23]. However, they found that the designed fractured wells do not usually result in the expected performance, which is attributed to the existence of nonlinear flows. From the viewpoint of multi-physics, the flowing in the underground porous media is the result of fluid–thermal–solid coupling [24]. The reservoir production is generally regarded as the isothermal process; therefore, the thermal is not considered in the production simulation. Strictly speaking, the rock deformation may have a significant effect on fluid flow owing to fluid-particle interactions. According to the conclusions presented by Zhang and Tahmasebi [25,26], the effect of solid deformation on fluid flow, especially in hydraulic fractures, could be incorporated by regarding the porosity and permeability as a function of stress changes. In other words, the conductivity of hydraulic fracture changes because of a variety of mechanisms, including fracture closure under high effective stress and permeability loss as a result of proppant embedment and crushing. Besides, non-Darcy flow may occur in the fracture when the well is produced at a high flow rate [27]. Therefore, the nonlinear flow effects have to be taken into account to investigate the effect of nonlinear behaviors on the performance of a vertical well and fractured well in the infinite-acting reservoir. Limited efforts have been made to examine the effects of non-Darcy flow and pressure-sensitive conductivity on the performance of horizontal wells with multiple fractures. Wu et al. [28] used a numerical approach to investigate the non-Darcy effect in an unconventional gas reservoir. Lin and Zhu [29] developed a slab source method to evaluate the performance of horizontal wells with or without fractures using an analytical approach. Besides the non-Darcy effect, pressure sensitivity in the fracture may also exert an important influence on the transient response. A transient pressure analysis for a vertical well intercepting a dynamic-conductivity fracture has been conducted using analytical and numerical methods. Pedroso and Correa [30] developed a new model to investigate the effects of pressure-sensitive permeability on the flow behavior of a fractured well. Zhang et al. [31] presented the results of a study that analyzed the build-up of pressure in a vertically fractured well with a stress-sensitive conductivity in the fracture. Yao et al. [32] studies the effect of dynamic conductivity on the transient performance of multiple fractured horizontal well.
To the best of our knowledge, few papers are published to study the productivity index of a multiple-fractured horizontal well (MFHW) in a finite drainage area with an elliptical boundary under nonlinear flow conditions. The objective of this work is to quantify the effect of non-Darcy flow and pressure sensitivity on the inflow performance of a MFHW. First, we establish a hybrid model that describes Darcy flow in the reservoirs toward hydraulic fractures, as well as non-Darcy flow along the fracture with pressure-dependent conductivity under steady-state conditions. First, a mathematical model is established to describe the fluid flow depleted by a MFHW in a reservoir with elliptical-shaped boundary, with the flexibility of accounting for fracture properties (i.e., fracture dimensions and configuration). Next, nonlinear flow mechanisms in the fracture, including stress-sensitivity effect and non-Darcy flow condition, are taken into account, which results in the nonlinearity in partial differential equation. Then, the technique of dimension transformation is used to render the resulting nonlinear equations amenable to linear treatment, and an effective and efficient algorithm is developed to solve the model. Finally, a study for sensitivity analysis is introduced to investigate the effect of fracture properties and nonlinear flow mechanisms on the productivity-index behavior.

2. Model Assumption

2.1. Physical Model

In this physical model, a horizontal well is located in the center of a homogeneous reservoir with an elongated elliptical boundary. The whole drainage area consists of several different flow patterns in four contiguous flow regions, as sketched in Figure 2, including (1) the hydraulic fractures, (2) the inner stimulated rock volume (SRV) region impacted by hydraulic fracturing, (3) the outer linear region beyond the extent of the SRV, and (4) the outer elliptical region beyond the tips of the horizontal wellbore.
Different from the anisotropic characteristics in a naturally fractured reservoir, where a large number of nature fractures are uniformly and globally distributed [33,34,35], this paper assumes that the reservoir is stimulated by several hydraulic fractures, so the resulting hydraulically fractured reservoir is regarded as isotropic. In addition, there are some necessary assumptions, as follows:
  • The reservoir is horizontal and of uniform thickness, with impermeable lower and upper boundaries. The pressure on the outer boundary in the xy plane keeps constant.
  • Along the horizontal wellbore, multiple hydraulic fractures are evenly distributed. Hydraulic fractures are considered to be a finite conductivity.
  • The wellbore is produced under constant-rate condition.
  • Flow in the reservoir is assumed to be single-phase fluid (i.e., pure gas) with slight compressibility and constant viscosity.
  • Fluid flowing in the fracture is assumed to be nonlinear.
  • Fluids from the reservoir enter the wellbore only through hydraulic fractures, and the pressure loss in the wellbore is neglected.

2.2. Variables Definitions

To introduce the mathematical model conveniently and concisely, the dimensionless variables are defined, respectively, by
p D = 2 π k m h ( p i p ) q μ B ,   ξ D = ξ L r e f ,   q D = q q r e f ,   C f D = k f w f k m L r e f
When the flow in the fracture satisfies non-Darcy effect owing to high velocity, the corresponding dimensionless variables are defined as
( q D N D ) f = k f β ρ μ w f h q r e f ,   and   q c D = q c q r e f ,
When the fracture conductivity behaves pressure sensitivity, the corresponding dimensionless variables are defined as
γ f D = q μ B 2 π k m h γ f
where pm is the pressure in the reservoir, Pa; pf is the pressure in the fracture, Pa; pi is the pressure on the elliptical boundary, Pa; pex is the pressure on the width of inner SRV region, Pa; pey is the pressure on the length of inner SRV region, Pa; (xofm, yofm) is the location of m-th fracture tip, m; xe is the width of inner SRV region, m; ye is the length of inner SRV region, m; xR is the width of the whole drainage region, m3; h is the reservoir thickness, m; wf is the fracture width, m; Lf is the fracture length, m; rw is the radius of wellbore, m; Lref is the reference length, m; qfm is the function of inflow distribution along the m-th fracture, m2/s; qwm is the production rate, m3/s; qc is the influx rate, m3/s; μ is the viscosity, Pa s; B is the volume factor; γf is the permeability modulus, Pa−1; and β is the Beta factor, m−1. Here, the subscript D represents dimensionless variables.

3. Mathematical Model

3.1. Fluid Flow in the Inner SRV Region

The SRV region is regarded as a rectangular reservoir with constant-pressure boundary containing a set of hydraulic fractures, as shown in Figure 3. The diffusion equation governing fluid flow is given as the following dimensionless form:
2 p m D x D 2 + 2 p m D y D 2 + 2 π m = 1 N f 0 L f D m q f D m ( u D ) δ ( x D x o f D m u D ) δ ( y D y o f D m ) d u D = 0
The outer boundary conditions are satisfied by
{ p m D ( x D = x e D , y D ) = p m D ( x D = 0 , y D ) = p e D x p m D ( x D , y D = 0 ) = p m D ( x D , y D = y e D ) = p e D y
where there are Nf fractures in the SRV region, the dimensionless length of m-th fracture is LfDm, and the dimensionless coordinate of m-th fracture tip is (xofDm, yofDm). The pressure on the outer boundary parallel to the fracture is peDx, and the pressure on the outer boundary perpendicular to the fracture is peDy.
Using Fourier transformation and inverse transformation (Appendix A), the dimensionless pressure distribution caused by Nf fracture is given by
p m D ( x D , y D ) = p e D x I B ( x D , y D ) + p e D y I D ( x D , y D ) + 2 m = 1 N f 0 L f D m q f D m ( u D ) δ p u D ( x D , y D , u D ) d u D

3.2. Fluid Flow in the Outer Region

The fluids stored in the outer region are depleted and then flow into the inner SRV region through the interface (denoted by the yellow block in Figure 4). The streamlines satisfy the hyperbolic shape (orange solid line in Figure 4), which are the counterpart of elliptical-shaped isopotential line (black dashed line). Once the fluids enter the inner region, the streamlines are in a linear shape around the interface. Therefore, the coordinate systems are divided into two subsets: 1D elliptical coordinate in the outer region (blue region), and 1D linear coordinate (yellow region).
The outer boundary satisfies the constant pressure of peD.
In the outer elliptic region, as shown in Figure 4, there are two sets of coordinates: the Cartesian coordinates (x’,y’) and the elliptic coordinates (ξ,η). The relationship is given by
x = ( x e / 2 ) cosh ξ cos η ,   y = ( x e / 2 ) sinh ξ sin η
We used the elliptic coordinates to simulate the flowing process. The outer and inner boundaries of the elliptic region, represented by ξ in the elliptic coordinates, could be expressed using spatial variables in the Cartesian coordinates:
{ ξ e = ξ ( x = x R / 2 , y = 0 ) = ln [ ( x R + x R 2 x e 2 ) / x e ] ξ w = ξ ( x = 0 , y = 0 ) = ln 1 = 0
The corresponding pressure gradient is written as the following dimensionless expression,
( p m D ξ D ) ξ D = 0 = p i p e y ξ e ξ w = p e y D ln [ ( x R D + x R D 2 x e D 2 ) / x e D ]
Note that the flux inflow in Equation (9) is continuous on the interface between the outer elliptical region and inner SRV region. Meanwhile, the flowing is modeled using linear coordinates. As a result, the continuous condition is written in the form of the average integral, which is given by
1 x e D 0 π p m D ξ D | ξ D = 0 d η D e l l i p t i c   f l o w = 1 x e D 0 x e D p m D y D | y D = y e D d x D l i n e a r   f l o w
It could be expressed as the function of peDx and peDy (Appendix B), which is written as follows:
A 2 p e D x + B 2 p e y D = m = 1 N f q w D m C 2 , m
Likewise, the continuity condition of flux inflow between outer linear region and inner SRV region contributes to
A 1 p e D x + B 1 p e y D = m = 1 N f q w D m C 1 , m
We can obtain the dimensionless pressure on the interface by combing Equations (11) and (12),
{ p e D x = m = 1 N f q w D m ( B 2 C 1 , m B 1 C 2 , m A 1 B 2 A 2 B 1 ) p e D y = m = 1 N f q w D m ( A 1 C 2 , m A 2 C 1 , m A 1 B 2 A 2 B 1 )
where the relationship between influx distribution and production rate of m-th fracture is given as.
Substituting Equation (13) into Equation (6) could yield the pressure distribution caused by Nf fractures, which is the function with regard to influx distribution along the fracture.

3.3. Fluid Flow in the Fracture

For hydraulic fracture, the flow pattern is widely assumed to be 1D linear and incompressible [36]. As shown in Figure 5, the governing equation of fluid flow on the m-th fracture is described as the following dimensionless form:
2 p f D n x D m 2 2 π C f D m q f D m ( x D m ) + 2 π C f D m [ q w f D m δ ( x D m x w D m ) ] = 0
The corresponding boundary conditions are described as
p f D m x D m | x D m = 0 = p f D x D m | x D m = L f D m = 0
Integrating Equation (14) from 0 to xDm with respect to xDm based on Equation (15) twice, the result is given by
p f D m ( x D m ) = p f D m ( 0 ) + 2 π C f D m 0 x D m d ζ 0 ζ q f D m ( ξ ) d ξ 2 π C f D m q w D m H ( x D m , x w D m )
where δ ( x , x 0 ) = { ,   x = x 0 0 ,   x x 0 is the Dirac function, H ( x , x 0 ) = { 1 ,   x x 0 0 ,   x < x 0 is the Heaviside function, and G ( x , x 0 ) = { x x 0 ,   x x 0 0 ,   x < x 0 is the integral of Heaviside function.
As shown in Figure 5, the fluid flow is regarded as a 1D linear pattern in the region far from wellbore, and becomes a 1D radial pattern around the wellbore. It is different from the vertical fracture, where the vertical fracture is in lateral contact with the wellbore and only linear flow is considered. The additional pressure loss owing to the convergence of fluids into the horizontal wellbore should be considered. The effect of flow convergence is regarded as a skin factor, namely convergence flow skin:
S c , m = h D L f D m 1 C f D m [ ln ( h D 2 r w D ) π 2 ]
Equation (16) is rewritten as the final expression by
p f D m ( x D m ) = p f D m ( 0 ) + q w D m S c , m + 2 π C f D m 0 x D m d ζ 0 ζ q f D m ( ξ ) d ξ 2 π C f D m q w D m G ( x D m , x w D m )

3.3.1. Model of a Conductivity Fracture with Non-Darcy Flow

Non-Darcy flow may occur in the fracture when the inertial forces may no longer be neglected compared with viscous forces. It is very common near production wells where local velocities can be very high. Extra pressure drop is required to overcome the inertial forces.
To account for the effect of extra pressure loss caused by high-velocity flow, which is proportional to the square of the velocity, the Forchheimer flow equation is presented [27]:
p = μ k v + β ρ v | v |
The apparent fracture permeability in Equation (19) is defined as
k f , a p p = k f 1 1 + k f β ( ρ | v | / μ )
The corresponding apparent permeability is rewritten by
k f , a p p = k f 1 1 + ( q D N D ) f | q c D ( x D ) |
Substituting Equation (21) into Equation (18) yields the following form:
C f D , a p p ( q c D m ) p f D m x D m + 2 π q c D m + 2 π q w D m H ( x D m , x w D m ) = 0
where the apparent conductivity is written as C f D , a p p ( q c D m ) = C f D m ( x D ) 1 + ( q D N D ) f | q c D m | , and the flux on the cross section is given as q c D m = C f D m 2 π p f D m x D m . The corresponding integral form is given by
q c D m = q w D m H ( x D m , x w D m ) 0 x D m q f D m ( x ) d x
Note that Equation (22) presents a general model for a uniform-conductivity under the non-Darcy-flow condition.

3.3.2. Model of a Conductivity Fracture with Pressure Sensitivity Effect

The variation of permeability caused by a stress change can be expressed as a function of pore pressure. For hydraulic fracture, a general model presented by Zhang et al. [31] is introduced to describe the relationship between fracture permeability and pore pressure next:
k f ( p f ) = k min + [ k f ( p i ) k f min ] exp [ γ f ( p i p f ) ]
According to the definitions of the dimensionless parameters, the pressure-sensitive conductivity of a fracture is
C f D ( p f D ) C f D i = ( 1 C f D min C f D i ) exp [ γ f D p f D ] + C f D min C f D i
Substituting Equation (25) into Equation (14) yields the following equation governing flow in the m-th fracture flowing:
x D m ( C f D m ( p f D m ) p f D m x D m ) 2 π q f D m ( x D m ) + 2 π q w D m δ ( x D m , x w D m ) = 0
Note that Equation (26) presents a general model for a uniform-conductivity with the pressure sensitivity effect.

4. Semi-Analytical Solution for Coupled Model

4.1. Dimension Transformation

Compared with the equation of uniform conductivity fracture described by Equations (14), (22) and (26) are subject to nonlinear equations. In the case of considering non-Darcy flow, the fracture conductivity is a function of flowing rate; in the case of considering pressure-sensitivity effect, the fracture conductivity is a function of pressure. In essence, both flow rate and pressure chang with the spatial variable. Here, we attempt to solve the problem in the same way as the linear treatment in Equation (14) by a dimension transformation, which is defined as
{ ξ D m ( x D m ) = C ^ f D m 0 x D m 1 C f D m ( x D ) d x D C ^ f D m = L f D m / 0 L f D m 1 C f D n ( x D ) d x D
Then, the fracture-flow equations (Equations (22) and (26)) can be written as
p f D m ( ξ D m ) = p f D m ( ξ D m = 0 ) + q w D m S c , m + 2 π C ^ f D m ( 0 ξ D m d ζ 0 ζ q ˜ f D m ( ς ) d ς q w D m G ( ξ D m , ξ w D m ) )
Equation (28) has the same form as the fracture-flow equations for a uniform-conductivity fracture with constant conductivity given by Equation (18). The difference is that the variable ξD is a function of the variable xD and depends on the distribution of conductivity CfD(xD).

4.2. Discretization

In this study, the productivity index is achieved by discretizing the fracture panel into M segments with equal length ∆xDi, as shown in Figure 6a. Note that the equal-length fracture segment would be transformed into an unequal-length segment after using dimension transformation, as shown in Figure 6b.
In the reservoirs, the dimensionless pressure drop on the j-th segment of the n-th fracture caused by the i-th segment of the m-th fracture is given by
p m D n , j = m = 1 N f q w D m B D m n , j + m = 1 N f i = 1 M m q f D m , i Δ p u D m , i n , j Δ ξ D m , i Δ x D m , i
where B D m n , j = I B n , j ( B 2 C 1 , m B 1 C 2 , m A 1 B 2 A 2 B 1 ) + I D n , j ( A 1 C 2 , m A 2 C 1 , m A 1 B 2 A 2 B 1 ) , and Δ p u D m , i n , j = x o f D m , i 1 x o f D m , i δ p u D ( u ) d u .
In the fracture, the dimensionless pressure on the j-th segment of the n-th fracture is expressed as
p w D n p f D n , j = 2 π C ^ f D n q w D n G ( ξ D n , j , ξ w D n ) [ I ( ξ D n , j ) I ( ξ w D ) ] q w D n S c , n
where I ( ξ D ) = 0 ξ D d ζ 0 ζ q ˜ f D n ( ς ) d ς .
According to the continuity condition that the pressure must be continuous along the fracture surface, the following conditions must hold along the fracture plane:
p m D ( x o f D n + x D n , j , y o f D n ) = p f D ( x D n , j )
Moreover, the sum of all the segments’ flow rate equals to the total flow rate at the m-th fracture:
q w f D m = i = 1 M m ( q ˜ f D m , i Δ ξ D m , i Δ x D m , i )
Combining Equations (29)–(32), we can obtain the fundamental coupled solutions by computing the apparent linear equations [36,37].

4.3. Computation Consideration

According to the previous statement, the dimensionless conductivity of CfDn is the function of pressure pfDn in the fracture with pressure-dependent conductivity, and the dimensionless conductivity of CfDn is the function of flowing rate qcDn in the fracture under non-Darcy flow. In a mathematical context, CfDn is a function with regard to the spatial variable. At a given k-th step, if the distribution of pressure pfDn or flowing rate qcDn was known, CfDn along the fracture would be obtained. Thus, the nonlinear governing equation for the (k + 1)-th step could be linearized on the assumption of known conductivity distribution on the k-th step.
In the case of considering the non-Darcy flow effect (Case 1), the formation for the k-th step yields
C f D n 1 + ( q D N D ) f n q c D n k p f D n k + 1 x D n + 2 π q c D n k + 1 ( x D n ) + 2 π q w D n k + 1 H ( x D n , x w D n ) = 0
In the case of considering the pressure-dependent conductivity effect (Case 2), the formation for the k-th step yields
C f D n ( p f D n k , x D ) p f D n k + 1 x D n + 2 π [ q w D n k + 1 H ( x D m , x w D m ) 0 x D m q f D m k + 1 ( x ) d x ] + 2 π q w D n k + 1 H ( x D n , x w D n ) = 0
The flow distribution qfDn and pressure distribution pfDn along the fracture at the k-th time step would be achieved based on coupled solution. Next, the calculated pfDn was used to update the fracture conductivity (Case 1); the calculated qcDn was used to update the fracture conductivity (Case 2). The iterative procedure is repeated until the wellbore-pressure was converged. The calculation procedure is given as follows:
  • Step 1: Initial calculation, with k = 0, the fracture is assumed to be uniform (Case 1 and Case 2). By combing through Equation (29) to Equation (32), we can obtain (qfDn)k. The fracture pressure (pfDn)k (Case 1) and flowing rate (qfDn)k (Case 2) would be then achieved from Equation (28).
  • Step 2: Calculating the pressure-sensitive conductivity CfDn[(pfDn)k] (Case 1) and CfDn[(qcDn)k] (Case 2) along fracture, and then transforming xDn into ξDn based on Equation (26).
  • Step 3: Solving Equation (33) with the updated CfDn[(pfDn)k] to achieve (pfDn)k+1 and (pwD)k+1 in Case 1; solving Equation (34) with the updated CfDn[(qcDn)k] to achieve (qcDn)k+1 and (pwD)k+1 in Case 2.
  • Step 4: Terminate the iterative procedure if |(pwD)k+1 − (pwD)k|/(pwD)k < ε; otherwise, update pressure distribution along fracture by setting (pfDn)k = (pfDn)k+1 (Case 1) and flow distribution along fracture by setting (qcDn)k = (qcDn)k+1 (Case 2), and k = k + 1 back to step 2 until the convergence is achieved.

5. Results and Sensitivity Analysis

The productivity index (PI) is defined as the ratio of the production rate to pressure drawdown, as follows:
J D = μ B 2 π k m h q p i p w = 1 p w D = ( 1 2 ln 4 A e γ C A r w e 2 ) 1
Put another way, the dimensionless PI is the reciprocal of flow resistivity. Here, A is the drainage area, γ is the Eular constant, and CA is the shape factor. Specially, rwe is the effective wellbore radius determined by the geometry of drainage area, well configuration, and nonlinear flow mechanism.

5.1. Influence of Fracture Properties on PI

In this subsection, the nonlinear flow mechanisms are firstly not taken into account, and the focus is put on the impacts of the dimensions and configuration of MFHW on PI behavior. For the purpose of discussion, homogeneous fracture properties and MFHW configuration are assumed; the fractures are evenly spaced along horizontal wellbore, and all fractures are of identical properties. As a result, the PI is mainly determined by five factors, including dimensionless conductivity (CfD), fracture number (Nf), penetration ratio of fracture length to inner SRV width (Ix = Lf/xe), penetration ratio of fractured horizontal wellbore to inner SRV length (Iy = Df/ye) (Df is defined as the distance between two outmost fractures), and the drainage ratio of inner drainage to the whole drainage (Ie = xe/xR).
Figure 7 shows the effect of dimensionless conductivity on PI in the case of different penetration ratios. We further assume that the fracture is symmetrical (i.e., the fracture length on the right hand side equals to the left hand side with regard to wellbore), and all fractures are located in the center of the drainage. The fracture conductivity changes from 0.1 to 1000 (i.e., CfD = 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000). As shown in Figure 7a, at a given penetration ratio (Ix), PI is increased with the increase in dimensionless conductivity. However, after a certain dimensionless conductivity, the increase could be miniscule. In other words, beyond a certain dimensionless conductivity, the PI increase with conductivity essentially equals the PI decrease caused by the inner interference within the fracture. Generally, the threshold value is regarded as CfD,threshold = 300. When the CfD is larger than 300, there is no significant pressure drop within the fracture; the PI reaches the maximum, and does not increase with the conductivity. Figure 7b shows the PI derivative with regard to ln(CfD) in the semi-log plot. For each given Ix, the conductivity, corresponding to the maximum PI derivative, is defined as certain conductivity. Taking Ix = 1, for example (JDmax = 7.25), when the derivative is at a maximum, and that the certain value is CfD,certain = 7.5, the corresponding value of PI is JD = 4.54. This indicates that the ratio of JD/JDmax equals 62.6. This means that PI of MFHW at CfD = 7.5 could reach 62.6% of the maximum, but PI only increases by 37.4% when the CfD is further increased from 7.5 to 300. Therefore, the optimum dimensionless conductivity is defined to as the certain value. The optimum conductivity indicates that the inflow from the reservoir could match the outflow of the fracture. As analyzed from Figure 7b, the larger the Ix, the larger the optimum CfD.
Figure 8 shows the effect of the number of fractures on PI behavior. As expected, the cases with a larger number exhibit a higher value of PI. That is, increasing the number contributes to an increased connected area of MHFW with the reservoirs, and an increase in the fracture–fracture interference. However, as the number of fractures increases beyond a certain number, the increase of PI could be miniscule. For example, PI increases by 8.97% when the number is increased from Nf = 2 to Nf = 3. In comparison, incremental PI would decline to 0.39% when the fracture number is further increased from Nf = 9 to Nf = 10. Overall, the effect of the increased connected area could offset the effect of the increase of fracture interference, which results from the increased fracture number. Besides, the optimum conductivity is decreased with the increase in fracture number, as shown in Figure 8b.
Figure 9 presents the effect of the ratio of inner SRV area to the whole area, i.e., Ie. Here the value that Ie = 1 means that the area of inner SRV equals to the whole drainage area. Compared with the PI that Ie = 1, the more approaching unity the value is, the larger the PI is, which is caused by the decrease of the proportion of inner SRV. Moreover, the tendency is more noticeable in the condition of high conductivity. When CfD = 0.1, the PI at Ie = 0.5 is that JD@Ie=0.5 = 0.45, while the PI at Ie = 1 is that JD@Ie=1 = 0.65. Hence, the ratio of Ie = 1 to Ie = 0.5 in PI value is 1.44. As a comparison, when CfD = 1000, the PI at Ie = 0.5 is that JD@Ie=0.5 = 0.72, while the PI at Ie = 1 is that JD@Ie=1 = 1.95. The ratio in PI is 2.7.

5.2. Influence of Non-Darcy Flow on PI

In this subsection, the non-Darcy flow mechanism is taken into account, and we investigate the effect of non-Darcy flow on PI behavior. Figure 10 shows the flow distribution along the fracture, where CfD = 1. In this figure, the conductivity is set to be relatively low. It is shown that the influx density, which is defined as the influx rate per length along fracture face, is concentrated around the wellbore and fracture tips. It is noted that the concentration region of the influx density is controlled by a larger pressure drop/depletion. The variable of (qDND)f, defined in Equation (2), is the Forchheimer number for a given inertial factor β. When the (qDND)f is increased, the influx density is increasingly more concentrated towards wellbore. Put another way, an extra pressure drop is required to offset the effect of non-Darcy flow caused by the inertial force. Figure 11 presents the effect of fracture conductivity on influx-density distribution along the fracture. Without the non-Darcy effect (Forchheimer number equals zero) shown in Figure 11a, the increase in conductivity makes the distribution of fluid influx more gentle. When CfD > 300 (i.e., infinite conductivity), which indicates that the pressure on any point in fracture panel is same as the pressure on the wellbore, the inflow flux is concentrated on the fracture tips. With the non-Darcy effect (Figure 11b), when CfD > 1, more volume of influx fluid is concentrated on the nearby region of wellbore. The phenomenon is consistent with the characteristics given in Figure 11a. As the conductivity increases, the effect of non-Darcy flow on the distribution of inflow flux becomes weak, until the distribution of inflow flux is independent of the non-Darcy flow effect.
Figure 12 shows the effect of non-Darcy flow on PI behavior. As shown in Figure 12a, at any given conductivity, PI is the highest in the condition of (qDND)f = 0. PI increases with the increase of CfD, until it reaches the maximum value (JDmax) at CfD = 300. However, the increasing rate of PI is gradually slowed down with CfD; likewise, the maximum PI is independent of (qDND)f. In other words, non-Darcy flow plays a negative role in determining PI for the range of low- and intermediate-conductivity (CfD < 300), but has no effect on the maximum PI, which corresponds to the infinite conductivity. The relationship between PI and the Forchheimer number is further shown in Figure 12b. It is shown that, although PI declined with the Forchheimer number, the decreasing rate slows down.
Figure 13 shows the effect of configuration and dimension of fractures on PI behavior with the consideration of the non-Darcy flow effect. As shown in Figure 13a, PI loss in the case of a lower fracture number is relatively more noticeable than that in the case of a larger number. For example, when Nf = 2 and (qDND)f = 0, PI is that JD = 0.96; when Nf = 2 and (qDND)f = 10, PI is that JD = 0.75. The relative loss in PI is up to 22%, which is defined as 1 − JD@(qDND)f=10/JD@(qDND)f=0. In comparison, the relative loss in PI is only 9%. Figure 13b shows the penetration ratio of fracture length with respect to the inner SRV width (Ix). The effect of Forchheimer number is weaker when the Ix is relatively small, but its effect would be amplified when the Ix is relatively small. Figure 13c shows the penetration ratio of the inner SRV region with regard to whole drainage (Ie). In the small range (Ie < 0.9), the relationship between Ie and PI exhibits an approximate linear behavior. When Ie > 0.9, PI is increased rapidly with the Ie.

5.3. Influence of the Pressure-Sensitivity Effect on PI

The characteristic of the PI behavior is similar to the behavior of non-Darcy flow. First, Figure 14 shows the influx-distribution along the fracture with the consideration of the pressure-sensitivity effect. Compared with the nonsensitive case (Figure 14a), more influx density is distributed around the wellbore in the sensitive case (Figure 14b), and the tendency would be more remarkable with the decrease of the conductivity. The reason is explained in that the apparent conductivity is degraded and the flow resistance is consequently increased, which is caused by the closing of the fracture under the influence of rock compaction.
Figure 15 shows the effect of pressure sensitivity on PI behavior. The intensity is measured by the dimensionless variable γfD. As shown in Figure 15a, the effect of pressure sensitivity on PI becomes weak with the increase of the conductivity. When the magnitude of the conductivity reaches the level of infinite conductivity, the effect of pressure sensitivity is almost neglected. As expected, as presented in Figure 15b, PI declines until the minimum with the increasing γfD, and corresponding γfD is defined as the threshold. Meanwhile, the decreasing rate slows down. For example, if CfD = 1, PI reaches the minimum JDmin when γfDthreshold = 1.2; if CfD = 10, PI reaches the minimum JDmin when γfDthreshold = 2.7. Thus, a larger conductivity contributes to a larger γfDthreshold.
Figure 16 shows the effect of configuration and dimension of fractures on PI in the consideration of pressure sensitivity in the fracture. The overall features are similar to the case of considering non-Darcy flow, so we will not rewrite again.

6. Conclusions

We derived a new steady-state productivity model for a multi-fractured horizontal well in a reservoir with pseudo-elliptical boundary. We summarize some important findings as follows:
  • PI increases with the increase of CfD, until it reaches the maximum (JDmax) at CfD = 300.
  • PI is deteriorated under the influence of nonlinear flow mechanisms. With the consideration of non-Darcy flow, for the small range of the penetration ratio of the inner SRV region with regard to whole drainage (Ie < 0.9), the relationship between Ie and PI exhibits an approximately linear behavior. When Ie > 0.9, PI is increased rapidly with the Ie.
  • With the consideration of pressure sensitivity, the apparent fracture is degraded from the initial conductivity to the minimal conductivity, which is caused by fracture closure. As a result, an extra pressure drop is acquired to offset the conductivity degradation, and PI would be declined to the minimal PI.
  • The disadvantage dimensions of fracture (such as small conductivity, penetration ratio, and less fractures) contribute to severe pressure depletion, while, in turn, the severe pressure depletion will strengthen the effect of the nonlinear flow mechanism on PI behavior.
  • If the conductivity in the fracture reaches the level of infinite conductivity, the influence of nonlinear flow mechanism on PI could be neglected.

Author Contributions

M.C. established the model and wrote the manuscript, and H.X. supervised this research and revised the manuscript. C.W. validated the model and conducted the sensitivity analysis. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Science and Technology Major Project (No. 2017ZX05019005-006) and Heilongjiang Natural Science Foundation (LH2019F004). Special thanks to the anonymous reviewers and editors for their valuable comments.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Analytical Solution for the Inner SRV Region

The appendix is an optional section that can contain details and data supplemental to the main text. For example, explanations of experimental details that would disrupt the flow of the main text, but nonetheless remain crucial to understanding and reproducing the research shown; figures of replicates for experiments of which representative data is shown in the main text can be added here if brief, or as Supplementary data. Mathematical proofs of results not central to the paper can be added as an appendix.
According to the governing equation and boundary conditions, Fourier sine transformation are used to deal with Equation (4), which is defined as
F x [ f ( x D ) ] = 0 x e D f ( x D ) sin ( β υ x D ) d x D = f ¯ ( β υ )
F y [ f ( y D ) ] = 0 y e D f ( y D ) sin ( γ ν y D ) d y D = f ¯ ( γ ν )
β υ = υ π x e D ,   υ = 0 , 1 , 2 , 3 ,   γ ν = ν π y e D ,   ν = 0 , 1 , 2 , 3 ,
Taking the Fourier transformation of Equations (4) and (5), the first term on the left hand side (LHS) of the original equation is given as
F y [ F x ( 2 p m D x D 2 ) ] = β υ p e D x [ 1 ( 1 ) υ ] γ ν [ 1 cos ( γ ν y e D ) ] β υ 2 p ¯ ^ m D
Likewise, the second term on the LHS is given as
F x [ F y ( 2 p m D y D 2 ) ] = γ ν p e D y [ 1 ( 1 ) ν ] β υ [ 1 cos ( β υ x e D ) ] γ ν 2 p ¯ ^ m D
The third term on the LHS is given as
0 L f D m q f D m ( u D ) δ ( x D x o f D m u D ) δ ( y D y o f D m ) d u D = 0 L f D q f D ( u D ) [ sin [ β n ( x o f D + u D ) ] sin ( γ m y o f D ) ] d u D
Finally, the Fourier solution for Equation (4) is given as
β υ p e D x [ 1 ( 1 ) υ ] [ 1 cos ( γ ν y e D ) ] γ ν β υ 2 p ¯ ^ m D + γ ν p e D y [ 1 ( 1 ) ν ] [ 1 cos ( β υ x e D ) ] β υ γ ν 2 p ¯ ^ m D + 2 π m = 1 N f 0 L f D m q f D m ( u D ) { sin [ β υ ( x o f D m + u D ) ] sin ( γ ν y o f D m ) } d u D = 0
The inversion Fourier sine transformation is written as
p ¯ ^ m D = p e D x β υ [ 1 ( 1 ) υ ] [ 1 cos ( γ ν y e D ) ] γ ν ( β υ 2 + γ ν 2 ) + p e D y γ ν [ 1 ( 1 ) ν ] [ 1 cos ( β υ x e D ) ] β υ ( β υ 2 + γ ν 2 ) + 2 π m = 1 N f 0 L f D m q f D m ( u D ) [ sin [ β υ ( x o f D m + u D ) ] sin ( γ ν y o f D m ) β υ 2 + γ ν 2 ] d u D
where
N ( β υ ) = 0 x e D sin 2 ( β υ x D ) d x D = x e D 2 ( 1 sin ( β υ x e D ) cos ( β υ x e D ) β υ x e D )
N ( γ ν ) = 0 y e D sin 2 ( γ ν y D ) d y D = y e D 2 ( 1 sin ( γ ν y e D ) cos ( γ ν y e D ) γ ν y e D )
Accordingly, the original solution is expressed in forms of Fourier solution,
p m D = υ = 0 sin ( β υ x D ) N ( β υ ) [ ν = 0 sin ( γ ν y D ) N ( γ ν ) p ¯ ^ m D ( β υ , γ ν ) ]
Equation (A11) could be further expressed as
p m D = 4 m e D x x e D y e D ν = 1 sin ( γ ν y D ) [ 1 cos ( γ ν y e D ) ] γ ν υ = 1 sin ( β n x D ) β υ [ 1 ( 1 ) υ ] ( β υ 2 + γ ν 2 ) V 1 + 4 m e D y x e D y e D υ = 1 sin ( β υ x D ) [ 1 cos ( β υ x e D ) ] β υ ν = 1 sin ( γ ν y D ) γ ν [ 1 ( 1 ) ν ] ( γ ν 2 + β υ 2 ) V 2 + 8 π x e D y e D m = 1 N f 0 L f D m q f D m ( u D ) ( υ = 1 sin ( β υ x D ) sin [ β υ ( x o f D m + u D ) ] ν = 1 sin ( γ ν y D ) sin ( γ ν y o f D m ) β υ 2 + γ ν 2 V 3 ) d u D
where
V 1 = υ = 1 β υ sin ( β υ x D ) β υ 2 + γ ν 2 + υ = 1 ( 1 ) υ 1 β υ sin ( β υ x D ) β υ 2 + γ ν 2 = x e D 2 sinh [ ν π ( x e D x D ) / y e D ] + sinh ( ν π x D / y e D ) sinh ( ν π x e D / y e D )
V 2 = y e D 2 sinh [ υ π ( y e D y D ) / x e D ] + sinh ( υ π y D / x e D ) sinh ( υ π y e D / x e D )
V 3 = x e D y e D 4 υ π cosh [ υ π ( y e D | y D y o f D | ) / x e D ] cosh { υ π [ y e D ( y D + y o f D ) ] / x e D } sinh ( υ π y e D / x e D )
Correspondingly, Equation (A12) is finalized by
p m D ( x D , y D ) = p e D x I B ( x D , y D ) + p e D y I D ( x D , y D ) + 2 m = 1 N f 0 L f D m q f D m ( u D ) δ p u D ( x D , y D , u D ) d u D
where
{ I B = 2 π υ = 1 1 υ sin ( υ π y D y e D ) [ 1 + ( 1 ) υ 1 ] sinh ( υ π x e D / y e D ) ( sinh υ π ( x e D x D ) y e D + sinh υ π x D y e D ) I D = 2 π υ = 1 1 υ sin ( υ π x D x e D ) [ 1 + ( 1 ) υ 1 ] sinh ( υ π y e D / x e D ) ( sinh υ π ( y e D y D ) x e D + sinh υ π y D x e D ) δ p u D = υ = 1 1 υ sin ( β υ x D ) sin [ β υ ( x o f D m + u D ) ] sinh ( υ π y e D / x e D ) cosh n π y e D | y D ± y o f D m | x e D

Appendix B. Analytical Solution for the Outer Region

On the interface between the inner SRV region and outer elliptical region, the continuity condition is given as
0 π p m D ξ D | ξ D = 0 d η D e l l i p t i c   f l o w = 0 x e D p m D y D | y D = y e D d x D l i n e a r   f l o w
In detail, the right hand side (RHS) of Equation (A18) is expressed as
R H S = p e D x A 2 p e D y ( 1 x e D ln [ x R D + x R D 2 x e D 2 / x e D ] + B 2 ) + m = 1 N f q w D m C 2 , m )
where
{ A 2 = 4 π x e D υ = 1 1 υ [ 1 + ( 1 ) υ 1 ] sinh ( υ π x e D / y e D ) ( cosh υ π x e D y e D 1 ) B 2 = 1 x e D ln [ ( x R D + x R D 2 x e D 2 ) / x e D ] 2 π x e D υ = 1 1 υ [ 1 + ( 1 ) υ 1 ] 2 sinh ( υ π y e D / x e D ) ( cosh υ π y e D x e D 1 ) C 2 , m = 8 π L f D m υ = 1 1 υ 2 [ 1 ( 1 ) υ ] sinh ( υ π y e D / x e D ) sin ( υ π L f D m 2 x e D ) sin ( υ π ( x o f D m + 0.5 L f D m ) x e D ) sinh υ π y o f D m x e D
Thus, Equation (A18) is written as
A 2 p e D x + B 2 p e y D = m = 1 N f q w D m C 2 , m
Likewise, on the interface between the inner SRV region and outer linear region, the continuity condition is given as
1 y e D 0 y e D p m D x D | x D = x e D d y D = 1 y e D 0 y e D p m D x D | x D = x e D d y D
In detail, the right hand side (RHS) of Equation (A22) is expressed as
R H S = p e D x ( 1 x R D x e D A 1 ) p e D y B 1 m = 1 N f q w D m C 1 , m
where
{ A 1 = 1 ( x R D x e D ) 1 π y e D υ = 1 1 υ [ 1 + ( 1 ) υ 1 ] 2 sinh ( υ π x e D / y e D ) ( cosh υ π x e D y e D 1 ) B 1 = 2 π y e D υ = 1 1 υ [ 1 + ( 1 ) υ 1 ] sinh ( υ π y e D / x e D ) ( cosh υ π y e D x e D 1 ) C 1 , m = 4 x e D π y e D L f D m υ = 1 1 υ 2 ( 1 ) υ sinh ( υ π y e D / x e D ) sin ( υ π L f D m 2 x e D ) sin ( υ π ( x o f D m + 0.5 L f D m ) x e D ) sinh υ π ( y e D y o f D m ) x e D
Thus, Equation (A22) is written as
A 2 p e D x + B 2 p e y D = m = 1 N f q w D m C 2 , m

References

  1. Clarkson, C.R.; Bustin, R.M.; Seidle, J.P. Production-data analysis of single-phase (gas) coalbed-methane wells. SPE Reserv. Eval. Eng. 2007, 3, 312–341. [Google Scholar] [CrossRef]
  2. Wang, J.H.; Wang, X.D.; Dong, W.X. Rate decline curves analysis of multiple-fractured horizontal wells in heterogeneous reservoirs. J. Hydrol. 2017, 553, 527–539. [Google Scholar] [CrossRef]
  3. Deng, J.; Zhu, W.Y.; Ma, Q. A new seepage model for shale gas reservoir and productivity analysis of fractured well. Fuel 2014, 124, 232–240. [Google Scholar] [CrossRef]
  4. Clarkson, C.R. Production data analysis of unconventional gas wells: Review of theory and best practices. Int. J. Coal Geol. 2013, 109–110, 101–146. [Google Scholar] [CrossRef]
  5. Zhang, Q.; Su, Y.L.; Wang, W.D.; Lu, M.; Ren, L. Performance analysis of fractured wells with elliptical SRV in shale reservoirs. J. Nat. Gas Sci. Eng. 2017, 45, 380–390. [Google Scholar] [CrossRef]
  6. Wang, J.L.; Wei, Y.S.; Luo, W.J. A unified approach to optimize fracture design of a horizontal well intercepted by primary- and secondary-fracture networks. SPE J. 2019, 24, 1270–1287. [Google Scholar] [CrossRef]
  7. Al-Rbeawi, S. Productivity-index behavior for hydraulically fractured reservoirs depleted by constant production rate considering transient-state and semisteady-state conditions. SPE Prod. Oper. 2018, 33. [Google Scholar] [CrossRef]
  8. Warren, J.E.; Root, P.J. The behavior of naturally fractured reservoirs. SPE J. 1963, 3, 245–255. [Google Scholar] [CrossRef] [Green Version]
  9. Obinna, E.D.; Hassan, D. A model for simultaneous matrix depletion into natural and hydraulic fracture networks. J. Nat. Gas Sci. Eng. 2014, 16, 57–69. [Google Scholar]
  10. Stalgorova, E.; Mattar, L. Analytical model for unconventional multifractured composite systems. SPE Reserv. Eval. Eng. 2013, 16, 246–256. [Google Scholar] [CrossRef]
  11. Al-Ghamdi, A.; Chen, B.; Behmanesh, H.; Qanbari, F.; Aguilera, R. An improved triple-porosity model for evaluation of naturally fractured reservoirs. SPE Reserv. Eval. Eng. 2011, 8, 397–404. [Google Scholar] [CrossRef]
  12. Karimi-Fard, M.; Durlofsky, L.J. A general gridding, discretization, and coarsening methodology for modeling flow in porous formations with discrete geological features. Adv. Water Resour. 2016, 96, 354–372. [Google Scholar] [CrossRef]
  13. Xu, Y.; Cavalcante Filho, J.S.A.; Yu, W. Discrete-fracture modeling of complex hydraulic-fracture geometries in reservoir simulators. SPE Reserv. Eval. Eng. 2017, 20, 403–422. [Google Scholar] [CrossRef]
  14. Li, L.Y.; Lee, S.H. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reserv. Eval. Eng. 2008, 8, 750–758. [Google Scholar] [CrossRef]
  15. Cinco-Ley, H.; Samaniego, F.; Dominguez, N. Transient pressure behavior for a well with a finite-conductivity vertical fracture. Soc. Pet. Eng. J. 1978, 8, 254–265. [Google Scholar]
  16. Wang, J.L.; Jia, A.L.; Wei, Y.S.; Qi, Y.D. Approximate semi-analytical modeling of transient behavior of horizontal well intercepted by multiple pressure-dependent conductivity fractures in pressure-sensitive reservoir. J. Pet. Sci. Eng. 2017, 153, 157–177. [Google Scholar] [CrossRef]
  17. Wang, J.L.; Jia, A.L.; Wei, Y.S. A general framework model for simulating transient response of a well with complex fracture network by use of source and Green’s function. J. Nat. Gas Sci. Eng. 2018, 55, 254–275. [Google Scholar] [CrossRef]
  18. Zhou, W.T.; Banerjee, R.; Poe, B. Semianalytical production simulation of complex hydraulic-fracture networks. SPE J. 2013, 19, 6–18. [Google Scholar] [CrossRef]
  19. Kuchuk, F.; Biryukov, D. Pressure-transient behavior of continuously and discretely fractured reservoirs. SPE Reserv. Eval. Eng. 2014, 17, 82–97. [Google Scholar] [CrossRef]
  20. Al-Rbeawi, S. The optimal reservoir configuration for maximum productivity index of gas reservoirs depleted by horizontal wells under Darcy and non-Darcy flow conditions. J. Nat. Gas Sci. Eng. 2018, 49, 179–193. [Google Scholar] [CrossRef]
  21. Chen, Z.M.; Liao, X.W.; Zhao, X.L.; Lyu, S.; Zhu, L. A comprehensive productivity equation for multiple fractured vertical well with non-linear effects under steady-state flow. J. Pet. Sci. Eng. 2017, 149, 9–24. [Google Scholar] [CrossRef]
  22. Li, J.; Guo, B.; Gao, D.; Ai, C. The effect of fracture-face matrix damage on productivity of fractures with infinite and finite conductivities in shale-gas reservoirs. SPE Drill. Complet. 2012, 9, 347–353. [Google Scholar] [CrossRef]
  23. Prats, M.; Raghavan, R. Finite horizontal well in a uniform-thickness reservoir crossing a natural fracture normally. SPE J. 2013, 10, 982–992. [Google Scholar] [CrossRef]
  24. Tahmasebi, P.; Kamrava, S. A pore-scale mathematical modeling of fluid-particle interactions: Thermo-hydro-mechanical coupling. Int. J. Greenh. Gas Control 2019, 83, 245–255. [Google Scholar] [CrossRef]
  25. Zhang, X.M.; Tahmasebi, P. Effects of grain size on deformation in porous media. Transp. Porous Media 2019, 129, 321–341. [Google Scholar] [CrossRef]
  26. Zhang, X.M.; Tahmasebi, P. Micromechanical evaluation of rock and fluid interactions. Int. J. Greenh. Gas Control 2018, 76, 266–277. [Google Scholar] [CrossRef]
  27. Guppy, K.H.; Cinco-Ley, H.; Ramey, H.J. Effect of non-Darcy flow on the constant-pressure production of fractured wells. Soc. Pet. Eng. J. 1981, 12, 390–400. [Google Scholar] [CrossRef]
  28. Wu, Y.S.; Li, J.F.; Ding, D.; Wang, C.; Di, Y. A generalized framework model for simulation of gas production in unconventional gas reservoir. In Proceedings of the SPE Reservoir Simulation Symposium, The Woodlands, TX, USA, 18–20 February 2013. [Google Scholar]
  29. Lin, J.; Zhu, D. Predicting well performance in complex fracture systems by slab source method. In Proceedings of the SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, 6–8 February 2012. [Google Scholar]
  30. Pedroso, C.A.; Correa, A.C. A new model of a pressure-dependent-conductivity hydraulic fracture in a finite reservoir: Constant rate production case. In Proceedings of the Latin American and Caribbean Petroleum Engineering Conference, Rio de Janeiro, Brazil, 30 August–3 September 1997. [Google Scholar]
  31. Zhang, Z.; He, S.L.; Liu, G.; Guo, X.; Mo, S. Pressure buildup behavior of vertically fractured wells with stress-sensitive conductivity. J. Pet. Sci. Eng. 2014, 122, 48–55. [Google Scholar] [CrossRef]
  32. Yao, S.S.; Zeng, F.H.; Liu, H. A semi-analytical model for hydraulically fractured wells with stress-sensitive conductivities. In Proceedings of the SPE Unconventional Resources Conference, Calgary, AB, Canada, 5–7 November 2012. [Google Scholar]
  33. Biryukov, D.; Kuchuk, F.J. Transient pressure behavior of reservoirs with discrete conductive faults and fractures. Transp. Porous Media 2012, 95, 239–268. [Google Scholar] [CrossRef] [Green Version]
  34. Rasoulzadeh, M.; Kuchuk, F.J. Pressure transient behavior of high-fracture-density reservoirs (dual-porosity models). Transp. Porous Media 2019, 127, 222–243. [Google Scholar] [CrossRef]
  35. Kuchuk, F.J.; Habashy, T. Pressure behavior of laterally composite reservoirs. SPE Form. Eval. 1997, 3, 47–56. [Google Scholar] [CrossRef]
  36. Luo, W.J.; Tang, C.F. A semianalytical solution of a vertical fractured well with varying conductivity under non-Darcy-flow condition. SPE J. 2015, 20, 1028–1040. [Google Scholar] [CrossRef]
  37. Luo, W.J.; Tang, C.F.; Feng, Y. Mechanism of fluid flow along a dynamic conductivity fracture with pressure-dependent permeability under constant wellbore pressure. J. Pet. Sci. Eng. 2018, 166, 465–475. [Google Scholar] [CrossRef]
Figure 1. The flow regime characteristics of MHFW, where (a) sequence of flow-regimes for a MFHW, completed in a tight sand gas reservoir modified from Clarkson [4]. Lines represent isopotential lines; arrows represent streamlines. Flow-regime 1 corresponds to linear flow; flow-regime 2 corresponds to elliptical flow; flow-regime 3 corresponds to fracture interference; flow-regime 4 corresponds to compound linear, and flow-regime 5 corresponds to compound elliptical flow. (b) simulation of flow regimes for a stimulated MFHW in a tight sand gas reservoir. (simulated through numerical simulator).
Figure 1. The flow regime characteristics of MHFW, where (a) sequence of flow-regimes for a MFHW, completed in a tight sand gas reservoir modified from Clarkson [4]. Lines represent isopotential lines; arrows represent streamlines. Flow-regime 1 corresponds to linear flow; flow-regime 2 corresponds to elliptical flow; flow-regime 3 corresponds to fracture interference; flow-regime 4 corresponds to compound linear, and flow-regime 5 corresponds to compound elliptical flow. (b) simulation of flow regimes for a stimulated MFHW in a tight sand gas reservoir. (simulated through numerical simulator).
Energies 13 02015 g001
Figure 2. Schematic of the compound elliptical-flow model representing four contiguous regions for a multiple-fractured horizontal well (MFHW). SRV, stimulated rock volume.
Figure 2. Schematic of the compound elliptical-flow model representing four contiguous regions for a multiple-fractured horizontal well (MFHW). SRV, stimulated rock volume.
Energies 13 02015 g002
Figure 3. Schematic of MFHW in the inner region with constant-pressure outer boundary.
Figure 3. Schematic of MFHW in the inner region with constant-pressure outer boundary.
Energies 13 02015 g003
Figure 4. Elliptical coordinate used in the outer elliptical region.
Figure 4. Elliptical coordinate used in the outer elliptical region.
Energies 13 02015 g004
Figure 5. Schematic of the flow pattern in the transverse fracture.
Figure 5. Schematic of the flow pattern in the transverse fracture.
Energies 13 02015 g005
Figure 6. Illustration of dimension transformation and fracture discretization.
Figure 6. Illustration of dimension transformation and fracture discretization.
Energies 13 02015 g006
Figure 7. Productivity-index behavior of multiple fractures, including (a) the productivity index vs. dimensionless conductivity, and (b) the derivative of productivity index vs. dimensionless conductivity.
Figure 7. Productivity-index behavior of multiple fractures, including (a) the productivity index vs. dimensionless conductivity, and (b) the derivative of productivity index vs. dimensionless conductivity.
Energies 13 02015 g007
Figure 8. Influence of number of fractures on PI behavior, including (a) the productivity index vs. dimensionless conductivity, and (b) the derivative of productivity index vs. dimensionless conductivity.
Figure 8. Influence of number of fractures on PI behavior, including (a) the productivity index vs. dimensionless conductivity, and (b) the derivative of productivity index vs. dimensionless conductivity.
Energies 13 02015 g008
Figure 9. Influence of penetration ratio of inner SRV region with regard to the whole region on productivity-index (PI) behavior.
Figure 9. Influence of penetration ratio of inner SRV region with regard to the whole region on productivity-index (PI) behavior.
Energies 13 02015 g009
Figure 10. Influx-flow distribution along the fracture under the influence of non-Darcy flow.
Figure 10. Influx-flow distribution along the fracture under the influence of non-Darcy flow.
Energies 13 02015 g010
Figure 11. Influx-flow distribution along the fracture under the influence of conductivity (a) without and (b) with non-Darcy flow effect.
Figure 11. Influx-flow distribution along the fracture under the influence of conductivity (a) without and (b) with non-Darcy flow effect.
Energies 13 02015 g011
Figure 12. (a) Influence of fracture conductivity on PI behavior under the effect of non-Darcy flow, (b) influence of Forchheimer number on PI behavior with the consideration of finite conductivity.
Figure 12. (a) Influence of fracture conductivity on PI behavior under the effect of non-Darcy flow, (b) influence of Forchheimer number on PI behavior with the consideration of finite conductivity.
Energies 13 02015 g012
Figure 13. Influence of dimensions of fracture on PI behavior in the consideration of non-Darcy flow, including (a) the influence of the number of fractures, (b) the influence of penetration ratio of fracture, and (c) the influence of penetration ratio of horizontal well.
Figure 13. Influence of dimensions of fracture on PI behavior in the consideration of non-Darcy flow, including (a) the influence of the number of fractures, (b) the influence of penetration ratio of fracture, and (c) the influence of penetration ratio of horizontal well.
Energies 13 02015 g013
Figure 14. Influence of dimensionless conductivity on influx-flow distribution (a) with and (b) without the pressure-sensitivity effect.
Figure 14. Influence of dimensionless conductivity on influx-flow distribution (a) with and (b) without the pressure-sensitivity effect.
Energies 13 02015 g014
Figure 15. (a) Influence of fracture conductivity on PI behavior with the pressure-sensitivity effect, and (b) influence of permeability modulus on PI behavior with consideration of finite conductivity.
Figure 15. (a) Influence of fracture conductivity on PI behavior with the pressure-sensitivity effect, and (b) influence of permeability modulus on PI behavior with consideration of finite conductivity.
Energies 13 02015 g015
Figure 16. Influence of dimensions of the fracture on PI behavior in the consideration of pressure sensitivity, including (a) number of fracture, (b) penetration index of fracture, and (c) penetration of horizontal wellbore.
Figure 16. Influence of dimensions of the fracture on PI behavior in the consideration of pressure sensitivity, including (a) number of fracture, (b) penetration index of fracture, and (c) penetration of horizontal wellbore.
Energies 13 02015 g016

Share and Cite

MDPI and ACS Style

Cao, M.; Xiao, H.; Wang, C. Productivity-Index Behavior for a Horizontal Well Intercepted by Multiple Finite-Conductivity Fractures Considering Nonlinear Flow Mechanisms under Steady-State Condition. Energies 2020, 13, 2015. https://doi.org/10.3390/en13082015

AMA Style

Cao M, Xiao H, Wang C. Productivity-Index Behavior for a Horizontal Well Intercepted by Multiple Finite-Conductivity Fractures Considering Nonlinear Flow Mechanisms under Steady-State Condition. Energies. 2020; 13(8):2015. https://doi.org/10.3390/en13082015

Chicago/Turabian Style

Cao, Maojun, Hong Xiao, and Caizhi Wang. 2020. "Productivity-Index Behavior for a Horizontal Well Intercepted by Multiple Finite-Conductivity Fractures Considering Nonlinear Flow Mechanisms under Steady-State Condition" Energies 13, no. 8: 2015. https://doi.org/10.3390/en13082015

APA Style

Cao, M., Xiao, H., & Wang, C. (2020). Productivity-Index Behavior for a Horizontal Well Intercepted by Multiple Finite-Conductivity Fractures Considering Nonlinear Flow Mechanisms under Steady-State Condition. Energies, 13(8), 2015. https://doi.org/10.3390/en13082015

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop