A Semianalytical Coupling Model between Reservoir and Horizontal Well with Different Well Completions

Although currently, large-scale and multilateral horizontal wells are an important way to improve the oil recovery in the unconventional reservoirs, the ﬂ ow behavior of ﬂ uid from the reservoir into the horizontal wellbore becomes more challenged compared to the single small-scale horizontal well. One of the main challenges is that pressure loss from the well completion section and wellbore cannot be ignored in the coupling process between the reservoir and the horizontal well. In this paper, a new method is presented to solve the coupling ﬂ ow between the reservoir and the horizontal well with di ﬀ erent well completions. The new coupling model is compared with Ouyang ’ s model (1998) and Penmatcha ’ s model (1997), and the predicted data are consistent with each other at both early and late times. Meanwhile, four di ﬀ erent cases have been proposed to verify the application of the new coupling model with di ﬀ erent well completions, and the results indicate that the uneven in ﬂ ow pro ﬁ le can be e ﬀ ectively alleviated via reasonable completion parameters and di ﬀ erent well completions. Based on two types of ﬂ ow-node units, it can quickly model and solve the coupling problem between the reservoir and the horizontal well with complex completion cases. It can also depict the in ﬂ ow pro ﬁ le of the horizontal well with di ﬀ erent well completions, which is conducive to understand the coupling process. The new coupling model can provide theoretical support for further optimization of completion parameters and well completions and ﬁ nally improve oil recovery.


Introduction
The horizontal well technology has been applied as early as the 20th century. Nowadays, it has gradually developed into an effective technology for the oil field to create huge economic benefits. Joshi is the first one to propose the analytical model on horizontal well productivity [1]. Later, Giger and Francois, Renard and Dupuy, and Hongen improved the analytical model based on his model [2][3][4]. Although it is easy to solve the analytical models, they still suffer some limitations. One of the limitations is that they cannot describe the influences of the transient flow and boundary effect, especially in the complex reservoirs. Therefore, a method that can solve this problem is urgently needed.
The numerical simulation method is such a compelling method that offers several advantages over the analytical models. Vicente et al. developed a fully implicit and threedimensional simulator to simulate the transient pressure and flow rate behavior using the numerical simulation method [5]. Subsequently, Oliveira studied the coupling progress between the reservoir and the horizontal well via the numerical simulation method, and a three-dimensional model was established [6]. Gui et al., Akim et al., and Wilson used the numerical simulation method to study the coupling problem in unconventional oil-gas recovery [7][8][9]. Wang and Leung studied the effects of capillarity and geomechanics on water loss in the fracture-matrix system, and the coupling models were constructed to simulate the multiphase flow and fluid distribution during shut-in and flowback processes [10].
Although the numerical simulation method provides various advantages to solve the flow characteristic of the reservoir, its accuracy and solution speed are directly restricted by the size of the grid. Therefore, semianalytical models are a more practical method that can solve as quickly as the analytical model and also simulate the instantaneous flow and boundary effect. Dikken established the first semianalytical coupling model between the reservoir and the horizontal well on calculating productivity [11]. Islam and Chakma, Folefac et al., Ozkan and Sarica, and Seines et al. set up different coupling models between the reservoir and horizontal well [12][13][14][15]. However, these semianalytical models only considered the influence of friction pressure loss in the horizontal wellbore but ignored the pressure loss caused by well completions and other factors such as wall friction.
Later, Ouyang et al. studied the flow rule of single-phase or multiphase in the horizontal wellbore and improved the model to solve the pressure-loss issue in the horizontal wellbore [16][17][18]. Penmatcha and Aziz established a semianalytical model and considered the pressure loss caused by wall friction and acceleration [19]. Based on Green's functions, Valvatne et al. established a semianalytical model on complex well configurations with intelligent completions and solved the pressure loss in the annulus, tubing, and downhole chokes [20]. Fokker and Verga considered the well interference effects and established a semianalytical model to calculate the productivity of vertical, horizontal, or multilateral wells [21]. Valko and Amini [22,23] Jiang et al. [24], and Sheng et al. [25] proposed different semianalytical models to predict production performance of the horizontal wells located in a fractured reservoir, and high accuracy and computational efficiency were obtained. Li et al. presented a semianalytical model to predict the inflow profile of horizontal wells in a bottom-water gas reservoir. The near-well skin is used to simulated the effect of completion on the coupling process in his model [26]. Luo et al. established a semianalytical model to calculate the productivity index of a horizontal well with pressure loss along the wellbore [27]. Yue et al. prospered a coupling model for estimating the water breakthrough time of complex structured multilateral horizontal wells. The pressure loss caused by well completion is not considered in the model [28]. However, all the mentioned semianalytical models oversimplify the horizontal well as a single channel that cannot describe the characteristics of fluid flow in different well completions. Luo et al. developed a semianalytical model for predicting the performance of horizontal wells completed by inflow control devices in bottom-water reservoirs. The model considers pressure drops due to fluid flow in the tubing, ICDs, and annulus, but it is only suitable for the single completion (just ICDs) and cannot solve the coupling problems with different well completions [29]. Therefore, in this paper, a new semianalytical coupling model between the reservoir and the horizontal well with different well completions is proposed using the idea of node analysis. Based on two types of flow-node units, a flownode network with different well completions is constructed. This method is suitable for the coupling problem between the reservoir and the horizontal well with different well completions and also can be solved efficiently. Meanwhile, it can more accurately simulate the coupling process between the reservoir and the horizontal well and improve the accuracy of production prediction in the horizontal well.

Description of the Coupling Mechanism between Reservoir and Horizontal Well
The coupling is a process of interaction and energy transmission between the reservoir and the horizontal well. The coupling process includes two different types of flow forms: seepage in the reservoir and flow in the wellbore. Both of them can lead to different response speeds in their flow fields and finally make changes in the entire system. For the coupling system, any change in production parameters will cause dramatic changes in the flow field of the wellbore, while the reservoir is the opposite. Moreover, the well completion also adds the pressure loss to the coupling process and changes the fluid flow behavior in the wellbore, which makes the coupling process more complex. Meanwhile, the coupling process cannot be solved directly, because this solving process will increase not only the difficulty but also the complexity of the problem. Therefore, the entire coupling system needs to be discretized for simplification. The node analysis method is adopted in this paper. As shown in Figures 1 and 2, the coupling system is divided into different flow sections according to the flow characteristics. The inflow from the well toe is not considered. The box represents the pressure loss between the nodes, and the node's location represents different flow sections. The framework of the coupling process is built up based on the flow process. And it is equivalent to change the whole coupling system into a series and parallel "pipe network." That is, the nodes in the same fluid layer are connected in series with each other and parallel with nodes in other flow layers. This can show the entire system more clearly and keep the results undistorted.

Geofluids
Currently, the completion technologies used in horizontal wells are open-hole, perforation, prepacked screen, slotting liner, etc. They can be divided into two major categories according to with/without the casing. For completion without the casing, according to its flow characteristics, the corresponding flow-node unit is extracted and recorded as the flow-node unit of type I in Figure 1. The fluid flow in the annulus between the wellbore and the completion string is not considered. Similarly, the flowing-node unit of completion with the casing is recorded as the flow-node unit of type II in Figure 2. Due to the casing separation, the fluid will first flow into the casing-tubing annulus, then enter the tubing through the well completion. Therefore, the flow-node unit of type II has an additional flow layer of casing-tubing annulus compared with that of type I. The coupling model between the horizontal well and the reservoir is established based on the two types of flow-node units.

Mathematical Models for Reservoir
As shown in Figure 3, there is a horizontal well in the boxed reservoir with coordinates from ðx 0 , y 1 , z 0 Þ to ðx 0 , y N+1 , z 0 Þ. The line source model based on Babu et al. (1989) has been used to evaluate pressure loss in the reservoir. The superposition principle and Duhamel's principle have been used to deal with the effects of internal and external boundaries in the model [30]. Assumptions: (1) The length, width, and thickness of the boxed reservoir are b, a, and h, respectively (2) The reservoir is isothermal, and temperature distribution is uniform throughout the reservoir The porous medium is anisotropic, and permeability in the three directions is K x , K y , and K z (4) Flow in the porous medium is governed by Darcy's law (5) The fluid is slightly compressible and single-phase oil Then, pressure loss at any position ðx, y, zÞ in the reservoir with all six closed boundaries in a horizontal wellbore at time t is as follows:

Geofluids
Discretizing Equation (1) at the kth time layer, the time step is set to Δt.
Therefore, Equation (6) can be deformed into the following equation.
Based on Van Everdingen and Hurst's method, the additional pressure caused by completion, formation damage, and rock crushing around each wellbore segment in the near-wellbore area is incorporated into the line source model using the expression below [31]: Therefore, the pressure loss in the reservoir can be calculated by

Mathematical Model for Wellbore
In the actual production process, except for the wellbore with open hole completion, the fluid in the reservoir flows to the casing-tubing annulus firstly and then enters the tubing through well completions or flows to the tubing through well completions directly. Therefore, the pressure loss of the wellbore includes pressure loss of well completions, casingtubing annulus, and tubing. The pressure loss of wellbore based on Ouyang has been used to calculate the pressure loss of tubing by [17] There is a flow area of casing-tubing annulus for completion wellbore with a casing, then pressure loss of casingtubing annulus is as follows: Ouyang (1998) has discovered that the parabolic velocity profile is inappropriate for describing the velocity distribution over the cross-section with mass transfer due to inflow and outflow. The wall friction increases for the inflow case but decreases for the outflow case when fluid flows in laminar flow, and the wall friction in turbulent flow is the opposite [17]. 4 Geofluids The wall friction in laminar flow is calculated by Each well completion has its characteristics, which are reflected not only in different internal structures but also in different flow forms.
The commonly used characteristic relationship of well completions is shown as follows. The perforation pressure loss is calculated by EI-Rabba et al. in [32] Δp c = 8ρ The pressure loss of fluids through the screen is calculated by Azizi in [33] The pressure loss of fluids through the slotting liner is calculated by Huang and Wang in [34] For simplicity, the flow rule can be described by the characteristic relationship of well completions, that is, Δp c = f ðQÞ .

Mathematical Models for the Coupling Progress between Horizontal Well and Reservoir with Different Well Completions
The coupling progress between the reservoir and the horizontal well with different well completions can be divided into the horizontal well, well completion, and reservoir. These three sections can be calculated separately and then coupled together in the near-well area according to the continuity of pressure and flow. If the horizontal wellbore is taken as the research object, the reservoir section provides an external boundary condition for it and vice versa. To study the influence of different well completions on the coupling process between the reservoir and the horizontal well, this paper takes the wellbore as the research object and elaborates the flow law in the completion section of the coupling process between the reservoir and the horizontal well.

Mathematical Models for the Two Types of Flow-Node Unit
(1) The mathematical model for the flow-node unit of type I Figure 4 shows the flow-node unit of type I, which is naturally sorted according to the flow direction. The sorting rules of the rest layer are the same, except that the subscripts are different. The subscript r indicates the reservoir node, the subscript c indicates the completion node, the subscript t indicates the tubing node, and the subscripts i, j, and k indicate the serial numbers, respectively.
The nonlinear equations of the flow-node unit of type I are as follows: (2) The mathematical model for the flow-node unit of type II Figure 5 shows the flow-node unit of type II, and the sorting rule is the same as that of type I. The main differences are as follows: first, the type II has one more flow layer (annulus) than type I, which is represented by the subscript a; second, since the fluid in the annulus flows not only to the next annulus node but also into the pipe node through well completions, there is a progress of flow distribution. Therefore, the Figure 4: Flow-node unit of type I. 5 Geofluids coefficients of flow distribution α n and β n are introduced, and α n + β n = 1, n = i, i + 1, ⋯, j.
The nonlinear equations of the flow-node unit of type II are as follows: First is discretization. According to the length of the complete section and the well completions at different positions, the horizontal section is divided into several microelement segments, and the number of nodes is set. The node position can be arbitrarily selected, and different nodes can make node refinements or nonrefinements according to the research when performing node division.
Second, we select the flow-node units. We can select the corresponding types of flow-node units according to the well completions and draw the coupling structure diagram. Then, the nodes are sorted, and the natural ordering along the flow direction is convenient to solve. That is, the node of the horizontal well at the toe is recorded as node 1, and the nodes are sorted till the heel of the horizontal well. Other sorting methods are also available, as long as the solution variables and sorting numbers are kept consistent.
Third, we establish the coupling model. The nonlinear equations are connected with the coupling structure diagram, that is, the coupling model, as shown in Equation (27). For a single well completion, we can simply select the nonlinear equations and boundary conditions of the corresponding flow-node units.
The first well completion, , The junction between the first and the second well completion, The mth well completion, p aj p ai p a (j-1) Figure 5: Flow-node unit of type II. 6 Geofluids where It can be seen from the steps of establishing the coupling model between the reservoir and the horizontal well that equation (27) is a set of nonlinear equations composed by equations of different types of flow nodes, junction equations between different completion modes, and boundary conditions. Simply substituting the pressure loss of different flow layers, the coupling model of corresponding well completions can be obtained.

The Solution of Mathematical Coupling
Model. The coupling progress between the reservoir and the horizontal well is a typical "grey-box" problem, that is, the flow rule of the coupling progress is known, but the flow field distribution of the coupling progress is unknown. The coupling model established above is exactly based on the flow rule of the coupling process, and it is the coupling relationship among the flow variables at different positions in time. To solve the "grey-box" problem, the essence is to use the flow law to solve its distribution law, which is to use the obtained conditions to solve the established nonlinear equations. It can be known from the coupling model that it is a function of the pressure loss of different flow layers and the inner and outer boundary conditions of the near-wellbore area. The Newton-Raphson method is adopted to solve the problem in this paper. The detailed solution of the nonlinear equations (Equation (27)) is as shown in the following equation (Equations (29) and 30). We can set the time interval and loop to get the solution of the coupling model at different times.  [17,19]. A horizontal well is centered in a boxed reservoir. Other parameters of the reservoir and the fluid are listed in Table 1. Figure 6 displays a comparison of inflow rate distribution by the new coupling model, Penmatcha's model (1997), and Ouyang's model (1998) at different times. It indicates that those three coupling models give the identical inflow rate Δp ti − p wf = 0, Under given flowing bottom-hole pressure, Under given production rate: distribution at both early times (t = 0:1 days) and late times (t = 1000 days), demonstrating the accuracy of our model. Although both Penmatcha's (1997) and Ouyang's (1998) models can give an accurate prediction on single well completion, they cannot handle different well completions. To verify the application of our new coupling model in different well completions, we will discuss four different cases later.

Application
Four different cases have been proposed to verify the application of the new coupling model: (1) completion with uniform perforation density, (2) completion with different perforation densities, (3) segregated well completion, and (4) different well completions. A horizontal well is centered in a boxed reservoir, and the data listed in Table 1 is used for all four examples.

Case 1.
Completion with uniform perforation density: in this case, a horizontal is perforated with uniform perforation density thought the whole wellbore. Specific well completion parameters are listed in Table 2. Figure 7 displays the inflow rate distribution of completion with uniform perforation density at different times. It indicates that the inflow profile presents a strong asymmetry due to mass flow in the tubing. In the early stage, the influence of the near-well area plays a major role, and the inflow profile changes obviously at these times. With continuous production, the pressure wave propagates to the boundary. The influence of the boundary is gradually enhanced and plays a leading role. When the production well enters the pseudo-steady state, the inflow profile tends to be gentle over time.  Table 3. Figure 8 displays the inflow rate distribution of completion with different perforation densities at different times. It indicates that the magnitude and range of the inflow profile are smaller due to the reduction of the perforation density     Figure 7, and the asymmetry of the curve is relieved. The reason is that the closer to the heel, the greater the resistance of the fluid which will enter the wellbore with a smaller perforation density. The toe is exactly the opposite.
Case 3. Segregated well completion: in this case, a horizontal well is divided into three completion areas with a slotted liner and connected with a blank pipe. Fluid in the reservoir can only flow into the wellbore through the slot. Specific well completion parameters are listed in Table 4. Figure 9 displays the inflow rate distribution of segregated well completions at different times. It indicates that inflow rate distribution forms a "U-shape" in each completion area, and the curve moves up over time, because the completion area divided by the blind pipe acts independently as "horizontal wells" and is coupled under the effect of well interference with each other at the same time.
Case 4. Different well completions: in this case, a horizontal well is completed alternately by perforating and slotted liner completions from heel to toe. Perforation parameters are listed in Table 2, and the slotted liner parameter is listed in Table 3. Divided into four completion areas in the horizontal direction, specific well completion parameters are listed in Table 5. Figure 10 displays the inflow rate distribution of different well completions at different times. Compared with perforation inflow, horizontal slotting makes the contact surface of the wellbore and the reservoir larger per unit length. So, there is a significant difference between the inflow rates of the two well completions. It indicates that the inflow profile changes stepwise, which can clearly distinguish different completion methods.
Comparison of the four cases: the inflow rate distribution of the four cases is compared to evaluate the effect of alleviating uneven inflow rate distribution at both early times (t = 0:1 days) and late times (t = 1000 days). Figure 11 displays the inflow rate distribution of the four cases at both early times (t = 0:1 days) and late times (t = 1000 days). It indicates that the inflow rate distribution of completion with uniform perforation density is the most uneven, in which the inflow rate varies greatly from the heel to the toe. The inflow rate distribution of segregated well completion is effectively relieved, but the inflow rate changes significantly during each completion section at late times. The mitigation effect of completion with different perforation densities is slightly better than that of different well completions, but they are inferior to segregated well completion.

Conclusions
The coupling problem between the horizontal well and the reservoir is transformed into the interaction relationship between different flow sections using node discretization. According to the flow characteristics of different well completions, the corresponding flow-node type is extracted. And the new coupling model is established based on the two types of flow-node units. The model can solve the coupling problem on not only single well completion but also different well completions.
Comparison among the predicted data obtained from the new coupling model, Ouyang's model (1998), and Penmatcha's model (1997) reveals that the new model has a potential practical application in predicting the performance of horizontal wells with different well completions.
Case studies of four well completions are carried out in this paper, and the inflow rate distribution of the four cases are compared to evaluate the effect of alleviating uneven inflow rate distribution. This work can provide a guideline on optimizing completion parameters, which are beneficial to alleviate the uneven inflow distribution, and finally can improve oil recovery. The friction factors f a : The friction factor of the casing-tubing annulus f t : The friction factor of tubing h: The thickness of the boxed reservoir (m) i: Serial number j: Serial number K x : Permeability in the x-axis direction (D) K y : Permeability in the y-axis direction (D) K z : Permeability in the z-axis direction (D) K s : Pressure loss coefficient of screen (D) k: Serial number k w : Effective reservoir permeability = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K x K y K z 3 p (D) L: The length of the horizontal well (m) L s : The length of slotting (m)