Introduction

The critical production of oil wells is one of the important parameters in the development process of bottom water reservoirs (Dupuit 1863; Folefac et al. 1991; Ouyang et al. 1998). When there are various impermeable or weakly permeable interlayers in bottom water reservoirs due to geological causes, the existing horizontal well productivity model of bottom water reservoirs is used to evaluate the reservoir productivity (Al-Fatlawi et al. 2019; Dheyauldeen et al. 2022; Li et al. 2018b). It is found that the calculation results of the model are quite different from the actual productivity. In addition, in order to control bottom water coning, prolong the water-free production period of oil wells, and improve oil and gas recovery, artificial interlayer technology is developed near the oil–water contact surface (Jin et al. 2010; Qin et al. 2014). However, the final yield increase effect and main influencing factors of the artificial interlayer are still fuzzy and uncertain. The main purpose of this study is to effectively predict the productivity of horizontal wells in bottom water reservoirs, establish a new productivity formula for horizontal wells in bottom water reservoirs, and determine reasonable production parameters for the development of bottom water by horizontal wells.

Horizontal well technology has been widely used around the world. The calculation methods for horizontal well productivity mainly include analytical methods and numerical simulation technology (Aziz and Ouyang 2001; Bahadori et al. 2013; Al-Fatlawi et al. 2019; Economides et al. 1991). At present, it is recognized that the Borisov formula, the Giger formula, and the Joshi formula have strong applicability. Borisov proposed a series of assumptions that horizontal wells are located in the reservoir block and established the steady-state productivity model of horizontal wells (Borisov 1984). On the basis of the Borisov formula, Giger uses an electrical model to study the reservoir problem of horizontal wells and puts forward the productivity formula at the center of the reservoir (Giger 1985). Joshi uses the superposition theory of potential to approximate the seepage problem of the horizontal well to the two-dimensional plane seepage problem in the vertical planes and horizontal planes, and obtains the steady-state productivity formula for the horizontal well (Joshi 1988). Later, many scholars put forward the horizontal well productivity formula, and the above formula is very similar.

Some scholars have shown that reservoir barriers can significantly increase the critical rate and delay the water breakthrough time, thereby increasing productivity (Huang et al. 2015; Xue and Cheng 2013). The rationale is that the main problem in developing bottom water-thick oil reservoirs is the ridge of bottom water, and the presence of an interlayer can significantly slow down the rise of bottom water (Kim and Shin 2020; Li et al. 2011). However, at the same time, it is found that the traditional horizontal well productivity model of bottom water reservoirs is used to evaluate the productivity of interlayer reservoirs, and the calculation results of this model are quite different from the actual productivity. A common point affecting these wells is that there are many extended horizontal interlayers in the reservoir. The main reason why the traditional bottom water reservoir productivity model cannot be effectively applied to the field is that the influence of interlayers on critical productivity is ignored. Hou (2007) demonstrates the genesis and classification criteria of interlayer deposition and predicts and analyzes the distribution of interlayers in different types of sedimentary sand bodies. Hou et al. (2011) used numerical simulation to discuss the influence of factors on the incremental distribution of polymer flooding, including the location, area, and permeability of the thin, low-permeability interlayer and the perforation position relative to the interlayer. Based on the theory of water coning and non-Darcy flow, Xinzhi and Shengju (2012) used the material balance principle to derive the prediction equation of the water breakthrough time of a low-permeability bottom water reservoir with a barrier under the hemispherical radial flow and the plane radial flow over the barrier. Yue et al. (2015) deduced the horizontal well model of a bottom water drive reservoir with a semi-permeable barrier considering thickness, which can be used to predict the key production parameters of horizontal wells in bottom water reservoirs with different sizes, thicknesses, and semi-permeable barrier permeability. And it is concluded that the larger the barrier size, thickness, or barrier position, the higher the key parameters of horizontal wells. Fu (2019) proposed a new equation for the critical productivity of horizontal wells by using the equivalent flow resistance method, obtained the productivity equation of the external seepage field by using Dupuit (1863) critical productivity equation, and obtained the productivity equation of the internal field by using the imaging principle. And the breakthrough time of the formation of water is predicted, and it is concluded that a low-permeability interlayer can effectively inhibit the rising speed of water and delay its breakthrough. Most studies have explored the state of reservoir bottom water based on the distribution of interlayer, and the summary of the impact on reservoir productivity is inadequate.

The advantage of the wellbore-reservoir coupling model is that it can analyze the coupling mode, mechanism, and influencing factors of wellbore pipe flow and reservoir seepage. At present, the most common coupling model is that the wellbore and reservoir are boundary conditions and interact with each other to form a complex wellbore-reservoir coupling system (Johansen and Khoriakov 2007; Penmatcha and Aziz 1999). Some of the articles focus on the study of the wellbore flow model, which is generally based on continuity and flow equations to predict wellbore pressure, flow distribution, and fluid state (Johansen and Khoriakov 2007; Khoriakov et al. 2012; Li et al. 2018a). Whether the time variable is introduced is also one of the classification conditions of the coupling model. The researchers studied a steady-state coupling model based on a one-dimensional grid connection network (Johansen and Khoriakov 2007; Ma et al. 2020). Some scholars have also introduced time variables through auxiliary equations such as continuity and energy equations to establish transient coupling models (Galvao et al. 2019; Khoriakov et al. 2012; Tang et al. 2018). In addition, the productivity index is a key parameter in determining the accuracy of the coupled model. Some scholars choose the appropriate productivity formula according to the reservoir characteristics or the functional requirements of the coupling model (Al-Kabbawi 2022; Luo et al. 2015; Tan et al. 2018). However, the productivity index is a key parameter that determines the accuracy of the coupling model. However, there is no productivity formula that can be used to calculate the productivity of horizontal wells in bottom water reservoirs with interlayers, and the optimization of completion facilities cannot be established. Therefore, the focus of this paper is to fully couple the wellbore model with the reservoir productivity model with interlayers, which is a steady-state model.

Compared with the traditional reservoir capacity equation, the innovation and reliability of this paper are reflected in the use of the two-dimensional capacity equation for horizontal wells, the mirror reflection method, and the superposition principle of potential to resolve the physical model of horizontal wells with interlayer bottom water reservoirs and establish a new model for the critical production of horizontal wells in interlayer bottom water reservoirs. Based on the node model, the wellbore is coupled with the interlayer bottom water reservoir. The influence of interlayer range and horizontal well length on oil well productivity is analyzed, and the action law of the interlayer on the horizontal well in the bottom water reservoir is summarized.

Analytical model and solution

The model, which focuses on the horizontal wells in bottom water reservoirs with an impermeable interlayer, is based on the following assumptions:

  1. 1.

    The upper and lower boundaries are closed, the horizontal direction is infinitely extended, and there is no boundary around them. The schematic diagram of a horizontal well is shown in Fig. 1.

  2. 2.

    In the wellbore, the flow is infinite, and the horizontal well is parallel to the top–bottom boundary.

  3. 3.

    Single-phase and steady-state flow are considered; the production process is isothermal and governed by Darcy’s law.

Fig. 1
figure 1

The inner seepage field of the horizontal well

Full size image

The radial flow velocity accords with Darcy’s radial flow equation. The relationship between fluid velocity and potential can be expressed as:

$${V = - \frac{K}{\mu }\nabla \Phi }$$
(1)

The relationship between potential and pressure is as follows:

$${\Phi = p + \rho gz}$$
(2)

p—Pressure of fluid; z—Vertical coordinate

Based on the fundamental principle of seepage flow mechanics, the continuity equation of a fluid can be expressed as:

$$\begin{array}{c}\nabla \cdot V=0\end{array}$$
(3)

Combination Eqs. (2) and (3), further inferred as:

$${\nabla \cdot \left( {\nabla \Phi } \right) = \nabla ^{2} \Phi = \Delta \Phi = 0}$$
(4)

At the same time, if the critical production of horizontal wells is taken as the production system, bottom water pressure equals hydrostatic column pressure \({p}_{{\text{w}}}\); the potential of bottom water is constant, and the relational formula is as follows:

$${\Phi _{{\text{w}}} = p_{{\text{w}}} + \rho _{{\text{w}}} g}$$
(5)

\({p}_{{\text{w}}}\)—Water phase pressure.

Since the effect of capillary pressure is not considered, the pressure on both sides of the oil–water interface should be equal. For any height Z, the oil–water interface satisfies the following relationship:

$${p_{{\text{o}}} \left( z \right) = p_{{\text{w}}} \left( z \right) = \Phi _{{\text{w}}} - \rho _{{\text{w}}} gz}$$
(6)

\({p}_{{\text{o}}}\)—oil phase pressure.

Therefore, the oil phase potential at the oil–water interface at any height is as follows:

$${\Phi _{{\text{o}}} \left( z \right) = p_{{\text{o}}} \left( z \right) + \rho _{{\text{o}}} gz = \Phi _{{\text{w}}} - \rho _{{\text{w}}} gz + \rho _{{\text{o}}} gz = \Phi _{{\text{w}}} - \Delta \rho _{{{\text{wo}}}} gz}$$
(7)

It can be concluded that:

Point 1, at the well point: \(z=H-{h}_{\text{p}}\), \(R={R}_{w}\), then:

$${\Phi _{{\text{o}}} = \Phi _{{\text{w}}} - \Delta \rho _{{{\text{wo}}}} g\left( {H - h_{p} } \right)}$$
(8)

Point 2, at the boundary of the diaphragm: \(z={h}_{b}\), \(R={R}_{b}\), then:

$${\Phi _{{\text{o}}} = \Phi _{{\text{w}}} - \Delta \rho _{{{\text{wo}}}} g\left( {h_{b} } \right)}$$
(9)

Point 3, outside the boundary: \(z=0\), \(R={R}_{e}\), then:

$${\Phi _{{\text{o}}} = \Phi _{{\text{w}}} }$$
(10)

The second formula of GREEN is used for the reservoir volume \(V\) occupied by oil phase. For two functions \(U\) and \(W\) defined in volume \(V\) with first-order and second-order continuous derivatives, the following is:

$${\iiint\limits_{V} {\left( {W\Delta U - U\Delta W} \right){\text{d}}V} = \iint\limits_{S} {\left( {W\frac{{\partial U}}{{\partial \mathop n\limits^{ \rightharpoonup } }} - W\frac{{\partial W}}{{\partial \mathop n\limits^{ \rightharpoonup } }}} \right){\text{d}}S}}$$
(11)

\(\mathop{n}\limits^{\rightharpoonup}\)—Unit vector; \(S\)—Surface of volume \(V\); \(\Delta\)—Laplace’s operator;

As shown in Fig. 1, along line AB, we have that:

$$\begin{array}{c}{\Phi }_{{\text{o}}}={\Phi }_{{\text{w}}}\end{array}$$
(12)
$$\begin{array}{c}{\text{ln}}\left(\frac{R}{{R}_{b}}\right)={\text{ln}}\left(\frac{{R}_{c}}{{R}_{b}}\right)\end{array}$$
(13)
$$\frac{{\partial \Phi_{{\text{o}}} }}{{\partial \mathop{n}\limits^{\rightharpoonup} }} = \frac{{Q_{{{\text{ch}}}} \mu_{{\text{o}}} }}{{2\pi K_{{\text{h}}} HR_{c} }}$$
(14)
$$\frac{\partial }{{\partial \mathop{n}\limits^{\rightharpoonup} }}\left( {\ln \left( {\frac{R}{{R_{b} }}} \right)} \right) = \frac{1}{{R_{c} }}$$
(15)
$${{\text{d}}S = 2\pi R_{c} {\text{d}}z}$$
(16)

Along line BC, we have that:

$$\frac{{\partial \Phi_{o} }}{{\partial \mathop{n}\limits^{\rightharpoonup} }} = 0$$
(17)
$$\frac{\partial }{{\partial \mathop{n}\limits^{\rightharpoonup} }}\left( {\ln \left( {\frac{R}{{R_{b} }}} \right)} \right) = - \frac{{{\text{sin }}\alpha }}{R}$$
(18)
$$\begin{array}{c}{\rm d}S=\frac{2\pi R{\text{d}}z}{\mathrm{sin }\alpha }\end{array}$$
(19)

Along line CD, we have that:

$$\begin{array}{c}{\text{ln}}\left(\frac{R}{{R}_{b}}\right)={\text{ln}}\left(\frac{{R}_{b}}{{R}_{b}}\right)=0\end{array}$$
(20)
$$\frac{\partial }{{\partial \mathop{n}\limits^{\rightharpoonup} }}\left( {\ln \left( {\frac{R}{{R_{b} }}} \right)} \right) = - \frac{1}{{R_{b} }}$$
(21)
$${{\text{d}}S = 2\pi R_{b} {\text{d}}z}$$
(22)

Along line DA, we have that:

$$\frac{{\partial \Phi_{{\text{o}}} }}{{\partial \mathop{n}\limits^{\rightharpoonup} }} = 0$$
(23)
$$\frac{\partial }{{\partial \mathop{n}\limits^{\rightharpoonup} }}\left( {\ln \left( {\frac{R}{{R_{b} }}} \right)} \right) = 0$$
(24)

Bring the integral term conversion of four surface boundaries into Eq. (12), we then have:

$${\frac{{Q_{{{\text{ch}}}} \mu _{{\text{o}}} }}{{2\pi K_{{\text{h}}} }}{\text{ln}}\left( {\frac{{R_{c} }}{{R_{b} }}} \right) = \int\limits_{0}^{H} \Phi _{{\text{o}}} \left( {R_{c} } \right){\text{d}}z - \int\limits_{0}^{{h_{b} }} \Phi _{{\text{o}}} \left( {z_{{{\text{wo}}}} } \right){\text{d}}z - \int\limits_{{h_{b} }}^{h} \Phi _{{\text{o}}} \left( {R_{b} } \right){\text{d}}z}$$
(25)

Bring the boundary conditions into Eq. (25), the yield formula of external seepage field can be obtained:

$$\begin{array}{c}{Q}_{{\text{ch}}}=\frac{\pi {K}_{{\text{h}}}\Delta {\rho }_{{\text{wo}}}g\left(2H\cdot {h}_{b}-{h}_{b}^{2}\right)}{{\mu }_{{\text{o}}}{\text{ln}}\left(\frac{{R}_{c}}{{R}_{b}}\right)}\end{array}$$
(26)

Setting oil discharge area of horizontal well as \(M\), according to Joshi horizontal well theory, the oil discharge radius is as follows:

$$R_{c} = \sqrt {{M \mathord{\left/ {\vphantom {M \pi }} \right. \kern-\nulldelimiterspace} \pi }} = \frac{{\left[ {a + \sqrt {a^{2} - {{L^{2} } \mathord{\left/ {\vphantom {{L^{2} } 4}} \right. \kern-\nulldelimiterspace} 4}} } \right]}}{2}$$
(27)

External seepage resistance is as follows:

$$\begin{array}{c}{R}_{{\text{ch}}}=\frac{{\mu }_{{\text{o}}}{\text{ln}}\left(\frac{{R}_{c}}{{R}_{b}}\right)}{2\pi {K}_{{\text{h}}}H}\end{array}$$
(28)

The length of a horizontal well is regarded as reservoir thickness. Using the principles of mirror reflection and superposition of potential to calculate the output of a borehole. It can be regarded as a production well in the upper and lower closed formations to be mapped into a producing well array. The coordinates are (\(0\), \(2nh+a\)) and (\(0\), \(2nh-a\)), where \(n\) = 0, ± 1, ….

According to the superposition principle of potential, the potential in the formation can be expressed by the following equation:

$$\begin{array}{c}\Phi \left(y,z\right)=\frac{q}{4\pi }\sum_{n=-\infty }^{+\infty }{\text{ln}}\left[{y}^{2}+{\left(z-2nh-a\right)}^{2}\right]\left[{y}^{2}+{\left(z-2nh+a\right)}^{2}\right]+C\end{array}$$
(29)

Reference formula:

$$\begin{array}{c}\sum_{-\infty }^{+\infty }{\text{ln}}\left[{\left(X-x\right)}^{2}+{\left(Y-2nh-y\right)}^{2}\right]={\text{ln}}\left[{\text{ch}}\frac{\pi \left(X-x\right)}{h}-{\text{cos}}\frac{\pi \left(Y-y\right)}{h}\right]\end{array}$$
(30)

Simplify Eq. (30):

$${\Phi \left( {y,z} \right) = \frac{q}{{4\pi }}{\text{ln}}\left[ {{\text{ch}}\frac{{\pi y}}{h} - {\text{cos}}\frac{{\pi \left( {z - a} \right)}}{h}} \right]\left[ {{\text{ch}}\frac{{\pi y}}{h} - {\text{cos}}\frac{{\pi \left( {z + a} \right)}}{h}} \right] + C}$$
(31)

Internal boundary conditions for: When \(y=0\), \(z=a-{r}_{w}\), \(\Phi \left(y,z\right)=\Phi \left(y,a-{r}_{w}\right)={\Phi }_{{\text{w}}}\). Here is the potential at the shaft wall. It is considered that \({r}_{w}\ll a\) can be calculated as:

$${\Phi _{{\text{w}}} = \frac{q}{{4\pi }}{\text{ln}}\left[ {1 - {\text{cos}}\frac{{\pi r_{w} }}{h}} \right]\left[ {1 - {\text{cos}}\frac{{2\pi a}}{h}} \right] + C}$$
(32)

By using the trigonometric relation of complex numbers, we can use these formulas:

$${1 - {\text{cos}}\frac{{\pi r_{w} }}{{2h}} \approx \frac{1}{2}\left( {\frac{{\pi r_{w} }}{{2h}}} \right)^{2} }$$
(33)
$${1 - {\text{cos}}\frac{{\pi r_{w} }}{h} \approx 2{\text{sin}}^{2} \left( {\frac{{\pi r_{w} }}{{2h}}} \right)}$$
(34)

According to Eqs. (33) and (34), above Eq. (32) can be simplified as:

$$\begin{array}{c}{\Phi }_{{\text{w}}}\approx \frac{q}{2\pi }{\text{ln}}\left(\frac{\pi {r}_{w}}{h}{\text{sin}}\frac{\pi a}{h}\right)+C\end{array}$$
(35)

External boundary conditions for: When \({r}_{e}=y\), \(\Phi ={\Phi }_{y}\), the equation is as follows:

$${\Phi _{y} \approx \frac{q}{{2\pi }}{\text{ln}}\left( {{\text{ch}}\frac{{\pi y}}{h}} \right) + C}$$
(36)

At this time, the potential difference can be expressed as:

$${\Delta \Phi = \Phi _{y} - \Phi _{{\text{w}}} = \frac{q}{{2\pi }}{\text{ln}}\left( {\frac{{{\text{ch}}\frac{{\pi y}}{h}}}{{\frac{{\pi r_{w} }}{h}{\text{sin}}\frac{{\pi a}}{h}}}} \right)}$$
(37)

Furthermore:

$$\begin{array}{c}{Q}_{{\text{cv}}}=\frac{{P}_{y}-{P}_{{\text{w}}}}{\frac{\mu y}{2{K}_{{\text{v}}}hL}+\frac{u}{2\pi {K}_{{\text{v}}}L}{\text{ln}}\frac{h}{{2\pi r}_{w}{\text{sin}}\frac{\pi a}{h}}}\end{array}$$
(38)

For the external seepage field, if the oil well is produced at a critical production rate, the production pressure difference of the oil well is as follows:

$$\begin{array}{c}\Delta p=\frac{\Delta {\rho }_{{\text{wo}}}g\left(2H\cdot {h}_{b}-{h}_{b}^{2}\right)}{2H}\end{array}$$
(39)

According to the Darcy-Seepage formula, production is calculated by linear supply. We can get the following equation:

$$\begin{array}{c}{Q}_{{\text{cv}}}=\frac{{P}_{e}-{P}_{{\text{w}}}}{\frac{uy}{2{K}_{{\text{v}}}hL}}\end{array}$$
(40)

Comparisons Eqs. (39) and (40) can be drawn: Fig. 1 can be seen as a closed reservoir. \(\frac{u}{2\pi {K}_{{\text{v}}}L}{\text{ln}}\frac{h}{{2\pi r}_{w}{\text{sin}}\frac{\pi a}{h}}\) can be regarded as internal seepage resistance.

$${R_{{{\text{cv}}}} = \frac{u}{{2\pi K_{{\text{v}}} L}}{\text{ln}}\frac{h}{{2\pi r_{w} {\text{sin}}\frac{{\pi a}}{h}}}}$$
(41)

According to the law of equivalent percolation resistance, the total seepage resistance can be expressed as:

$$\begin{aligned} R_{t} & = R_{{{\text{ch}}}} + R_{{{\text{cv}}}} = \frac{\mu }{{2\pi K_{{\text{h}}} H}}\ln \left( {\frac{{R_{c} }}{{R_{b} }}} \right) + \frac{u}{{2\uppi K_{{\text{v}}} L}}\ln \frac{h}{{2\pi r_{w} \sin \frac{\pi a}{h}}} \\ & = \frac{\mu }{2\pi KH}\left[ {\ln \left( {\frac{{R_{c} }}{{R_{b} }}} \right) + \frac{H}{L}\ln \frac{h}{{2\pi r_{w} \sin \frac{\pi a}{h}}}} \right] \\ \end{aligned}$$
(42)

The anisotropic parameters can be corrected by the Besson correction method, thus deriving a critical production equation for bottom water reservoirs with barriers.

$$\begin{array}{c}{Q}_{{\text{cL}}}=\frac{\pi K\Delta \rho g\left(2H\cdot {h}_{b}-{h}_{b}^{2}\right)}{{\text{ln}}\left(\frac{{R}_{c}}{{R}_{b}}\right)+\frac{H}{L}{\text{ln}}\frac{h}{{2\pi r}_{w}{\text{sin}}\frac{\pi a}{h}}}\end{array}$$
(43)
$$R_{c} = \sqrt {{M \mathord{\left/ {\vphantom {M \pi }} \right. \kern-\nulldelimiterspace} \pi }} = \frac{1}{2}\left[ {a + \sqrt {a^{2} - {{L^{2} } \mathord{\left/ {\vphantom {{L^{2} } 4}} \right. \kern-\nulldelimiterspace} 4}} } \right]$$
(44)

Basic model

Network structure

Flow in a horizontal wellbore is not equal mass flow, but variable mass flow. The flow of fluid from the reservoir continues to flow into each section of the horizontal well, which increases the flow rate and quality of the fluid along the direction of the horizontal well. Therefore, the flow in the horizontal wellbore is more complex than the conventional pipe flow, and various calculations are more complex.

By coupling the flow of the reservoir and horizontal well together and establishing a unified coupling equation, the influence of the flow rate of fluid flowing into the horizontal well and the pressure in the horizontal well on reservoir pressure can be correctly analyzed. The horizontal segment of the completion model system is composed of nodes and flow bridge segments. As shown in Fig. 2, the uppermost row of nodes represents the reservoir, the middle row represents the annulus, and the lower most row represents the tubing. Standing, Vasquez-Beggs, Marhoun, and other methods can be selected to calculate the viscosity, density, volume coefficient, and other parameters. So, the specific calculation formulas for these parameters are no longer detailed.

Fig. 2
figure 2

Node network model system

Full size image

Mass conservation equation and momentum conservation equation

The basic principle of the model is the discretization method. Fluid flows from the reservoir radially into the horizontal wellbore and then flows from the toe to the heel. The fluid mass flow in the wellbore increases gradually, which is a wellbore variable mass flow. Restriction relations are formed by mass conservation law and momentum conservation law between nodes. The flow diagram is shown in Fig. 3.

Fig. 3
figure 3

Micro-volume diagram

Full size image

Mass conservation equation

Fluid flow is regarded as steady flow, and there is no mass aggregation at nodes. Requirement of boundary conditions that satisfy the principle of mass conservation. So, the total amount of fluid input into the node equals the total amount of fluid output from the node. The outflow and inflow of each node can be expressed as:

$$\begin{array}{*{20}c} {\sum\limits_{{i = 1}}^{n} {m_{i} = 0} } \\ \end{array}$$
(45)
Momentum conservation equation

The momentum conservation equation between nodes is the calculation model of the pressure drop between pipes.

  1. 1.

    Productivity formulas

The flow of oil and gas from the reservoir to the annulus is governed by an equation of the form:

$${q_{{r,i}} = {\text{PI}}*\Delta p_{{r,i}} }$$
(46)

Therefore, the pressure loss in the reservoir can be calculated by:

$${\Delta p_{{r,i}} = {{{\text{PI}}} \mathord{\left/ {\vphantom {{{\text{PI}}} {q_{{r,i}} }}} \right. \kern-\nulldelimiterspace} {q_{{r,i}} }}}$$
(47)
  1. 2.

    Frictional pressure drop

The friction factor \(f\) can be expressed by using the Reynolds number under laminar flow and turbulent flow:

Here, the continuous-phase Reynolds number is defined through:

$$R_{e} = {{\rho vd} \mathord{\left/ {\vphantom {{\rho vd} \mu }} \right. \kern-\nulldelimiterspace} \mu }$$
(48)

The friction factor calculation method for radial inflow generally adopts the recent research of Ouyang et al. The calculation process is as follows:

If \({R}_{e}\le 2300\), the wall friction in laminar flow is calculated by:

For the inflow case:

$$\begin{array}{c}f=\frac{16}{{R}_{e}}\left(1+0.04304{R}_{ew}^{0.6142}\right)\end{array}$$
(49)

For the outflow case:

$$\begin{array}{c}f=\frac{16}{{R}_{e}}\left[1-0.0625\frac{{\left(-{R}_{ew}\right)}^{1.3056}}{{\left({R}_{ew}+4.626\right)}^{-0.2724}}\right]\end{array}$$
(50)

If \({R}_{e}\ge 2300\), the wall friction in turbulent flow is calculated by:

For the inflow case:

$$\begin{array}{c}f=\frac{16}{{R}_{e}}\left(1-0.0153{R}_{ew}^{0.3978}\right)\end{array}$$
(51)

For the outflow case:

$$\begin{array}{c}f=\frac{16}{{R}_{e}}\left[1-17.5\frac{{R}_{ew}}{{R}_{e}^{0.75}}\right]\end{array}$$
(52)

The pressure loss of wellbore based on Ouyang has been used to calculate the pressure loss of tubing by:

$$\begin{array}{c}{\Delta p}_{f}=\frac{{f}_{t}\rho L}{2{D}_{{\text{in}}}}{\left(\frac{{q}_{t,i}}{{A}_{t}}\right)}^{2}\end{array}$$
(53)

There is a flow area of casing-tubing annulus for completion wellbore with a casing, and the pressure loss of casing-tubing annulus is as follows:

$$\begin{array}{c}\Delta {p}_{f}=\frac{{f}_{a}\rho L}{2{(D}_{{\text{out}}}-{D}_{{\text{in}}})}{\left(\frac{{q}_{a,i}}{{A}_{a}}\right)}^{2}\end{array}$$
(54)
  1. 3.

    Choke pressure drop

Each well completion has its characteristics, which are reflected not only in different internal structures but also in different flow forms. Fluid flow from the annulus to the tubing is considered a radial flow model. The pressure difference formula of the nozzle throttling model is as follows:

$$\begin{array}{c}\Delta {p}_{a-t,i}=255.16\frac{\rho {q}_{a-t,i}^{2}}{{C}_{d}^{2}{d}_{e}^{4}}\end{array}$$
(55)

The calculation method for multiple nozzles with different diameters:

$$\begin{array}{c}{d}_{e}=\sqrt{{d}_{1}^{2}+{d}_{2}^{2}\dots +{d}_{n}^{2}}\end{array}$$
(56)

For slotted tubular, the formula for throttling pressure drop can be calculated as follows:

$$\begin{array}{c}\Delta {p}_{a-t,i}=0.0036\frac{\rho {q}_{a-t,i}^{2}}{{2C}_{d}^{2}{A}_{s}^{2}}\end{array}$$
(57)

Solution of mathematical model

The complete base solver model is now given by Eqs. (43) (45) (47) (54) and (57), and the equations associated with this as described above. The model control equation is a system of nonlinear equations. The iterative method is used to solve the nonlinear equations with high convergence speed, so this method is selected to solve the network model. The technical roadmap for solving the model is as follows (Fig. 4):

Fig. 4
figure 4

Workflow and method for the resolution of the coupling model

Full size image

Applications

Model validation

The purpose of the case is to solve the pressure field and flow field. Analyzing the influence of different parameters on productivity can provide technical support for the study of downhole oil and gas dynamics. This study verifies the calculation results of the developed model through a case study of the Tarim oilfield. The data in Table 1 are used for model calculation to obtain relevant data, such as inflow rate distribution along the well, wellbore pressure distribution, and capacity prediction, and compare them with other models, as shown in Figs. 5, 6 and 7.

Table 1 Basic reservoir and fluid parameters
Full size table
Fig. 5
figure 5

Comparison of pressure along the well for various productivity models

Full size image
Fig. 6
figure 6

Comparison of inflow rate distribution along the well for various productivity models

Full size image
Fig. 7
figure 7

Comparison of pipe flow along the well for various productivity models

Full size image

The Joshi's model and the new capacity model belong to the mechanism model. Their advantage is that they can clearly reveal the physical processes and mechanisms within the system. Eclipse is a numerical simulation software that can handle complex physical phenomena and geometric structures. By comparing these two methods, the changes in pressure field and flow field after considering interlayer factors are analyzed.

As shown in Figs. 5, 6 and 7, in order to better compare the accuracy of the established model, compare and analyze it with eclipse and traditional productivity models, and observe the flow and pressure along the well. As can be seen from Fig. 5, the pressures of the traditional capacity model, eclipse simulation, and the new capacity model at the toe end are 34.0632 MPa, 34.0507 MPa, and 34.0454 MPa, respectively. Take the simulation results of commercial software Eclipse as the standard. The relative error of the traditional productivity model is 4.35–20.39%, and the relative error of the capacity model established in this paper is 2.35–11.81%. Figure 6 shows the comparison of inflow distribution along the well under different productivity models; Fig. 7 shows the comparison of pipe flow along the well for various productivity models. The inflow of the traditional capacity model, Eclipse simulation, and the new capacity model at the heel end are 1189.74 m3/d, 1753.41 m3/d, and 1604.77 m3/d, respectively. The relative error of the traditional productivity model is 36.67%, and the relative error of the capacity model established in this paper is 8.56%. These results show that the model established in this paper is more compatible with the commercial software Eclipse simulation.

Effect of horizontal well length

Based on the above data, change the single variable combined with the reservoir wellbore coupling model to explore the influencing factors. Figure 8 shows a linear relationship between the length of a horizontal well and its production. With the increase in the length of horizontal wells, the yield increases, but the yield per unit length decreases gradually. It indicates that for a certain reservoir, the horizontal section length of the horizontal well has an optimal length. These objective laws should be used as important reference materials in pursuit of economic benefits.

Fig. 8
figure 8

Variation curves of horizontal well production under different wellbore lengths

Full size image

It can be seen from Fig. 8 that there are significant differences in the total inflow of horizontal wells with different lengths. When the length of the horizontal well is constant, the productivity index of the new model is greater than that of Joshi's model, which is in line with the view of scholars that the interlayer can effectively block the rise of bottom water and improve productivity. It can be seen that the productivity index PI of horizontal wells increases with the length of the wells. At the same time, it can be observed in the new PI model that when the horizontal well length is less than 300 m, the productivity index PI increases significantly with the increase in horizontal well length; when the horizontal well length is more than 300 m, the productivity index PI increases slowly with the increase in horizontal well length. Therefore, in the process of horizontal well design, it is important to consider the reasonable length of the horizontal section for the efficient application of horizontal well technology.

The total pressure drop in a horizontal wellbore is the pressure loss caused by friction, acceleration, gravity, and other factors. In the long horizontal well section, the pressure drop in the wellbore is relatively large, resulting in a large production pressure difference and flow rate at the heel of the horizontal well and a small production pressure difference and flow rate at the toe of the horizontal well, which makes the long horizontal well section not effectively used. To further explore the relationship between horizontal well length and productivity, the wellbore heel pressure is set to 35 MPa, and the length of horizontal wells is set to 200 m, 300 m, 400 m, 500 m, and 600 m, respectively.

The flow distribution in the reservoir inflow, the flow distribution in the pipeline, and the pressure distribution in horizontal wells of different lengths are plotted as shown in Figs. 9, 10, and 11, respectively. It can be observed from the following figures: (1) The pressure values along the horizontal well obtained by the two models are close by coupling calculation; (2) the inflow profiles of both models are asymmetric in the calculation results. The Joshi's model has inflows of approximately 29.88 m3/d and 70.51 m3/d at the toe and heel of the well, respectively, while the new model has inflows of approximately 4.35 m3/d and 245.74 m3/d at the toe and heel of the well, respectively. It can be seen from the above that the use of different productivity models will lead to a large difference in the flow distribution of reservoir inflow, resulting in a difference in the cumulative flow along the tubing. Therefore, it is very important to adopt a productivity model suitable for an interlayer-reservoir.

Fig. 9
figure 9

Comparison of inflow rate distribution along the well for various well length

Full size image
Fig. 10
figure 10

Comparison of pipe flow along the well for various well length

Full size image
Fig. 11
figure 11

Comparison of pressure along the well for various well length

Full size image

Effect of interlayer radius

Through the parameters given in the above case, observe the influence of interlayer radius on reservoir productivity, and the dynamic curve is shown in Figs. 12, 13 and 14. As shown in Fig. 13, when the partition lengths are 50 m, 100 m, 150 m, and 200 m, respectively, the flow in the pipe at the heel end of the well is 617.16 m3/d, 895.29 m3/d, 1211.29 m3/d, and 1604.77 m3/d in sequence. The results show that the production of horizontal wells increases with the increase in interlayer radius under the same production pressure difference. There are two main reasons: (1) The existence of an interlayer inhibits bottom water coning; (2) the non-permeable interlayer used in this paper makes the bottom water flow around the interlayer edge to form edge water flooding or secondary bottom water drive.

Fig. 12
figure 12

Comparison of inflow rate distribution along the well for various interlayer radii

Full size image
Fig. 13
figure 13

Comparison of pipe flow along the well for various interlayer radii

Full size image
Fig. 14
figure 14

Comparison of pressure along the well for various interlayer radii

Full size image

Figure 12 describes the inflow along the well with different interlayer radii. As can be seen from Fig. 12, when the partition lengths are 50 m, 100 m, 150 m, and 200 m, respectively, the inflow flows of the reservoir at the heel end of the well are 198.95 m3/d, 220.44 m3/d, 234.91 m3/d, and 245.74 m3/d in sequence, and those at 15 m away from the heel end are 107.71 m3/d, 140.71 m3/d, 167.27 m3/d, and 190.19 m3/d in sequence. At the same time, Fig. 14 depicts the effect of different interlayer radii on pressure along the well. When the interlayer radii are 50 m, 100 m, 150 m, and 200 m, respectively, the wellbore pressure at the toe end is 34.0040 MPa, 34.0107 MPa, 34.0238 MPa, and 34.0504 MPa, respectively. Affected by the wellbore pressure loss, the wellbore radial flow decreases gradually from heel end to toe end, and the greater the loss of wellbore pressure, the greater the difference in flow rate between heel and toe.

Conclusions

The use of appropriate capacity models for particular reservoirs is important in evaluating pressure and flow along horizontal wells. Taking the interlayer into account, this work derives the reservoir productivity model and the wellbore-reservoir coupling model. The effects of horizontal well length and interlayer radius on productivity were evaluated. Some condensed conclusions are as follows:

  1. 1.

    This study deduces the bottom water reservoir productivity model with an interlayer. The productivity of this model is quite different from that of the traditional production model, indicating that the interlayer factors in the reservoir cannot be ignored. Compared with the simulation results of the commercial simulator (ECLIPSE), the results are in good agreement, indicating that the model has potential practical application value.

  2. 2.

    In this paper, an analytical model coupling the bottom water reservoir with an interlayer and wellbore is established. The sensitivity analysis of productivity coefficient, horizontal well length, and interlayer length is carried out using the coupling model. The analytical model can not only be used as a basis for predicting the inflow profile of horizontal wells in bottom water reservoirs but also as a basis for horizontal well completion design and production prediction.

  3. 3.

    This study only considers the case of a single interlayer. In fact, the number of interlayers and their distribution in the reservoir will affect the reservoir productivity and the flow dynamics of fluid from the reservoir to the horizontal well. We will continue to introduce this part in future research.