A Volume Fracturing Percolation Model for Tight Reservoir Vertical Wells

Abstract

Based on the non-linear seepage characteristics of tight reservoirs and the reconstruction mode of vertical wells with actual volume fracturing, a seven-area percolation model for volume fracturing vertical wells in tight reservoirs is established. Laplace transform and Pedrosa transform are applied to obtain analytical solutions of bottom hole pressure and vertical well production under a constant production regime. After verifying the correctness of the model, the influence of the fracture network parameters on the pressure and production is studied. The research results indicate that as the permeability modulus increases, the production of volume fracturing vertical wells decreases. The penetration ratio of the main crack and the half-length of the main crack have a small impact on production, while the diversion capacity of the main crack has a significant impact on the initial production, but it is ultimately limited by the effective volume of the transformation. Under constant pressure conditions, the greater the width and permeability of the ESRV region, the higher the vertical well production rate is. The smaller the aspect ratio of the ESRV region, the higher the mid-term yield and the faster the yield decrease. The research results show guiding significance for the design of vertical well volume fracturing in tight reservoirs.

1. Introduction

The development potential of tight oil reservoirs is enormous. According to incomplete statistics, the geological reserves of tight oil resources in China are 7.4–8 billion tons, which is of great significance to the Chinese petroleum industry. Due to the low permeability and strong heterogeneity of tight reservoirs, the production of tight oil is relatively low [1,2,3,4,5]. Volume fracturing is the main technical means for the transformation and development of tight oil reservoirs [6,7,8,9]. Unlike conventional fracturing techniques, volume fracturing can crush the reservoir, creating a secondary fracture network near the main fracture that has been fractured, increasing the volume of the transformation, shortening the fluid migration distance in the area, reducing seepage resistance, and improving oil recovery [10,11,12,13,14].

At present, the research on analytical models for volume fracturing of horizontal wells in tight reservoirs is relatively in-depth, and the types of models are diverse [15,16,17,18]. Brown et al. [19] established a three-area composite model for fracturing horizontal wells, with reservoirs divided into three areas: main fracture area, reformed area (SRV), and unreformed area, and through analytical solutions, the pressure and productivity dynamics of conventional oil reservoir fracturing wells were studied. Brohi et al. [20] established a three-area composite model for unconventional gas reservoirs based on the Brown model by characterizing the modified area with dual pore media and the unmodified area with single pore media. In addition, Stalgorova et al. [21] established a five-area composite model for horizontal wells to describe the situation where the main fractures of volume fracturing were not fully modified, and analyzed the problems in the theoretical derivation and practical application of linear flow models. On the basis of previous research, Zeng et al. [22] established a seven-area composite model for volume fracturing of horizontal wells in tight oil reservoirs. This model considers the heterogeneity of geological parameters in different main fracture areas and the situation where longitudinal fractures do not penetrate the reservoir. At the same time, the wellbore pressure drop loss between different fracture sections was calculated, making the model closer to the actual situation of the mine. Unlike volume fracturing horizontal wells, the analytical model for volume fracturing vertical wells is mainly based on the “composite radial flow” model. Zhu et al. [23] established a three-area composite model for volume fracturing vertical wells in tight oil reservoirs, which is divided into three regions: the main fracture area, the elliptical modified area, and the unreformed area. The transformed area adopts a fractal medium model to characterize the fracture network, while the unreformed area takes into account the starting pressure gradient of the matrix. The model applies the equivalent seepage resistance method to obtain a steady-state production capacity solution. Zhu et al. [24] found that the reconstruction area of vertical wells with volume fracturing is closer to a rectangle according to the microseismic cloud map, applied the linear flow model of horizontal wells with volume fracturing to vertical wells, and established a five-area model of vertical wells with volume fracturing that considers the reservoir stress sensitivity and has a rectangular reconstruction volume. The regional division is consistent with the five-area model of horizontal wells with volume fracturing. Although scholars have conducted extensive research on the seepage law of volume fracturing, there is relatively little research on multi-stage volume fracturing in vertical wells [25,26,27,28]. Currently, most of the seepage models for tight oil reservoirs with volume fracturing in vertical wells are radial composite models, and the assumption that the reconstruction area is a cylinder does not match the results of on-site microseismic monitoring [29,30,31,32]. Therefore, the establishment of a seven-area seepage model for volume fracturing of tight oil reservoirs is of great significance.

Tight reservoirs have dense lithology, porosity less than 10%, permeability less than 0.1 × 10⁻³ μm², and pore throat diameter at the micro-nanometer level, with obvious microscale effect, which leads to tight reservoirs presenting nonlinear seepage characteristics. This article divides the renovation area into seven seepage areas based on the nonlinear seepage characteristics of tight oil reservoirs and the renovation mode of actual volume fracturing vertical wells. An analytical modeling approach was used to establish a seepage model for seven areas of volume-fractured straight wells in tight reservoirs after reasonable simplification of the fracture by certain physical assumptions, and the pressure and production solutions under a constant production regime were obtained through Laplace transformation. The flow stages were divided and the influence of fracture network parameters on the seepage law was studied. In the actual exploration and development process, it is generally multi-stage combined production after segmented volume fracturing. For multi-stage combined production, the dynamic splitting method of multi-stage combined production and the Blasingame curve algorithm fitting technology can be used to invert the ESRV and other fracture network parameters. Due to the limitation of space, this paper does not cover the issue of multi-stage combined production, which will be discussed separately in the future.

2. Physical Model

Vertical well volume fracturing technology can “break” the reservoir, forming a transformation area composed of primary and secondary fracture networks. The complex fracture network generated increases the effective transformation volume (ESRV) and increases the fracture conductivity. The complex fracture network in the ESRV area of volume fracturing wells in tight oil reservoirs results in different seepage characteristics from porous media, requiring regional analysis. The author adopts a hypothetical method, assuming that the vertical well is located at the center of a rectangular sealed reservoir and that there is a finite conductivity fracture on the vertical well section. The model is divided into the following areas: the unreformed areas (areas 3, 4, 5, and 6) are single medium; the effective renovation area (area 1), is characterized by a dual medium model; the main crack area F. According to symmetry, one-eighth of the seepage area is taken for research and calculation, as shown in Figure 1.

Based on the volume fracturing characteristics of tight oil reservoirs, the following assumptions are made: ① The influence of starting pressure gradient is considered in the unreformed area, and the influence of permeability stress sensitivity of cracks is considered in the transformed area, and the main fracture area. ② Sealing of the outer boundary of the oil reservoir. ③ The main crack is a vertical crack with symmetrical wings, considering that the crack network does not fully penetrate the reservoir vertically. ④ The liquid is single-phase and slightly compressible. ⑤ The seepage process is isothermal seepage. ⑥ Neglecting the influence of gravity and capillary force. ⑦ The fluid flow direction, boundary, and coupling conditions in each region of the oil reservoir are:

Unreformed area: The outer boundary of area 6, which is not penetrated by the main crack and not compressed longitudinally, is closed, the inner boundary pressure is continuous, and the fluid flows to areas 4 and 2. The boundary conditions of area 5 and area 6 are the same, with fluid flowing towards areas 3 and 1. The outer boundary of area 4 is closed, and the inner boundary pressure is continuous, which is supplemented by the fluid in area 6 and flows towards area 2. The boundary conditions of area 3 are the same as those of area 4, supplemented by fluid from area 5 and flowing toward area 1. The outer boundary of area 2 is closed, and the inner boundary pressure is continuous, supplemented by fluids from areas 6 and 4 and flows towards area 1.

Effective transformation area: The flow rate at the outer boundary of area 1 is continuous, the pressure at the inner boundary is continuous, and the fracture system is supplemented by the matrix system and fluids from areas 5, 3, and 2, flowing towards the main fracture.

Main fracture area F: The fluid in the main fracture flows into the wellbore, with the outer boundary closed, and the inner boundary conditions determined by the production system.

3. Mathematical Model

3.1. Definition of Dimensionless Quantity

Dimensionless pressure (under fixed production conditions):

$$p_{j\mathrm{D}}=\frac{K_{\mathrm{ref}}h\left(p_i-p_j\right)}{1.842\times {10}^{-3}{q}_{\mathrm{F}}\mu B}$$

(1)

where K_ref is the reference permeability, μm², the subscript ref represents the reference value; h is the thickness of the formation, m; p_j is the pressure, MPa, the subscript j represents different partitions; q_F is the main fracture yield, m³·d⁻¹; and B is the volume coefficient of crude oil, dimensionless.

Dimensionless production:

$$q_{\mathrm{FD}}=\frac{q_{\mathrm{F}}\mu B}{2\mathsf{\pi }{K}_{\mathrm{ref}}{h}_{\mathrm{ref}}\left(p_i-p_{\mathrm{wf}}\right)}$$

(2)

where p_wf is the bottom hole flow pressure, MPa; h_ref is the reference thickness, m.

Dimensionless conductivity coefficient:

$${\eta }_{i\mathrm{D}}=\frac{{\eta }_i}{{\eta }_{\mathrm{ref}}}$$

(3)

$${\eta }_{\mathrm{ref}}=\frac{K_{\mathrm{ref}}}{u_{\mathrm{ref}}{\left(\varphi C_{\mathrm{t}}\right)}_{\mathrm{ref}}}$$

(4)

where μ_ref is the reference viscosity, mPa·s; ϕ is porosity, dimensionless; C_t is the comprehensive compression coefficient, MPa⁻¹.

Dimensionless time:

$$t_{\mathrm{D}}=\frac{3.6{\eta }_{\mathrm{ref}}t}{L_{\mathrm{ref}}^2}$$

(5)

where L_ref is the reference length, m; t is time, s.

Dimensionless starting pressure gradient:

$$G_{\mathrm{D}}=\frac{K_{\mathrm{ref}}{L}_{\mathrm{ref}}h{\lambda }_{\mathrm{m}}}{1.842\times {10}^{-3}{q}_{\mathrm{F}}\mu B}$$

(6)

where λ_m is the starting pressure gradient, MPa·m⁻¹.

Dimensionless distance:

$$\begin{array}{l}{x}_{\mathrm{FD}}=x_{\mathrm{f}}/L_{\mathrm{ref}}\hfill \\ x_{\mathrm{eD}}=x_{\mathrm{e}}/L_{\mathrm{ref}}\hfill \\ y_{1\mathrm{D}}=y_1/L_{\mathrm{ref}}\hfill \\ y_{\mathrm{eD}}=y_{\mathrm{e}}/L_{\mathrm{ref}}\hfill \\ z_{1\mathrm{D}}=z_1/L_{\mathrm{ref}}=h_{\mathrm{F}}/\left(2L_{\mathrm{ref}}\right)\hfill \\ z_{2\mathrm{D}}=z_2/L_{\mathrm{ref}}=h/\left(2L_{\mathrm{ref}}\right)\hfill \\ w_{\mathrm{D}}=w_{\mathrm{F}}/L_{\mathrm{ref}}\hfill \end{array}$$

(7)

where x_f is the half-length of the main crack, m; x_e is the distance from the well in the x-direction to the reservoir boundary, m; y₁ is the half-width of the effective renovation area, m; y_e is the distance from the well in the y-direction to the reservoir boundary, m; z₁ is the half-height of the main crack, m; h_F is the main crack height, m, the subscript F represents the main crack; z₂ is the half-height of the reservoir, m; and w_F is the half-height of the reservoir, m.

Dimensionless fracture conductivity:

$$F_{\mathrm{CD}}=\frac{K_{\mathrm{F}i}{w}_{\mathrm{F}}}{K_{\mathrm{ref}}{L}_{\mathrm{ref}}}$$

(8)

Crossflow coefficient:

$$\lambda =\alpha \frac{K_{1\mathrm{m}}}{K_{1\mathrm{fi}}}{L}_{\mathrm{ref}}^2$$

(9)

where K_1m is the matrix permeability of area 1, μm²; α is the form factor, dimensionless; K_1fi is the initial permeability of the secondary crack in the area 1, μm².

Elastic storage capacity ratio:

$$\omega =\frac{{(\varphi C_{\mathrm{t}})}_{1\mathrm{f}}}{{(\varphi C_{\mathrm{t}})}_{1\mathrm{f}}+{(\varphi C_{\mathrm{t}})}_{1\mathrm{m}}}$$

(10)

Dimensionless starting pressure gradient:

$${\lambda }_{\mathrm{D}}=C_{\mathrm{L}}{\lambda }_{\mathrm{m}}{L}_{\mathrm{ref}}$$

(11)

where C_L is the compressibility coefficient of the liquid volume, MPa⁻¹.

3.2. Mathematical Model Establishment and Solution

3.2.1. Area 6

This area is characterized by fractures that do not penetrate the reservoir, assuming a one-dimensional flow in the z-direction with a starting pressure gradient.

According to the principle of material balance, the continuity equation is obtained as follows:

$$\frac{\partial ({\rho }_6{v}_6)}{\partial z}=-\frac{\partial ({\rho }_6{\varphi }_6)}{\partial t}$$

(12)

where ρ is the fluid density, kg·m⁻³; v is the seepage rate m·h⁻¹.

The motion equation considering the starting pressure gradient is:

$$v_6=-\frac{K_6}{\mu }(\frac{\partial {\rho }_6}{\partial z}-{\lambda }_{\mathrm{m}})$$

(13)

The equation of state for rocks and fluids is:

$${\rho }_6={\rho }_0[1+C_{\mathrm{L}}(p_6-p_0)]$$

(14)

$${\varphi }_6={\varphi }_0+C_{\mathrm{f}}(p_6-p_0)$$

(15)

By substituting the motion equation and rock fluid state equation into the continuity equation, the seepage control equation for this area can be obtained as follows:

$$\frac{{\partial }^2{p}_6}{\partial z^2}-{\lambda }_{\mathrm{m}}{C}_{\mathrm{L}}\frac{\partial p_6}{\partial z}=\frac{\varphi \mu C_{\mathrm{t}6}}{K_6}\frac{\partial p_6}{\partial t}$$

(16)

The outer boundary conditions of a finite enclosed oil reservoir are:

$${\left.\frac{\partial p_6}{\partial z}\right|}_{z_{\mathrm{D}}=z_{2\mathrm{D}}}=0$$

(17)

The pressure at the interface between area 6 and areas 2 and 4 is continuous, and the internal boundary conditions are:

$${\left.p_6\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.p_2\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.p_4\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(18)

According to the definition of dimensionless variables, the equation is dimensionless and subjected to a Laplace transformation based on dimensionless time t_D, resulting in:

$$\frac{{\partial }^2{\overline{p}}_{6\mathrm{D}}}{\partial z_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{6\mathrm{D}}}{\partial z_{\mathrm{D}}}=\frac{s}{{\eta }_{6\mathrm{D}}}{\overline{p}}_{6\mathrm{D}}$$

(19)

$${\left.\frac{\partial {\overline{p}}_{6\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{2\mathrm{D}}}=0$$

(20)

$${\left.{\overline{p}}_{6\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.{\overline{p}}_{2\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.{\overline{p}}_{4\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(21)

where s is the Laplace space variable.

The solution of the seepage differential equation can be obtained by combining the above equation:

$${\overline{p}}_{6\mathrm{D}}={\left.c_6{\overline{p}}_{2\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.c_6{\overline{p}}_{4\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(22)

Among them:

$$c_6=\frac{\left(-\frac{r_{62}}{r_{61}}\right){\mathrm{e}}^{r_{61}\left(z_{\mathrm{D}}-z_{2\mathrm{D}}\right)}+{\mathrm{e}}^{r_{62}\left(z_{\mathrm{D}}-z_{2\mathrm{D}}\right)}}{\left(-\frac{r_{62}}{r_{61}}\right){\mathrm{e}}^{r_{61}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}+{\mathrm{e}}^{r_{62}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}}$$

$$r_{61}=\frac{{\lambda }_{\mathrm{D}}+\sqrt{{\lambda }_{\mathrm{D}}^2+\frac{4s}{{\eta }_{6\mathrm{D}}}}}{2},r_{62}=\frac{{\lambda }_{\mathrm{D}}-\sqrt{{\lambda }_{\mathrm{D}}^2+\frac{4s}{{\eta }_{6\mathrm{D}}}}}{2}$$

Taking the partial derivative of ${\overline{p}}_{6\mathrm{D}}$ at $z_{\mathrm{D}}=z_{1\mathrm{D}}$ yields:

$${\left.\frac{\partial {\overline{p}}_{6\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.-{\beta }_6{\overline{p}}_{2\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.-{\beta }_6{\overline{p}}_{4\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(23)

Among them:

$${\beta }_6=-\frac{-r_{62}{\mathrm{e}}^{r_{61}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}+r_{62}{\mathrm{e}}^{r_{62}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}}{\left(-\frac{r_{62}}{r_{61}}\right){\mathrm{e}}^{r_{61}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}+{\mathrm{e}}^{r_{62}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}D}\right)}}$$

3.2.2. Area 5

The seepage control equation for region 5 is:

$$\frac{{\partial }^2{\overline{p}}_{5\mathrm{D}}}{\partial z_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{5\mathrm{D}}}{\partial z_{\mathrm{D}}}=\frac{s}{{\eta }_{5\mathrm{D}}}{\overline{p}}_{5\mathrm{D}}$$

(24)

$${\left.\frac{\partial {\overline{p}}_{5\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{2\mathrm{D}}}=0$$

(25)

$${\left.{\overline{p}}_{5\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.{\overline{p}}_{1\mathrm{fD}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.{\overline{p}}_{3\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(26)

The solution of the seepage differential equation is:

$${\overline{p}}_{5\mathrm{D}}={\left.c_5{\overline{p}}_{1\mathrm{fD}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.c_5{\overline{p}}_{3\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(27)

Among them:

$$c_5=\frac{\left(-\frac{r_{52}}{r_{51}}\right){\mathrm{e}}^{r_{51}\left(z_{\mathrm{D}}-z_{2\mathrm{D}}\right)}+{\mathrm{e}}^{r_{52}\left(z_{\mathrm{D}}-z_{2\mathrm{D}}\right)}}{\left(-\frac{r_{52}}{r_{51}}\right){\mathrm{e}}^{r_{51}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}+{\mathrm{e}}^{r_{52}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}}$$

$$r_{51}=\frac{{\lambda }_{\mathrm{D}}+\sqrt{{\lambda }_{\mathrm{D}}^2+\frac{4s}{{\eta }_{5\mathrm{D}}}}}{2},r_{52}=\frac{{\lambda }_{\mathrm{D}}-\sqrt{{\lambda }_{\mathrm{D}}^2+\frac{4s}{{\eta }_{5\mathrm{D}}}}}{2}$$

Taking the partial derivative of ${\overline{p}}_{5\mathrm{D}}$ at $z_{\mathrm{D}}=z_{1\mathrm{D}}$ yields:

$${\left.\frac{\partial {\overline{p}}_{5\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.-{\beta }_5{\overline{p}}_{1\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}={\left.-{\beta }_5{\overline{p}}_{3\mathrm{D}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}$$

(28)

Among them:

$${\beta }_5=-\frac{-r_{52}{\mathrm{e}}^{r_{51}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}+r_{52}{\mathrm{e}}^{r_{52}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}}{\left(-\frac{r_{52}}{r_{51}}\right){\mathrm{e}}^{r_{51}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}+{\mathrm{e}}^{r_{52}\left(z_{1\mathrm{D}}-z_{2\mathrm{D}}\right)}}$$

3.2.3. Area 4

The fluid in area 4 flows towards area 2 in the x-direction and is supplemented by fluid from area 6 in the z-direction. The seepage control equation is:

$$\frac{{\partial }^2{\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}}+\frac{K_6}{K_4{z}_{1\mathrm{D}}}\frac{\partial {\overline{p}}_{6\mathrm{D}}}{\partial z_{\mathrm{D}}}{|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}=\frac{s}{{\eta }_{4\mathrm{D}}}{\overline{p}}_{4\mathrm{D}}$$

(29)

order ${\alpha }_4=\frac{K_6{\beta }_6}{K_4{z}_{1\mathrm{D}}}+\frac{s}{{\eta }_{4\mathrm{D}}}$, and by simplifying it into the above equation, we can obtain:

$$\frac{{\partial }^2{\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}}={\alpha }_4{\overline{p}}_{4\mathrm{D}}$$

(30)

Closed oil reservoir, with external boundary conditions as follows:

$${\left.\frac{\partial {\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{e\mathrm{D}}}=0$$

(31)

The pressure at the interface between area 4 and area 2 is continuous, and the internal boundary conditions are:

$${\left.{\overline{p}}_{4\mathrm{D}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}={\left.{\overline{p}}_{2\mathrm{D}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}$$

(32)

By combining the internal and external boundary conditions and the seepage control equation, it is obtained that:

$${\overline{p}}_{4\mathrm{D}}={\left.c_4{\overline{p}}_{2\mathrm{D}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}$$

(33)

Among them:

$$c_4=\frac{\left(-\frac{r_{42}}{r_{41}}\right){\mathrm{e}}^{r_{41}\left(x_{\mathrm{D}}-x_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{42}\left(x_{\mathrm{D}}-x_{\mathrm{eD}}\right)}}{\left(-\frac{r_{42}}{r_{41}}\right){\mathrm{e}}^{r_{41}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{42}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}}$$

$$r_{41}=\frac{{\lambda }_{\mathrm{D}}+\sqrt{{\lambda }_{\mathrm{D}}^2+4{\alpha }_4}}{2},r_{42}=\frac{{\lambda }_{\mathrm{D}}-\sqrt{{\lambda }_{\mathrm{D}}^2+4{\alpha }_4}}{2}$$

Taking the partial derivative of ${\overline{p}}_{4\mathrm{D}}$ at $x_{\mathrm{D}}=x_{\mathrm{FD}}$ yields:

$${\left.\frac{\partial {\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}={\left.-{\beta }_4{\overline{p}}_{2\mathrm{D}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}$$

(34)

Among them:

$${\beta }_4=-\frac{-r_{42}{\mathrm{e}}^{r_{41}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}+r_{42}{\mathrm{e}}^{r_{42}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}}{\left(-\frac{r_{42}}{r_{41}}\right){\mathrm{e}}^{r_{41}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{42}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}}$$

3.2.4. Area 3

The fluid flows in the x-direction towards area 1 and is supplemented by fluid from area 5 in the z-direction. The seepage control equation is:

$$\frac{{\partial }^2{\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}}+\frac{K_5}{K_3{z}_{1\mathrm{D}}}\frac{\partial {\overline{p}}_{5\mathrm{D}}}{\partial z_{\mathrm{D}}}{|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}=\frac{s}{{\eta }_{3\mathrm{D}}}{\overline{p}}_{3\mathrm{D}}$$

(35)

Order ${\alpha }_3=\frac{K_5{\beta }_5}{K_3{z}_{1\mathrm{D}}}+\frac{s}{{\eta }_{3\mathrm{D}}}$, and by simplifying it into the above equation, we can obtain:

$$\frac{{\partial }^2{\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}}={\alpha }_3{\overline{p}}_{3\mathrm{D}}$$

(36)

Closed oil reservoir, with external boundary conditions as follows:

$${\left.\frac{\partial {\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{eD}}}=0$$

(37)

The pressure at the interface between area 3 and area 1 is continuous, and the internal boundary conditions are:

$${\left.{\overline{p}}_{3\mathrm{D}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}={\left.{\overline{p}}_{1\mathrm{fD}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}$$

(38)

By combining the internal and external boundary conditions and the seepage control equation, it is obtained that:

$${\overline{p}}_{3\mathrm{D}}={\left.c_3{\overline{p}}_{1\mathrm{fD}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}$$

(39)

Among them:

$$c_3=\frac{\left(-\frac{r_{32}}{r_{31}}\right){\mathrm{e}}^{r_{31}\left(x_{\mathrm{D}}-x_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{32}\left(x_{\mathrm{D}}-x_{\mathrm{eD}}\right)}}{\left(-\frac{r_{32}}{r_{31}}\right){\mathrm{e}}^{r_{31}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{32}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}}$$

$$r_{31}=\frac{{\lambda }_{\mathrm{D}}+\sqrt{{\lambda }_{\mathrm{D}}^2+4{\alpha }_3}}{2},r_{32}=\frac{{\lambda }_{\mathrm{D}}-\sqrt{{\lambda }_{\mathrm{D}}^2+4{\alpha }_3}}{2}$$

Taking the partial derivative of ${\overline{p}}_{3\mathrm{D}}$ at $x_{\mathrm{D}}=x_{\mathrm{FD}}$ yields:

$${\left.\frac{\partial {\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}={\left.-{\beta }_3{\overline{p}}_{1\mathrm{fD}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}$$

(40)

Among them:

$${\beta }_3=-\frac{-r_{32}{\mathrm{e}}^{r_{31}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}+r_{32}{\mathrm{e}}^{r_{32}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}}{\left(-\frac{r_{32}}{r_{31}}\right){\mathrm{e}}^{r_{31}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{32}\left(x_{\mathrm{FD}}-x_{\mathrm{eD}}\right)}}$$

3.2.5. Area 2

The fluid flows in the y-direction towards area 1 and is supplemented by fluid from area 6 in the z-direction, and area 4 in the x-direction. The seepage control equation is:

$$\frac{{\partial }^2{\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}+{\left.\frac{K_6}{K_2{z}_{1\mathrm{D}}}\frac{\partial {\overline{p}}_{6\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}+{\left.\frac{K_4}{K_2{x}_{\mathrm{FD}}}\frac{\partial {\overline{p}}_{4\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}=\frac{s}{{\eta }_{2\mathrm{D}}}{\overline{p}}_{2\mathrm{D}}$$

(41)

Order ${\alpha }_2=\frac{K_6{\beta }_6}{K_2{z}_{1\mathrm{D}}}+\frac{K_4{\beta }_4}{K_2{x}_{\mathrm{FD}}}+\frac{s}{{\eta }_{2\mathrm{D}}}$, and by simplifying it into the above equation, we can obtain:

$$\frac{{\partial }^2{\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}{}^2}-{\lambda }_{\mathrm{D}}\frac{\partial {\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}={\alpha }_2{\overline{p}}_{2\mathrm{D}}$$

(42)

Closed oil reservoir, with external boundary conditions as follows:

$${\left.\frac{\partial {\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{\mathrm{eD}}}=0$$

(43)

The pressure at the interface between area 2 and area 1 is continuous, and the internal boundary conditions are:

$${\left.{\overline{p}}_{2\mathrm{D}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}={\left.{\overline{p}}_{1\mathrm{fD}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}$$

(44)

By combining the internal and external boundary conditions and the seepage control equation, it is obtained that:

$${\overline{p}}_{2\mathrm{D}}={\left.c_2{\overline{p}}_{1\mathrm{fD}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}$$

(45)

Among them:

$$c_2=\frac{\left(-\frac{r_{22}}{r_{21}}\right){\mathrm{e}}^{r_{21}\left(y_{\mathrm{D}}-y_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{22}\left(y_{\mathrm{D}}-y_{\mathrm{eD}}\right)}}{\left(-\frac{r_{22}}{r_{21}}\right){\mathrm{e}}^{r_{21}\left(y_{1\mathrm{D}}-y_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{22}\left(y_{1\mathrm{D}}-y_{\mathrm{eD}}\right)}}$$

$$r_{21}=\frac{{\lambda }_{\mathrm{D}}+\sqrt{{\lambda }_{\mathrm{D}}^2+4{\alpha }_2}}{2},r_{22}=\frac{{\lambda }_{\mathrm{D}}-\sqrt{{\lambda }_{\mathrm{D}}^2+4{\alpha }_2}}{2}$$

Taking the partial derivative of ${\overline{p}}_{2\mathrm{D}}$ at $y_{\mathrm{D}}=y_{1\mathrm{D}}$ yields:

$${\left.\frac{\partial {\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}={\left.-{\beta }_2{\overline{p}}_{1\mathrm{fD}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}$$

(46)

Among them:

$${\beta }_2=-\frac{-r_{22}{\mathrm{e}}^{r_{21}\left(y_{1\mathrm{D}}-y_{\mathrm{eD}}\right)}+r_{22}{\mathrm{e}}^{r_{32}\left(y_{1\mathrm{D}}-y_{\mathrm{eD}}\right)}}{\left(-\frac{r_{22}}{r_{21}}\right){\mathrm{e}}^{r_{21}\left(y_{1\mathrm{D}}-y_{\mathrm{eD}}\right)}+{\mathrm{e}}^{r_{32}\left(y_{1\mathrm{D}}-y_{\mathrm{eD}}\right)}}$$

3.2.6. Area 1

This area is effectively modified with a developed secondary fracture network, where we assume a dual medium and use the W R quasi steady state model. The fluid in the fracture system of area 1 flows in the y-direction towards the main fracture area F, while being supplemented by fluids from area 5, area 3, area 2, and the matrix system. The fracture system is sensitive to permeability.

The continuity equation of the crack system in this area is:

$$-\frac{\partial ({\rho }_1{v}_{1\mathrm{f}})}{\partial y}+{\rho }_1{q}_{21}+{\rho }_1{q}_{31}+{\rho }_1{q}_{51}+{\rho }_1{q}_{1\mathrm{mf}}=\frac{\partial ({\rho }_1{\varphi }_{1\mathrm{f}})}{\partial t}$$

(47)

The motion equation considering permeability sensitivity is:

$$v_{1\mathrm{f}}=-\frac{K_{1\mathrm{fi}}{\mathrm{e}}^{\gamma (p_{1\mathrm{f}}-p_0)}}{\mu }\frac{\partial p_{1\mathrm{f}}}{\partial y}$$

(48)

where γ is the permeability modulus, MPa⁻¹.

The equation of state for rocks and fluids is:

$${\rho }_1={\rho }_0[1+C_{\mathrm{L}}(p_{1\mathrm{f}}-p_0)]$$

(49)

$${\varphi }_{1\mathrm{f}}={\varphi }_0+C_{\mathrm{f}}(p_{1\mathrm{f}}-p_0)$$

(50)

The supplementary items of source and exchange:

$$q_{21}=\frac{K_2}{\mu y_{1\mathrm{D}}}{\left.\frac{\partial p_2}{\partial y}\right|}_{y=y_{1\mathrm{D}}}$$

(51)

$$q_{31}=\frac{K_3}{\mu x_{\mathrm{FD}}}{\left.\frac{\partial p_3}{\partial x}\right|}_{x=x_{\mathrm{FD}}}$$

(52)

$$q_{51}=\frac{K_5}{\mu z_{1\mathrm{D}}}{\left.\frac{\partial p_5}{\partial z}\right|}_{z=z_{1\mathrm{D}}}$$

(53)

$$q_{1\mathrm{mf}}=\alpha \frac{K_{1\mathrm{m}}}{\mu }(p_{1\mathrm{m}}-p_{1\mathrm{f}})$$

(54)

The continuity equation of the matrix system is:

$$q_{1\mathrm{mf}}=-{\varphi }_{1\mathrm{m}}{C}_{\mathrm{t}1\mathrm{m}}\frac{\partial p_{1\mathrm{m}}}{\partial t}$$

(55)

The motion, state equation and supplementary equation are introduced into the continuity equation and dimensionless to obtain the seepage control equation:

Matrix system:

$$\lambda (p_{1\mathrm{fD}}-p_{1\mathrm{mD}})=\frac{1-\omega }{{\eta }_{1\mathrm{D}}}\frac{\partial p_{1\mathrm{mD}}}{\partial t_{\mathrm{D}}}$$

(56)

Crack system:

$$\begin{array}{l}{\mathrm{e}}^{-{\gamma }_{\mathrm{D}}{p}_{1\mathrm{fD}}}[\frac{{\partial }^2{p}_{1\mathrm{fD}}}{\partial y_{\mathrm{D}}{}^2}-{\gamma }_{\mathrm{D}}{(\frac{\partial p_{1\mathrm{fD}}}{\partial y_{\mathrm{D}}})}^2]+\frac{K_2}{K_{\mathrm{fi}}{y}_{1\mathrm{D}}}{\left.\frac{\partial p_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}+\frac{K_3}{K_{1\mathrm{fi}}{x}_{\mathrm{FD}}}{\left.\frac{\partial p_{3\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}+\\ \frac{K_5}{K_{1\mathrm{fi}}{z}_{1\mathrm{D}}}{\left.\frac{\partial p_{5\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}-\lambda (p_{1\mathrm{fD}}-p_{1\mathrm{mD}})=\frac{\omega }{{\eta }_{1\mathrm{D}}}\frac{\partial p_{1\mathrm{fD}}}{\partial t_{\mathrm{D}}}\end{array}$$

(57)

The flow at the interface between area 1 and area 2 is continuous, and the outer boundary conditions are:

$${\left.{\mathrm{e}}^{{\gamma }_{\mathrm{D}}{p}_{1\mathrm{fD}}}\frac{\partial p_{1\mathrm{fD}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}=\frac{K_2}{K_{1\mathrm{fi}}}{\left.\frac{\partial p_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}$$

(58)

The pressure at the interface between area 1 and the main crack area F is continuous, and the internal boundary conditions are:

$${\left.p_{1\mathrm{fD}}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}={\left.p_{\mathrm{FD}}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}$$

(59)

Due to the presence of crack stress sensitivity, the equation is nonlinear and difficult to solve directly. By applying the Pedrosa transformation, it can be transformed into a linear equation.

Transforming pressure p_1fD into functions of perturbation τ_1fD:

$$p_{1\mathrm{fD}}=-\frac{1}{{\gamma }_{\mathrm{D}}}\ln (1-{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}})$$

(60)

When τ_1fD is expanded to n-th order, there are:

$$\frac{1}{1-{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}}}=1+{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}}+{\gamma }_{\mathrm{D}}^2{\tau }_{1\mathrm{fD}}^2+\cdots {\gamma }_{\mathrm{D}}^n{\tau }_{1\mathrm{fD}}^n$$

(61)

$$-\frac{1}{{\gamma }_{\mathrm{D}}}\ln (1-{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}})={\tau }_{1\mathrm{fD}}+\frac{1}{2}{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}}^2+\frac{1}{3}{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}}^3+\cdots \frac{1}{n}{\gamma }_{\mathrm{D}}{\tau }_{1\mathrm{fD}}^n$$

(62)

Considering that γ_D is generally a small value, only the zero-order expansion to obtain τ_1fD0 can satisfy the accuracy requirement.

After performing zero-order Pedrosa transformation on p_1fD and Laplace transformation on t_D, the seepage control equation is obtained as follows:

$$\lambda ({\overline{\tau }}_{1\mathrm{fD}0}-{\overline{p}}_{1\mathrm{mD}})=\frac{(1-\omega )s}{{\eta }_{1\mathrm{D}}}{\overline{p}}_{1\mathrm{mD}}$$

(63)

$$\begin{array}{l}\frac{{\partial }^2{\overline{\tau }}_{1\mathrm{fD}0}}{\partial y_{\mathrm{D}}{}^2}+\frac{K_2}{K_{1\mathrm{fi}}{y}_{1\mathrm{D}}}{\left.\frac{\partial {\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}+\frac{K_3}{K_{1\mathrm{fi}}{x}_{\mathrm{FD}}}{\left.\frac{\partial {\overline{p}}_{3\mathrm{D}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}+\frac{K_5}{K_{1\mathrm{fi}}{z}_{1\mathrm{D}}}{\left.\frac{\partial {\overline{p}}_{5\mathrm{D}}}{\partial z_{\mathrm{D}}}\right|}_{z_{\mathrm{D}}=z_{1\mathrm{D}}}\\ -\lambda ({\overline{\tau }}_{1\mathrm{fD}0}-{\overline{p}}_{1\mathrm{mD}})=\frac{\omega s}{{\eta }_{1\mathrm{D}}}{\overline{\tau }}_{1\mathrm{fD}0}\end{array}$$

(64)

$${\left.\frac{\partial {\overline{\tau }}_{1\mathrm{fD}0}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}=\frac{K_2}{K_{1\mathrm{fi}}}{\left.\frac{\partial {\overline{p}}_{2\mathrm{D}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=y_{1\mathrm{D}}}$$

(65)

$${\left.{\overline{\tau }}_{1\mathrm{fD}0}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}={\left.{\overline{\tau }}_{\mathrm{FD}0}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}$$

(66)

Order ${\alpha }_1=\frac{K_2{\beta }_2}{K_{1\mathrm{fi}}{y}_{1\mathrm{D}}}+\frac{K_3{\beta }_3}{K_{1\mathrm{fi}}{x}_{\mathrm{FD}}}+\frac{K_5{\beta }_5}{K_{1\mathrm{fi}}{z}_{1\mathrm{D}}}+\frac{\omega s}{{\eta }_{1\mathrm{D}}}+\frac{\lambda (1-\omega )s}{\lambda {\eta }_{1\mathrm{D}}+(1-\omega )s}$, the above equation can be simplified as:

$$\frac{{\partial }^2{\overline{\tau }}_{1\mathrm{fD}0}}{\partial y_{\mathrm{D}}{}^2}={\alpha }_1{\overline{\tau }}_{1\mathrm{fD}0}$$

(67)

The solution for area 1 is obtained by combining the internal and external boundary conditions as follows:

$${\overline{\tau }}_{1\mathrm{D}0}=A_1\cosh [(y_{\mathrm{D}}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]+B_1\sinh [(y_{\mathrm{D}}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]$$

(68)

Among them:

$$A_1=\frac{{\left.{\overline{\tau }}_{\mathrm{FD}0}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}}{\cosh [(\frac{w_{\mathrm{D}}}{2}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]+b\sinh [(\frac{w_{\mathrm{D}}}{2}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]}$$

$$B_1=b{A}_1,b=-\frac{K_2{\beta }_2}{K_{1\mathrm{fi}}\sqrt{{\alpha }_1}}$$

Taking the partial derivative of ${\overline{\tau }}_{1\mathrm{D}0}$ at $y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}$ yields:

$${\left.\frac{\partial {\overline{\tau }}_{1\mathrm{D}0}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}=e_1\sqrt{{\alpha }_1}{\left.{\overline{\tau }}_{\mathrm{FD}0}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}$$

(69)

Among them:

$$e_1=\frac{b\cosh [(\frac{w_{\mathrm{D}}}{2}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]+\sinh [(\frac{w_{\mathrm{D}}}{2}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]}{\cosh [(\frac{w_{\mathrm{D}}}{2}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]+b\sinh [(\frac{w_{\mathrm{D}}}{2}-y_{1\mathrm{D}})\sqrt{{\alpha }_1}]}$$

3.2.7. Area F

The main fracture area, where fluid flows towards the wellbore in the x-direction and is supplemented by fluid from area 1 in the y-direction, has permeability sensitivity. The seepage control equation in this area is:

$${\mathrm{e}}^{-{\gamma }_{\mathrm{D}}{p}_{\mathrm{FD}}}[\frac{{\partial }^2{p}_{\mathrm{FD}}}{\partial x_{\mathrm{D}}{}^2}-{\gamma }_{\mathrm{D}}(\frac{\partial p_{\mathrm{FD}}}{\partial x_{\mathrm{D}}}){}^2]+\frac{2K_{1\mathrm{f}i}{\mathrm{e}}^{-{\gamma }_{\mathrm{D}}{p}_{1\mathrm{fD}}}}{K_{\mathrm{Fi}}{w}_{\mathrm{D}}}{\left.\frac{\partial p_{1\mathrm{fD}}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}=\frac{1}{{\eta }_{\mathrm{FD}}}\frac{\partial p_{\mathrm{FD}}}{\partial t_{\mathrm{D}}}$$

(70)

Assuming the crack tip is closed, the outer boundary condition is:

$${\left.{\mathrm{e}}^{{\gamma }_{\mathrm{D}}{p}_{\mathrm{FD}}}\frac{\partial p_{\mathrm{FD}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}=0$$

(71)

The oil well is produced at a fixed production rate, and the internal boundary condition is:

$${\mathrm{e}}^{-{\gamma }_{\mathrm{D}}{p}_{\mathrm{FD}}}{\left.\frac{\partial p_{\mathrm{FD}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=0}=-\frac{\mathsf{\pi }}{F_{\mathrm{CD}}}$$

(72)

where F_CD is the dimensionless conductivity of the main crack.

The internal boundary condition for constant pressure production in oil wells is:

$${\left.p_{\mathrm{FD}}\right|}_{x_{\mathrm{D}}=0}=1,{\mathrm{e}}^{-{\gamma }_{\mathrm{D}}{p}_{\mathrm{FD}}}{\left.\frac{\partial p_{\mathrm{FD}}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=0}=-\frac{\mathsf{\pi }{q}_{\mathrm{wfD}}}{F_{\mathrm{CD}}}$$

(73)

Similarly, the zero-order Pedrosa transformation is introduced to transform it into a linear equation and Laplace transformation is performed, resulting in:

$$\frac{{\partial }^2{\overline{\tau }}_{\mathrm{FD}0}}{\partial x_{\mathrm{D}}{}^2}+\frac{2K_{\mathrm{fi}i}}{K_{\mathrm{Fi}}{w}_{\mathrm{D}}}{\left.\frac{\partial {\overline{\tau }}_{1\mathrm{fD}0}}{\partial y_{\mathrm{D}}}\right|}_{y_{\mathrm{D}}=\frac{w_{\mathrm{D}}}{2}}=\frac{s}{{\eta }_{\mathrm{FD}}}{\overline{\tau }}_{\mathrm{FD}0}$$

(74)

$${\left.\frac{\partial {\overline{\tau }}_{\mathrm{FD}0}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=x_{\mathrm{FD}}}=0$$

(75)

Fixed production system:

$${\left.\frac{\partial {\overline{\tau }}_{1\mathrm{fD}0}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=0}=-\frac{\mathsf{\pi }}{F_{\mathrm{CD}}s}$$

(76)

Fixed pressure production system:

$${\left.{\overline{\tau }}_{\mathrm{FD}0}\right|}_{x_{\mathrm{D}}=0}=\frac{1-{\mathrm{e}}^{-{\gamma }_{\mathrm{D}}}}{s{\gamma }_{\mathrm{D}}},{\left.\frac{\partial {\overline{\tau }}_{\mathrm{FD}0}}{\partial x_{\mathrm{D}}}\right|}_{x_{\mathrm{D}}=0}=-\frac{\mathsf{\pi }{\overline{q}}_{\mathrm{wfD}}}{F_{\mathrm{CD}}}$$

(77)

Order ${\alpha }_{\mathrm{F}}=\frac{s}{{\eta }_{\mathrm{FD}}}-\frac{2e_1{K}_{1\mathrm{fi}}\sqrt{{\alpha }_1}}{K_{\mathrm{Fi}}{w}_{\mathrm{D}}}$, the above equation can be simplified as:

$$\frac{{\partial }^2{\overline{\tau }}_{\mathrm{FD}0}}{\partial x_{\mathrm{D}}{}^2}={\alpha }_{\mathrm{F}}{\overline{\tau }}_{\mathrm{FD}0}$$

(78)

The boundary conditions for simultaneous production can be solved as follows:

$${\overline{\tau }}_{\mathrm{FD}0}=\frac{\mathsf{\pi }\cosh [(x_{\mathrm{D}}-x_{\mathrm{FD}})\sqrt{{\alpha }_{\mathrm{F}}}]}{F_{\mathrm{CD}}s\sqrt{{\alpha }_{\mathrm{F}}}\sinh (x_{\mathrm{FD}}\sqrt{{\alpha }_{\mathrm{F}}})}$$

(79)

When $x_{\mathrm{D}}=0$, the dimensionless bottom hole pressure solution in the Lagrangian space for constant production rate is obtained as:

$${\overline{\tau }}_{\mathrm{wfD}}=\frac{\mathsf{\pi }}{F_{\mathrm{CD}}s\sqrt{{\alpha }_{\mathrm{F}}}\tanh (x_{\mathrm{FD}}\sqrt{{\alpha }_{\mathrm{F}}})}$$

(80)

By combining the boundary conditions of constant pressure, the dimensionless bottom hole production in the Lagrangian space can be solved as:

$${\overline{q}}_{\mathrm{wfD}}=\frac{(1-{\mathrm{e}}^{-{\gamma }_{\mathrm{D}}})F_{\mathrm{CD}}\sqrt{{\alpha }_{\mathrm{F}}}}{{\gamma }_{\mathrm{D}}\mathsf{\pi }s}\tanh (x_{\mathrm{FD}}\sqrt{{\alpha }_{\mathrm{F}}})$$

(81)

Perform Stehfest numerical inversion on the pressure solution in Laplace space to obtain the dimensionless perturbation pressure solution in real space:

$${\tau }_{\mathrm{wfD}}(t_{\mathrm{D}})=\frac{\ln 2}{t_{\mathrm{D}}}{\displaystyle \sum _{i=1}^N{V}_i{\overline{\tau }}_{\mathrm{wfD}}(\frac{\ln 2}{t_{\mathrm{D}}}i)}$$

(82)

Among them:

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

In engineering applications, N is an even number between 4 and 12, and is generally selected based on actual debugging.

Finally, the perturbed pressure solution is subjected to inverse Pedrosa transformation to obtain a dimensionless pressure solution under constant production conditions:

$$p_{\mathrm{wfD}}(t_{\mathrm{D}})=-\frac{\ln \{1-{\gamma }_{\mathrm{D}}{\tau }_{\mathrm{wfD}}(t_{\mathrm{D}})\}}{{\gamma }_{\mathrm{D}}}$$

(83)

Similarly, the dimensionless production solution under constant pressure conditions can be obtained:

$$q_{\mathrm{wfD}}(t_{\mathrm{D}})=\frac{\ln 2}{t_{\mathrm{D}}}{\displaystyle \sum _{i=1}^N{V}_i{\overline{q}}_{\mathrm{wfD}}(\frac{\ln 2}{t_{\mathrm{D}}}i)}$$

(84)

4. Model Validation

The established volume fracturing seven-area composite flow model was validated using analytical and numerical simulation methods.

4.1. Comparative Validation of Analytical Models

When the fractures in the volume fracturing seven-area composite flow model completely penetrate the reservoir, and the permeability sensitivity of main fractures and the secondary fracture network in the effectively modified area are not considered, neither is the starting pressure gradient of the unmodified area, and the model degenerates into the dual medium five-area linear flow model proposed by Stalgorova et al. Here, we compare the degraded seven-area composite flow model with the five-area linear flow model.

As shown in Figure 2, the pressure and pressure derivative results for both coincide exactly, proving the correctness of the pressure solution of the model.

4.2. Numerical Simulation Comparison Verification

Compare the seven-area model with numerical simulation software using numerical simulation methods. A numerical model was established using CMG, which can reflect the characteristics of seepage in actual volume fractured vertical wells. The numerical simulation model is shown in Figure 3. The effective renovation area adopts a dual medium model, while the unmodified area adopts a single pore medium to simulate. By using the zoning function, the stress sensitivity effect of cracks is set in the effective renovation area and the main crack area, and the starting pressure gradient is set in the original unmodified area. The other model data are shown in

Table 1. Table of basic parameters of numerical model.

Parameter	Numerical Value	Parameter	Numerical Value
Reservoir length/m	370	Initial permeability of ESRV fracture system/μm2	5 × 10−3
Reservoir width/m	150	Fracture porosity	0.2
Reservoir thickness/m	30	Initial permeability of main fracture/μm2	1
Main crack seam height/m	20	Half-length of main crack/m	75
Matrix porosity	0.1	Main crack width/m	0.01
Matrix permeability/μm2	1 × 10−3	ESRV half-width/m	35
Comprehensive compressibility coefficient of matrix	0.0001	Original formation pressure/MPa	20
Comprehensive compression coefficient of cracks/MPa−1	0.001	Bottom hole flowing pressure/MPa	15
Oil volume factor	1.1	Starting pressure gradient/(MPa·m−1)	0.02
Crude oil viscosity/(mPa·s)	5	Permeability modulus/MPa−1	0.01

, and Figure 4 shows the comparison between the theoretical calculation and numerical simulation yield results of the volume fracturing seven-area composite flow model (without considering skin factors). It can be seen that the yield curve is basically consistent, further verifying the correctness of the model yield solution.

4.3. Division of Flow Stages

The dimensionless bottom hole pressure and pressure derivative curve is drawn under the fixed production system of volume fracturing vertical wells in double logarithmic coordinates, as shown in Figure 5. It can be divided into six flow stages: ① In the linear flow stage of the main fracture and ESRV regional fracture network, the slope of the dimensionless pressure and pressure derivative curves is 0.25. ② Matrix fracture scouring phase in the ESRV area, where the pressure and pressure derivative curves begin to be non-parallel at the beginning of the phase and the pressure derivative becomes depressed. ③ In the linear flow stage of the ESRV area, the pressure and pressure derivatives are once again parallel to each other, with slopes of 0.5. ④ During the boundary flow stage of the ESRV area, the flow continues to extend along the crack network of the renovation area and reaches the boundary of the ESRV area. ⑤ Linear flow in unmodified areas, liquid supply to ESRV area in unmodified areas. ⑥ During the stage of reservoir boundary influence, the slope is related to the type of reservoir boundary, and the slope of a closed reservoir is 1.

5. The Influence of Mesh Parameters on Flow Patterns

We use Matlab to program and calculate the analytical solution of the model, analyzing the impact of different parameters on pressure and production. In the pressure curve graph, the relevant parameters with the symbol in the legend represent the pressure derivative curve.

5.1. Starting Pressure Gradient

From Figure 6, it can be seen that the starting pressure gradient has a significant impact on production capacity from the boundary flow stage of the ESRV region in the middle to late stages. In the early stage, the fracture network of the effective modification area reduced the seepage resistance and weakened the impact of the starting pressure gradient on dimensionless pressure, pressure reciprocal, and dimensionless production. As the starting pressure gradient increases in the middle and later stages, the pressure consumption in the unmodified area increases to maintain production under fixed production conditions. As a result, the dimensionless pressure rises faster, and the pressure derivative also increases accordingly. Under constant pressure conditions, a larger starting pressure gradient leads to greater seepage resistance of fluid flowing into the effective modified area in the unmodified area, and the dimensionless production is lower.

5.2. Permeability Modulus

From Figure 7, it can be seen that the permeability modulus measures the sensitivity of fractures to pressure changes. The stress sensitivity effect causes the permeability of primary and secondary fractures to decrease with the decrease in formation pressure during the production process. From the production curve, it can be seen that the stress sensitivity effect affects the entire reservoir development stage. As the permeability modulus increases, the degree of fracture closure increases, the permeability decreases, the production decreases, and the pressure derivative curve rises faster under constant pressure conditions.

5.3. Main Crack Penetration Ratio

From Figure 8, it can be seen that the fracture penetration ratio affects the pressure and productivity of volume fracturing wells starting from the boundary flow stage of the effective transformation area. When the length of the main crack and the width of the renovation area are fixed, the crack penetration ratio increases and the effective renovation volume increases. In the renovation area, the crack network develops, and the pressure conductivity is stronger than in the unreformed area. Therefore, under fixed production conditions, as the crack penetration ratio increases, the rise of dimensionless pressure and pressure derivative curves slows down, and the start time of upward warping becomes later. For constant pressure conditions, as the fracture penetration ratio increases, the degree of transformation increases in the reservoir, and there are more flow channels in the fracture network, resulting in higher oil well production.

5.4. Half-Length of Main Crack

From Figure 9, it can be seen that the half-length of the main fracture is similar to the penetration ratio of the main fracture, which directly determines the size of the effective reconstruction area. Starting from the linear flow stage of the reconstruction area, it affects the pressure and productivity of the volume fracturing well. When the height of the main crack and the width of the renovation area are constant, the half-length of the main crack increases, the rise of the dimensionless pressure and pressure derivative curve is slower for fixed production conditions, and the start time of upward warping is later. For constant pressure conditions, as the half-length of the main fracture increases, the oil drainage area becomes larger, and the production of the oil well increases. Unlike the limitation of reservoir thickness on fracture height, the variation range of fracture length is relatively large, and its impact on dimensionless pressure and production curve is more significant.

5.5. ESRV Area Width

From Figure 10, it can be seen that the width of the ESRV area has an impact starting from the matrix crack flow stage of the effective transformation area. As the width of the ESRV region increases, the dimensionless pressure and pressure derivative decrease under constant production conditions, while the dimensionless production increases under constant pressure conditions. From the production curve, it can be seen that the width of the secondary fracture network has a significant impact on production. When the ESRV width is 0, it is conventional fracturing. From the production curve, it can be seen that volume fracturing is significantly better than conventional fracturing, which can significantly improve the production capacity of tight oil reservoirs.

5.6. Main Fracture Conductivity

From Figure 11, it can be seen that the conductivity of the main fracture mainly affects the early and middle flow stages covered by the main fracture and the effective renovation area. As the fracture conductivity increases, the dimensionless pressure curve is lower in the early and middle stages and corresponds to higher production in the early and middle stages. From the production curve, it can be seen that the high conductivity of the main fracture can quickly increase the production capacity of the effective transformation area, and shorten the development cycle of tight oil reservoirs, but it is ultimately limited by the size of the effective transformation volume and the permeability of the formation matrix, indicating that the matching relationship between the three should be optimized during the development process.

5.7. ESRV Regional Permeability

From Figure 12, it can be seen that the permeability of the fracture network in the ESRV region has a significant impact on the linear flow and matrix flow of cracks in the effectively modified area characterized by dual media. The permeability of the fracture network represents the degree of transformation in the volume fracturing ESRV area. As the permeability of the fracture network increases, the liquid supply capacity of the ESRV area becomes stronger, the dimensionless pressure and pressure derivative become smaller, the dimensionless production becomes larger, and entering the boundary flow stage of the reformed area happens earlier, which is ultimately limited by the effective modification volume and formation matrix permeability.

5.8. ESRV Area Aspect Ratio

From Figure 13, it can be seen that under the same effective renovation volume, the aspect ratio of the ESRV area determines the starting time of the linear flow stage in the renovation area. The aspect ratio of the ESRV area is the ratio of the length of the main crack to the renovation bandwidth. With the increase in the aspect ratio, the effective renovation area is narrower and longer, the bandwidth is smaller, the pressure wave propagates to the ESRV bandwidth boundary earlier, contacts the unmodified area earlier, and enters the flow phase at the boundary of the modified area, which is reflected in the mid-term dimensionless pressure and pressure derivative curves start to rise earlier. From the production curve, it can be seen that as the aspect ratio of the ESRV region increases, the mid-term production becomes higher and the production declines faster. This is because the aspect ratio of the ESRV region is smaller (when the aspect ratio is maintained to be larger than 1), the area of direct communication between the ESRV region and the matrix of the unmodified region, i.e., the oil drainage area is larger, and the production capacity is released more quickly.

6. Conclusions

• Based on the reservoir transformation characteristics of actual volume fracturing vertical wells, the volume fracturing transformation area was finely divided, and a seven-area seepage mathematical model for volume fracturing vertical wells in tight oil reservoirs was established. The pressure solution under constant production conditions and production solution under constant pressure conditions were solved using methods such as Laplace transform and Pedrosa transform. The accuracy of the model was verified using analytical and numerical simulation methods, respectively.
• Based on the dimensionless pressure and pressure derivative curve, the flow of vertical well volume fracturing wells in tight oil reservoirs is divided into six stages: main fractures, linear network flow in the ESRV region, matrix fracture channeling flow in the ESRV region, linear flow in the ESRV region, boundary flow in the ESRV region, linear flow in the unmodified region, and reservoir boundary influence.
• The production of volume fracturing vertical wells is influenced by permeability modulus, main fracture conductivity, ESRV area width, and ESRV area permeability. The production decreases with the increase in permeability modulus. The high conductivity of the main fracture can quickly release the production capacity of the effectively transformed area, and shorten the development cycle of tight oil reservoirs, but it is ultimately limited by the size of the effectively transformed volume. Under constant pressure conditions, the yield increases with increasing ESRV area width and ESRV area permeability.
• By using well testing parameters, the seven-area seepage model for volume fracturing in tight reservoirs can be used to reverse calculate fracture parameters, providing guidance for reservoir development design.