Study on a New Transient Productivity Model of Horizontal Well Coupled with Seepage and Wellbore Flow

Abstract

Digital transformation has become one of the major themes of the development of the global oil industry today. With the development of digital transformation, on-site production will surely achieve further automated management, that is, on-site production data automatic collection, real-time tracking, diagnosis and optimization, and remote control of on-site automatic adjustment devices. In this process, the realization of real-time optimization work based on massive data collection needs to be carried out combined with oil and gas well transient simulation. Therefore, research of the horizontal well capacity prediction transient model is one of the important basic works in the work of oil and gas digital transformation. In this paper, the method and process of establihing the transient calculation model of single-phase flow in horizontal wells are introduced in detail from three aspects: reservoir seepage, horizontal wellbore flow (taking one kind of flow as an example), and the coupling model of two flows. The model is more reliable through the verification of pressure recovery data from multiple field logs. The transient model of single-phase seepage in horizontal wells will lay the foundation for the establishment of transient models of oil-gas two-phase seepage and oil-gas-water three-phase seepage.

1. Introduction

Digital transformation has become one of the major themes of the development of the global oil industry today. As a traditional industrial industry which is an industry with a relatively low degree of automation, the oil and gas industry faces the new situation and new trends of the accelerated energy revolution and energy transformation. It must effectively utilize digital technologies to drive business model reconstruction, management model reform, business model innovation, and core competence enhancement. Digital technologies are related technologies represented by cloud computing, Internet of Things, 5G, big data, artificial intelligence, etc. Ultimately, the transformation and upgrading of the industry and value growth will be realized. The transformation and upgrading of oil and gas production involves automatic collection of on-site production data, real-time tracking, diagnosis and optimization, and remote control of on-site automatic adjustment devices. Through transformation and upgrading, on-site production automation management is finally realized, management efficiency is improved, and production and operation costs are saved. It is the need of industry transformation and development, and it is also an urgent need on site. Many aspects need to be involved in this upgrade process, not only to achieve breakthroughs in hardware, but also to form a supporting package with this hardware in software. For example, the realization of real-time optimization based on the collection of massive data needs to be combined with the transient simulation of oil and gas wells. Horizontal well and branch well development technology is a revolutionary technology in the petroleum industry today. As the economic benefits of using this technology for oilfield development become more and more obvious, its application becomes more and more extensive. It is not only applied to difficult-to-explore oil reservoirs such as low permeability and heavy oil, but also more and more oil reservoirs with better properties. Its main function is to increase the contact area with the formation and expand the oil flow channel, so that crude oil and natural gas can be more easily produced on the ground. Therefore, the development of the transient model of horizontal well productivity prediction is one of the important basic tasks in the digital transformation of oil and gas. This is the need to accelerate the development of digital transformation of oilfields. This is not only of great significance for enriching the basic theory of horizontal well production technology, but also of obvious value for the actual application of oilfield production.

In the development of horizontal wells, the continuous inflow of reservoir fluids leads to changes in the inflow and pressure of the horizontal wellbore. The changes in the inflow and pressure of the horizontal wellbore affect the total production of the horizontal well and the design of the horizontal wellbore structure parameters. Therefore, the development of horizontal well transient productivity prediction research involves the fluid flow state in the horizontal wellbore and the interaction between it and the reservoir. Many scholars at home and abroad have carried out related research.

In 1991, Ozkan et al. [1,2,3] derived the point source solution in the Laplace transform domain and obtained the general solution after considering a variety of well structure (vertical well, horizontal well. and fractured well) and reservoir boundary conditions. Subsequently, in 1995, Ozkan et al. [4] assumed that the single-phase adiabatic weakly compressible fluid had constant compressibility coefficient and viscosity and its flow state in the horizontal wellbore in the homogeneous reservoir was laminar or turbulent and regarded the flow in the horizontal wellbore as one-dimensional flow. They proposed a semi analytical model coupling the horizontal wellbore and reservoir seepage. The flow rate and pressure distribution along the horizontal wellbore were studied, and the effectiveness of the assumption of infinite conductivity was discussed.

In 1998, Penmatcha et al. [5] established a three-dimensional, anisotropic semi analytical solution transient model for the coupling of reservoir and horizontal wellbore with infinite conductivity and limited conductivity to calculate the productivity of horizontal wells in rectangular reservoirs. The finite conductivity model considers the friction pressure drop, acceleration pressure drop, and pressure drop caused by reservoir fluid inflow in the wellbore.

In 2004, Duan et al. [6] proposed an unsteady mathematical model of a coupling reservoir and wellbore based on the unsteady seepage theory and derived the numerical solution of the model by using the boundary integral method. The model regards the seepage in the reservoir and the flow in the wellbore as an interactive whole and considers many influencing factors, such as fluid friction resistance, momentum change, mixed interference of horizontal wellbore wall inflow, and so on. Huang Cheng et al. [7] put forward the productivity prediction model of horizontal wells on the basis of the former and according to the calculation method of Chino et al. [8]. In order to simplify the calculation, a calculation model in Laplace space is further proposed. Through iterative solution, the distribution of pressure drop and production along the length of a horizontal wellbore under constant production can be obtained. Chen Wei et al. [9] verified the capacity prediction model proposed by Duan Yonggang and Huang Cheng by using the published production data. The results show that the model is true and reliable.

In 2005, Wang et al. [10] proposed a pressure drop model coupled with a horizontal wellbore based on the elastic unsteady seepage model in a horizontal equal thickness infinite homogeneous reservoir. The distribution of pressure drop and production in the length direction and the model can optimize the design of perforation parameters for perforated completion horizontal wells and extend the waterless recovery period.

In 2013, Tajer et al. [11] proposed an analytical model of three-phase flow in multi-stage fractured horizontal wells in low permeability reservoirs based on the time-dependent trilinear flow model and the dissolved gas drive model. Compared with numerical simulation, the production forecast using this model is much faster.

In 2017, Chen et al. [12] established a transient pressure model of a multi-fracture horizontal well with limited conductivity in consideration of wellbore pressure loss and used the Stehfest method to perform Laplace transform numerical inversion to obtain the transient pressure solution.

In 2020, Chu et al. [13] considered that the transient study of horizontal well productivity is an important method for unconventional oil and gas reservoir productivity, reserves prediction, and completion evaluation. Combing the Laplace transform and finite difference method to establish a semi-analytical solution, a transient model of horizontal well productivity is established. This model is suitable for multi-fracture horizontal wells with different bottom hole pressures, different hydraulic fracture properties, and different start-up times.

Although there is a lot of research on transient models for productivity prediction of horizontal wells in recent years, most of them are too complex or lack of example verification or are based on some assumptions and simplification, as shown in

Table 1. Statistical analysis of research results.

Literatures	Main Research Content	There Are Problems or Shortcomings
Ozkan et al. [1,2,3,4]	A well test model for the pressure recovery of a dimensionless horizontal well is established.	The horizontal well is simplified as a straight line, the model is too complicated, and there are many constraints on the assumptions.
Penmatcha et al. [5]	Derived horizontal well productivity prediction model based on Babu and Odeh.	Horizontal wells are simplified as straight lines, aiming at closed rectangular oil reservoirs.
Duan et al. [6,7]	Based on the calculation method of Chino et al. [8], an unsteady mathematical model of the reservoir and wellbore coupling is established, and the numerical solution of the model is derived.	Horizontal wells are simplified as straight lines, aiming at closed rectangular oil reservoirs, and the model is complex.
Wang et al. [10]	Based on the transient pressure model of vertical wells, the unsteady seepage model of horizontal wells in infinitely large homogeneous reservoirs of horizontal thickness is derived, and a pressure drop model coupled with horizontal wellbore is proposed.	Horizontal wells are simplified to straight lines. There are limitations in assuming that the horizontal well segment is a vertical well model because the oil layer height is not infinite.
Tajer et al. [11]	An analytical model of three-phase flow in multi-stage fractured horizontal wells in low permeability reservoirs is proposed.	Horizontal wells and fractures are simplified as straight lines, and it is impossible to analyze the influence of well trajectory on productivity. The productivity model of fractured wells is derived from the permeability characteristics of shale tight oil reservoirs.
Chen et al. [12]	A transient pressure model for multi-fracture horizontal wells with limited conductivity is established, and the Stehfest method is used to perform Laplace transform numerical inversion to obtain the transient pressure solution.
Chu et al. [13]	The Laplace transform and the micro-element method are combined to establish a semi-analytical multi-stage fractured horizontal well pressure recovery model.

below. Therefore, this paper carries out transient model research on the basis of the previously established reliable steady-state model [14,15]. The corresponding reservoir transient seepage model and horizontal well production prediction model are established.

2. Establishment of Coupled Mathematical Model

The transient productivity prediction model of horizontal wells should consider not only the influence of horizontal wellbore pressure drop, but also the important factor of reservoir production time. As an important basic part of the horizontal well transient productivity prediction model, the horizontal well single-phase seepage transient model exists in all kinds of reservoirs when the pressure is high in the early stage of oilfield development. It is one of the important seepage modes in the process of reservoir development. Horizontal well production involves two flow processes: the flow of fluid in the formation section and the flow of fluid in the wellbore section. In this study, models for fluid flow in the reservoir and in the wellbore are developed. The two models are then coupled to establish a prediction model for the production rate of horizontal wells.

2.1. Calculation of the Spatial Potential of the Uniform Inflow into the Horizontal Segment

2.1.1. Spatial Steady-State Point Sink Function

Assuming that a point in the space sinks M, according to the seepage theory, with point M as the center, the output is q, and the seepage velocity of a spherical surface with an arbitrary radius r is:

$$V=\frac{q}{4\pi r^2}$$

(1)

In addition, according to the definition of the potential and Darcy’s law:

$$V=\frac{d\varphi }{dr}$$

(2)

The above two formulas are equal:

$$\frac{q}{4\pi r^2}=\frac{d\varphi }{dr}$$

(3)

The expression of spatial potential obtained by separating the two equations and integrating them is:

$$\varphi =-\frac{q}{4\pi r}+C$$

(4)

2.1.2. Spatial Instantaneous Point Source Function

It is common knowledge that in the three-dimensional unbounded homogeneous layer, the pressure distribution is uniform at the initial moment and the diffusion coefficient is η_r. Considering a spatial instantaneous point sink with constant intensity q₀, one-dimensional unstable radial seepage occurs at the moment t = 0. Taking spherical coordinates, the spatial point sink is represented as a small ball with a very small radius and the governing equation of unsteady seepage is as follows:

$$\frac{{\partial }^2(r\Delta p)}{\partial r^2}=\frac{1}{{\eta }_r}\frac{\partial (r\Delta p)}{\partial t}$$

(5)

$$r^2={(x-x_w)}^2+{(y-y_w)}^2+{(z-z_w)}^2$$

(6)

Condition of definite solution:

$$\Delta p(r,0) = 0, \Delta p(\infty ,t) = 0; \underset{\epsilon \to 0^{+}}{{\displaystyle \lim }}\frac{4\pi k}{\mu }{(r^2\frac{\partial \Delta p}{\partial r})}_{r=\epsilon }=-q_0\delta (t)$$

(7)

Obtained by Laplace transformation:

$$\Delta p(r,t)=\frac{q_0\mu }{8\pi k_r\sqrt{\pi {\eta }_r{t}^3}}\exp (-\frac{r^2}{4{\eta }_{r}t})$$

(8)

When q₀ = ϕc_t, there is an instantaneous Green source function:

$$G_s(r,t)=\frac{1}{\sqrt{{\left(4\pi {\eta }_{r}t\right)}^3}}\exp (-\frac{r^2}{4{\eta }_{r}t})$$

(9)

Under isotropic conditions (anisotropy can be processed by coordinate transformation), the above formula can be directly decomposed:

$$\begin{array}{c}{G}_s(r,t)=\frac{1}{\sqrt{{\left(4\pi {\eta }_{r}t\right)}^3}}\exp (-\frac{r^2}{4{\eta }_{r}t})\\ =\frac{1}{2\sqrt{\pi {\eta }_{x}t}}\exp (-\frac{{(x-x_w)}^2}{4{\eta }_{x}t})\cdot \frac{1}{2\sqrt{\pi {\eta }_{y}t}}\exp (-\frac{{(y-y_w)}^2}{4{\eta }_{y}t})\cdot \frac{1}{2\sqrt{\pi {\eta }_{z}t}}\exp (-\frac{{(z-z_w)}^2}{4{\eta }_{z}t})\\ =G_x(x,t)\cdot G_y(y,t)\cdot G_z(z,t)\end{array}$$

(10)

It can be seen that under certain conditions, the three-dimensional seepage problem can be decomposed into three one-dimensional seepage problems. Therefore, we can realize the way to solve the multi-dimensional seepage problem with one-dimensional point source function.

For the problem of continuous point source in space, it can also be obtained by integration of q(t) = q:

$$\Delta p(r,t)=\frac{1}{\varphi c_t}{\displaystyle {\int }_0^{t}q(t-\tau )G_s}(r,\tau )\mathrm{d}\tau =\frac{1}{\varphi c_t}q{\displaystyle {\int }_0^t\frac{1}{\sqrt{{\left(4\pi {\eta }_r\tau \right)}^3}}\exp (-\frac{r^2}{4{\eta }_r\tau })}\mathrm{d}\tau$$

(11)

Order $\eta =\frac{r}{2\sqrt{{\eta }_r\tau }}$, then:

$$\mathrm{d}\eta =\frac{r}{2\sqrt{{\eta }_r}}-\frac{1}{2}{\tau }^{-\frac{1}{2}-1}\mathrm{d}\tau =-\frac{r}{4\sqrt{{\eta }_r}}{\tau }^{-\frac{3}{2}}d\tau =-\frac{r}{4\sqrt{{\eta }_r}}\frac{1}{\sqrt{{\tau }^3}}d\tau$$

(12)

Then Equation (11) can become,

$$\begin{array}{ll}\Delta p(r,t)\hfill & =\frac{1}{\varphi c_t}q{\displaystyle {\int }_0^t\frac{1}{\sqrt{{\left(4\pi {\eta }_r\tau \right)}^3}}\exp (-\frac{r^2}{4{\eta }_r\tau })}\mathrm{d}\tau =\frac{1}{\varphi c_t}q{\displaystyle {\int }_0^t\frac{1}{4\pi {\eta }_r\sqrt{4\pi {\eta }_r{\tau }^3}}\exp (-\frac{r^2}{4{\eta }_r\tau })}\mathrm{d}\tau \hfill \\ & =\frac{1}{\varphi c_t}q{\displaystyle {\int }_0^t\frac{1}{2\pi {\eta }_r\sqrt{\pi }4\sqrt{{\eta }_r{\tau }^3}}\exp (-\frac{r^2}{4{\eta }_r\tau })}\mathrm{d}\tau =-\frac{1}{\varphi c_t}\frac{q}{2\pi {\eta }_{r}r\sqrt{\pi }}{\displaystyle {\int }_0^t-\frac{r}{4\sqrt{{\eta }_r{\tau }^3}}\exp (-\frac{r^2}{4{\eta }_r\tau })}\mathrm{d}\tau \hfill \\ & =-\frac{1}{\varphi c_t}\frac{q}{2\pi {\eta }_{r}r\sqrt{\pi }}{\displaystyle {\int }_{\infty }^{\frac{r}{2\sqrt{{\eta }_{r}t}}}\exp (-{\eta }^2)}\mathrm{d}\eta =\frac{q\mu 2}{4\pi kr\sqrt{\pi }}(-{\displaystyle {\int }_{\infty }^{\frac{r}{2\sqrt{{\eta }_{r}t}}}\exp (-{\eta }^2)}\mathrm{d}\eta )\hfill \\ & =\frac{q\mu }{4\pi kr}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\hfill \end{array}$$

$${\eta }_r=\frac{k}{\mu \varphi c_t}$$

(13)

$$p(r,t)-p_{wf}=\frac{q\mu }{4\pi kr}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})$$

(14)

Then, for the spatial continuity point sink, the total differential on both sides of the formula becomes:

$$p(r,t)-p_{wf}=-\frac{q\mu }{4\pi kr}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})$$

(15)

$$d(p(r,t)-p_{wf})=d\left[\frac{q\mu }{4\pi kr}erfc\left(\frac{r}{2\sqrt{{\eta }_{r}t}}\right)\right]$$

(16)

$$\mathrm{d}p=\frac{\mu }{k}\mathrm{d}\varphi =\frac{\mu }{4\pi kr}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\mathrm{d}q$$

(17)

2.1.3. Spatial Instantaneous Line Sink Function

This paper mainly carries out a spatial instantaneous line sink model on the basis of the previously established reliable steady-state model [14,15]. In the unbounded three-dimensional formation, there is a horizontal well with a measured length L (as shown in Figure 1). When producing with production q, the coordinates of the heel end and toe end are (x₁, y₁, z₁), (x₂, y₂, z₂). Assuming that single phase crude oil is flowing through the formation, the horizontal well is a linear sink with uniform inflow.

The horizontal well is equally divided into m segments according to its length. It can be seen that when m is large enough, each segment can be approximated as a straight line segment with the length of $\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$m$}\right.$. The start coordinate of each segment is (x_si, y_si, z_si) and the end coordinate is (x_ei, y_ei, z_ei), where i = 1, 2, 3... m.

Take any point in one segment with coordinates (x, y, z) as the end point, and the distance from the beginning of the segment is:

$$s=\sqrt{{(x-x_{si})}^2+{(y-y_{si})}^2+{(z-z_{si})}^2}$$

(18)

By obtaining the total differential of both sides of Equation (18), the infinitesimal part ds meets the following formula.

$$\mathrm{d}s=\frac{1}{s}[(x-x_{si})\mathrm{d}x+(y-y_{si})\mathrm{d}y+(z-z_{si})\mathrm{d}z]$$

(19)

For ds, the production of this part is $dq=\frac{q}{L}ds$. For the continuous point sink problem in space, the potential generated in space point (X, Y, Z) is

$$\mathrm{d}\varphi =-\frac{\mathrm{d}q}{4\pi r}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})$$

(20)

$$\mathrm{d}\varphi =-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\mathrm{d}s$$

(21)

$$\mathrm{d}\varphi =-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}[(x-x_{si})\mathrm{d}x+(y-y_{si})\mathrm{d}y+(z-z_{si})\mathrm{d}z]$$

(22)

We supposed that f(x, y, z, t), g(x, y, z, t), and h(x, y, z, t) are as follows:

$$f(x,y,z,t)=-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(x-x_{si})$$

(23)

$$g(x,y,z,t)=-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(y-y_{si})$$

(24)

$$h(x,y,z,t)=-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(z-z_{si})$$

(25)

Thus, as the spatial region belongs to 3D single connected open region G and f(x, y, z, t), g(x, y, z, t), and h(x, y, z, t) have the first-order partial derivative in this region (For infinitesimals, r and t are constants), the following formula is met.

$$\frac{\partial g}{\partial z}=\frac{\partial h}{\partial y}, \frac{\partial g}{\partial z}=\frac{\partial h}{\partial y}, \frac{\partial f}{\partial z}=\frac{\partial h}{\partial x}$$

(26)

Therefore, the potential generated by this segment on space (X, Y, Z) can be found by Equation (28):

$${\varphi }_i={\displaystyle {\int }_{(x_{si},y_{si},z_{si})}^{(x_{ei},y_{ei},z_{ei})}-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\mathrm{d}s}+C$$

(27)

$${\varphi }_i=-\frac{q}{4\pi L}({\displaystyle {\int }_{x_{si}}^{x_{ei}}f(x,y_{si},z_{si},t)\mathrm{d}x+{\displaystyle {\int }_{y_{si}}^{y_{ei}}g(x,y,z_{si},t)\mathrm{d}y}+{\displaystyle {\int }_{z_{si}}^{z_{ei}}h(x,y,z,t)\mathrm{d}z}})+C$$

(28)

That is:

$$\begin{array}{l}{\varphi }_i={\displaystyle {\int }_{x_{si}}^{x_{ei}}-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(x-x_{si})\mathrm{d}x}\\ +{\displaystyle {\int }_{y_{si}}^{y_{ei}}-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(y-y_{si})\mathrm{d}y}\\ +{\displaystyle {\int }_{z_{si}}^{z_{ei}}-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(z-z_{si})\mathrm{d}z}+C\end{array}$$

(29)

In the above formula, the first item on the right has x as the integration variable, and y and z are constant. The rest are integrated similarly.

The first item on the right is integrated as follows:

$$\begin{array}{c}{\displaystyle {\int }_{x_{si}}^{x_{ei}}-\frac{q}{4\pi rL}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(x-x_{si})\mathrm{d}x}=-\frac{q}{4\pi L}{\displaystyle {\int }_{x_{si}}^{x_{ei}}\frac{1}{r}erfc(\frac{r}{2\sqrt{{\eta }_{r}t}})\frac{1}{s}(x-x_{si})\mathrm{d}x}\\ =-\frac{q}{4\pi L}{\displaystyle {\int }_{x_{si}}^{x_{ei}}\left[\begin{array}{l}\frac{1}{\sqrt{{(x-X)}^2+{(y-Y)}^2+{(z-Z)}^2}}\cdot \\ erfc(\frac{\sqrt{{(x-X)}^2+{(y-Y)}^2+{(z-Z)}^2}}{2\sqrt{{\eta }_{r}t}})\cdot \\ \frac{1}{\sqrt{{(x-x_{si})}^2+{(y-y_{si})}^2+{(z-z_{si})}^2}}(x-x_{si})\end{array}\right]\mathrm{d}x}\end{array}$$

(30)

After simplifying the formula and supposing that a = (y − Y)² + (z − Z)², b = (y − y_si)² + (z − z_si)²,

$$=-\frac{q}{4\pi L}{\displaystyle {\int }_{x_{si}}^{x_{ei}}\frac{1}{\sqrt{{(x-X)}^2+a}}erfc(\frac{\sqrt{{(x-X)}^2+a}}{2\sqrt{{\eta }_{r}t}})\frac{1}{\sqrt{{(x-x_{si})}^2+b}}(x-x_{si})\mathrm{d}x}$$

(31)

Assuming that the function satisfies $f(x,y_{si},z_{si},t)=\frac{1}{\sqrt{{(x-X)}^2+a}}erfc(\frac{\sqrt{{(x-X)}^2+a}}{2\sqrt{{\eta }_{r}t}})\frac{1}{\sqrt{{(x-x_{si})}^2+b}}(x-x_{si})$, then substitute into Equation (30), and Equation (30) equals the integration of f(x, y_si, z_si, t) in the domain of [x_si, x_ei].

Finally, the entire potential in point (X, Y, Z) from the horizontal well can be:

$$\varphi ={\displaystyle \sum _{i=1}^m{\varphi }_i=-\frac{q}{4\pi L}{\displaystyle \sum _{i=1}^m({\displaystyle {\int }_{x_{si}}^{x_{ei}}f(x,y_{si},z_{si},t)dx+{\displaystyle {\int }_{y_{si}}^{y_{ei}}g(x,y,z_{si},t)dy}+{\displaystyle {\int }_{z_{si}}^{z_{ei}}h(x,y,z,t)dz}})}}$$

(32)

Due to the positional relationship, there is a difference between the confluence of fluids at both ends of the horizontal well in the reservoir and the confluence of fluids in the middle part of the reservoir. Due to the interference between the segments of the wellbore and the pressure drop of the fluid flow in the wellbore, the flow rate from the oil layer into the horizontal wellbore is different. For this reason, a horizontal well is divided into many segments of line sinks. Since the length of each segment is very short, assuming that the fluid flows in uniformly from all places along the segment of the oil layer, the potential generated by each segment is equivalent to the potential of a horizontal well in Equation (32).

Oil well productivity prediction is also closely related to the types of reservoirs. Generally, four types of reservoirs can be distributed: top closed bottom water reservoirs, gas cap bottom water reservoirs, upper and lower closed edge water reservoirs, and upper and lower closed boundary reservoirs. Taking the reservoir type as an example of the upper and lower closed boundary reservoirs, the calculation method of its potential is introduced below.

2.2. Calculation of Spatial Potential of Uniform Inflow into Horizontal Section in Upper and Lower Closed Reservoir

The process of building the model in this part is similar to that in the previously established reliable steady-state model [14,15]. For the upper and lower closed boundary reservoir as shown in Figure 2, a horizontal well with length L is divided into N segments. According to the mirror reflection principle, as shown in Figure 3, we can get:

$${\varphi }_j(X,Y,Z,t)=-\frac{q_j}{4\pi }\{{\displaystyle \sum _{n=-\infty }^{\infty }[\xi (x,y,2nh+z,X,Y,Z,t)+\xi (x,y,2nh-z,X,Y,Z,t)}]\}+C_j$$

(33)

where ϕ_j is the potential generated by the j-th line sink at any point in the oil layer, q_j is the flow rate of the j-th line sink, h is the oil thickness, z is the distance between each part of the well and the bottom of the oil layer, C_j is a constant, and ξ is a function defined by the following formula:

$$\begin{array}{l}{\xi }_j(x,y,2nh+z,X,Y,Z,t)=\frac{1}{L_j}{\displaystyle \sum _{i=1}^m({\displaystyle {\int }_{x_{si}}^{x_{ei}}f(x,y_{si},z_{si},t)dx+{\displaystyle {\int }_{y_{si}}^{y_{ei}}g(x,y,z_{si},t)dy}}}\\ +{\displaystyle {\int }_{4nh+z_{si}}^{4nh+z_{ei}}h(x,y,z,t)dz)}\end{array}$$

(34)

where L_j are the length of the j-th segment line sink, x_s1 and x_em are the start and end abscissas of the j-th segment line in the x-axis direction, and the other parameters are y and z direction coordinates.

3. Horizontal Well Flow Relationship

The process of building the model in this part is similar to that in the previously established reliable steady-state model [14,15]. According to the potential superposition principle, the potential generated by the whole horizontal well in the oil layer is

$$\varphi (X,Y,Z,t)={\displaystyle \sum _{j=1}^N{\varphi }_j(X,Y,Z,t)}+C=-{\displaystyle \sum _{j=1}^N\frac{q_j}{4\pi }}{\phi }_j+C$$

(35)

For different types of reservoirs, the formula φ_j is respectively equal to the formulas in braces for different reservoir types in Equation (33).

It can be obtained from Equation (35).

$${\varphi }_e={\displaystyle \sum _{j=1}^N{\varphi }_{je}}+C$$

(36)

The formula ϕ_e is the potential function at the constant pressure boundary or oil drain boundary, and ϕ_je is the potential generated by the j-segment line sink at the constant pressure boundary or the oil drain boundary.

It is obtained by Equation (35) and Equation (36).

$$\varphi (X,Y,Z,t)={\varphi }_e+{\displaystyle \sum _{j=1}^N[{\varphi }_j(X,Y,Z,t)-{\varphi }_{je}]}$$

(37)

According to the relationship between pressure and potential, we can obtain

$$p(X,Y,Z,t)=\frac{\mu }{k}\varphi (X,Y,Z,t)-\rho gh$$

(38)

where p is the pressure at any point in the reservoir, k is reservoir permeability, μ is viscosity, ρ is density, and g is the acceleration of gravity.

Substitute Equation (37) into Equation (38) to obtain

$$p(X,Y,Z,t)=p_e+\frac{\mu }{k}{\displaystyle \sum _{j=1}^N[{\varphi }_j(X,Y,Z,t)-{\varphi }_{je}]}-\rho g(Z-z_e)$$

(39)

In the formula, p_e and z_e are the pressure and z coordinates at the corresponding boundary respectively.

Equation (39) can reflect the seepage law of the wellbore in the formation, that is, the relationship between the external pressure of the wellbore and the output flowing into the wellbore, and the variable mass flow law in the wellbore needs to be considered in establishing the coupling model.

4. Variable Mass Flow of Horizontal Wellbore with Different Completions

Due to the complexity of the problem, in actual calculations, the numerical calculation method is generally used to obtain the pressure drop of each section in sections. Suppose the horizontal wellbore of length L is divided into N sections, and the pressure drop calculation method of single-phase variable mass flow is obtained according to the principle of conservation of mass and the principle of conservation of momentum. That is, the calculation formula for the pressure drop of the oil single-phase variable mass flow at any micro-element section △x(dx) length under different well completion conditions. They are as follows.

4.1. Open Hole Completion

$$\Delta p_w=\frac{2f_{oh}\rho }{{\pi }^2{D}^5}{\left(2Q+q\right)}^2\Delta x+\frac{16\rho q}{{\pi }^2{D}^4}\left(2Q+q\right)$$

(40)

where f_oh is the frictional resistance coefficient of the wall with radial inflow, z_e f_oh = C_ohf, and f is the wall friction coefficient of common horizontal circular pipe flow. The wall friction coefficient calculation method for the following three completion modes is the same as here. C_oh is the correction coefficient considering the effect of radial inflow on wall friction in a horizontal wellbore of open hole completion, and its value can be obtained via the test; Q is the mainstream flow at the upstream section of this section; q is the total flow from the reservoir into this section of the horizontal wellbore; ρ is fluid density; D is the horizontal wellbore diameter.

4.2. Perforation Completion

The pressure drop of the j-th section in the wellbore of the perforation completion is:

$$\Delta p_{wj}=\frac{8f_{cp}\rho Q_j^2\Delta x}{{\pi }^2{D}^5}\left[1+\frac{q_j}{Q_j}+\left(\frac{1}{3}+\frac{1}{6n^2}\right){\left(\frac{q_j}{Q_j}\right)}^2\right]+\frac{32\rho Q_j{q}_j}{{\pi }^2{D}^4}\left(1+\frac{q_j}{2Q_j}\right)$$

(41)

where n is the number of perforation in this section.

4.3. Screen Completion

The process of building the model in this part is similar to that in the previously established reliable steady-state model [15]. As shown in Figure 4, the pressure loss along the center wellbore according to momentum conservation is:

$$-\frac{\mathrm{d}{p}_{w,i}}{\mathrm{d}x}=\rho g\sin {\theta }_i+\frac{f_{f,i}\rho }{2D}{\left(\frac{V_{1,i}+V_{2,i}}{2}\right)}^2+\frac{8\rho V_{1,i}{V}_{r,i}}{D}+\frac{16\rho V_{r,i}^2}{D^2}\mathrm{d}x+\frac{\mathrm{d}{p}_{mix,i}}{\mathrm{d}x}$$

(42)

According to the conservation of momentum, the pressure loss along the annulus is:

$$-\frac{\mathrm{d}{p}_{aw,i}}{\mathrm{d}x}=\rho g\sin {\theta }_i+\frac{f_{mf,i}\rho }{2\overline{D}}{\left(\frac{V_{m1,i}+V_{m2,i}}{2}\right)}^2+\rho \frac{V_{m2,i}^2-V_{m1,i}^2}{\mathrm{d}x}+\frac{\mathrm{d}{p}_{mmix,i}}{\mathrm{d}x}$$

(43)

4.4. ICD Completion

The nozzle type ICD uses nozzles to create pressure resistance. The number and diameter of the nozzles can be selected to produce the required pressure drop when the fluid flows through the nozzles at a certain flow rate. It mainly includes the pressure drop caused by the contraction effect (the restrictor valve) and the friction pressure drop of the fluid passing through the pipe system. The calculation method is as follows:

$$\Delta p=\Delta p_{cons}+\Delta p_{fric}$$

(44)

The long tube type ICD uses many tubes of a certain diameter and length to apply a specific homeopathic pressure to a specific flow rate of fluid. This equipment combines the pressure drop produced by the fluid flowing through the throttling chamber and the pressure drop produced by the straight pipe. The calculation method is as follows:

$$\Delta p=\left(H_1+H_2\right)\rho$$

(45)

$$H_1=f\frac{L_{tube}}{D_{tube}}\frac{V_{tube}^2}{2},H_2=K\frac{V_{tube}^2}{2}$$

(46)

The spiral ICD uses many spiral channels with preset diameters and lengths to apply a homeopathic pressure to the fluid passing through at a certain flow rate. The calculation method is as follows:

$$\Delta p={\left(\frac{\rho }{{\rho }_{mix}}\frac{{\mu }_{mix}}{\mu }\right)}^{1/4}\frac{{\rho }_{mix}}{\rho }{a}_{ICD}{q}^2$$

(47)

The above-mentioned pressure drop is the radial pressure drop of the ICD tool, that is, the pressure drop loss that the fluid needs to overcome when flowing from the pipe wall into the central wellbore. It is coupled with the ICD central wellbore (axial) pressure drop along the path to perform the ICD wellbore variable mass flow pressure drop. The calculation method of ICD central wellbore (axial) pressure drop along the way is the same as the calculation method of screen completion pressure drop.

5. Coupling Model of Inflow Performance and Flow in Wellbore and Its Solution

The process of building the model in this part is similar to that in the previously established reliable steady-state model [14,15]. According to the flow in the wellbore, a coupling equation is established to solve the flow in the formation, and the coordinated production in accordance with the two flow laws is obtained, that is, the coordinated production of the oil well.

The three-dimensional transient seepage of fluid in the oil layer and the flow in the wellbore are both interrelated and affect each other. Suppose the pressure at the midpoint of the j section line sink on the horizontal well is p_w,j, and the potential generated by the i section line sink at the midpoint of the j section line sink is ϕ_ij, then according to Equation (39), obtain

$$p_{w,j}=p_e+\frac{\mu }{k}{\displaystyle \sum _{i=1}^N({\varphi }_{ij}-{\varphi }_{ie})}+\rho g(z_e-z_w) (j=1,2,\phi ,\mathrm{N})$$

(48)

The above formula is deformed as

$${\displaystyle \sum _{i=1}^N\lambda q_i}({\phi }_{ij}-{\phi }_{ie})=p_e-p_{w,j}+\rho g(z_e-z_w) (j=1,2,\phi ,\mathrm{N})$$

(49)

where $\lambda =\frac{\mu }{4\pi k}$.

The pressure at the midpoint of section j in the wellbore is

$$p_{w,j}=p_{1,j}-0.5d{p}_{w,j} (j=1,2,\phi ,\mathrm{N})$$

(50)

where p_2,N = p_wf, and p_wf is the flowing pressure at the heel end of the wellbore.

$$p_{1,j+1}=p_{2,j}=p_{1,j}-\Delta p_{w,j} (j=1,2,\phi ,\mathrm{N})$$

(51)

The total production of the whole well is

$$Q_o=\frac{(q_1+q_2+q_3+\cdot \cdot \cdot +q_N)}{B_o}$$

(52)

where B_o is the volume factor of the formation of crude oil.

In the above-mentioned coupled model, p_w and q are unknowns, which can be solved using the iterative method. Assuming a set of p_w values, q array is solved using Equation (49). Then, q array is substituted into the pressure drop formulas in the above-mentioned variable mass flow of a horizontal wellbore with different completions, and Equation (50) is used to update p_w array from heel to toe. Then use Equation (49) to solve the q array, and so on, until q array and p_w array reach a certain computational accuracy. Finally, the total well production can be obtained from Equation (52). The flowchart is shown in Figure 5 below.

6. Samples Calculation and Verification

6.1. Samples Calculation

This paper carries out transient model research on the basis of the previously established reliable steady-state model [14,15]. The prediction accuracy of the basic model has been described in these documents. The newly established model can consider the real flow in the annulus and can simulate the formation energy and output, which changes with time as the production proceeds, and the prediction is more in line with the actual production situation. For example, with constant bottom hole pressure production, the simulated production changes over time are as follows. The basic parameters of oil reservoirs and horizontal wells are shown in

Table 2. Parameters of oil reservoir and horizontal well.

Name	Value	Unit
Horizontal permeability	10	mD
Vertical permeability	10	mD
Crude oil volume factor	1.281
Porosity	0.3
Total compressibility coefficient of formation	0.0002	1/MPa
Crude oil viscosity	0.6365	mPa.s
Reservoir thickness	30	m
Horizontal well length	400	m
Original formation pressure	30	MPa
Well bottom flowing pressure	28	MPa
Completion method	Open hole completion
Oil reservoir type	the upper and lower closed boundary reservoir

.

The software algorithm is compiled using VC++, and the interface design uses C#. The calculation results are shown in

Table 3. Production prediction with constant bottom hole flowing pressure.

Development Time (Year)	Average Formation Pressure (MPa)	Recovery Degree (%)	Average Production (m3/d)
1	29.65	1.14	204.52
2	29.37	2.1	170.25
3	29.13	2.91	141.49
4	28.93	3.59	117.42
5	28.77	4.154	97.32
6	28.64	4.624	80.58
7	28.53	5.014	66.66
8	28.44	5.334	55.10
9	28.36	5.604	45.52
10	28.30	5.824	37.59

. It can be seen from the table that with the increase of the production time, crude oil is continuously being produced, the degree of production gradually increases, and the energy contained in the formation gradually decreases. The formation pressure gradually decreases. As the formation pressure decreases, the bottom hole pressure is the same. The output of oil wells gradually decreases. The dynamic process predicted by the simulation is the same as the change of the actual oil reservoir production process, which is in line with the actual production situation.

The fluid production and pressure distribution of a horizontal well in the 10th year with a constant bottom hole flowing pressure are shown in Figure 6.

6.2. Samples Verification

6.2.1. Pressure Recovery Verification of Iran’s MIS Oilfield

The basic parameters and test data of the two wells A and B of Missan oilfield are shown in

Table 4. Basic parameters of oil reservoir and two horizontal wells.

Well Parameters\Well Number	MIS A	MIS B
Oil layer thickness/m	157.79	137.16
Porosity/decimal	0.116	0.116
total compression coefficient/(1/MPa)	0.001276264	0.001276264
Absolute permeability (in the X, Y direction)/mD	43	578
Absolute permeability (Z direction)/mD	43	578
Reservoir temperature/°C
Initial saturation of formation crude oil/decimal
Crude oil volume factor	1.09	1.09
Crude oil viscosity/mPa.s	1.8	1.8
Relative density of crude oil	0.832	0.832
Well type	Horizontal well	Horizontal well
Oil reservoir type	Infinite homogeneous reservoir	Infinite homogeneous reservoir
Diameter of oil drain area/m		1524
Horizontal wellbore length/m	394.1	422.34
Wellbore diameter/in	6 1/8″	6 1/8″
Skin factor	10	−3
Wellbore storage factor/(m3/MPa)	7.093	7.093
Reservoir pressure/MPa	2.91	3.059
Test production time before shutting in/h	18	18
Shut-in bottom hole flowing pressure after test production/MPa	2.21	2.834

,

Table 5. Test production data of MIS A.

No.	Duration	ESP (Hz)	Pwf (psi)	Δp (psi)	Rate (bbl/d)	PI (bbl/d/psi)
1	22:00–2:00	35	351	72	255	3.86
2	2:04–6:15	40	346	77	330
3	6:19–10:15	45	333	90	398
4	10:21–16:00	50	252	171	596

, and

Table 6. Test production data of MIS 322 CN-H2 wells.

No.	Oil Rate (bbl/d)	ΔP (psi)	Pwf (psi)	Watercut (%)	Choke Size (1/64″)	ESP (Hz)	PI (bbl/d/psi)
1	3392	25	419	0.2	42	45	136.57
2	3820	28	416	0.2	52	47
3	4260	31	413	0.1	50	50

, respectively. After the oil well is completed, the well test is started before production. To produce with a certain output, test the bottom hole flow pressure. When the oil well is stable, that is, when the bottom hole flow pressure and production do not change, shut in the well and test the bottom hole pressure changes to obtain the pressure recovery data after shut-in.

The verification results of MIS A are shown in Figure 7 and

Table 7. Error analysis of the results of pressure recovery simulation calculation and test data.

Project	Average Relative Error (%)	Absolute Average Relative Error (%)	The Relative Standard Deviation (%)
Simulation calculation	0.2821	6.7422	1.1514

.

The verification results of MIS B are shown in Figure 8 and

Table 8. Error analysis of the results of pressure recovery simulation calculation and test data.

Project	Average Relative Error (%)	Absolute Average Relative Error (%)	The Relative Standard Deviation (%)
Simulation calculation	1.2948	1.356	0.0989

.

6.2.2. Pressure Recovery Verification of Hafaya Oilfield

The basic parameters of Hafaya well A are shown in

Table 9. Basic parameters of oil reservoir and horizontal well A.

Well Parameters\Well Number	Hafaya Well A
Oil layer thickness/m	70
Porosity/decimal	0.177
Comprehensive compression coefficient/(1/MPa)Total compression coefficient/(1/MPa)	0.0013154
Absolute permeability (in the X, Y direction)/mD	0.076
Absolute permeability (Z direction/mD)	0.076*0.36
Reservoir temperature/°C	81.2
Initial saturation of formation crude oil/decimal
Crude oil volume factor	1.237
Crude oil viscosity/mPa.s	1.33
Relative density of crude oil	0.904
Well type	Horizontal well
Reservoir type	Infinite homogeneous oil reservoir
Horizontal wellbore length/m	532.1
Wellbore diameter/m	0.15
Skin factor	−5.81
Wellbore storage factor/(m3/MPa)	2.49
Reservoir pressure/MPa	30.38
Test production time before shutting in/h	51.46
Shut-in bottom hole flowing pressure after test production/MPa	15.10

.

The production index is 0.093 bbl/d/psi. The verification results of Hafaya well A are shown in Figure 9 and

Table 10. Error analysis of the results of pressure recovery simulation calculation and test data.

Project	Average Relative Error (%)	Absolute Average Relative Error (%)	The Relative Standard Deviation (%)
Simulation calculation	−5.0319	5.1707	0.1733

below.

The pressure recovery simulation of three oil wells in the above two oil fields is compared with the actual pressure recovery test data. It can be seen that the new model established can well simulate the production of oil wells under different bottom hole flow pressure conditions. From the comparison, it can be seen that after the oil well is shut in, different bottom hole flow pressure conditions will produce different oil production rates. The different oil production rates cause the bottom hole flow pressure to rise and change. The law of these changes is the same as the actual test results, which verifies that the new transient productivity prediction model established is reliable and feasible.

7. Conclusions

(1) Considering the influence of formation single-phase seepage, horizontal section pressure drop, completion method, screen and formation annulus flow, establishing a single-phase horizontal wellbore variable-mass flow and formation seepage coupling coupled productivity prediction transient model, the prediction results are more in line with the actual situation. The transient model of single-phase seepage in horizontal wells is the basis for the establishment of the transient models of oil-gas two-phase seepage and oil-gas-water three-phase seepage, because we can still build a horizontal well potential model based on the principle of superposition of transient potential, the principle of mirror reflection, and the law of two-phase and three-phase permeability.

(2) The simulation calculation shows that the single-phase seepage transient model of horizontal wells can simulate and reflect the changes of horizontal well productivity with production time. That is to say, it can simultaneously reflect the two important factors that need to be considered in the horizontal well productivity prediction: horizontal wellbore pressure drop and oil reservoir production time.

(3) The model is more reliable through the verification of pressure recovery data from multiple field logs.