Numerical Simulation of Multifracture Growth under Extremely Limited Entry Fracturing of Horizontal Well

Abstract

The multifracture competitive growth from a horizontal well is an essential issue in multi-cluster fracturing design. In recent years, extremely limited entry (ELE) fracturing has been implemented to promote uniform multifracture growth. However, the mechanism of multifracture growth and ELE design remain unclear. Based on the planar three-dimensional multifracture propagation model, a multi-cluster horizontal well fracturing model that considers ELE design has been developed. The model considers flow in the wellbore and fluid filtration loss in the fracture. The simulator enables the simulation and analysis of non-uniform in situ stress, filtration loss, and fracture properties. Using this program, we simulated the propagation process of multiple clusters of fractures in ELE fracturing of horizontal wells. The results show the following: The perforation friction in the ELE fracturing can counteract the difference in fluid allocation caused by stress interference, allowing all clusters of perforations to have even fluid allocation but to differ significantly in fracture geometry. The in situ stress profile and 3D fracture stress interference determine the fracture geometry, and the fracture of the middle cluster could cross through the layer with relatively higher in situ stress, resulting in a decrease in effective fracture area in the pay zone. Furthermore, an increase in perforation diameter causes the flow-limiting effect of the perforations to decrease. The fluid volumes entering different clusters of perforations become less uniform. The difference in fracture toughness within a perforated stage has a minor influence on the fluid allocation between different clusters, while the in situ stress distribution within a perforated stage has a significant impact on the fluid allocation between different perforation clusters in the stage. Fractures preferentially propagate at the perforation points with lower in situ stress and stress interference. This study can be helpful to understand multifracture competitive growth and the optimization of ELE fracturing design.

1. Introduction

Multi-cluster fracturing technology for petroleum and geothermal energies is the primary technology of reservoir stimulation under the current low oil prices [1,2]. Temperature and acoustic monitoring with optical fiber and perforation imaging after the multi-cluster fracturing of horizontal wells have revealed non-uniform fluid allocation between different clusters of fractures [3,4,5]. Non-uniform fluid allocation between different perforation clusters may lead to a fast extension of some fractures and a premature stop in others, and, thus, a decrease in stimulated reservoir volume, making fracturing design challenging [6].

The multifracture propagation mechanism is currently mainly analyzed theoretically by numerical simulation. The main simulation methods include the finite element, extended finite element, discrete element, and boundary element method [7]. The finite element, extended finite element, and discrete element methods are full-field discretization techniques with high storage and computational costs. In contrast, the boundary element method only requires discretization at the fractures and, thus, is less computationally intensive and is better for engineering applications, particularly for fracture simulation [8]. It has been found that reservoir heterogeneity and stress interference between multiple fractures are the primary factors that affect the uneven propagation of multiple fractures [9,10,11,12]. To promote the simultaneous propagation of multiple fractures, Shell proposed an extremely limited entry (ELE) operation in 2018 [12]. In this process, the perforation friction is significantly increased and the number of perforations is limited to 2–3/clusters to achieve uniform fracture initiation and fluid allocation between different fracture clusters. To date, the effectiveness of limited entry fracturing is mainly based on simulation results from two-dimensional (2D) and pseudo-three-dimensional (P3D) fracture propagation models, and the growth pattern of multi-cluster fractures in ELE fracturing in three-dimensional fracture morphology has scarcely been studied [13,14]. Meanwhile, the ELE fracturing process was tested in China in the fracturing operations of Shengli shale oil and Fuling shale gas reservoirs [14]. Analysis methods for this fracturing process must be urgently developed to demonstrate their effectiveness.

The hydraulic fracture simulation model is an effective way to optimize fracturing design [15]. It contains 2D, P3D, planar 3D (PL3D), and fully 3D models. The 2D and P3D models are typically used in the industry due to their high efficiency, while the PL3D and fully 3D models have increased accuracy but are less efficient. A comprehensive review of the four kinds of models can be found in [16] and [17]. Recently, drilling cores and distributed acoustic sensing from adjacent wells show that fractures are mostly planar in the reservoir scale [18,19,20,21]. These results demonstrate that using PL3D is the best choice for optimizing fracturing design due to its accuracy and efficiency [22].

Therefore, a planar 3D multifracture propagation model was proposed in this work based on our previous work [22,23,24,25], the extension of multiple fractures in extremely limited entry fracturing was simulated numerically using parameters of a shale oil well to demonstrate the effectiveness of this fracturing method in the hope of providing a reference basis for fracturing design.

The paper is organized as follows. The mathematical model of the planar 3D multifracture is given in Section 2. Section 3 presents the numerical scheme for the mathematical model. A case study from the Shengli shale oil formation is then conducted and analyzed in Section 3. A parametric sensitivity is conducted in Section 4. Finally, conclusions are drawn in Section 5.

2. Mathematical Model

Fracture models are simplifications of actual fractures, so recognizing real fracture patterns is the first step in fracture modeling. Many observations and interpretations have been performed to find patterns of fractures generated in fracturing. The most reliable method to find fracture patterns is the observation of drilled cores containing hydraulic fractures. Many planar fractures were found in cores taken from shale gas wells after fracturing in the United States [20], complex fractures were only seen in local parts, so it is reasonable to simplify the fractures as planar fractures on a macroscopic level [21]. Meanwhile, unconventional oil and gas reservoirs in China generally have large horizontal stress differences [13], so it is difficult for them to be deflected theoretically, which further confirms that simplifying fractures into planar fractures is reasonable. Therefore, in this work, the fractures were assumed to extend along a fixed plane, i.e., the assumption of a planar three-dimensional fracture model was adopted.

2.1. Governing Equations

Hydraulic fracturing is a fracture-extension process driven by high-pressure viscous fluid and consists of four main physical processes: rock deformation, fluid flow in the fracture, fluid filtration loss to the formation, and fracture propagation at the fracture tip [26]. For multi-cluster fracturing of a horizontal well, it also includes fluid distribution through the wellbore. The geometric model is shown in Figure 1.

2.1.1. Elasticity

A three-dimensional displacement discontinuity method (DDM) was used to calculate the rock deformation. The advantage of DDM is that it reduces the problem by one dimension, so the computation burden is less than FEM or extended FEM. However, the disadvantage of the DDM used in our model is based on 3D Green’s function of dislocation dipoles in a linear elastic media with uniform elastic properties, so the heterogeneous elastic properties cannot be considered. In contrast, such properties can be easily handled by FEM or extended FEM. In our model, we use the boundary element method to solve solid elasticity, as it is suitable for the fracture problems in an infinite domain.

The amount of shear displacement discontinuity for planar fracture extension is zero. Therefore, the relationship between pressure in a fracture and the fracture width can be described by Equation (1) [27]:

$$p_{\mathrm{f}}-{\sigma }_{\mathrm{h}}={\displaystyle {\int }_{A\left(t\right)}C\left(x-x^{\prime },y-y^{\prime },z-z^{\prime }\right)w\left(x^{\prime },y^{\prime },z^{\prime }\right)dA}$$

(1)

where p_f is the fluid pressure, MPa; σ_h is the minimum horizontal principle stress, MPa; C is the Green’s function; t is the time, s; A(t) is the area of opened fracture at moment t, m²; w is the fracture width, m; (x′, y′, z′) are the coordinates of the source point; and (x, y, z) are the coordinates of the field point.

The specific form of Green’s function C is [27]

$$C=\frac{E}{8\pi \left(1-v^2\right)}\frac{1}{{\left[{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z-z^{\prime }\right)}^2\right]}^{1.5}}$$

(2)

where E is Young’s modulus, MPa, and v is Poisson’s ratio, dimensionless.

2.1.2. Fluid Flow in a Fracture

Assuming that the fluid flow inside the fracture conforms to Poiseuille flow [28], the constitutive equation of fluid flow inside the fracture is

$$\mathit{q}=-\frac{w^3}{12\mu }\nabla p_{\mathrm{f}}$$

(3)

where q is the volume flow rate, m²/s, and μ is the fluid viscosity, Pa·s.

The component form of Equation (2) is

$$\{\begin{array}{l}{q}_x=-\frac{w^3}{12\mu }\frac{\partial p_{\mathrm{f}}}{\partial x}\\ q_y=-\frac{w^3}{12\mu }\frac{\partial p_{\mathrm{f}}}{\partial y}\end{array}$$

(4)

where q_x and q_y are the volume flow rates in the x direction and the y direction, respectively.

The continuity equation of fluid flow in the fracture is

$$\frac{\partial w}{\partial t}+\nabla \cdot \mathit{q}+\frac{2C_l}{\sqrt{t-t_0}}=\delta \left(x-x_{\mathrm{in},k},y-y_{\mathrm{in},k},z-z_{\mathrm{in},k}\right)Q_k, k=1,2,\dots ,n_{\mathrm{f}}$$

(5)

where Q_k is the flow rate at the inlet of fracture k, m³/s; C_l is the filtration loss coefficient, m/min^0.5; t₀ is the moment when filtration loss starts, min; n_f is the number of fractures; and (x_in,k, y_in,k, z_in,k) are the coordinates of the inlet of fracture k.

Substituting Equation (3) into Equation (5), the following transport equation is obtained:

$$\frac{\partial w}{\partial t}=\nabla \cdot \left(\frac{w^3}{12\mu }\nabla p_{\mathrm{f}}\right)+\frac{2C_l}{\sqrt{t-t_0}}+\delta \left(x-x_{\mathrm{in},k},y-y_{\mathrm{in},k},z-z_{\mathrm{in},k}\right)Q_k, k=1,2,\dots ,n_{\mathrm{f}}$$

(6)

2.1.3. Wellbore Conditions

In multi-cluster fracturing, the injected fracturing fluid passes through the wellbore and perforations into different clusters of fractures. The amount of fluid entering each cluster is controlled by the “wellbore–perforation–fracture” system. The pressure at the inlet of each fracture cluster satisfies [29]

$$p_w=p_{\mathrm{p},k}+p_{\mathrm{t},k}+p_{\mathrm{in},k}, k=1,2,\dots ,n_{\mathrm{f}}$$

(7)

where p_t,k is the wellbore flow friction of fracture k, MPa; p_p,k is the friction of perforation corresponding to fracture k, MPa; and p_in,k is the pressure at the inlet of fracture k, MPa.

The friction of perforation is calculated by [30]

$$p_{\mathrm{p},k}=\frac{0.807\rho Q_k^2}{n_k^2{d}_k^4{K}_k^2}$$

(8)

The flow friction in the wellbore is calculated by Churchill’s entire flow equation [31].

At the same time, the flow rates of all clusters satisfy the mass conservation law, i.e.,

$$Q_{\mathrm{T}}={\displaystyle \sum _{I=1}^{n_{\mathrm{f}}}{Q}_{\mathrm{in},k}}$$

(9)

2.2. Initial Boundary Condition

The flow at inlet of each fracture satisfies

$${Q|}_{\mathrm{in}}=-\frac{w^3}{12\mu }\frac{\partial p}{\partial n}=Q_{\mathrm{in},k}$$

(10)

The flow at the fracture boundary is zero, i.e.,

$${Q|}_{\partial A}=-\frac{w^3}{12\mu }\frac{\partial p}{\partial n}=0$$

(11)

When the stress intensity factor at the fracture tip reaches the rock fracture toughness, the fracture starts to grow.

$$K_{\mathrm{tip}}=K_{\mathrm{Ic}}$$

(12)

where K_Ic is the rock type I fracture toughness, MPa·m^0.5.

The calculation equation of the tip stress intensity factor is [32]

$$K_{\mathrm{tip}}=\frac{0.806E\sqrt{\pi }{w}_{\mathrm{tip}}}{4\left(1-v^2\right)\sqrt{d}}$$

(13)

where w_tip is the width of the fracture tip cell, m, and d is the length of the tip cell, m.

3. Numerical Scheme

A fixed grid (rectangular cell-structured grid) was used to calculate the fracture-extension process. The cell was labeled as (i,j,k), and the corresponding position was (x_i,y_j,z_k). There were four types of cells: channel cells, tip cells, cells to be activated, and inactive cells (Figure 2). The cell type of the grid system was updated each time by determining whether the tip cell met the propagation condition. The cell’s center point was the unknown quantity (width and pressure) solution point, and the cell boundary was the position of the flow rate.

All open cells and tip cells were labeled sequentially as I, and the number of cells at the current moment was Ne, I = 1, 2, ..., N_e. Equation (1) was discretized by the constant displacement discontinuity method.

$$p_I={\displaystyle \sum _{J=1}^{N_e}{C}_{IJ}{w}_J}+{\sigma }_h$$

(14)

where N_e is the number of cells, and the capital I and J cell number labels were used to distinguish from the overall cell serial number (i,j,k). C_IJ was derived from Equation (2).

The matrix form of Equation (14) is

$$\mathit{p}=\mathit{C}\mathit{w}+{\mathsf{\sigma }}_h$$

(15)

The discrete form of the flow Equation (6) was established by the finite volume method:

$$\begin{array}{l}\Delta w_{i,j,k}=\frac{\Delta t}{12\mu dx}\left[\left(w_{i+\frac{1}{2},j,k}^3\frac{p_{i+1,j,k}-p_{i,j,k}}{dx}-w_{i-\frac{1}{2},j,k}^3\frac{p_{i,j,k}-p_{i-1,j,k}}{dx}\right)\right]\\ +\frac{\Delta t}{12\mu dy}\left[\left(w_{i,j+\frac{1}{2},k}^3\frac{p_{i,j+1,k}-p_{i,j,k}}{dy}-w_{i,j-\frac{1}{2},k}^3\frac{p_{i,j,k}-p_{i,j-1,k}}{dy}\right)\right]-4C_l\sqrt{t-t_0}+\Delta t{Q}_k\end{array}$$

(16)

To calculate the wellbore flow distribution according to the wellbore conditions, we use

$$\mathit{F}\left(Q_1,Q_2,\dots ,Q_{n_{\mathrm{f}}},p_w\right)=\mathbf{0}$$

(17)

The specific components of Equation (17) are

$$\{\begin{array}{l}{F}_1\left(Q_1,Q_2,\dots ,Q_{n_{\mathrm{f}}},p_w\right)=p_w-\frac{0.807\rho Q_1^2}{n_1^2{d}_1^4{K}_1^2}-p_{t,1}-p_{in,1}\left(Q_1,Q_2,\dots ,Q_{n_{\mathrm{f}}},p_w\right)\\ \dots \\ F_k\left(Q_1,Q_2,\dots ,Q_{n_{\mathrm{f}}},p_w\right)=p_w-\frac{0.807\rho Q_k^2}{n_k^2{d}_k^4{K}_k^2}-p_{t,1}-p_{in,k}\left(Q_1,Q_2,\dots ,Q_{n_{\mathrm{f}}},p_w\right)\\ \dots \\ F_{n_{\mathrm{f}}}\left(Q_1,Q_2,\dots ,Q_N,p_w\right)=p_w-\frac{0.807\rho Q_N^2}{n_{n_{\mathrm{f}}}^2{d}_{n_{\mathrm{f}}}^4{K}_{n_{\mathrm{f}}}^2}-p_{t,1}-p_{in,n_{\mathrm{f}}}\left(Q_1,Q_2,\dots ,Q_{n_{\mathrm{f}}},p_w\right)\\ F_{n_{\mathrm{f}}+1}\left(Q_1,Q_2,\dots ,Q_N,p_w\right)=Q_T-{\displaystyle \sum _{k=1}^{n_{\mathrm{f}}}{Q}_k}\end{array}$$

(18)

The Newton–Raphson method was used to solve Equation (17), and the solution was iterated by coupling Equations (15) and (16) to calculate the flow rate of each cluster. The coupling equations were solved by the second-order Rough–Kutta–Legendre method. It uses strategically devised explicit RK stages to stabilize the solution, with larger time-steps than that of the explicit Euler finite difference scheme. The timestep was determined adaptively based on the fracture propagation velocity. The accuracy of the model was fully validated by penny fracture analytical solutions, hydraulic fracturing physical modeling experiments, and published literature results; see the details in References [22,23,24,25]. The model is solved programmatically using MATLAB (2020a), and the calculations are fully vectorized with high computational efficiency [33,34].

4. Results and Analysis of a Field Well

Parameters of a shale oil well X in Shengli Oilfield were taken for fracture extension calculation and analysis. The well was treated by ELE fracturing with slick water as the fracturing fluid. Each fracturing stage had six clusters of fractures at an average fracture spacing of 15 m. The locations of the horizontal well and fracture clusters are shown in Figure 3. The horizontal in situ stress difference at the wellbore is 12 MPa, so the fractures are mainly planar fractures. The perforation parameters were three perforations with 12 mm diameters in one cluster, discharge coefficients of 0.8, a fracturing fluid density of 1030 kg/m³, and an injection rate of 10 m³/min. Microfractures and natural fractures in the formation can increase the fluid filtration loss, so a larger filtration loss coefficient, 2 × 10⁻⁴ m/min^0.5, was taken.

4.1. Base Case

Figure 4 shows patterns of the fractures of the base case. To clearly demonstrate the patterns of the fractures, three-dimensional (Figure 4a) and planar (Figure 4b) graphs of all the fractures, and planar graphs (Figure 4c) of each fracture were drawn separately. The plane with y = 0 in Figure 3 was the perforation point (injection point). Since multiple fractures extend along the paths with the lowest energy dissipation of the system, the fractures have large differences in length and height. For example, the HF2, HF4, and HF5 fractures have larger heights, and nearly half of the fracture areas in the low-stress layer lie above the perforation section.

The variation coefficient of fluid allocation is defined as

$$V_d=\frac{V_{\max }-V_{\min }}{\overline{V}}\times 100\%$$

(19)

where $V_{\max } ,V_{min}, \mathrm{and} \overline{V}$ are the cluster with the maximum fluid volume allocation, the cluster with minimum fluid volume allocation, and the average fluid volume allocation of the clusters in m³, respectively. Compared with measuring the variation coefficient of fluid allocation with a standard deviation, this parameter is more intuitive.

Figure 5 shows the flow and fluid volume allocation of each cluster. The heel-side fracture had the maximum fluid volume, while the middle fracture had the lowest fluid volume. This observation results from the stress interference among multiple fractures. The middle fracture takes more interaction stresses induced by other fractures. According to Equation (19), the variation coefficient of fluid allocation was calculated at 9.1%.

Figure 6 shows the inlet width and area of each cluster of fractures. The results show that the fracture width generated by ELE fracturing ranged from 1.8 to 2.3 mm, with a range of 0.5 mm. However, the fractures had small differences in total area, as shown in Figure 6b.

Figure 7 shows the friction of each cluster perforation and bottom-hole pressure. The results show that the different clusters of perforations (three perforations with 12 mm diameters in one cluster) had friction values of 5–6 MPa at an injection rate of 10 m³/min. This friction could counteract the differences in fluid resistance caused by stress interference between the fractures. Therefore, the fluid volumes entering the fractures were relatively uniform. Despite the uniform fluid allocation, the different clusters of fractures vary widely in shape because fracture extension morphology is controlled by both in situ stress and stress interference between fractures [29]. Multiple fractures take the propagation paths with the lowest energy consumption (the embodiment of the principle of least action).

4.2. Effect of Perforation Diameter

In this section, the influence of perforation diameter was analyzed. Due to erosion of the fluid and proppant, perforations are subject to abrasion during the fracturing process and increase in diameter to varying degrees. Perforation abrasion models are mostly developed from the results of surface experiments, which are quite different from the monitoring of downhole abrasion. In this section, the effect of perforation diameter was analyzed by changing the perforation diameter from 10–= to 16 mm.

Figure 8 shows the extension patterns of six clusters of fractures at different perforation diameters. Figure 9 shows the patterns of fractures (HF1) at the heel at different perforation diameters.

Figure 10 shows the fluid allocations at different perforation diameters. Figure 11 shows the relationship between fluid allocation and the perforation diameter variation coefficient. Figure 10 and Figure 11 show that as the perforation diameter increases, the variance coefficient of fluid allocation increases. When the perforation diameter increases from 10 mm to 16 mm, the variation coefficient of fluid allocation increases from 5% to 30%. Perforations may increase by 50–120% in diameter after fracturing according to field monitoring, so even in extremely limited entry fracturing, fluid volumes entering different clusters of fractures may still be uneven.

4.3. Fracture Toughness Heterogeneity of Perforated Sections

Four cases were simulated to determine the effect of fracture toughness heterogeneity in a perforated stage. The heterogeneous fracture toughness was set to increase linearly at 0.25 MPa·m^0.5, 0.5 MPa·m^0.5, 0.75 MPa·m^0.5, and 1 MPa·m^0.5 per cluster from HF1 to HF6, respectively.

Figure 12 shows the patterns of fractures and fluid allocations in each of the four cases. Figure 13 shows the relationship between the heterogeneous fracture toughness and the variation coefficient of fluid allocation. Figure 13 shows that the variation coefficient of fluid allocation decreases with the increase in heterogeneous fracture toughness. However, the variation amplitude is small, so the influence of heterogeneous fracture toughness between clusters on fluid allocation between different clusters is also small. This is because the bottom-hole pressure is high in ELE fracturing. The high bottom-hole pressure is sufficient to counteract the fracture pressure differences caused by the heterogeneous fracture toughness of different clusters.

4.4. Heterogeneous In Situ Stress of the Perforated Section

Four different cases of heterogeneous in situ stress were simulated, in which the perforated section containing HF1 to HF6 fractures with heterogeneous fracture toughness was set to linearly increase the minimum horizontal principal stress. The stresses were set to 0.25 MPa, 0.5 MPa, 0.75 MPa, and 1 MPa per cluster from HF1 to HF6, respectively. Figure 14 shows the fracture geometry and fluid volume distribution among multiple fractures for different in-stage in-situ stress distributions. It can be observed that the multifracture geometry differs minorly among different in-situ stress distributions; whereas with the increase of in-situ stress in a perforation cluster, the fluid volume decreases significantly. Figure 15 shows that, for the well in this study, the variation coefficient of fluid allocation does not vary monotonically with the stress differences between clusters. When the stress difference between two clusters is less than 0.5 MPa, the variation coefficient of fluid allocation decreases with an increase in stress difference. In contrast, when the stress difference between two clusters is greater than 0.5 MPa, the variation coefficient of fluid allocation increases with an increase in stress difference between two clusters. This variation pattern is caused by an extremely limited entry fracturing operation mode and the profile of interlayer stress. When the stress difference between two clusters is less than 0.5 MPa, the friction of extremely limited entry flow can counteract the in situ stress difference between clusters, allowing the fluid allocation between clusters to be relatively uniform. In contrast, when the stress difference between clusters is higher, at 0.5 MPa for example, the fractures at the toe end section with higher in situ stress tend to propagate to the layer above the perforated section with lower in situ stress. When the fractures at the toe-end enter the layer with lower in situ stress, the stress interference between the fracture clusters reduces due to uneven widths of fractures longitudinally. This further reduces the difference in fluid allocation between the clusters. When the stress difference between the clusters is higher than 0.5 MPa, the difference in fluid allocation between the clusters is dominated by the stress difference between the clusters. Thus, the variation coefficient of fluid allocation increased with the increase in stress difference between the clusters.

4.5. Therotiecal Analysis of Fluid Allocation among Multiple Fractures

The parameter analysis above reveals that the fluid distribution among multiple fractures for ELE fracturing is primarily uniform due to high perforation friction. In this section, we provide a theoretical method to analyze this observation. The flux distribution of fluid flow is similar to current flow in an electrical circuit network, as shown in Figure 16. Without considering the wellbore friction, each fracture has the same pressure drop. The flow rate in each fracture is thus dependent on the flow resistance in the deformable fracture and perforations. The total pressure drop is calculated as

$$\Delta p=p_w-p_0=f_{p,k}{Q}_k^2+f_{f,k}{Q}_k$$

(20)

where $\mathsf{\Delta }p$ is the pressure drop, MPa; p_w is the wellbore pressure, MPa; p₀ is the pore pressure, MPa; $f_{p,k}$ is the coefficient of perforation friction for the k^th fracture; and $f_{hf,k}$ is the coefficient of fluid flow friction in the k^th fracture.

The resistance of flow in a hydraulic fracture should be solved by the solid-fluid coupling equations, where the stress shadowing effect and spatial variation of in situ stresses play an important role. The resistance of flow through perforations is calculated by Equation (8).

For the coupled system, the perforation friction, spatial variation of in situ stress, and flow friction in each fracture (dependent on stress shadowing effect, etc.) control the flow distribution among each fracture. For ELE fracturing, the perforation friction is much higher than the interaction stresses among multiple fractures. Thus, Equation (20) shows that, when all the fractures have the same perforation parameters, fluid allocation is relatively uniform.

5. Conclusions

In this work, a planar three-dimensional multifracture extension model coupling wellbore flow was developed, and the fracture extension dynamics in extremely limited entry fracturing were analyzed with this model. Through the study, the following conclusions have been reached:

(1)

The perforation friction in extremely limited entry fracturing of up to 5–6 MPa can counteract the difference in fluid allocation between fracture clusters caused by stress interference, resulting in an even fluid allocation between different clusters of perforations.

(2)

The path of hydraulic fracture reflects the principle of least action. In other words, fractures automatically grow along paths with the least resistance. Although fluid allocation between different fracture clusters in extremely limited entry fracturing is quite even, fractures of different clusters vary widely in geometry. The in situ stress profile and 3D fracture stress interference are the major factors affecting the patterns of fractures; fractures of the middle clusters could cross layers, causing the effective fracture area within the pay zone to decrease.

(3)

With an increase in perforation diameter, the flow limiting effect decreases, and the fluid allocation between clusters of fractures becomes less uniform.

(4)

The fracture toughness difference in a perforated section has little effect on the fluid allocation between different clusters of fractures in the section. In contrast, the stress distribution in a perforated section significantly affects the fluid allocation between the clusters of perforations. Fractures preferentially extend at perforation clusters with lower in situ stress and stress interference.

(5)

The reasonable perforation erosion model should be investigated to further analyze the effectiveness of ELE fracturing in future studies.