Study on Deliverability Evaluation of Staged Fractured Horizontal Wells in Tight Oil Reservoirs

Abstract

At present, the existing deliverability evaluation models mainly consider the impact of specific factors on production, and the description of the complex fracture network structure primarily remains at the stage of an ideal dual-pore medium with uniform distribution. However, this cannot reflect the actual fracture network structure and fluid flow law of fractured horizontal wells. Thus, in this paper, a non-uniform fracture network structure is proposed considering the influence of the threshold pressure gradient and stress sensitivity characteristics on the production performance of horizontal wells. The stress sensitivity and the fractal theory are combined to characterize the permeability of the complex fracture network, and a three-zone compound unsteady deliverability model for staged fractured horizontal wells in tight oil reservoirs is successfully developed. Laplace transformation, perturbation theory, and numerical inversion are applied to obtain the semi-analytical solution of the proposed deliverability model. The reliability and accuracy of the analytical solution are verified by the classical tri-linear flow model and an oil field example. The effects of related influential parameters on the production of horizontal wells are analyzed. The deliverability evaluation method proposed in this paper can provide a theoretical basis for formulating rational development technology policies for tight oil reservoirs.

1. Introduction

Tight oil reserves are abundant and the development of tight oil reservoirs has long been a popular study topic in petroleum research. At present, the technology of “long horizontal section horizontal well + large-scale volume fracturing” is often employed to develop tight oil reservoirs [1,2]. A complex fracture network is formed after large-scale volume fracturing of horizontal wells, which shortens the seepage distance of oil and gas from the matrix to the fracture [3,4], and changes the fluid seepage pattern of the area around the horizontal well, thus improving the development effect of tight oil reservoirs.

Numerous studies have been undertaken globally on the productivity evaluation of fractured horizontal wells in tight oil reservoirs, including steady-state and unsteady-state productivity models. Numerical simulation, analytical, and semi-analytical methods are extensively applied in relevant studies [5,6,7,8,9,10,11,12,13,14]. Li Longlong et al. [15] applied potential theory and the superposition principle to derive the steady-state productivity linear equations of staged multi-cluster fractured horizontal wells in the two cases of infinite and finite conductivity fractures. However, the model lacked the description of the reservoir reformation area around the primary fracture. Hu Jinghong et al. [16] established a three-zone flow numerical model and a power-rate declining model based on the theory of unsteady seepage mechanics and the principle of superposition to calculate the productivity of fractured horizontal wells. The productivity model described the primary fracture formed by fracturing and the surrounding matrix as an elliptical seepage area, without considering the complex fracture network formed around the primary fracture. In reality, because of the development of natural fractures in tight reservoirs, the fractures formed by fracturing can connect with the natural fractures to form secondary fractures, thus generating a complex fracture network. This leads to the complexity of the seepage law of horizontal wells, the changes in the productivity of horizontal wells, and the increase in the difficulty of deliverability evaluation. Therefore, experts and scholars have further studied the productivity of horizontal wells, aiming at the complex fracture network formed by fracturing. According to the fundamental physical characteristics of multi-fractured horizontal wells in unconventional reservoirs, M. Brown et al. [17] established a trilinear unsteady seepage model and obtained an analytical solution by employing idealized dual-porosity media to simulate natural fractured porous media in the reservoir. Considering the non-Darcy flow and stress sensitivity characteristics of fluids in tight reservoirs, Ji Jinghao et al. [18] established a five-zone flow productivity model for fractured horizontal wells. The reservoir reconstruction area in the model (in which the reservoir volume area containing the complex fracture network can improve the production of horizontal wells) was an idealized dual-porous medium area. Fang Wenchao et al. [19] established a multi-scale flow discrete fracture numerical simulation model considering the compressibility of rocks and fluids in the reservoir and the nonlinear flow characteristics of the matrix. The model can be employed to simulate tight reservoirs with a large-scale complex fracture network system. Sun ruofan et al. [20] established an elliptical seepage model for fractured horizontal wells and obtained a semi-analytical solution to the problem of micro-fractures in tight reservoirs affecting the fluid flow and oil well production. These models have limitations in describing the morphology of complex fracture networks. In addition, the fracture networks considered in the models were homogeneous, and were thus not in accordance with the heterogeneous characteristics of actual reservoirs after fracturing. Considering the heterogeneity of dual-porosity reservoirs and pressure sensitivity, Jiang Ruizhong et al. [21] established a dual-porosity fractal reservoir nonlinear flow model. Li Xianwen et al. [22] obtained the distribution of secondary fractures based on the inversion of micro-seismic data, used fractal theory to describe the complex fracture network, and established a five-zone seepage model for fractured horizontal wells. Xie Bin et al. [23] established a three-zone composite seepage model by applying a fractal dual-media model to characterize the pattern of the fracture network rebuilding area. However, the model lacked the description of the shape of the fractal fracture network area and ignored the stress-sensitive characteristics of the medium. Among these studies, although the fractal theory introduced to characterize the heterogeneous morphological characteristics of the complex fracture network formed by fracturing is more suitable for the actual fracturing production than the idealized homogeneous dual-pore medium model, the morphological description of the heterogeneous distribution of the fracture network remains unclear. Therefore, the productivity model considering the complex morphology of the fracture network needs to be further studied.

Based on previous research, in this study, considering the changes in the fluid flow capacity in the volume area between the primary fractures, and combining the stress sensitivity and fractal theory, a non-uniform fracture network structure model of the heterogeneous characteristics of the fracture network was established and a new deliverability evaluation method for fractured horizontal wells in tight oil reservoirs was proposed. In this method, the fractured horizontal well in the tight reservoir is divided into three seepage areas. The fracture network in the volume area between primary fractures is unevenly distributed along the horizontal wellbore, which is different from the classical tri-linear flow model. Laplace transformation, perturbation theory, and numerical inversion are applied to solve the established mathematical model for fractured horizontal wells to obtain the semi-analytical productivity formula. The accuracy of the analytical solution of the model was verified by comparing the analytical results of the simplified model and the classic model. Finally, the reliability of the model was further proved by matching the model calculation results with the actual production data of an oil well in the X Oilfield (in China).

2. Materials and Methods

2.1. Physical Model and Assumptions

2.1.1. Non-Uniform Fracture Network Structure Model

A cluster of primary fractures is first formed during the process of fracturing horizontal wells, and then the primary fractures connect with the natural fractures in the reservoir to generate secondary fractures, thus forming a fracture network around the primary fractures. It has been demonstrated in relevant studies that the density of the fracture network decreases away from the primary fracture. The Fractal Random Fracture Network Generation Algorithm (FRFNA) proposed by Sheng Guanglong et al. [24] based on micro-seismic data inverted the complex fracture network for horizontal wells in Eagle Ford shows that the further away from hydraulic fractures, the fewer the number of secondary fractures and the lower their density.

Based on this, our paper proposes a non-uniform fracture network structure model, as depicted in Figure 1 (the primary fracture is located at the left side of the model, and the secondary fracture extends to the right direction). As shown in Figure 1, the matrix rock block is unevenly divided by the fracture network, and the fracture network is non-uniformly distributed in the horizontal direction. The farther from the primary fracture, the smaller the fracture density. The gradual change in fracture density along the horizontal well leads to the uneven change in fluid flow capacity. Therefore, to consider the stress sensitivity of the fracture network, the fractal theory is introduced to characterize the complex fluid seepage law under the structure of a non-uniform fracture network.

2.1.2. Seepage Physical Model and Condition Assumptions for Fractured Horizontal Well

According to the fluid seepage characteristics of fractured horizontal wells in tight reservoirs, we established an unsteady productivity seepage physical model for a fractured horizontal well, as shown in Figure 2. The quarter reservoir controlled by a single fracture is viewed as the research object, and the seepage physical model is composed of three parts containing the main fracture area (region I), the volume area between primary fractures (region II), and the unstimulated volume area at the reservoir boundary (region III). Region I represents the primary fracture, which is a linear flow in this area, and the fluid flows directly into the horizontal wellbore along the primary fracture. Region II is the volume area between the primary fractures, which is a dual-porous medium area. There are two kinds of media in zone II: the matrix and the fracture network. The fluid in region II flows in the direction perpendicular to the primary fracture into region I. The complex fracture network in this area is composed of secondary fractures and natural fractures. Fractures are unevenly distributed in zone II. The density of secondary fractures near the primary fracture is the largest and gradually decreases in the y-direction. There are almost no secondary fractures and only a small number natural fractures at the position distant from the primary fracture. Therefore, the reservoir volume distant from the primary fracture can be regarded as the unstimulated reservoir volume. Region III is an unstimulated volume area at the reservoir boundary with a small number of natural fractures. The threshold pressure gradient is introduced to characterize the non-Darcy seepage of the fluid in this area. The fluid in region III flows in the direction parallel to the primary fracture into region II.

The basic assumptions of the physical model are as follows: (1) The tight oil reservoir is a homogeneous box-shaped formation, of which the upper and lower boundaries are closed. The horizontal wellbore passes through the center of the reservoir. The half-width of the reservoir is $x_e$ and the thickness of the reservoir is h. (2) The half distance between primary fractures is $y_e$. (3) The width of the hydraulic fracture is $w_f$ and the half-length of the fracture is $x_f$. The hydraulic fractures are symmetrically distributed perpendicular to the horizontal wellbore, and the height of the fracture is equal to the thickness of the reservoir. (4) There is a single-phase isothermal unsteady seepage in the reservoir, without considering the influence of gravity and capillary force. (5) Both rock and fluid are slightly compressible and have a constant compression coefficient. (6) The fracture network in region II is a fractal reservoir with a fractal dimension of D, which is embedded in a two-dimensional Euclidean rock block. (7) The pressure loss in the horizontal section is negligible. (8) The stress sensitivity of the fracture system and the non-Darcy seepage characteristics of the fluid in region III are taken into consideration. (9) The skin and wellbore storage effects are not considered in the model. (10) The horizontal well produces with the fixed bottom hole pressure.

2.2. Establishment of the Mathematical Model

2.2.1. Characterization of Porosity and Permeability of the Heterogeneous Fracture Network

Considering the stress-sensitive characteristics of the fracture network, the fractal theory is introduced to characterize the permeability and porosity of the complex non-uniform fracture network. It is assumed that the fluid storage space in the fractal dual- medium is at the node with the unit volume of $V_s$, each node has the same volume, the radial distance between the node and the wellbore is x, and the density of the node is $N\left(x\right)$.

(1)$N(x)=\alpha \cdot x^{D-1}$

where $\alpha$—A parameter related to fracture porosity; D—Fractal dimension.

There are many calculation methods for the fractal dimension, among which the most commonly used is the box-counting dimension method [25]. This method obtains the box-counting dimension by processing the image and then using the least-squares fitting function. Therefore, the fractal dimension of the fracture network can be obtained by the box-counting dimension method after acquiring the pattern of the fracture network based on the inversion of micro-seismic data [24].

The fracture fractal porosity is defined as ${\varphi }_{2f}$. G represents the corresponding symmetry relationship. The equation at $\left(x+dx\right)$ is obtained as follows:

(2)${\varphi }_{2f}=\frac{V_{\varphi }}{V}=\frac{V_s\cdot N(x)dx}{G{x}^{d-1}dx}$

Substituting Equation (1) into Equation (2), we can get:

(3)${\varphi }_{2f}=\frac{V_s\cdot \alpha \cdot x^{D-1}dx}{G{x}^{d-1}dx}$

Rearranging the formula, we can obtain:

(4)${\varphi }_{2f}=\frac{V_s\cdot \alpha }{G}{x}^{D-d}$

where $d$—dimension of Euclidean space embedded by the fractal.

When $x=x_w$, the equation is given by:

(5)${\varphi }_{2fw}=\frac{V_s\cdot \alpha }{G}{x}_w{}^{D-d}$

Combining Equations (4) and (5), the fracture fractal porosity is obtained as follows (d equals 2 in two-dimensional space):

(6)${\varphi }_{2f}={\varphi }_{2fw}{\left(\frac{x}{x_w}\right)}^{D-2}$

Considering the pressure sensitivity of the fracture system, its permeability changes with the effective confining pressure and has an exponential relationship with its effective overlying pressure. Referring to the definition of fractal porosity and the abnormal diffusion coefficient, the fracture fractal permeability considering stress sensitivity can be obtained as:

(7)$k_{2f}=k_{2fw}{\left(\frac{x}{x_w}\right)}^{D-2-\theta }{e}^{-\gamma (p_0-p_{2f})}$

The non-uniform fracture network structure model shows that the closer to the primary fracture, the denser the distribution of secondary fractures. The boundary of the hydraulic fracture is taken as the reference point in this paper (the location with the largest fracture density, i.e., $x_w=w/2$), and the value of fractal dimension D is generally 0–2 [22].

2.2.2. Definition of Dimensionless Parameters

The dimensionless seepage mathematical model is obtained by dimensionless treatment of the relevant parameters.

Dimensionless pressure:

(8)$p_{iD}=\frac{p_0-p_i}{p_0-p_w},i=1,2,3$

Dimensionless time:

(9)$t_D=\frac{k_{2f0}t}{{\left(\varphi C_t\right)}_2{\mu }_2{x}_f^2}$

In above equation, ${\left(\varphi C_t\right)}_2={\left(\varphi C_t\right)}_{2f}+{\left(\varphi C_t\right)}_{2m}$.

Dimensionless threshold pressure gradient:

(10)$G_{3D}=C_{l2f}{x}_f{G}_3$

(11)$g_{3D}=-\frac{x_f{G}_3}{p_0-p_w}$

Dimensionless productivity:

(12)$q_D=\frac{{\mu }_{2}B}{k_{2f0}h(p_0-p_w)}q$

Dimensionless distance:

(13)$x_D=\frac{x}{x_f}$

(14)$y_D=\frac{y}{x_f}$

(15)$w_{fD}=\frac{w_f}{x_f}$

Dimensionless reservoir pressure sensitivity coefficient:

(16)${\gamma }_D=\gamma \left(p_0-p_w\right)$

2.2.3. Establishment of the Dimensionless Seepage Mathematical Model

(1) Mathematical Model for Region III

Region III is the unstimulated volume area at the reservoir boundary. The fluid in this area flows into the volume area between primary fractures (region II) in the x-direction, which is a one-dimensional linear flow. Considering the influence of the threshold pressure gradient, the continuity equation is obtained as:

(17)$\frac{k_3}{{\mu }_3}\left(\frac{{\partial }^2{p}_3}{\partial x^2}-C_{l3}{G}_3\frac{\partial p_3}{\partial x}\right)={\varphi }_3{C}_{t3}\frac{\partial p_3}{\partial t}$

Incorporating the boundary conditions and initial conditions, we can obtain the dimensionless seepage mathematical model for region III after dimensionless transformation:

(18)$\{\begin{array}{l}\frac{{\partial }^2{p}_{3D}}{\partial x_D^2}-G_{3D}\frac{C_{l3}}{C_{l2f}}\frac{\partial p_{3D}}{\partial x_D}=\frac{{\varphi }_3{\mu }_3{C}_{t3}{k}_{2f0}}{{\left(\varphi C_t\right)}_2{\mu }_2{k}_3}\frac{\partial p_{3D}}{\partial t_D}\\ p_{3D}|{}_{x_D=1}=p_{2fD}|{}_{x_D=1}\\ \frac{\partial p_{3D}}{\partial x_D}|{}_{x_D=x_{eD}}=0\\ p_{3D}|{}_{t_D=0}=0\end{array}$

(2) Mathematical Model for Region II

Region II is the volume area between primary fractures. The two kinds of media in this area are the matrix and the fracture network, and the fluid seepage law under the complex fracture network is very complicated. Fluid exchange occurs due to the pressure difference between the matrix and the fracture. A portion of the fluid in this area flows from the matrix to the fracture network, which is regarded as a stable process. The mass of fluid discharged from the matrix to the fracture per unit time is described as follows:

(19)$q^{\prime }=\frac{\alpha {\rho }_0{k}_{2m}}{{\mu }_2}\left(p_{2m}-p_{2f}\right)$

where $q^{\prime }$—Mass of fluid flowing out in per unit rock volume and per unit time, $\mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$, ${\rho }_0$—Density of crude oil, $\mathrm{kg}/{\mathrm{m}}^3$.

Therefore, the governing equation for the matrix system is given by:

(20)${\varphi }_{2m}{C}_{t2m}\frac{\partial p_{2m}}{\partial t}+\frac{\alpha k_{2m}}{{\mu }_2}\left(p_{2m}-p_{2f}\right)=0$

Secondary fractures are unevenly distributed in the y-direction in this area. The farther away from the primary fracture, the smaller the density of secondary fractures, thus leading to the change in fluid flow capacity in region II in the y-axes. Considering the stress sensitivity, fractal theory based on the self-similarity and the scale invariance of the fracture network are introduced to obtain the expression of the porosity and permeability of the complex fracture network. The change in fluid flow capacity is characterized by permeability. With the exception that part of the fluid flows from the matrix system to the fracture system, the fluid in region III flows into region II. Therefore, the governing equation of the fluid in the fracture system is given by:

(21)$\begin{array}{l}{\varphi }_{2f0}{C}_{l2f}{\mu }_2{\left(\frac{y}{w_f/2}\right)}^{D-2}\frac{\partial p_{2f}}{\partial t}\\ -k_{2f0}{e}^{-\gamma \left(p_0-p_{2f}\right)}{\left(\frac{y}{w_f/2}\right)}^{D-2-\theta }\left(\frac{{\partial }^2{p}_{2f}}{\partial y^2}+\left(C_{l2f}+\gamma \right){\left(\frac{\partial p_{2f}}{\partial y}\right)}^2+\frac{D-2-\theta }{y}\frac{\partial p_{2f}}{\partial y}\right)\\ -\alpha k_{2m}\left(p_{2m}-p_{2f}\right)-\frac{{\mu }_2{k}_3}{{\mu }_3{x}_f}\left(\frac{\partial p_3}{\partial x}-G_3\right)|{}_{x=x_f}=0\end{array}$

Combining the boundary conditions and initial conditions, the dimensionless seepage mathematical model for region II can be obtained after dimensionless transformation:

(22)$\{\begin{array}{l}\left(1-{\omega }_2\right)\frac{\partial p_{2mD}}{\partial t_D}+{\lambda }_2\left(p_{2mD}-p_{2fD}\right)=0\\ {\lambda }_2(p_{2mD}-p_{2fD})-\frac{C_{l2f}{\varphi }_{2f0}}{{\left(\varphi C_t\right)}_2}{\left(\frac{y_D}{w_{fD}/2}\right)}^{D-2}\frac{\partial p_{2fD}}{\partial t_D}\\ -e^{-{\gamma }_D{p}_{2fD}}{\left(\frac{y_D}{w_{fD}/2}\right)}^{D-2-\theta }\left[-\frac{D-2-\theta }{y_D}\frac{\partial p_{2fD}}{\partial y_D}+{\gamma }_D{\left(\frac{\partial p_{2fD}}{\partial y_D}\right)}^2-\frac{{\partial }^2{p}_{2fD}}{\partial y_D{}^2}\right]+\frac{k_3{\mu }_2}{k_{2f0}{\mu }_3}(\frac{\partial p_{3D}}{\partial x_D}-g_{3D})|{}_{x_D=1}=0\\ p_{2fD}|{}_{y_D=w_{fD}/2}=p_{1D}|{}_{y_D=w_{fD}/2}\\ e^{-{\gamma }_D{p}_{2fD}}\frac{\partial p_{2fD}}{\partial y_D}|{}_{y_D=y_{eD}}=0\\ p_{2fD}|{}_{t_D=0}=p_{2mD}|{}_{t_D=0}=0\end{array}$

In the above mathematical model, ${\omega }_2$ is the elastic storage capacity ratio, and its expression is as follows:

(23)${\omega }_2=\frac{{\left(\varphi C_t\right)}_{2f}}{{\left(\varphi C_t\right)}_{2f}+{\left(\varphi C_t\right)}_{2m}}$

${\lambda }_2$ is the dimensionless interporosity flow coefficient, and its expression is given by:

(24)${\lambda }_2=\alpha \frac{k_{2m}}{k_{2f0}}{x}_f^2$

(3) Mathematical Model for Region I

Region I is the area of the primary fracture. The fluid in this area flows into the horizontal wellbore along the fracture, which is a one-dimensional linear flow. Considering the influence of stress sensitivity, the permeability of the main fracture changes with the effective confining pressure, and has an exponential function relationship with the overlying pressure. We can acquire the expression of the permeability of the main fracture as:

(25)$k_1=k_{10}{e}^{-\gamma (p_0-p_1)}$

The tip of the fracture is closed, and the fluid in region II flows into region I. Without considering the skin effect of the primary fracture and the wellbore storage effect, the governing equation of fluid in this area can be expressed as:

(26)$\frac{k_{10}}{{\mu }_1}{e}^{-\gamma (p_0-p_1)}\left[\frac{{\partial }^2{p}_1}{\partial x^2}+\left(C_{l1}+\gamma \right){\left(\frac{\partial p_1}{\partial x}\right)}^2\right]+\frac{2k_{2f0}{e}^{-\gamma \left(p_0-p_{2f}\right)}}{{\mu }_2{w}_f}\frac{\partial p_{2f}}{\partial y}|{}_{y=\frac{w_f}{2}}={\varphi }_1{C}_{t1}\frac{\partial p_1}{\partial t}$

Incorporating the boundary and initial conditions, the dimensionless seepage mathematical model for region I can be obtained after dimensionless treatment:

(27)$\{\begin{array}{l}{e}^{-{\gamma }_D{p}_{1D}}\left[\frac{{\partial }^2{p}_{1D}}{\partial x_D^2}-{\gamma }_D{\left(\frac{\partial p_{1D}}{\partial x_D}\right)}^2\right]+\frac{2{\mu }_1{k}_{2f0}}{{\mu }_2{k}_1{w}_{fD}}{e}^{-{\gamma }_D{p}_{2fD}}\frac{\partial p_{2fD}}{\partial y_D}|{}_{y_D=\frac{w_{fD}}{2}}=\frac{{\varphi }_1{\mu }_1{C}_{t1}{k}_{2f0}}{{\left(\varphi C_t\right)}_2{\mu }_2{k}_{10}}\frac{\partial p_{1D}}{\partial t_D}\\ p_{1D}|{}_{x_D=0}=1\\ e^{-{\gamma }_D{p}_{1D}}\frac{\partial p_{1D}}{\partial x_D}|{}_{x_D=1}=0\\ p_{1D}|{}_{t_D=0}=0\end{array}$

3. Results

3.1. Solution of the Mathematical Model

The mathematical equations for region I and region II have strong nonlinear characteristics due to consideration of the influence of the stress sensitivity of the fracture system. The transformation method proposed by Pedrosa [26] can eliminate the nonlinearity of the equation, and its transformation formula is given by:

(28)$p_{iD}=-\frac{1}{{\gamma }_D}\ln (1-{\gamma }_D{\zeta }_{iD}),i=1,2$

The value of ${\gamma }_D$ is small in general. According to the perturbation theory, the following equations are defined as follows:

(29)${\zeta }_{iD}={\zeta }_{iD0}+{\gamma }_D{\zeta }_{iD1}+{\gamma }_D^2{\zeta }_{iD2}+\cdots ,i=1,2$

(30)$\frac{1}{1-{\gamma }_D{\zeta }_{iD}}=1+{\gamma }_D{\zeta }_{iD}+{\gamma }_D^2{\zeta }_{{}_{iD}}^2+\cdots ,i=1,2$

(31)$-\frac{1}{{\gamma }_D}\ln (1-{\gamma }_D{\zeta }_{iD})={\zeta }_{iD}+\frac{1}{2}{\gamma }_D{\zeta }_{iD}^2+\cdots ,i=1,2$

Considering that ${\gamma }_D$ is very small, zero-order perturbation can provide the required accuracy of the solution. The above method is employed to solve the seepage mathematical model for area I and area II.

3.1.1. Solution of the Seepage Mathematical Model for Region III

Laplace transformation is applied to the dimensionless seepage mathematical model for region III to obtain the corresponding mathematical model in Laplace space:

(32)$\{\begin{array}{l}\frac{{\partial }^2{\overline{p}}_{3D}}{\partial x_D{}^2}-b_3\frac{\partial {\overline{p}}_{3D}}{\partial x_D}={\eta }_{3}s{\overline{p}}_{3D}\\ {\overline{p}}_{3D}|{}_{x_D=1}={\overline{p}}_{2fD}|{}_{x_D=1}\\ \frac{\partial {\overline{p}}_{3D}}{\partial x_D}|{}_{x_D=x_{eD}}=0\\ {\overline{p}}_{3D}|{}_{t_D=0}=0\end{array}$

In the above equation, $b_3=G_{3D}\frac{C_{l3}}{C_{l2f}}$, ${\eta }_3=\frac{C_{t3}{\varphi }_3{k}_{2f0}{\mu }_3}{{(\varphi C_t)}_2{k}_3{\mu }_2}$.

Solving the mathematical model, we can get:

(33)$\frac{\partial {\overline{p}}_{3D}}{\partial x_D}|{}_{x_D=1}=F_{32}{\overline{p}}_{2fD}|{}_{x_D=1}$

where $F_{32}=\frac{r_{31}{r}_{32}\left[e^{r_{31}}-e^{(r_{31}-r_{32})x_{eD}}{e}^{r_{32}}\right]}{r_{32}{e}^{r_{31}}-r_{31}{e}^{(r_{31}-r_{32})x_{eD}}{e}^{r_{32}}}$, $r_{31}=\frac{b_3-\sqrt{b_3^2+4s{\eta }_3}}{2}$, $r_{32}=\frac{b_3+\sqrt{b_3^2+4s{\eta }_3}}{2}$.

3.1.2. Solution of the Seepage Mathematical Model for Region II

By performing zero-order perturbation transformation and Laplace transformation on the dimensionless seepage mathematical model for region II, the mathematical model in Laplace space is obtained as follows:

(34)$\{\begin{array}{l}{\lambda }_2\left[\frac{{\lambda }_2}{s(1-{\omega }_2)+{\lambda }_2}-1\right]{\overline{\zeta }}_{2fD0}-\frac{C_{l2f}{\varphi }_{2f0}}{{\left(\varphi C_t\right)}_2}{\left(\frac{y_D}{w_{fD}/2}\right)}^{D-2}s{\overline{\zeta }}_{2fD0}\\ +{\left(\frac{y_D}{w_{fD}/2}\right)}^{D-2-\theta }\left[\frac{D-2-\theta }{y_D}\frac{\partial {\overline{\zeta }}_{2fD0}}{\partial y_D}+\frac{{\partial }^2{\overline{\zeta }}_{2fD0}}{\partial y_D{}^2}\right]+a_2(\frac{\partial {\overline{p}}_{3D}}{\partial x_D}-\frac{g_{3D}}{s})|{}_{x_D=1}=0\\ {\overline{\zeta }}_{2fD0}|{}_{y_D=w_{fD}/2}={\overline{\zeta }}_{1D0}|{}_{y_D=w_{fD}/2}\\ \frac{\partial {\overline{\zeta }}_{2fD0}}{\partial y_D}|{}_{y_D=y_{eD}}=0\\ {\overline{\zeta }}_{2fD0}|{}_{t_D=0}=0\end{array}$

Substituting Equation (33) into Equation (34) and simplifying the governing equations, we can acquire:

(35)$\frac{{\partial }^2{\overline{\zeta }}_{2fD0}}{\partial y_D{}^2}+\frac{D-2-\theta }{y_D}\frac{\partial {\overline{\zeta }}_{2fD0}}{\partial y_D}-C_1{C}_2\left({\overline{\zeta }}_{2fD0}+a_2\frac{g_{3D}}{s{C}_1}\right)y_D^{\theta }=0$

where $C_1=W_2-a_2{F}_{32}+{\eta }_{2}s,C_2={\left(\frac{2}{w_{fD}}\right)}^{\theta }$, $a_2=\frac{k_3{\mu }_2}{k_{2f0}{\mu }_3}$, $W_2=\frac{{\lambda }_{2}s(1-{\omega }_2)}{s(1-{\omega }_2)+{\lambda }_2},{\eta }_2=\frac{C_{l2f}{\varphi }_{2f0}}{{\left(\varphi C_t\right)}_2}$.

Combined with the boundary conditions, the model is solved according to the Bessel function, and the result is obtained as follows:

(36)$\frac{\partial {\overline{\zeta }}_{2fD0}}{\partial y_D}|{}_{y_D=\frac{w_{fD}}{2}}=F_{21}{\overline{\zeta }}_{1D0}|{}_{y_D=\frac{w_{fD}}{2}}+Z_2$

In above equation,

$F_{21}=\frac{\frac{(\theta +2)}{2}{E}_2{\left(\frac{w_{fD}}{2}\right)}^{\frac{\theta }{2}}\left\{I_{\frac{1-D}{2+\theta }}\left[E_2{\left(\frac{w_{fD}}{2}\right)}^{\frac{2+\theta }{2}}\right]-\frac{h_1}{h_2}{K}_{\frac{1-D}{2+\theta }}\left[E_2{\left(\frac{w_{fD}}{2}\right)}^{\frac{2+\theta }{2}}\right]\right\}}{I_{\frac{3+\theta -D}{2+\theta }}\left[E_2{\left(\frac{w_{fD}}{2}\right)}^{\frac{2+\theta }{2}}\right]+\frac{h_1}{h_2}{K}_{\frac{3+\theta -D}{2+\theta }}\left[E_2{\left(\frac{w_{fD}}{2}\right)}^{\frac{2+\theta }{2}}\right]}$

$Z_2=\frac{a_2{g}_{3D}}{s{C}_1}{F}_{21}$, $h_1=I_{\frac{1-D}{2+\theta }}\left(E_2{y}_{eD}{}^{\frac{2+\theta }{2}}\right)$, $h_2=K_{\frac{1-D}{2+\theta }}\left(E_2{y}_{eD}{}^{\frac{2+\theta }{2}}\right)$, $E_2=\frac{2\sqrt{C_1{C}_2}}{2+\theta }$.

3.1.3. Solution of the Seepage Mathematical Model for Region I

Laplace transformation is applied to the dimensionless seepage mathematical model for region I to obtain the corresponding mathematical model in Laplace space:

(37)$\{\begin{array}{l}\frac{{\partial }^2{\overline{\zeta }}_{1D0}}{\partial x_D^2}+a_1\frac{\partial {\overline{\zeta }}_{2fD0}}{\partial y_D}|{}_{y_D=\frac{w_{fD}}{2}}={\eta }_{1}s{\overline{\zeta }}_{1D0}\\ {\overline{\zeta }}_{1D0}|{}_{x_D=0}=\frac{1-e^{-{\gamma }_D}}{s{\gamma }_D}\\ \frac{\partial {\overline{\zeta }}_{1D0}}{\partial x_D}|{}_{x_D=1}=0\\ {\overline{\zeta }}_{1D0}|{}_{t_D=0}=0\end{array}$

where $a_1=\frac{2{\mu }_1{k}_{2f0}}{{\mu }_2{k}_1{w}_{fD}}$, ${\eta }_1=\frac{{\varphi }_1{\mu }_1{C}_{t1}{k}_{2f0}}{{\left(\varphi C_t\right)}_2{\mu }_2{k}_{10}}$.

Substituting Formula (36) into the above equation and solving the mathematical model, we can get:

(38)$\frac{\partial {\overline{\zeta }}_{1D0}}{\partial x_D}|{}_{x_D=0}=-\frac{\left(r_{11}{r}_{12}{e}^{r_{11}-r_{12}}-r_{11}{r}_{12}\right)}{s{\gamma }_D(e^{r_{11}-r_{12}}{r}_{11}-r_{12})}(e^{-{\gamma }_D}-1+d_{1}s{\gamma }_D)$

In above equation, $r_{11}=\sqrt{s{\eta }_1-a_1{F}_{21}},r_{12}=-\sqrt{s{\eta }_1-a_1{F}_{21}},d_1=\frac{a_1{Z}_2}{s{\eta }_1-a_1{F}_{21}}$.

The productivity at the inner boundary of the primary fracture is a quarter of the production of the bottom hole, and its relational expression is as follows:

(39)$\frac{\partial {\overline{\varsigma }}_{1D0}}{\partial x_D}|{}_{x_D=0}=-\frac{k_{2f0}{\mu }_1}{2k_{10}{w}_{fD}{\mu }_2}{\overline{q}}_D$

Substituting Equation (38) into the above equation, the production of a single fracture of a horizontal well under Laplace space is obtained as follows:

(40)${\overline{q}}_D=a_0\frac{\left(r_{11}{r}_{12}{e}^{r_{11}-r_{12}}-r_{11}{r}_{12}\right)}{s{\gamma }_D(e^{r_{11}-r_{12}}{r}_{11}-r_{12})}(e^{-{\gamma }_D}-1+d_{1}s{\gamma }_D)$

where $a_0=\frac{2k_{10}{w}_{fD}{\mu }_2}{k_{2f0}{\mu }_1}$.

3.1.4. Solution of the Production of a Single Primary Fracture for Horizontal Wells

We obtained the production expression of a single fracture of the horizontal well in Laplace space in Section 3.1.3. The production expression in real space can be obtained by the Stehfest [27] numerical inversion according to the productivity formula in Laplace space. Combining with Equation (40), the corresponding dimensionless production expression of a single fracture at the bottom hole for any dimensionless time point $t_D$ is represented as follows:

(41)$q_D(t_D)=\frac{\ln 2}{t_D}{\displaystyle \sum _{i=1}^N{V}_i{\overline{q}}_D(s)}$

In the above expression, the Laplace variable $s=\frac{\ln 2}{t_D}i$.

$V_i={(-1)}^{\frac{N}{2}+i}{\displaystyle \sum _{k=[\frac{i+1}{2}]}^{\min (i,\frac{N}{2})}\frac{k^{\frac{N}{2}}(2k)!}{(\frac{N}{2}-k)!k!(k-1)!(i-k)!(2k-i)!}}$

Therefore, the dimensionless production expression at the bottom hole of a single fracture in real space is given by:

(42)$q_D(t_D)=\frac{\ln 2}{t_D}{\displaystyle \sum _{i=1}^N\left\{\begin{array}{l}{\left(-1\right)}^{\frac{N}{2}+i}{\displaystyle \sum _{k=[\frac{i+1}{2}]}^{\min (i,\frac{N}{2})}\frac{k^{\frac{N}{2}}(2k)!}{(\frac{N}{2}-k)!k!(k-1)!(i-k)!(2k-i)!}\times }\\ a_0\frac{\left(r_{11}{r}_{12}{e}^{r_{11}-r_{12}}-r_{11}{r}_{12}\right)}{s{\gamma }_D(e^{r_{11}-r_{12}}{r}_{11}-r_{12})}(e^{-{\gamma }_D}-1+d_{1}s{\gamma }_D)\end{array}\right\}}$

3.1.5. Solution of the Production of Staged Multi-Cluster Fractured Horizontal Wells with Non-Uniformly Distributed Fractures

Taking into account the current situation of tight oil exploitation with bridge plug staged multi-cluster perforation fracturing [28], a productivity prediction model for staged multi-cluster fractured horizontal wells was established in this study, and the fracture distribution diagram is represented in the following Figure 3.

As shown in Figure 3, fractures of the horizontal well are unevenly arranged. The number of fractures in each fracturing section, the distance between fracturing sections, and the space between each fracture of fracturing sections are different. Taking the fracture at any position as the central axis, its control range is the sum of half the distance from the left and right adjacent fractures. (In particular, the range that can be controlled by fractures at both ends of the horizontal well is composed of two parts: the distance from the fractures at both ends of the horizontal well to the reservoir boundary and 1/2 of the distance between the fractures and their adjacent fractures.) We divide the fractures into three types based on their position at the horizontal well: the fractures at both ends of the horizontal well, the fractures at both ends of the fracturing segment, and the fractures inside the fracturing segment. The shape parameters of these are defined as follows.

The fractures at both ends of the horizontal well ($\Delta y_d$ is the distance between the fractures at both ends of the horizontal well and the reservoir boundary):

(43)${\lambda }_d=W/\Delta y_d$

The fractures at both ends of the fracturing segment ($\Delta y_{ni}$ is the distance between the last fracture of the ith fracturing segment and the first fracture of the (i + 1)th fracturing segment, i represents the sequence of the fracturing segment):

(44)${\lambda }_i=W/\Delta y_{\mathrm{n}i},i=1,2,3\cdots$

The fractures inside the fracturing segment ($\Delta y_{mij}$ is the distance between the jth fracture inside the ith fracturing segment and the (j + 1)th fracture inside the same fracturing segment, j is the sequence of the fracture):

(45)${\lambda }_j=W/\Delta y_{mij},j=1,2,3\cdots$

According to the fracture production superposition calculation method proposed by Mayer [29], the single well production formula of staged multi-cluster fractured horizontal wells can be obtained as follows:

(46)$Q_D(t_D)={\displaystyle \sum _{i=1}^{ns}{\displaystyle \sum _{j=1}^{nc-1}{q}_D({\lambda }_j,t_D)}}+{\displaystyle \sum _{i=1}^{ns-1}{q}_D({\lambda }_i,t_D)}+q_D({\lambda }_d,t_D)$

In the above equation, ns expresses the number of fracturing segments and nc means the number of fractures in the fracturing segment. The total number of fracturing clusters in all fracturing sections is the total number of fractures.

3.2. Description of Symbols

$p_0$ is initial formation pressure, ${\mathrm{P}}_{\mathrm{a}}$; $p_w$ is bottom hole pressure, ${\mathrm{P}}_{\mathrm{a}}$; $p$ is current fluid flow pressure, ${\mathrm{P}}_{\mathrm{a}}$; $N_F$ is the number of fracture; $x_f$ is the half-length of fracture, $\mathrm{m}$; $w_f$ is the width of fracture, $\mathrm{m}$; $x_e$ is the half-width of formation, $\mathrm{m}$; $y_e$ is width of region II, $\mathrm{m}$; $B$ is oil volume factor; h is formation thickness, $\mathrm{m}$; $\mu$ is oil viscosity, ${\mathrm{P}}_{\mathrm{a}}·\mathrm{s}; k$ is permeability, ${\mathrm{m}}^2$; $\varphi$ is porosity, $\%$; $G$ is threshold pressure gradient, ${\mathrm{P}}_{\mathrm{a}}/\mathrm{m}$; $\gamma$ is stress sensitivity coefficient, ${\mathrm{P}}_{\mathrm{a}}^{-1}$; $D$ is fractal dimension; $\theta$ is anomalous diffusion coefficient; $\alpha$ is shape factor of reservoir, ${\mathrm{m}}^{-2}$; $C_t$ is total compression coefficient, ${\mathrm{P}}_{\mathrm{a}}^{-1}$; $C_l$ is liquid phase compression coefficient, ${\mathrm{P}}_{\mathrm{a}}^{-1}$; $q$ is the production of a single fracture of horizontal wells, ${\mathrm{m}}^3/\mathrm{s}$; $Q$ is the productivity of a single staged multi-cluster fractured horizontal well, ${\mathrm{m}}^3/\mathrm{s}$; $t$ is the production time of horizontal well, $\mathrm{s}$; $I,K$ represent the Bessel function; the subscripts 1, 2, and 3, respectively, represent the regions I, II, and III; the subscripts m, f respectively represent the matrix system and fracture system, respectively; the subscript 0 represents the initial value; subscript D means dimensionless; the superscript – means Laplace space; $s$ is the Laplace variable.

4. Discussion

4.1. Validation of the Productivity Model

4.1.1. Verification of the Analytical Solution

Without considering the non-uniform structure characteristics of the fracture network, the influence of threshold pressure gradient on fluid flow, and the stress sensitivity of the fracture system, our simplified model is similar to the classical trilinear flow model proposed by M. Brown et al. [17]. According to the expression of the dimensionless bottom hole pressure of the classical model, the dimensionless production at the bottom hole of a single fracture can be obtained as [30]:

(47)${\overline{q}}_D=\frac{1}{s^2{\overline{p}}_{wfD}}$

Setting $G_3=0,\gamma =0,D=2,\theta =0$ in the original mathematical model, the results of the simplified model were compared with the results of the classical model to verify the accuracy of the analytical solution of our original model. The simulation parameters are listed in Table 1. As shown in Figure 4, the results of the simplified model are highly consistent with the results of the classical model, which proves the accuracy of the analytical solution of our original model.

4.1.2. An Example in an Oil Field

We can calculate the production of a single fractured horizontal well according to Equation (46). To verify the reliability of our model, the basic parameters of an oil well in the X Oilfield (in China) are substituted into the mathematical model for calculation, and the calculation results of the model fitted well with the production data of the well. The basic parameters and fracturing design parameters of the oil well, respectively, are shown in Table 2 and Table 3.

The fitting results of the model and production data are shown in Figure 5. It can be seen from the figure that the calculation results of this model are basically consistent with the actual production data of the oil field within the allowable error range, which confirms the reliability of the model. The calculation results of the model deviate from the actual production data without considering the influence of stress sensitivity or the threshold pressure gradient. In addition, the fitting effect is obviously inferior to the original model when simultaneously ignoring both stress sensitivity and the threshold pressure gradient, which further proves the accuracy of the original model.

4.2. Analysis of Factors Affecting Productivity of Fractured Horizontal Wells

To facilitate the analysis of the factors affecting productivity, it is assumed that the fractures are evenly distributed along the horizontal well, with a total of 36 fracturing segments, four fractures within each segment, 28 m between each fracturing segment, and 7 m between each fracture inside the fracturing segment. Equation (46) can be rewritten as follows:

(48)$Q_D(t_D)=ns(nc-1)q_D({\lambda }_j,t_D)+(ns-1)q_D({\lambda }_i,t_D)+q_D({\lambda }_d,t_D)$

Based on Equation (48) and the reference data in Table 2, the factors influencing the productivity of horizontal wells were investigated according to the deliverability evaluation method proposed in this paper, mainly including fracturing parameters (fracture half-length, fracture number) and reservoir characteristic parameters (including matrix permeability, threshold pressure gradient, and stress sensitive coefficient).

4.2.1. The Effect of Fracture Number on Productivity

The number of fracturing sections was set as 30, 32, 34, and 36, with other parameters in Table 2 unchanged. The impact of the number of fractures on the productivity of the horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as shown in Figure 6.

As depicted in Figure 6, the production of horizontal wells rises with the increase in the number of fractures. Fractures can improve the fluid seepage environment. The number of fractures increases, leading to the increase in fluid seepage channels and a significant improvement in the seepage environment, which results in an increase in the production of the horizontal well. At the same time, the phenomenon is reflected in Figure 6 that the production growth of the horizontal wells slows with the increase in the number of fractures. There is mutual interference between hydraulic fractures, and the increase in the number of fractures leads to the stronger interference between the fractures, thus resulting in the decrease in the production growth of horizontal wells. In addition, the increase in the number of fractures leads to the growth in the fracturing construction costs. Therefore, the number of fractures should be controlled within a reasonable range during the process of fracturing design to ensure the positive growth in production benefits.

4.2.2. The Effect of Fracture Half-Length on Productivity

The value of $x_f$ was assumed to be 60, 90, 120, and 150 m, with other parameters in Table 2 unchanged. The impact of the half-length of fractures on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal wells were obtained, as shown in Figure 7.

As depicted in Figure 7a, the productivity of horizontal wells rises with the increase in fracture length. At the initial stage during the exploitation of tight oil reservoirs, the half-length of fractures increases by 30 m, and the daily production of the horizontal well rises by about 10–50 m³. As presented in Figure 7b, the growth rate of final cumulative production decreases obviously with the increase in fracture length, and the fracturing cost expands with the increase in the hydraulic fracture length. Therefore, there is an optimal fracture half-length, which can not only improve the production of horizontal wells, but also save expenditure.

4.2.3. The Effect of Matrix Permeability on Productivity

The value of matrix permeability was set as 0.13, 0.15, 0.17, and 0.19 $\mathrm{mD}$, with other parameters in Table 2 unchanged. The influence of matrix permeability on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were acquired, as depicted in Figure 8.

It can be seen from Figure 8a that, in the early and middle stages of exploitation, the daily productivity of the horizontal well rises with the increase in matrix permeability, and the higher the matrix permeability, the bigger the deceleration rate of the daily production of horizontal wells. As depicted in Figure 8b, the trend of the rising cumulative productivity slows with the increase in matrix permeability. After the horizontal well has produced for a certain time, the cumulative productivity remains largely unchanged with the increase in matrix permeability, which indicates that increasing matrix permeability cannot effectively improve the final cumulative productivity of the horizontal well, but can shorten the time to reach the same cumulative productivity. In general, matrix permeability has a significant impact on the early stage and, in particular, the middle stage, of the whole producing process.

4.2.4. The Effect of the Stress Sensitivity Coefficient on Productivity

The value of the stress sensitivity coefficient was assumed to be 0.00, 0.01, 0.03, and 0.05 ${\mathrm{MP}}_{\mathrm{a}}^{-1}$, with other parameters in Table 2 unchanged. The effect of the stress sensitivity coefficient on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as represented in Figure 9.

As presented in Figure 9, the stress sensitivity coefficient has a significant impact on the productivity of horizontal wells. The production of the horizontal well decreases with the increase in the stress sensitivity coefficient, and the trend of the reduction in productivity slows. The fracture system is compressed because of the increase in the effective pressure in the process of development and exploration of the tight oil reservoir, thus resulting in the decrease in permeability, which further reduces the productivity of horizontal wells. As depicted in the figures, the productivity of the horizontal well when considering the stress sensitivity of the fracture system is significantly lower than that without stress sensitivity, which shows that stress sensitivity is an important factor that cannot be ignored in the deliverability evaluation of fractured horizontal wells in tight reservoirs.

4.2.5. The Effect of Threshold Pressure Gradient on Productivity

$G_3$ was assumed to be 0.00, 0.02, 0.04, and 0.06${\mathrm{MP}}_{\mathrm{a}}/\mathrm{m}$, and other parameters in Table 2 were unchanged. The influence of the threshold pressure gradient on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were obtained, as presented in Figure 10. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 10a.

It can be seen from Figure 10 that the threshold pressure gradient has an obvious influence on the productivity of horizontal wells. The greater the threshold pressure gradient, the bigger the reduction rate of the daily productivity of horizontal wells and the smaller the final cumulative productivity. This is because a larger threshold pressure gradient leads to stronger seepage resistance, and difficulty in obtaining the fluid flow results for the reduction in productivity. The threshold pressure gradient was used to characterize the non-Darcy seepage of fluid in this study. As shown in Figure 10, the threshold pressure gradient is a non-negligible factor leading to the reduction in productivity of horizontal wells. Therefore, the fluid seepage process cannot be simply viewed as a flow following Darcy’s law in the process of the deliverability evaluation for staged multi-cluster fractured horizontal wells in tight reservoirs.

4.2.6. The Effect of the Fractal Dimension on Productivity

D was set as 1.70, 1.72, 1.74, and 1.76, and other parameters in Table 2 were unchanged. The impact of the fractal dimension on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as shown in Figure 11. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 11a.

In this study, the complex fracture network considering the fractal characteristics is called the non-uniform fracture network. As presented in Figure 11, the fractal dimension has a significant influence during the whole process of exploration of horizontal wells. The productivity of horizontal wells rises with the increase in the fractal dimension of the non-uniform fracture network, and the trend of rising productivity is also more obvious.

5. Conclusions

This paper proposes a new method of deliverability evaluation for staged fractured horizontal wells. The solution of the productivity model is obtained using the semi-analytical method, and the reliability of the model was verified with a classic model and an oilfield example. According to the research into the deliverability evaluation method and factors affecting productivity, several conclusions of this paper are represented as follows.

• (1). The non-uniform distribution fracture network structure along the horizontal wellbore is proposed in this paper. Combining fractal theory with stress sensitivity to characterize the permeability and porosity, a new three-zone composite unsteady deliverability evaluation model for staged multi-cluster fractured horizontal wells is established.

• (2). The number of fractures, fracture half-length, and matrix permeability can affect the productivity of horizontal wells to a certain extent. The oil production rises with the increase in fracture number, fracture half-length, and matrix permeability, but the growth rate of productivity decreases at the same time.

• (3). The stress sensitivity coefficient, threshold pressure gradient, and fractal dimension have obvious effects on the productivity of horizontal wells. The smaller the stress sensitivity coefficient and threshold pressure gradient, the greater the oil production. A larger fractal dimension leads to greater productivity of horizontal wells.

• (4). The deliverability evaluation method for staged fractured horizontal wells in tight reservoirs based on the non-uniform fracture network structure can effectively predict the actual production of oil wells. Additionally, it can provide a theoretical basis for the optimization of fracturing parameters of horizontal wells and the formulation of technical policies for rational development of tight reservoirs, which has great significance for realizing the effective development of tight reservoirs.

Although some innovations were achieved in this study, the productivity model does not consider the imbibition phenomenon in the development of tight oil reservoirs and the pressure loss along the horizontal wellbore. Therefore, a more comprehensive deliverability evaluation model will be established by considering both the imbibition effect and the pressure drop effect in our future research work.

Author Contributions

Conceptualization, S.L., S.X. and D.W.; Formal analysis, S.L. and S.X.; Funding acquisition, S.X.; Investigation, S.L.; Methodology, S.L. and R.C.; Project administration, S.X.; Resources, D.W., P.S., Y.H. and S.C.; Supervision, S.X.; Validation, S.L.; Writing—original draft, S.L.; Writing—review & editing, S.L., S.X. and X.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by PetroChina Major Scientific and Technological Projects, China National Petroleum Corporation Science and Technology research project, grant number 2016E-0505, KT2020-13-01.

Data Availability Statement

Restrictions apply to the availability of these data. Data was partly obtained from PetroChina Changqing Oilfield Company and partly generated during the process of calculations. Oilfield data are available from the corresponding author with the permission of PetroChina Changqing Oilfield Company and calculation data can be acquired from the lead author.

Acknowledgments

Authors would like to thank Zhiyuan Wang for solving some problems in the research.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Non-uniform fracture network structure model.

Figure 2. Physical model of seepage flow in fractured horizontal wells in tight oil reservoirs.

Figure 3. Schematic diagram of non-uniform fracture distribution of staged multi-cluster fractured horizontal wells.

Figure 4. Comparison curve of results between our model and the model [17] in the literature.

Figure 5. Model results fitted with production data.

Figure 6. The effect of fracture number on productivity: (a) effect of fracture number on well production rate; (b) effect of fracture number on cumulative production rate.

Figure 7. The effect of fracture half-length on productivity: (a) effect of fracture half-length on well production rate; (b) effect of fracture half-length on cumulative production rate.

Figure 8. The effect of matrix permeability on productivity: (a) effect of matrix permeability on well production rate; (b) effect of matrix permeability on cumulative production rate.

Figure 9. The effect of the stress sensitivity coefficient on productivity: (a) effect of the stress sensitivity coefficient on well production rate; (b) effect of the stress sensitivity coefficient on cumulative production rate.

Figure 10. The effect of the threshold pressure gradient on productivity: (a) effect of the threshold pressure gradient on well production rate; (b) effect of the threshold pressure gradient on cumulative production rate.

Figure 11. The effect of the fractal dimension on productivity: (a) effect of the fractal dimension on well production rate; (b) effect of the fractal dimension on cumulative production rate.

Table 1

Simulation parameters.

Parameters	Value	Parameters	Value
Initial formation pressure, $p_0/{\mathrm{MP}}_{\mathrm{a}}$	15	Oil viscosity, ${\mu }_1,{\mu }_2,{\mu }_3/{\mathrm{mP}}_{\mathrm{a}}·s$	1.27
Bottom hole pressure, $p_w/{\mathrm{MP}}_{\mathrm{a}}$	8	Permeability of HF, $k_{10}/\mathrm{D}$	12
Shape factor, $\alpha \times {10}^{-3}/{\mathrm{m}}^{-2}$	$1.2$	Permeability of fracture network, $k_{2f0}/\mathrm{mD}$	45
HF width, $w_f\times {10}^{-2}/\mathrm{m}$	1	Permeability of matrix, $k_{2m},k_3/\mathrm{mD}$	0.13
Reservoir half width, $x_e/\mathrm{m}$	300	Formation thickness, $h/\mathrm{m}$	8
Porosity of matrix system, ${\varphi }_{2\mathrm{m}},{\varphi }_3/\%$	9.87	Oil volume factor, $B$	1.261
Porosity of fracture system,${\varphi }_1,{\varphi }_{2f0}/\%$	12	HF number, $N_F$	1
Width of region II, $y_e/\mathrm{m}$	7.5	HF half length, $x_f/\mathrm{m}$	250

HF is the abbreviation of hydraulic structure (the same as Table 2 and Table 3).

Table 2

Basic parameters of an oil well in the X Oilfield.

Parameters	Value	Parameters	Value
Initial formation pressure, $p_0/{\mathrm{MP}}_{\mathrm{a}}$	15.8	Oil viscosity, ${\mu }_1,{\mu }_2,{\mu }_3/{\mathrm{mP}}_{\mathrm{a}}·\mathrm{s}$	1.27
Bottom hole pressure, $p_w/{\mathrm{MP}}_{\mathrm{a}}$	8.5	Permeability of HF, $k_{10}/\mathrm{D}$	12
Shape factor, $\alpha \times {10}^{-3}/{\mathrm{m}}^{-2}$	$1.2$	Permeability of fracture network, $k_{2f0}/\mathrm{mD}$	45
HF width, $w_f\times {10}^{-2}/\mathrm{m}$	1	Permeability of matrix, $k_{2\mathrm{m}},k_3/\mathrm{mD}$	0.13
Reservoir half width, $x_e/\mathrm{m}$	300	Formation thickness, $h/\mathrm{m}$	8
Porosity of matrix system, ${\varphi }_{2m},{\varphi }_3/\%$	9.87	Oil volume factor, $B$	1.261
Porosity of fracture system, ${\varphi }_1,{\varphi }_{2f0}/\%$	12	Fractal dimension, $\mathrm{D}$	1.7
Horizontal well length, $L/\mathrm{m}$	1850.5	Anomalous diffusion coefficient, $\theta$	0.1
Stress sensitivity coefficient, $\gamma /{10}^{-2}{\mathrm{MP}}_{\mathrm{a}}^{-1}$	3	Threshold pressure gradient, $G_3/{10}^{-2}{\mathrm{MP}}_{\mathrm{a}}·{\mathrm{m}}^{-1}$	6.67

Table 3

Fracture parameters of an oil well in the X Oilfield.

Fracturing Section Sequence	Cluster Number$\mathit{n}\mathit{c}$	Distance between Section$\Delta {\mathit{y}}_{\mathit{n}\mathit{i}}/\mathit{m}$	Distance between Cluster$\Delta {\mathit{y}}_{\mathit{m}\mathit{i}\mathit{j}}/\mathit{m}$	HF Half-Length${\mathit{x}}_{\mathit{f}}/\mathit{m}$
1	4	0	9,8,5	250
2	5	72.4	7,6.4,7.4,7.4	250
3	6	54.8	7.8,8,6.8,9,13.4	250
4	4	32	6.8,8,8.8	250
5	3	24	8,6	250
6	4	31.4	6.4,10.8,11	250
7	4	27.4	8.4,5.8,17	250
8	5	22.4	12,5.4,6.4,9.4	250
9	4	24.8	6.4,4.4,11.4	250
10	4	19.4	7.4,7.4,5.4	250
11	3	19	6,7	250
12	4	21	6.4,5.8,7.4	250
13	6	25.4	5.8,6.4,8,11.4,9.4	250
14	6	62.4	6.4,8.8,14,7.4,6.4	250
15	5	51.4	5.4,7.4,6.4,9.4	250
16	4	21.4	3.8,5.4,6	250
17	6	41	4.8,9,6.4,7.4,8.8	250
18	5	35.4	6.4,8.4,5.4,11.4	250
19	4	21.4	7,9.4,5.4	250
20	4	20.8	6.4,5.4,5.4	250
21	5	21.4	7.4,6.4,5.4,5.4	250
22	6	44.4	7.4,4.4,5.8,14,7.4	250
23	3	25	6,6	250
24	6	25.4	9.4,5.4,6.8,6,5.4	250
25	4	49.4	3.8,5,5.4	250
26	5	20.4	5.4,6.4,7.4,5.4	250
27	4	24.4	4.4,5,4.4	250
28	4	18.8	5.4,3.8,5	250
29	4	22	5.4,8,4.8	250
30	4	18.8	5.4,4.4,4.4	250
31	4	19.4	10.8,4,7.4	250
32	4	28	6.4,6.4,4.8	250
33	4	19.4	5.4,5.4,5.4	250