© M. Wang et al., published by IFP Energies nouvelles, 2020

Nomenclature

x : Distance in the x-axis, m

y : Distance in the y-axis, m

x_e : Reservoir length in the x-axis, m

y_e : Reservoir length in the y-axis, m

h : Reservoir thickness, m

L_f1 : Principal fracture length, m

L_f2, L_f3: Reoriented fracture length, m

L_R : Reference length, m

l : Discrete segment length, m

w_f : Fracture width, m

k_x : Permeability in the x-axis, 10⁻³ μm²

k_y : Permeability in the y-axis, 10⁻³ μm²

$\bar{k}=\sqrt{k_x{k}_y}$ : Equivalent system permeability, 10⁻³ μm²

k_f : Fracture permeability, 10⁻³ μm²

θ₁ : Principal fracture angle, degree

θ₂, θ₃: Reoriented fracture angle, degree

p_i : Initial reservoir pressure, MPa

p_avg : Average reservoir pressure, MPa

p_w : Wellbore pressure, MPa

φ : Reservoir porosity, fraction

μ : Fluid viscosity, mPa s

c_t : Total compressibility, 1/MPa

β : Permeability anisotropic factor

Q : Wellbore flow rate, m³/d

q_fw : Flow rate of per unit fracture length, m³/d

q_f : Flow rate strength, m²/d

p_D : Dimensionless pressure in real time domain

${\tilde{p}}_D$ : Dimensionless pressure in Laplace domain

p_avgD : Dimensionless average reservoir pressure

p_wD : Dimensionless wellbore pressure

q_fD : Dimensionless flow rate strength

q_D : Dimensionless point source flux

J_D : Dimensionless productivity index

J_Dmax : Maximum dimensionless productivity index

C_fD : Dimensionless fracture conductivity

l_D : Dimensionless discrete segment length

t_D : Dimensionless time

x_D : Dimensionless distance in the x-axis

y_D : Dimensionless distance in the y-axis

x_eD : Dimensionless reservoir length in the x-axis

y_eD : Dimensionless reservoir width in the y-axis

F_D : Dimensionless influence function

σ : Reoriented factor

s : Dimensionless time variable in Laplace domain

χ : Dimensionless integral variable

ψ_n : Dimensionless coefficient, n = 0, 1, 2……

N_i : Total discrete segments of the ith reorientation fracture, i = 1, 2, … N

N : Total number of reorientation fractures

Sum: Total numbers of discrete segments

△l_Di: Dimensionless length of the ith discrete segment, i = 1, 2, 3…

sinh: Hyperbolic sine function

cosh: Hyperbolic cosine function

1 Introduction

Large scale hydrocarbons produce from conventional reservoirs, but unconventional reservoirs still serve as the reliable formations to obtain oil and gas, such as tight reservoirs and shale reservoirs [1–3]. Hydraulic fracturing horizontal well is one of the most effective approaches to enhance well productivity in tight or low-permeability reservoirs, where non-fractured wells in these reservoirs are of no economic value [4]. Recently, most researchers have focused on pressure transient analysis of tight reservoirs intercepted by planar fractures, but they don’t pay much attention on the analysis of pseudo-steady-state productivity for such reservoirs. In general, hydraulic or volume fracturing can form reorientation fractures in tight reservoirs. However, there is a lack of a detailed research on calculating pseudo-steady-state productivity for a horizontal well with multiple reorientation fractures in a tight reservoir, which is exactly the original intention of this study.

For productivity analysis and evaluation, many mathematical models have been extensively reported to calculate the productivity index of non-fractured or fractured vertical wells in various reservoirs [5–13]. Prats [5] first introduced the concept of “Optimum Fracture Design” (OFD) to obtain the maximum productivity of fractured vertical wells, which laid a solid foundation for optimizing fracture conductivity. Raghavan and Joshi [6] presented a rigorous model to evaluate the productivity of fractured wells by taking productivity index as a function of reservoir properties and wellbore radius. Valko and Economides [7] first applied the method of “Unified Fracture Design” (UFD) to optimize fracture conductivity by taking proppant number into account. Based on Valko and Economides’ method, a series of research on optimization of fracture conductivity were conducted. Romero et al. [8] further employed the direct boundary method to estimate the effect of fracture face skin and choke fracture skin on productivity index, which was an extension of the work presented by Valko and Economides. Luo et al. [9] presented a semi-analytical solution of vertically fractured wells under non-Darcy flow based on the UFD method and suggested that productivity index is strongly affected by non-Darcy flow when Reynolds number is less than 5. Hagoort [10] provided a semi-steady-state productivity index of a vertical well in a rectangular reservoir in the form of analytical solutions. However, the influence of permeability anisotropy on dimensionless productivity index was not taken into consideration in Hagoort’s work. Lately, Johansen et al. [12] filled the deficiency and established a new model to calculate the productivity index of a well in anisotropic reservoirs by employing coordinate transformation, which provided a feasible method to evaluate the productivity in anisotropic reservoirs. On the basis of the superposition principle and mirror-image method, Shi et al. [13] developed a semi-analytical productivity index model for a vertically fractured well with arbitrary fracture length under various boundary conditions and insisted that their model can be employed to calculate the productivity index of vertically fractured wells or fractured horizontal wells in low-permeability and unconventional reservoirs. These models mainly focus on the productivity index of non-fractured or fractured vertical wells in various formations, but the related approaches can provide important ideas for calculating the productivity index of fractured horizontal wells.

Recently, there have been many mathematical models to analyze the productivity of fractured horizontal wells. Medeiros et al. [14] established a semi-analytical model to discuss the productivity of fractured horizontal wells in heterogeneous and tight gas formations, and suggested that this system’s productivity increased significantly if natural fractures were created around the well during fracturing process. Bhattacharya et al. [15] extended the UFD method to fractured horizontal wells in low permeability reservoirs and obtained the general correlation between the maximum dimensionless productivity index and optimal fracture conductivity. Based on the instantaneous source solutions, Al Rbeawi and Tiab [16] introduced the pseudo-steady-state productivity index into fractured horizontal wells. Wang and Jia [17] filled the gap in a multistage fractured horizontal well for simultaneously optimizing multiple fractures with different properties. They pointed out that the maximum productivity index can be obtained by optimizing fracture properties for a given proppant volume. Kaul et al. [18] proposed a new approach, based on flow rate, to calculate the productivity index and corresponding derivative for a horizontal well with multiple vertical fractures in double-porosity reservoirs. Al Rbeawi [19] studied the behavior of the productivity index of fractured horizontal wells under Darcy and non-Darcy flow respectively, and demonstrated that the productivity index is affected by non-Darcy flow at early and intermediate flow regimes. Sorek et al. [20] extended the UFD method to maximize the productivity index of the horizontal well with multiple acute-angle transverse fractures and concluded that fracture width and fracture length depend on fracture angle and proppant number. Wang et al. [21] applied the UFD method to calculate the productivity index of a horizontal well with fracture networks in a closed rectangular reservoir and stated that the productivity index reaches the peak when fracture networks penetrate throughout the seepage area under infinite conductivity. Asadi et al. [22] combined the UFD method with Direct Boundary-Element Method (DBEM) to evaluate the influence of fracture-fluid-leak off on the productivity index of multiple-fracture horizontal wells and found that the optimal fracture with leak off should be shorter and wider for a fixed proppant number than the case without leak off. Guk et al. [23] employed the UFD method to develop a rigorous and unified optimization technique with type curves and calculated the optimal fracture number for a given proppant number by using type curves. Asadi and Zendehboudi [24] applied the Extended Distributed Volumetric Sources (EDVS) method to evaluate the fractured horizontal wells’ productivity index in unconventional reservoirs. In terms of economic aspects, they elaborated that the EDVS method offered a better way to design multi-fractured horizontal wells compared to the traditional techniques.

The previous models were established on the assumption of the ideal planar fractures [5–24]. In fact, many complex fractures can be formed in tight reservoirs, such as multi-wing fractures, non-planar fractures and reorientation fractures [25–28]. Generally, the productivity index for a horizontal well with multiple reorientation fractures has not been well investigated. In terms of multi-wing fractures, Luo et al. [29] developed a new fracture-wing model to investigate the effect of complex fractures on productivity index and provided an effective approach to calculate the productivity index for multi-wing fractures or non-planar fractures. However, the fracture-wing model has an important assumption that the wellbore must be located at the intersection of two fracture wings. For reorientation fractures, the wellbore may be located at any position along the fractures, and these fractures may also have multiple reorientations, thus the fracture-wing model could not deal with reorientation fractures very well. Therefore, other methods should be explored to analyze the characteristic of productivity index for reorientation fractures. Zhou et al. [30] came up with the nodal analysis technique to investigate the production of complex planar fracture network. This technique provides a way to handle the flow problem at the reoriented section of reorientation fractures. Based on this technique, combining the material balance condition at discrete nodes and the flow rate equation at each segment presented in the fracture-wing method, we can study the productivity index of the horizontal well with multiple reorientation fractures in this work.

The objective of this paper is to analyse the productivity for a horizontal well intercepted by multiple reorientation fractures in an anisotropic reservoir. This paper is organized as follows: Firstly, a model for a horizontal well with multiple reorientation fractures was established to calculate its dimensionless productivity index. Secondly, the accuracy of this model’s solution was validated by comparing with a classic case in the literature. Thirdly, the effect of some key reorientation fracture parameters on the horizontal well’s productivity index was discussed in detail.

2 Physical model

Figure 1 presents the conceptual model used in this work. In this model, a fractured horizontal well consists of three reorientation fractures and each fracture can be divided into two wings by the horizontal wellbore. In order to obtain the productivity index of this fractured horizontal well, some basic assumptions are made as follows:

1. This tight reservoir is closed, anisotropic, rectangular and homogeneous. Its permeability in the x- and y-direction is k_x and k_y, respectively.

2. The reservoir has a constant thickness, h, and constant porosity, φ. The fluid in the reservoir is slightly compressible and the viscosity of the fluid is μ.

3. The flow in the reservoir and reorientation fractures is single-phase flow and obeys Darcy’s law. All the reorientation fractures have a finite conductivity.

4. The horizontal well is parallel to x-axis and the reservoir is fully penetrated by three finite-conductivity reorientation fractures. As illustrated in Figure 1, each fracture reorients twice due to stress change in the reservoir and can be divided into three sections: one principal fracture and two reoriented fractures.

5. The width of the reorientation fractures is assumed to be constant, w_f, and their permeability is k_f. For the ith reorientation fracture, the angle between the principal fracture and x-axis, and the angles between the reoriented fractures and x-axis, are θ_i,1, θ_i,2 and θ_i,3, respectively. Similarly, the length of the principal fracture and two reoriented fractures is L_{f i,1}, L_{f i,2} and L_{f i,3}, respectively.

6. The pressure loss along the horizontal wellbore is neglected. No fluid flows from the reservoir into the fracture tips. The rate of the reorientation fractures contributing to the total rate of the horizontal wellbore is Q.

3 Mathematical model

In this section, a semi-analytical solution was derived to obtain the productivity index of a horizontal well with multiple reorientation fractures in an anisotropic reservoir. For the sake of simplicity, the dimensionless variables used in this work were defined in Table 1.

3.1 Flow model in the anisotropic reservoir

In the horizontal well model, a horizontal well is fully penetrated by N reorientation fractures and the total rate of the horizontal wellbore is constant and equal to the sum rates of N reorientation fractures. According to the point source function theory, a fracture can be regarded as a line source, and thus the line source problem can be further converted into a point source problem and solved by integrating along the fracture against the point source [31–34]. In this work, all the reorientation fractures are discretized into a series of minute segments, and all the discrete segments are considered to be uniform-flux [35, 36]. By using the superposition principle, the dimensionless pressure in the pseudo-steady period for the ith discrete segment in a closed anisotropic rectangular reservoir can be expressed as [29]:$$p_{\mathrm{Di}}=\frac{2\pi }{x_{\mathrm{eD}}{y}_{\mathrm{eD}}}{t}_{\mathrm{D}}+\sum _{k=1}^N\mathrm{}\sum _{j=1}^{N_k}\mathrm{}{q}_{\mathrm{Dk},j}{F}_{\mathrm{Di},j}\left(x_{\mathrm{Di}},y_{\mathrm{Di}},x_{\mathrm{Dj}},y_{\mathrm{Dj}},\beta \right)\hspace{1em}k=\mathrm{1,2}...N,\hspace{1em}i=\mathrm{1,2}...N_k,$$(1)where N_k is the number of the discrete segments for the kth reorientation fracture, q_Dk,j represents the flow flux of the jth discrete segment in the kth reorientation fracture, and F_Di,j is the influence function describing the extra pseudo-state pressure drop at the ith discrete segment caused by the jth discrete segment, which has been derived in detail in Appendix.

According to the material balance for a closed rectangular anisotropic reservoir, the average reservoir pressure can be described as follows [37]:$$p_{\mathrm{avg}}=p_i-\frac{\mathrm{Q\mu }B}{2\pi \bar{k}h}\frac{2\pi \eta t}{x_e{y}_e}.$$(2)

Based on the dimensionless definitions in Table 1, the dimensionless average pressure of this closed anisotropic rectangular reservoir can be expressed as:$$p_{\mathrm{avgD}}=\frac{2\pi }{x_{\mathrm{eD}}{y}_{\mathrm{eD}}}{t}_D.$$(3)

Combining equations (1) and (3), it yields:$$p_{\mathrm{Di}}=p_{\mathrm{avgD}}+\sum _{k=1}^N\mathrm{}\sum _{j=1}^{N_k}\mathrm{}{q}_{\mathrm{Dk},j}{F}_{\mathrm{Di},j}\left(x_{\mathrm{Di}},y_{\mathrm{Di}},x_{\mathrm{Dj}},y_{\mathrm{Dj}},\beta \right),\hspace{1em}k=\mathrm{1,2}...N,\hspace{1em}i=\mathrm{1,2}...N_k.$$(4)

Further, for all the discrete segments, equation (4) can be rewritten in the matrix form:$${\mathit{p}}_{\mathit{D}}-{\mathit{p}}_{\mathbf{avgD}}-{\mathit{A}}_{\mathbf{sum}}{\mathit{q}}_{\mathit{D}}=\mathbf{0},$$(5)where:$$A_{\mathrm{sum}}=\left[\begin{array}{ccc}{F}_{D\mathrm{1,1}}\left(x_{D1},y_{D1},x_{D1},y_{D1},\beta \right)& F_{D\mathrm{1,2}}\left(x_{D1},y_{D1},x_{D2},y_{D2},\beta \right)& \begin{array}{cc}\cdots & F_{D1,\mathrm{sum}}\left(x_{D1},y_{D1},x_{D\mathrm{sum}},y_{D\mathrm{sum}},\beta \right)\end{array}\\ F_{D\mathrm{2,1}}\left(x_{D2},y_{D2},x_{D1},y_{D1},\beta \right)& F_{D\mathrm{2,2}}\left(x_{D2},y_{D2},x_{D2},y_{D2},\beta \right)& \begin{array}{cc}\cdots & F_{D2,\mathrm{sum}}\left(x_{D2},y_{D2},x_{D\mathrm{sum}},y_{D\mathrm{sum}},\beta \right)\end{array}\\ \begin{array}{c}⋮\\ F_{D\mathrm{sum},1}\left(x_{D\mathrm{sum}},y_{D\mathrm{sum}},x_{D1},y_{D1},\beta \right)\end{array}& \begin{array}{c}⋮\\ F_{D\mathrm{sum},2}\left(x_{D\mathrm{sum}},y_{D\mathrm{sum}},x_{D2},y_{D2},\beta \right)\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ F_{D\mathrm{sum},\mathrm{sum}}\left(x_{D\mathrm{sum}},y_{D\mathrm{sum}},x_{D\mathrm{sum}},y_{D\mathrm{sum}},\beta \right)\end{array}\end{array}\end{array}\right],\hspace{1em} \mathrm{Sum}=\sum _{k=1}^N\mathrm{}{N}_k.$$

3.2 Flow model in the reorientation fracture

Different from planar fractures, reorientation fractures include two reoriented sections. Similar to the boundary element method to handle the irregular boundary, we use a series of discrete planar segments to approximately describe the reoriented sections, as illustrated in Figure 2. It is worth noting that the number of the discrete segments mainly depends on the calculation accuracy. To improve the accuracy, we usually do the local refinement for the reoriented sections [38]. Meanwhile, regarding the midpoint of each discrete segment as the discrete node, we can take the pressure at the node as the average pressure of the corresponding discrete segment.

The flow in the reorientation fracture is one-dimensional Darcy’s flow. Based on the result of Zhou et al. [30], the pressure difference between two adjacent discrete nodes, the ith node and the (i-1)th node (shown in Fig. 3), is given by:$$p_{f,i}-p_{f,i-1}=\underset{l_i}{\overset{l_{i-1}}{\int }}\mathrm{}\frac{\mu B}{k_f{w}_{f}h}{q}_f\left(l,t\right)\mathrm{d}l.$$(6)

Note that the flow rate of each discrete segment may be different, equation (6) can be further rewritten as:$$p_{f,i}-p_{f,i-1}=\underset{l_i}{\overset{l_{i-0.5}}{\int }}\mathrm{}\frac{\mu B}{k_f{w}_{f}h}{q}_{\mathrm{fw}}\left(l,t\right)\mathrm{d}l+\underset{l_{i-0.5}}{\overset{l_{i-1}}{\int }}\mathrm{}\frac{\mu B}{k_f{w}_{f}h}{q}_{\mathrm{fw}}\left(l,t\right)\mathrm{d}l.$$(7)

According to the mass balance, the flow in each discrete segment should satisfy the following equation:$$q_{\mathrm{fw}}\left(l,t\right)=q_{\mathrm{fw},i}+q_{\mathrm{fi}}\left(l_i-l\right),\hspace{1em}{l}_{i-1}\le l\le l_i.$$(8)

Substituting equation (8) into equation (7), it yields:$$p_{f,i}-p_{f,i-1}=\underset{l_i}{\overset{l_{i-0.5}}{\int }}\mathrm{}\frac{\mu B}{k_f{w}_{f}h}\left[q_{\mathrm{fw},i}+q_{\mathrm{fi}}\left(l_i-l\right)\right]\mathrm{d}l+\underset{l_{i-0.5}}{\overset{l_{i-1}}{\int }}\mathrm{}\frac{\mu B}{k_f{w}_{f}h}\left[q_{\mathrm{fw},i-1}+q_{\mathrm{fi}-1}\left(l_{i-1}-l\right)\right]\mathrm{d}l.$$(9)

Meanwhile, the flow in every discrete node should satisfy the mass balance [39]:$$q_{\mathrm{fw},i}=\sum _{j=1}^{N_1}\mathrm{}\left(q_{f,j}\Delta l_j\right)-\sum _{j=1}^{i-1}\mathrm{}\left(q_{f,j}\Delta l_j\right).$$(10)

Thus, the pressure difference between the horizontal wellbore node and the ith discrete node can be obtained:$$\begin{array}{c}{p}_{f,i}-p_w=\underset{l_{1-0.5}}{\overset{l_{1+0.5}}{\int }}\mathrm{}{\left(\frac{\mu B}{k_f{w}_{f}h}\right)}_1\left[\sum _{j=1}^{N_1}\mathrm{}\left(q_{\mathrm{fj}}\Delta l_j\right)+q_{f1}\left(l_{1-0.5}-l\right)\right]\mathrm{d}l+\cdots \\ \hspace{1em}\hspace{1em}+\underset{l_{\left(i-1\right)-0.5}}{\overset{l_{\left(i-1\right)+0.5}}{\int }}\mathrm{}{\left(\frac{\mu B}{k_f{w}_{f}h}\right)}_{i-1}\left[\sum _{j=1}^{N_1}\mathrm{}\left(q_{\mathrm{fj}}\Delta l_j\right)-\sum _{j=1}^{i-2}\mathrm{}\left(q_{\mathrm{fj}}\Delta l_j\right)+q_{\mathrm{fi}-1}\left(l_{i-1-0.5}-l\right)\right]\mathrm{d}l\\ \hspace{1em}\hspace{1em}+\underset{l_{\left(i-1\right)+0.5}}{\overset{l_i}{\int }}\mathrm{}{\left(\frac{\mu B}{k_f{w}_{f}h}\right)}_i\left[\sum _{j=1}^{N_1}\mathrm{}\left(q_{\mathrm{fj}}\Delta l_j\right)-\sum _{j=1}^{i-1}\mathrm{}\left(q_{\mathrm{fj}}\Delta l_j\right)+q_{\mathrm{fi}}\left(l_{i-0.5}-l\right)\right]\mathrm{d}l.\end{array}$$(11)

Substituting the dimensionless definitions in Table 1 into equation (11), it can be rewritten as:$$p_{\mathrm{wD}}-p_{\mathrm{fD},i}=\frac{2\pi }{C_{\mathrm{fD}}}\left\{l_{\mathrm{Di}}\sum _{j=1}^{N_1}\mathrm{}{q}_{\mathrm{fDj}}-\frac{\Delta l_{\mathrm{Di}}}{8}{q}_{\mathrm{fDi}}-\sum _{j=1}^{i-1}\mathrm{}\left(\frac{\Delta l_{\mathrm{Dj}}}{2}+l_{\mathrm{Di}}-\sum _{n=1}^j\mathrm{}\Delta l_{\mathrm{Dn}}\right)q_{\mathrm{fDj}}\right\}.$$(12)

Equation (12) is a completely new equation to describe the flow inside a minute segment with arbitrary length, and can be employed to analyze the flow inside the reoriented sections which were approximately dealt with a series of minute planar segments. Thus, for one fracture wing, equation (12) can be rewritten in the matrix form:$${\mathit{p}}_{\mathbf{wD}}-{\mathit{p}}_{\mathbf{fD}}-{\mathit{g}}_{\mathit{i},\mathbf{1}}{\mathit{q}}_{\mathbf{fD}}=\mathbf{0}.$$(13)

The above matrix describes the flow in one wing of the ith reorientation fracture, and g_i,1 is a coefficient matrix which can be calculated from equation (12). Similarly, we can obtain the matrix describing the flow inside the other wing of the ith reorientation fracture and obtain another coefficient matrix g_i,2.

In general, the fractured horizontal well always consists of multiple reorientation fractures. Obviously, the flow in each reorientation fracture is independent and all the flow in these fractures satisfies equation (13). Therefore, for all the fracture wings, we can get the flow control matrix of the fractured horizontal well in this study:$${\mathit{p}}_{\mathrm{wD}}-{\mathit{p}}_{\mathrm{fD}}-{\mathit{G}}_{\mathrm{Sum}}{\mathit{q}}_{\mathrm{fD}}=\mathbf{0}.$$(14)

As the horizontal well is intercepted by N reorientation fractures, it naturally has 2N fracture wings. Each wing will correspond to a coefficient matrix, g_i,1 or g_i,2, which can also be derived from equation (12). Therefore, the coefficient matrix G_sum can be expressed as follows:

In the above matrix G_sum, two matrixes in one red box correspond to the coefficient matrixes of two wings of a reorientation fracture.

3.3 Productivity index for a horizontal well with multiple reorientation fractures

Due to the continuity of pressure and flux at all the discrete segments’ face, the following conditions must be satisfied along the fracture face:$$p_{\mathrm{Di}}=p_{\mathrm{fDi}}, q_{\mathrm{Di}}=q_{\mathrm{fDi}}.$$(16)

Considering the continuity conditions, coupling equations (5) and (14), we can obtain the following matrix to describe the flow at the face of the reorientation fractures:$${\mathit{p}}_{\mathbf{wD}}-{\mathit{p}}_{\mathbf{avgD}}-\left({\mathit{A}}_{\mathbf{sum}}+{\mathit{G}}_{\mathbf{sum}}\right){\mathit{q}}_{\mathbf{fD}}=\mathbf{0}.$$(17)

For a horizontal well with multiple fractures, the total dimensionless productivity index is expressed as:$$J_D=\frac{\sum _{k=1}^N\mathrm{}\sum _{j=1}^{N_k}\mathrm{}{q}_{\mathrm{fDj}}}{p_{\mathrm{wD}}-p_{\mathrm{avgD}}}=\sum _{k=1}^N\mathrm{}\sum _{j=1}^{N_k}\mathrm{}{J}_{\mathrm{Dj}}.$$(18)

Combining equations (17) and (18), we yield the following matrix equation:$$\left({\mathit{A}}_{\mathrm{sum}}+{\mathit{G}}_{\mathrm{sum}}\right){\mathit{J}}_D=\mathit{I}.$$(19)

The coefficient matrix (A_sum + G_sum) is not a sparse matrix, and thus equation (19) can be solved by using the Gauss elimination method. Therefore, all unknown dimensionless productivity indexes for the discrete segments can be obtained by solving equation (19). Further, the total dimensionless productivity index for the horizontal well with multiple reorientation fractures can be calculated through equation (18). Finally, we can use the productivity index to analyze the productivity of a horizontal well with multiple reorientation fractures.

4 Model verification

To the best of our knowledge, there is no report on discussion about the productivity index for a horizontal well intercepted by multiple reorientation fractures in an anisotropic reservoir. Luo et al. [29] investigated the productivity index of a horizontal well with multiple fractures with azimuth angle, which can be regarded as a special case of our model. Obviously, if principal fracture angle and reoriented fracture angle are equal in our proposed model, reorientation fractures can be reduced to planar fractures, which are exactly the same as that in Luo’s model. Analogously, when anisotropic factor β is equal to 1 in our model, tight reservoir becomes isotropic and is like Luo et al.’s reservoir model [29]. Figures 19 and 20 in reference [29] present the effect of fracture spacing and fracture rotation on horizontal well productivity index, respectively. Referring to these two figures [29], we take the same related parameters in our model and compare our calculation results with Luo et al.’s results to verify the accuracy of our model. Figures 4a and 4b show the comparison results under different fracture spacing and fracture rotation, respectively. It can be observed that the results of our work are in excellent agreement with the results of Luo et al. [29]. The verification indicates that the semi-analytical solution derived in this work is reliable.

5 Results and discussion

Kaul et al. [18] evaluated well productivity by calculating the productivity index and corresponding derivative. Similarly, in this study we also use these two parameters to analyze the productivity of the horizontal well with multiple reorientation fractures. In this section, we mainly focus on the effect of permeability anisotropy, fracture reorientation, fracture number, fracture spacing and fracture distribution on the productivity index and derivative. The investigation of fracture reorientation on well productivity is crucial and meaningful, and will be discussed from three aspects: principal fracture angle, reoriented fracture angle and reoriented factor. Their definitions have been given in the previous sections. Fracture distribution also has significant effect on well productivity, and we will select several kinds of typical fracture distributions to study their influence. Considering the importance of fracture conductivity on well productivity, it will be discussed throughout this section. Additionally, the reorientation fracture model established in this paper is aimed at horizontal wells, but it is easy to apply this model to vertical or inclined wells.

5.1 Effect of permeability anisotropy on horizontal well productivity

For an anisotropic reservoir, the permeability in the x- and y-axis is different, and anisotropic factor, β, is always used to assess the permeability anisotropy [40, 41]:$$\beta =\sqrt{\frac{k_x}{k_y}}.$$(20)

Figure 5 demonstrates the effect of permeability anisotropy on the dimensionless productivity index and its derivative under various fracture conductivities. For a fixed fracture configuration, with the increase of fracture conductivity, the productivity index first increases slowly and then almost remains constant. This characteristic can be explained as: when fracture conductivity is low, the ability of a fracture to provide fluid is limited and the productivity mainly depends on fracture conductivity; when fracture conductivity is high enough, the ability of a fracture to provide fluid can meet the requirements of fluid supply to the wellbore, and the productivity depends on other factors, such as principal fracture angle and fracture distributions. However, as fracture conductivity increases, the dimensionless productivity index derivative first increases and then decreases. The derivative curve shows a hump and has a maximum value (dash dot lines in Fig. 5), which just corresponds to the optimal fracture conductivity under this fracture configuration [42]. This feature also appears in the following discussions, implying that the pursuit of high fracture conductivity in fractured horizontal wells is meaningless and uneconomical during the oilfield stimulation.

When anisotropic factor increases, meaning that the flow in the x-axis becomes dominant, the productivity index gradually decreases and the corresponding curves move downward as shown in Figure 5. It indicates that for a given fracture configuration, permeability anisotropy is inversely correlated with the well productivity index. For example, as anisotropic factor increases from 2 to 4, the dimensionless productivity index decreases from 0.807 to 0.507 when C_fD is equal to 100. This suggests that strong permeability anisotropy is a negative factor to well production, which can be explained by the fact that when permeability anisotropy gets strong, the flow in the x-axis gradually dominates and the flow in the y-axis is inhibited, resulting in a decrease of fluid supply from the reservoir to the wellbore. Meanwhile, the peak on the derivative curve moves down to the left, meaning that the optimal fracture conductivity decreases slowly but limitedly as permeability anisotropy gets strong.

5.2 Effect of fracture reorientation on horizontal well productivity

5.2.1 Principal fracture angle

Figure 6 displays the effect of principal fracture angle (θ₁) on the well productivity under different fracture conductivities. For an intuitive comparison, an assumption was made that all the reoriented fracture angles are set to 30°. Similar to the last section, as fracture conductivity increases, all the productivity index curves first rise and then tend to a cluster of horizontal lines, while all the derivative curves show obvious inflection points (dash dot lines in Fig. 6).

In Figure 6, it is obviously shown that with the increase of the principal fracture angle, the dimensionless productivity index increases and the tendency of the increase for the productivity index varies slightly with the increase of fracture conductivity. For example, if we set C_fD = 100 and increase the principal fracture angle from 30° to 60°, the productivity index correspondingly increases from 0.744 to 0.807, nearly increasing by 8.5%. Large principal fracture angle corresponds to great reservoir contact area in the principal permeability direction, which is more conducive to fluid flow, and it suggests that we can increase the well’s productivity by changing the principal fracture angle (θ₁ ≤ 90°). However, in the case of C_fD = 0.1, the corresponding increasing rate is only 2.98% when the principal fracture angle increases from 30° to 60°, indicating that pseudo-steady-state productivity can be enhanced when fracture conductivity increases appropriately. Meanwhile, the optimal fracture conductivity slightly increases with the increase of the principal fracture angle. In total, the principal fracture angle has a limited effect on the optimal fracture conductivity.

5.2.2 Reoriented fracture angle

Similarly, in order to investigate the effect of reoriented fracture angle on the horizontal well productivity, all the principal fracture angles are fixed at 60°. Figure 7 illustrates the effect of reoriented fracture angle on the dimensionless productivity index and its derivative under different fracture conductivities.

Compared with Figure 6, the effect of reoriented fracture angle on the productivity index is not as strong as that of principal fracture angle, and all the productivity index curves for different reoriented fracture angles almost coincide and so do the derivative curves. For instance, when the reoriented fracture angle increases from 30° to 60°, the corresponding productivity index increases from 0.807 to 0.828 with an increasing rate of 2.6% (C_fD = 100). However, when the principal fracture angle increases from 30° to 60°, the productivity index increases by 8.48%, from 0.744 to 0.807 (Fig. 6). This phenomenon may be interpreted by that the length of the reoriented fractures is shorter than that of the principal fracture, and then the seepage area affected by the reoriented fractures is smaller than that affected by the principal fracture. It reminds us that changing the reoriented fracture angle to enhance the horizontal well’s productivity is not as effective as changing the principal fracture angle. Similar to the cases in Figure 6, the optimal fracture conductivities mildly increase with the increase of the reoriented fracture angle, and the reoriented fracture angle only has a slight effect on the optimal fracture conductivity.

5.2.3 Reoriented factor

For a horizontal well with multiple reorientation fractures, its seepage area depends on the length of the reoriented fractures to some extent and thus the length of the reoriented fractures is another significant factor affecting its productivity. In this study, we define a new parameter, reoriented factor, σ, to assess the extent of the reoriented fracture deviating from the principal fracture. According to its definition, reoriented factor is the ratio of the sum of the reoriented fractures’ length to the principal fracture’s length under a fixed total length of reorientation fracture. Letting σ = 1/3, 1/2, 1, 2 and 3, Figure 8 demonstrates the effect of reoriented factor on the well’s productivity under various fracture conductivities.

As shown in Figure 8, the effect of reoriented factor on the productivity index can be neglected when fracture conductivity is lower than 0.2. Under this circumstance, the well’s productivity mainly depends on fracture conductivity, and the influence of reoriented factor is relatively weak in contrast. As fracture conductivity increases (C_fD > 0.2), the influence of reoriented factor on the productivity index becomes stronger and a drop of the productivity index occurs as reoriented factor goes up. For instance, when C_fD is equal to 100, as the reoriented factor increases from 1/3 to 1, the productivity index reduces from 0.807 to 0.785, decreasing by 2.78%. It can be observed in Figure 8 that a large length of the reoriented fractures is not conducive to the horizontal well’s productivity when the reoriented fracture angle is smaller than the principal fracture angle, which suggests that the reoriented fracture angle should be greater than the principal fracture angle in the process of fracturing to achieve a higher well productivity. In addition, we can observe that the maximum productivity index derivative slightly decreases with the increase of the reoriented factor, but the reoriented factor has a weak influence on the optimal fracture conductivity.

5.3 Effect of reorientation fracture number on horizontal well productivity

To completely reveal the effect of reorientation Fracture Number (FN) on the horizontal well’s productivity, we set all the dimensionless fracture spacing equal to 2 and let FN = 3, 5, 7 and 9, respectively. Figure 9 presents the effect of reorientation fracture number on the well’s productivity. When fracture number is fixed, the influence of fracture conductivity on the shape and trend of the productivity index and derivative curves is almost the same as that discussed in previous cases (compared with Figs. 5–8).

In terms of reorientation fracture number, due to the enhanced seepage area near the horizontal wellbore, the productivity index increases significantly with the increase of fracture number and the productivity index curves are nearly parallel to each other in Figure 9. However, the increase gradually slows down when the fracture number increases to a certain degree. For instance, when the fracture number increases from 3 to 5 (C_fD = 100), the corresponding productivity index increases from 0.648 to 0.851, increasing by 31.3%, but the productivity index only increases by 9.61% when the fracture number increases from 7 to 9. It indicates that the fractured horizontal well has an optimal reorientation fracture number to obtain the economical productivity for a fixed seepage area. However, from the productivity index derivative curves, the optimal fracture conductivity has a little decrease with the increase of reorientation fracture number, demonstrating that pursing higher fracture conductivity is not necessary for fractured horizontal wells. For example, when reorientation fracture number increases from 3 to 5, the optimal fracture conductivity just decreases from 0.628 to 0.571, and even the optimal fracture conductivity tends to be the same for the cases of FN = 7 and FN = 9. This phenomenon can be explained by that for a fixed seepage area, the increase of fracture number creates additional flow interference among reorientation fractures, and when reorientation fracture number is big enough, the optimal fracture conductivity does not decrease with the increase of fracture number.

5.4 Effect of fracture spacing on horizontal well productivity

Considering reorientation fracture number is limited, we try to investigate the productivity behavior caused by the change of Fracture Spacing (FS) and further to pursue the optimal fracture spacing during hydraulic fracturing. Here we set FN equal to 3 and let FS = 2, 4, 6 and 8 respectively. Figure 10 demonstrates that the effect of fracture spacing on the well’s productivity under different fracture conductivities.

When fracture spacing is narrow (FS = 2), each reorientation fracture governs a small seepage area and these reorientation fractures are easy to interfere with each other, causing an extra decrease of productivity index (J_Dmax = 0.648). With the increase of fracture spacing, the dimensionless productivity index rises slowly. For FS = 4, the maximum productivity index increases to 0.807, with an increasing rate of 24.6%. However, the increase of the dimensionless productivity index slows down when fracture spacing increases to a certain degree. For example, the maximum increasing rate is only 10.5% when facture spacing increases from 4 to 6. Further, when fracture spacing continues to increase (FS = 8), the productivity index decreases instead. It may be explained by that the reorientation fractures near the ends of the horizontal wellbore are closer to the closed boundary when fracture spacing is wide enough, and thus the effect of the closed boundary on the fluid flow appears earlier and more severely. Therefore, the wider or narrower fracture spacing is not conducive to increase the productivity index, which implies that the optimal fracture spacing should be demonstrated to obtain the maximum productivity index. Additionally, the effect of fracture spacing on the optimal fracture conductivity also can be neglected.

5.5 Effect of fracture distribution on horizontal well productivity

There is an implicit condition in the previous discussion that all reorientation fractures are parallel to each other and also distribute symmetrically with respect to the corresponding intersection of the fracture and the horizontal wellbore. In this section, we focus on discussing the influence of several typical fracture distributions that do not satisfy this implicit condition on the horizontal well’s productivity, such as one or more reorientation fractures’ arbitrary rotation or asymmetrical distribution.

5.5.1 Arbitrary rotation of one reorientation fracture

Here we set fracture number and fracture spacing equal to 3 and 4, respectively. To illustrate the effect of arbitrary rotation of one reorientation fracture on productivity index, we discuss two kinds of rotations, the rotation of the central fracture and the rotation of the external fracture. Under this circumstance, we assume that both its principal fracture and reoriented fractures have the same rotation angle when the reorientation fracture occurs to rotate. Figures 11a and 11b present the effect of these two typical rotations on the well’s productivity index and its derivative, respectively.

As observed from Figure 11a, the arbitrary rotation of the central fracture has a weak impact on the productivity index. When the principal fracture angle of the central fracture increases from 30° to 75°, the corresponding productivity index goes up from 0.799 to 0.827 (C_fD = 100), with an increasing rate of 3.5%. Although the maximum productivity index derivative shows a little increase with the increase of the rotation angle of the central fracture, the optimal fracture conductivity remains the same when the rotation only occurs in the central fracture. Different from Figure 11a, Figure 11b shows that there is an obvious change of the productivity index when the rotation occurs in the external reorientation fracture. For example, the increasing rate of the productivity index is 5.3% when the principal fracture angle of the external reorientation fracture increases from 30° to 75° (C_fD = 100), larger than that caused by the same rotation angle of the central reorientation fracture. This phenomenon may be explained by that the external fracture has a larger seepage area and a stronger influence than the central fracture. Hence, in order to improve the stimulation performance, rotating the external reorientation fracture is more effective than rotating the central reorientation fracture.

5.5.2 Asymmetrical distribution of one reorientation fracture

To investigate the effect of the asymmetrical distribution of the reorientation fracture, we first discuss the simplest case, namely one reorientation fracture’s asymmetrical distribution. Here we assume that one of the fractures moves a certain distance toward a wing, making it asymmetrical about the intersection of the fracture and the horizontal wellbore. Similar to the last section, the reorientation fracture number and adjacent fracture spacing are set to 3 and 4, respectively. Figure 12 illustrates the effect of the asymmetrical distribution of one reorientation fracture on the well productivity, including three cases: (A) reference case: symmetrical distribution, (B) asymmetrical distribution of one external reorientation fracture, (C) asymmetrical distribution of the central reorientation fracture.

As shown in Figure 12, the asymmetrical distribution slightly affects the productivity index and its derivative. Obviously, all the productivity index curves almost overlap. Due to the fact that the area between the reorientation fractures is larger and the flow interference is weaker for Case B, it has a slightly larger productivity index when fracture conductivity is high enough. Considering the symmetry problem, no matter which wing the fracture moves toward, its influence on the productivity index and derivative is the same. Additionally, it can be observed from Figure 12 that the optimal conductivity almost keeps the same for different asymmetrical distributions of one reorientation fracture.

5.5.3 Asymmetrical distribution of two or more reorientation fractures

Based on the investigation of the effect of the asymmetrical distribution of one reorientation fracture on the productivity index and its derivative, we can further discuss the effect of two or more reorientation fractures’ asymmetrical distribution on fractured horizontal well productivity. The setting of fracture number and fracture spacing is the same as the last section.

Figure 13 shows the effect of asymmetrical distribution of two reorientation fractures on the productivity index, including four cases: (A) reference case: symmetrical distribution, (B) asymmetrical distribution of one external fracture and the central fracture, (C) asymmetrical distribution of two external fractures, (D) inverse asymmetrical distribution of one external fracture and the central fracture. Obviously, the seepage area affected by fractured stimulation in Case D is the largest. As shown in Figure 13, there is a slight difference in the productivity index when fracture conductivity is higher than 0.2. When 0.2 ≤ C_fD < 10, the productivity of Case A is slightly larger than the other three cases; when C_fD > 10, Case D has the largest productivity and it can be illustrated by the maximum seepage area affected by fractured stimulation in this case. Meanwhile, according to the derivative curves, the optimal fracture conductivity of Case A is the lowest and that of the other three cases is nearly the same, which can be explained by that the flow interference for the asymmetrical distribution of two reorientation fractures is more serious than that for the symmetrical distribution case.

Figure 14 demonstrates the effect of asymmetrical distribution of three reorientation fractures on well productivity, including four cases: (A) reference case: symmetrical distribution, (B) eccentric asymmetry distribution, (C) inverse asymmetry distribution I, D) inverse asymmetry distribution II. For a given fracture conductivity, the productivity index of Case C and Case D are almost the same and larger than that of another two cases, and the productivity index of Case B is the smallest. For example, when C_fD = 100, the productivity index drops from 0.902 to 0.859, a reduction of 4.8%, when the fracture configuration changes from Case C/D to Case B. Obviously, when three reorientation fractures distribute inversely and asymmetrically (Case C and Case D), these fractures are staggered with each other and the flow interference between them is the weakest. Thus, each fracture can take advantage of its flow ability under this condition and finally Case C and Case D obtain the maximum productivity index. Meanwhile, Case C and Case D have the same optimal fracture conductivity, which is equal to 1.9, but it’s rather different from another two cases’ optimal fracture conductivity, which is 0.72 for Case A and 1.3 for Case B, respectively. Therefore, different fracture distributions require different fracture conductivity to exploit their respective advantages.

6 Conclusion

Based on the above results and discussions in this work, several conclusions can be drawn as follows:

1. A semi-analytical model was established to calculate the dimensionless productivity index and derivative for a horizontal well with multiple reorientation fractures in an anisotropic rectangular reservoir. A new semi-analytical solution was obtained and its reliability was verified by a special case in the literature.

2. For a horizontal well with multiple reorientation fractures, strong permeability anisotropy is a negative factor for well production and the productivity index gradually decreases with the increase of anisotropic factor. The increase of principal fracture angle and reoriented fracture angle is beneficial to the productivity index. For a fixed fracture length, the length of reoriented fracture should be as short as possible to achieve the economic productivity when reoriented fracture angle is smaller than principal fracture angle.

3. Fracture number and fracture spacing have a strong influence on the well’s productivity, and the productivity index increases significantly with the increase of reorientation fracture number. A fractured horizontal well has an optimal reorientation fracture number for a fixed seepage area, and it’s necessary to optimize fracture number to obtain the maximum productivity. In addition, the external fractures should be kept away from the closed boundary as possible to obtain the economical productivity.

4. In an anisotropic reservoir, the effect of the rotation and distribution of the reorientation fractures on the well’s productivity should be taken into account. Compared with the rotation of the central fracture, the rotation of the external fracture has a greater influence on the productivity index. Meanwhile, the asymmetrical distribution of one or more reorientation fractures slightly affects the productivity index when fracture conductivity is large enough.

5. For a fixed fracture configuration, the optimal fracture conductivity can be observed from the derivative curves, and the pursuit of great fracture conductivity in hydraulic fracturing is meaningless and uneconomical for horizontal wells. In general, permeability anisotropy, principal fracture angle, reoriented fracture angle, reoriented factor, fracture number, fracture spacing and fracture distribution have a limited effect on the optimal fracture conductivity.

Acknowledgments

The authors gratefully acknowledge the financial support from the National Science and Technology Major Project of China (No. 2017ZX05030-002 & No. 2017ZX05013-004).

Author Contributions

All authors have contributed to this work. M.W established the mathematical model and wrote the main manuscript. As the corresponding author, G.X. made substantial contributions to the conception/design of the work. Z.F., L.Z., W.Z. and C.T. contributed to the modeling, programming, and results analysis and discussion.

Appendix

Derivation of Influence Function F_D

Based on the instantaneous point source solution, Wu et al. [41] obtained the pressure solution for a vertical well in a closed anisotropic rectangular reservoir:$$\begin{array}{c}s{\tilde{p}}_D=\frac{\pi \beta }{x_{\mathrm{eD}}}\left\{\frac{\cosh {\psi }_0\left(y_{\mathrm{eD}}-\left|y_D\pm y_{\mathrm{Di}}\right|\right)}{{\psi }_0\sinh {\psi }_0{y}_{\mathrm{eD}}}+2\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh {\psi }_n\left(y_{\mathrm{eD}}-\left|y_D\pm y_{\mathrm{Di}}\right|\right)}{{\psi }_n\sinh {\psi }_n{y}_{\mathrm{eD}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left(\frac{n\pi x_{\mathrm{Di}}}{x_{\mathrm{eD}}}\right)\right\}\\ {\psi }_n=\sqrt{\left(s\beta +{\left(\beta \frac{n\pi }{x_{\mathrm{eD}}}\right)}^2\right)},\hspace{1em}n=\mathrm{0,1},2\cdots \end{array}$$(A.1)

Note that if β is equal to 1, equation (A.1) is regressed to the pressure solution for a vertical well in an isotropic reservoir, which was presented by Ozkan [32].

Conducting the asymptotic analysis and Laplace inversion for equation (A.1), the pressure for a vertical well in a closed anisotropic rectangular reservoir in pseudo-state period can be obtained:$$p_D\approx \frac{2\pi }{x_{\mathrm{eD}}{y}_{\mathrm{eD}}}{t}_D+2\pi \beta \frac{y_{\mathrm{eD}}}{x_{\mathrm{eD}}}\left[\frac{1}{3}-\frac{\left|y_D\pm y_{\mathrm{Di}}\right|}{2y_{\mathrm{eD}}}+\frac{y_D^2+y_{\mathrm{Di}}^2}{2y_{\mathrm{eD}}^2}\right]+2\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left(y_{\mathrm{eD}}-\left|y_D\pm y_{\mathrm{Di}}\right|\right)}{n\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left(\frac{n\pi x_{\mathrm{Di}}}{x_{\mathrm{eD}}}\right).$$(A.2)

Based on the dimensionless definitions, the average pressure of this closed anisotropic rectangular reservoir can be expressed as:$$p_{\mathrm{avgD}}=\frac{2\pi }{x_{\mathrm{eD}}{y}_{\mathrm{eD}}}{t}_D.$$(A.3)

Combining equation (A.2) with (A.3), it yields:$$p_D-p_{\mathrm{avgD}}=2\pi \beta \frac{y_{\mathrm{eD}}}{x_{\mathrm{eD}}}\left[\frac{1}{3}-\frac{\left|y_D\pm y_{\mathrm{Di}}\right|}{2y_{\mathrm{eD}}}+\frac{y_D^2+y_{\mathrm{Di}}^2}{2y_{\mathrm{eD}}^2}\right]+2\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left(y_{\mathrm{eD}}-\left|y_D\pm y_{\mathrm{Di}}\right|\right)}{n\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left(\frac{n\pi x_{\mathrm{Di}}}{x_{\mathrm{eD}}}\right).$$(A.4)

By defining the following function, F_D, which represents the extra pressure drop caused by the closed boundary in the pseudo-state period,$$F_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right)=2\pi \beta \frac{y_{\mathrm{eD}}}{x_{\mathrm{eD}}}\left[\frac{1}{3}-\frac{\left|y_D\pm y_{\mathrm{Di}}\right|}{2y_{\mathrm{eD}}}+\frac{y_D^2+y_{\mathrm{Di}}^2}{2y_{\mathrm{eD}}^2}\right]+2\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left(y_{\mathrm{eD}}-\left|y_D\pm y_{\mathrm{Di}}\right|\right)}{n\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left(\frac{n\pi x_{\mathrm{Di}}}{x_{\mathrm{eD}}}\right).$$(A.5)

Equation (A.4) can be further rewritten as:$$p_D=p_{\mathrm{avgD}}+F_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right).$$(A.6)

In this work, we assume that the flow in every discrete segment is uniform-flux. For a uniform-flux discrete segment with an angle of θ_i, the following coordinate transformation relation can be employed:$${x\mathrm{\prime }}_{\mathrm{Di}}^{}=x_{\mathrm{Di}}+{\chi }_D\cos {\theta }_i,\hspace{0.5em}{y\mathrm{\prime }}_{\mathrm{Di}}^{}=y_{\mathrm{Di}}+{\chi }_D\sin {\theta }_i,\hspace{0.5em}-\Delta l_{\mathrm{Di}}/2{\chi }_D \Delta l_{\mathrm{Di}}/2.$$(A.7)

Integrating χ_D from −Δl_Di/2 to Δl_Di/2 in equation (A.5) gets the extra pseudo-state pressure drop of the point (x_Di, y_Di) caused by a uniform-flux discrete segment:$$\begin{array}{c}{F}_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right)=\frac{1}{\Delta l_{\mathrm{Di}}}\underset{-\frac{\Delta l_{\mathrm{Di}}}{2}}{\overset{\frac{\Delta l_{\mathrm{Di}}}{2}}{\int }}\mathrm{}2\pi \beta \frac{y_{\mathrm{eD}}}{x_{\mathrm{eD}}}\left[\frac{1}{3}-\frac{\left|y_D\pm \left(y_{\mathrm{Di}}+{\chi }_D\sin {\theta }_i\right)\right|}{2y_{\mathrm{eD}}}+\frac{y_D^2+{\left(y_{\mathrm{Di}}+{\chi }_D\sin {\theta }_i\right)}^2}{2y_{\mathrm{eD}}^2}\right]\mathrm{d}{\chi }_D\hspace{1em}\hspace{1em}\\ +\frac{2}{\Delta l_{\mathrm{Di}}}\underset{-\frac{\Delta l_{\mathrm{Di}}}{2}}{\overset{\frac{\Delta l_{\mathrm{Di}}}{2}}{\int }}\mathrm{}\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left|y_D\pm \left(y_{\mathrm{Di}}+{\chi }_D\sin {\theta }_i\right)\right|\right]}{n\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left[\frac{n\pi \left(x_{\mathrm{Di}}+{\chi }_D\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\right]\mathrm{d}{\chi }_D.\end{array}$$(A.8)

Equation (A.8) can be simplified as follows:$$F_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right)=\frac{2\pi \beta }{3}\frac{y_{\mathrm{eD}}}{x_{\mathrm{eD}}}-\frac{\beta \pi }{x_{\mathrm{eD}}\Delta l_{\mathrm{Di}}\sin {\theta }_i}\underset{y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\overset{y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\int }}\mathrm{}\left[\left|y_D+t\right|\right]\mathrm{d}t+\frac{\pi \beta }{x_{\mathrm{eD}}{y}_{\mathrm{eD}}\Delta l_{\mathrm{Di}}\sin {\theta }_i}\underset{y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\overset{y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\int }}\mathrm{}\left[y_D^2+t^2\right]\mathrm{d}t+\frac{2}{\Delta l_{\mathrm{Di}}\sin {\theta }_i}\underset{y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\overset{y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\int }}\mathrm{}\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left|y_D+t\right|\right]}{n\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left[\frac{n\pi \left(x_{\mathrm{Di}}+\left(t-y_{\mathrm{Di}}\right)\cot {\theta }_i\right)}{x_{\mathrm{eD}}}\right]\mathrm{d}t-\frac{\beta \pi }{x_{\mathrm{eD}}\Delta l_{\mathrm{Di}}\sin {\theta }_i}\underset{y_{\mathrm{Di}}-{\mathrm{y}}_D-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\overset{y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\int }}\mathrm{}\left|t\right|\mathrm{d}t+\frac{2}{\Delta l_{\mathrm{Di}}\sin {\theta }_i}\underset{y_{\mathrm{Di}}-y_D-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\overset{y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i}{\int }}\mathrm{}\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left|t\right|\right]}{n\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)\cos \left[\frac{n\pi \left(x_{\mathrm{Di}}+\left(t+y_D-y_{\mathrm{Di}}\right)\cot {\theta }_i\right)}{x_{\mathrm{eD}}}\right]\mathrm{d}t.$$(A.9)

Equation (A.9) is the extra pseudo-state pressure drop of any point in the closed anisotropic rectangular reservoir caused by a uniform-flux discrete segment, and we call it influence function.

In order to get the rapid speed of computation, equation (A.9) needs further treatment.

1. When $y_{\mathrm{Di}}-y_D-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\ge 0$, equation (A.9) can be rewritten as:

$$F_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right)=\frac{\beta \pi }{x_{\mathrm{eD}}}\left[\frac{2}{3}{y}_{\mathrm{eD}}-2y_{\mathrm{Di}}+\frac{\left(y_D^2+y_{\mathrm{Di}}^2\right)}{y_{\mathrm{eD}}}+\frac{{\Delta }^2{l}_{\mathrm{Di}}\mathrm{si}{\mathrm{n}}^2{\theta }_i}{12y_{\mathrm{eD}}}\right]+\frac{2x_{\mathrm{eD}}}{\pi \sin {\theta }_i\Delta l_{\mathrm{Di}}}\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)}{n^2\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\times \left[\begin{array}{c}\cot {\theta }_i\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\sin \frac{n\pi \left(x_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ -\cot {\theta }_i\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_{\mathrm{Di}}-y_D-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\sin \frac{n\pi \left(x_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ \begin{array}{c}+\beta \frac{\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_{\mathrm{Di}}-y_D-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\cos \frac{n\pi \left(x_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ -\beta \frac{\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\cos \frac{n\pi \left(x_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\end{array}\end{array}\right].$$(A.10)

1. When $y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_{i}0$, equation (A.9) can be rewritten as:

$$\begin{array}{c}{F}_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right)=\frac{\pi \beta }{x_{\mathrm{eD}}}\left[\frac{2y_{\mathrm{eD}}}{3}-2y_D+\frac{\left(y_D^2+y_{\mathrm{Di}}^2\right)}{y_{\mathrm{eD}}}+\frac{\Delta l_{\mathrm{Di}}^2\mathrm{si}{\mathrm{n}}^2{\theta }_i}{12y_{\mathrm{eD}}}\right]+\frac{2x_{\mathrm{eD}}}{\Delta l_{\mathrm{Di}}\pi \sin {\theta }_i}\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)}{n^2\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\cdot \\ \left[\begin{array}{c}\cot {\theta }_i\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D-y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\sin \frac{n\pi \left(x_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ -\cot {\theta }_i\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D-y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\sin \frac{n\pi \left(x_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ \begin{array}{c}+\beta \frac{\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D-y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]-\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\cos \frac{n\pi \left(x_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ +\beta \frac{\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]-\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D-y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\cos \frac{n\pi \left(x_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\end{array}\end{array}\right]\end{array}.$$(A.11)

1. When $y_{\mathrm{Di}}-y_D-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i<0<y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i$, equation (A.9) can be rewritten as:

$$\begin{array}{c}{F}_D\left(x_D,y_D,x_{\mathrm{Di}},y_{\mathrm{Di}},\beta \right)=\frac{\pi \beta }{x_{\mathrm{eD}}}\left[\frac{2y_{\mathrm{eD}}}{3}-\left(y_D+y_{\mathrm{Di}}\right)+\frac{\left(y_D^2+y_{\mathrm{Di}}^2\right)}{y_{\mathrm{eD}}}+\frac{\Delta l_{\mathrm{Di}}^2\mathrm{si}{\mathrm{n}}^2{\theta }_i}{12y_{\mathrm{eD}}}-\frac{{\left(y_{\mathrm{Di}}-y_D\right)}^2}{\Delta l_{\mathrm{Di}}\sin {\theta }_i}-\frac{\Delta l_{\mathrm{Di}}\sin {\theta }_i}{4}\right]+\frac{2x_{\mathrm{eD}}}{\Delta l_{\mathrm{Di}}\pi \sin {\theta }_i}\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\frac{\cos \left(\frac{n\pi x_D}{x_{\mathrm{eD}}}\right)}{n^2\sinh \beta \frac{n\pi y_{\mathrm{eD}}}{x_{\mathrm{eD}}}}\\ \left[\begin{array}{c}\cot {\theta }_i\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\sin \frac{n\pi \left(x_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ -\cot {\theta }_i\frac{\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\cosh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D-y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\sin \frac{n\pi \left(x_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ \begin{array}{c}-\beta \frac{\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_{\mathrm{Di}}-y_D+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]+\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\cos \frac{n\pi \left(x_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\\ -\beta \frac{\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D-y_{\mathrm{Di}}+\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]-\sinh \beta \frac{n\pi }{x_{\mathrm{eD}}}\left[y_{\mathrm{eD}}-\left(y_D+y_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\sin {\theta }_i\right)\right]}{{\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i}\cos \frac{n\pi \left(x_{\mathrm{Di}}-\frac{\Delta l_{\mathrm{Di}}}{2}\cos {\theta }_i\right)}{x_{\mathrm{eD}}}\end{array}\end{array}\right]+\frac{2\pi \beta x_{\mathrm{eD}}}{\Delta l_{\mathrm{Di}}\sin {\theta }_i\left({\beta }^2+\mathrm{co}{\mathrm{t}}^2{\theta }_i\right)}\\ \left[\frac{1}{3}-\frac{\left[\left|x_{\mathrm{Di}}+x_D+\left(y_D-y_{\mathrm{Di}}\right)\cot {\theta }_i\right|+\left|x_{\mathrm{Di}}-x_D+\left(y_D-y_{\mathrm{Di}}\right)\cot {\theta }_i\right|\right]}{2x_{\mathrm{eD}}}+\frac{{\left[\left(x_{\mathrm{Di}}+x_D+\left(y_D-y_{\mathrm{Di}}\right)\cot {\theta }_i\right)\right]}^2+{\left[\left(x_{\mathrm{Di}}-x_D+\left(y_D-y_{\mathrm{Di}}\right)\cot {\theta }_i\right)\right]}^2}{4x_{\mathrm{eD}}^2}\right]\end{array}.$$(A.12)

References

All Tables

All Figures