Numerical Simulation of Multi-Cluster Fracture Propagation In Horizontal Wells With Limited-Entry Designs

The effective propagation of multi-cluster fractures in horizontal wells is the key to the development of unconventional reservoirs. Due to the influence of pressure drops at perforating holes and the stress shadow effect, it is difficult to predict the fracturing fluid distribution and fracture dimensions in a fracturing stage. In this paper, a two-dimensional fluid-solid coupling model for simultaneous propagation of multiple fractures is established, and fluid distributions and dimensions of multiple fractures are studied with respect to different perforation designs. The model combines the User Amplitude Curve Subroutine (UAMP) in ABAQUS and the cohesive zone model (CZM) to calculate the perforating friction, fluid distribution and fracture propagation behaviors. After the accuracy of this model is verified by the analytical solution, a group of simulation is conducted to compare fracture propagations when the conventional limited-entry method (CLE) and extreme limited-entry method (less than 5 perforations per cluster, XLE) are used. Simulation results show that the edge and sub-central fractures in CLE cases almost get all the fluid and effectively propagate; central fractures receive little fluid and hardly propagate. In XLE cases, the fluid distribution of each fracture is relatively uniform, but the fracture lengths within one fracturing stage is still uneven; however, only reducing numbers or radii of perforation holes cannot achieve the uniform fracture propagation, where diverters might be further needed. Findings of this study provide a reference for the perforation optimization of multi-cluster horizontal wells in the field.


Introduction
Horizontal drilling and multi-cluster fracturing are the mainstream technology for developing unconventional oil and gas reservoirs. During hydraulic fracturing, the fracturing fluid is injected with a high pressure to break the rock and generate multiple fractures with high conductivity in the reservoir; this can greatly improve the seepage condition of the reservoir and make unconventional oil and gas reservoirs economically recoverable. As one of the main technologies to achieve the simultaneous propagation of multi-cluster fractures, plug and perf is widely used in domestic and foreign oilfields, where bridge plugs separate the horizontal well into different stages and allow a certain number of clusters to initiate fractures within one fracturing stage [1][2] . In the current environment of low oil prices, fewer bridge plugs and more clusters within one stage are applied to reduce the stimulation cost. However, numerous data from the field using the distributed optical sensing techniques (DTS and/or DAS) have shown that fractures seldom propagate evenly within one stage, where the stress shadow effect prevents fracture propagating when the shorter fracture spacing is used, and uneven distribution of proppants also leads to the heel-side biased fracture growth [3][4][5][6][7][8] . Production log even shows that 80% of production came from 20% fractures, and about 30% fractures had almost no production [9] . In order to achieve more uniform multi-fracture propagation, it is necessary to study the stress shadow effect and fluid distribution within one stage during hydraulic fracturing.
Limited-entry method can force the fracturing fluid to be more evenly distributed among perforations by limiting the sizes and numbers of perforations in one stage and thus enhancing the bottom-hole pressure. Initially, limited-entry method was mainly applied in vertical wells to simultaneously create fractures in interbedded layers with heterogeneous breakdown pressures [10] . Recently, this method has been applied to multi-cluster fracturing in horizontal wells, where fewer perforations and higher perforation frictions have been generated to achieve the uniform growth of multiple fractures. Although field observations indicate the enhancement of hydrocarbon production through the limited-entry method (or "extreme limited-entry" method) [11] , it remains unknown that how perforations are limited or optimized in cases with different fracture designs.
At present, numerical studies have been done to understand the uneven propagation of multicluster fractures, whose methods are mainly the finite element method, the extended finite element method, the discrete element method, and the boundary element method. Because the discrete element method uses microscopic rock particles to simulate fracture initiation and propagation in the reservoir, it requires a high computational capability [12] . Instead of presetting the global grids, the boundary element method only needs to set grids for fracture boundaries, and solves the induced stress through the approximate solution of Green's function, which enhance the calculation speed and thus being widely used [13] . However, the boundary element method cannot deal with problems with reservoir heterogeneity, which limits its application in the field. The finite element method with the cohesive force model can solve the problem of multi-cluster fracture propagation in the heterogenous reservoir with a moderate computation cost, and thus also being widely used in hydraulic fracturing simulations [14][15][16] . Wang et al. extended the finite element method with the cohesive zone model (CZM), and conducted a comprehensive study that showed the impact of primary fractures on the secondary fractures in different scenarios; however, the dynamic distribution of fracturing fluid and proppants among multiple fractures were ignored [15] . Shin et al. realized the automatic fluid distribution during fracture propagation using the cohesion model with User Amplitude Curve Subroutine (UAMP); however, impacts of perforation frictions no fluid distribution was not studied for optimizing the limited-entry designs [14] .
In this paper, the extended finite element method with the cohesive model and subroutine UAMP is applied, and a two-dimensional finite element model of multi-fracture propagation is established and studied. This model considers the perforating friction and flow distribution among multi-fractures, where the subroutine UAMP is used to calculate the dynamic distribution of perforating friction and fluid, and the CZM model is used to simulate fracture initiation and propagation. The nonlinear finite element simulator ABAQUS is used to simulate the simultaneous propagation of fractures. After the accuracy of the model is verified by comparing with results of the numerical solution, cases with five clusters and limited-entry designs are simulated to study the impacts of perforation friction on fracture dimensions and fluid distributions.

Numerical Model
Multi-cluster fracturing in a horizontal well is a complex fluid-solid coupling problem. According to the sequence of fracturing fluid flow, it can be divided into the following three processes [17] , as shown in Fig. 1. Firstly, the fracturing fluid flows into the horizontal wellbore by the high-pressure pump, and enters each fracture through perforating holes in the horizontal casing. Secondly, the fracturing fluid enters the fracture through the perforation holes, and the high-pressure fluid initiates fractures in the rock, followed by the fracturing fluid flowing into these fractures. Thirdly, the fracturing fluid leaks off into the reservoir rock adjacent to hydraulic fractures, which meanwhile causes the deformation of reservoir rocks.

Functions of rock deformation and fluid flow in porous media
Rock is a porous medium containing the solid skeleton and pores filled by the fluid. Assuming that rock units are isotropic and linear elastic, the stress-strain relationship satisfies the linear elastic equation [18] as follows.
, where ij  is the effective stress of the rock unit,  and G are Lamé parameters, C and M are elastic moduli of two phases, vol  is the volumetric strain,  is the strain used to describe the volumetric deformation of the fluid relative to the rock skeleton, p is the pore pressure.
Meanwhile, fluid flow in the porous media satisfies the conservation of mass as follows.

Cohesive zone model (CZM)
CZM is used to simulate the deformation and separation of fracture tips. Cohesion elements have the cohesive force, whose relationship with the displacement follows the tractionseparation law (TSL) [18] . Therefore, the separation of cohesive layers can represent the initiation and growth of hydraulic fractures, which follow the damage initiation criterion and damage evolution law respectively [15] . The cohesion element of a fracture consists of cohesive broken zone and cohesive unbroken zone, which represent the completely opened fracture body and the opening fracture tip respectively, as shown in

Damage Initiation Criterion
Based on the maximum nominal stress criterion, the fracture initiates when the maximum nominal stress ratio reaches 1, as described by Equation 4 below.

Damage Evolution Law
In order to describe the damage evolution of the cohesive element under tensile or shear stresses, the effective displacement is introduced and defined as follows.
, where m  is the effective displacement, n  is the nominal normal strain, s  and t  are the nominal shear strains of the cohesive element.
A damage variable D is also defined to show the property change of the element. After the damage starts, unloading is always assumed to occur linearly along the origin of the traction separation surface, as shown in Fig. 3. Reloading after unloading follows the same linear path until the softening envelope (AB line) is reached. Once the softening envelope is reached, further reloads follow the envelope shown by arrows in Fig. 3. In ABAQUS, D is defined by Equation 7 as below.
, where max m  is the maximum effective displacement during the loading process.

Fluid Flow in Fractures
Fluid flow in a fracture includes the flow along the fracture face and the flow perpendicular to the fracture face. Fluid flow along the fracture unit satisfies Darcy's law and can be described , where t k is the fracture permeability.
Fluid flow perpendicular to the fracture face is mainly the leak-off of the fracturing fluid. By defining the leak-off coefficients of the porous media, the leak-off rates can be described by Equation 9 as below. Within the fracture unit, the conservation of mass is followed, as described by Equation 11 as below.
, where w is the fracture width, and q is the flow rate of fluid in the fracture.

Fluid Distribution and UAMP of Perforation Units
During the simultaneous propagation of multi-cluster fractures, fluid flow through perforation holes are not only the function of pressure in each fracture, but also that of the extra friction exerted by the peroration hole (i.e., perforation friction). ABAQUS cannot simulate the dynamic distribution of fracturing fluid from the wellbore through perforation holes and among multiple fractures. However, this can be achieved by using UAMP in ABAQUS; sensors are preset at the inlet of each fracture to record the pressure change and send to UAMP, which then calculates the perforation pressure and automatically distributes the fracturing fluid among perforation holes within one fracturing stage. In UAMP, fluid distribution among multiple fractures should satisfy the following relationships. ;  is the slurry density, i n is the number of perforation holes in this cluster, i D is the diameter of each perforation hole, and C is the discharge coefficient.
In UAMP, initial pressure at each perforation cluster (fracture) are firstly calculated from the initial conditions; then, the fracture propagation is calculated using the CZM, from which displacements of all fracture elements, the stress and pressure fields can be determined; next, perforation frictions and fluids distribution can be determined through UAMP, during which Newton-Raphson iteration is applied. This is repeated till the end of fracture propagation, and provides fracture dimensions and fluid distribution among multi-clusters with consideration of stress shadow, fluid-solid coupling effects, and potential changes of perforation frictions.

Validation of UAMP Results
In order to validate the calculation results of perforation frictions and fluid distributions among multiple clusters, two cases with simple setting are simulated and compared with the analytical solutions, where the total injection flow rate t Q is 2 m 3 /min, the perforation number i n is 2 per cluster, the perforation diameter i D is 8 mm, and the initial discharge coefficient i C is 0.56. In Case #1, fracture inlet pressures are both 6 MPa at two clusters, the perforation number is 10 at the first cluster, while that of the second cluster changes from 6 to 14. In Case #2, perforation numbers are both 20 at both clusters, the fracture inlet pressure is 6 MPa at the first cluster, while that of the second cluster changes from 2 MPa to 10 MPa.
As shown in Fig. 4 and  5); reducing the perforation number can inhibit the flow of fracturing fluid and raise the bottomhole pressure, which can make the multiple fractures uniformly propagated, and this is also the physics behind the limited-entry method. Good match between analytical results and numerical results as shown in Fig. 4 and Fig. 5 indicates that the subroutine UAMP enables the extended finite element method to solve the problem of multiple fracture propagation with considerations of perforation friction change and fracturing fluid redistribution.

Simulation Results of Simultaneous Propagation of Multiple Fractures Using UAMP
In order to study the impact of different LE designs on dimensions of hydraulic fractures and distribution of fracturing fluid within one fracturing stage, a homogeneous 2D geological model is built, whose dimension is 100 m by 100 m. 5 clusters is initialized with an initial fracture spacing of 10 m; fractures are named HF1 to HF5 from the left to right, in which the central fracture is HF3. As listed in Table 1, CLE and XLE are compared, where the former design has 10 perforations per cluster and the latter design only has 5 perforations per cluster.  Fig. 6 shows the propagation of 5 fractures and the change of the stress field with time when CLE design is applied. As can be seen from Fig. 6 (a), 5 clusters of fractures are relatively uniform at the initial stage of fracturing, and the length of the central fracture (HF3) is slightly smaller than that of the sub-central fractures and edge fractures. This is the evidence showing the influence of stress shadow effect exerted by the neighboring fractures that inhibits the growth of the central fracture. As the fracturing time increases, the central fracture and subcentral fracture are more obviously inhibited by edge fracture rocks. Fig. 6 (b) and FIG. 6(c) show that the central fracture stops growing due to the compression of neighboring fractures, and the fracture width of the central fracture is even narrower than that of the initial fracture due to the lack of proppants. Fig. 6 (d) shows that at the end of fracturing (T=100 s in this simulation case), the edge fractures (HF1 and HF5) are the longest, followed by sub-central fractures (HF2 and HF4), and the central fracture (HF3) does not grow.  . 7 -Fig. 9 shows the changes of fracture inlet pressure, fluid distribution ratio, and fracture inlet width of 5 fractures within one fracturing stage with time when the conventional limitedentry method (CLE) is applied, where the perforation number is 10 per cluster. As shown in Fig. 7, pressures in all 5 fractures are same at the early stage of fracturing; however, the pressure of the central fracture (HF3) increases fast due to the stress shadow effect exerted by the neighboring fractures, and becomes the largest pressure among all fractures within this stage; the lowest fracture inlet pressure is observed at the edge fractures, and this is because they endure the compressive stress only from one side of the fracture during fracture propagation as can be observed in Fig. 6 (d). The difference in pressure results in an opposite trend in the fluid distribution among 5 fracture, where edge fractures HF1 and HF5 receive 60% of the total volume of fracturing fluid, while the central fracture HF3 barely receives any fluid (Fig. 8). In this simulation case, due to the coupling effects between the stress shadow effect and fluid redistribution, the central fracture barely grows, and the edge fractures forms long but narrow fractures; meanwhile, growth of the sub-central fractures are hindered by the edge fractures, and the injected fracturing fluid forces the fractures to increase in width ( Fig.  9). Simulation results show the extremely uneven propagation of fractures through CLE even though one stage only has 5 fractures with a fracture spacing of 10 m.  Fig. 10 shows the propagation of 5 fractures and the change of the stress field with time when XLE design is applied. Compared to the fracture propagation with the CLE design, fracture lengths of 5 fractures are more uniform with the XLE design. As can be seen from Fig. 10 (b), all 5 fractures can initiate and have similar fracture lengths. Because of the increase of perforation frictions, wellbore pressure increases and suppresses the heterogeneity among fractures, thus resulting in similar fracturing fluid distribution ratios and fracture lengths as can be seen in Fig. 11 -Fig. 13. With the propagation of fractures, the central fracture gradually becomes faster than the sub-central fractures, but the stress shadow effect generated is not large enough to limit the propagation of edge fractures ( Fig. 10 (b) and Fig.  10 (c)). At the late time of hydraulic fracturing, edge fractures become the largest fractures with this stage and limit the propagation of sub-central and central fractures, as can be observed in Fig. 10 (d). This simulation result indicates that reducing the perforation number can successfully balance the fracturing fluid flow into the multiple fractures, and thus resulting in a relatively uniform fluid distribution and fracture propagation.  Fig. 13 further shows the changes of fracture inlet pressure, fluid distribution ratio, and fracture inlet width of 5 fractures within one fracturing stage with time when the extreme limited-entry method (XLE) is applied, where the perforation number is reduced by half comparing to that in CLE. This reduction increases the wellbore pressure and the fracture inlet pressure. Comparing to the case using CLE, the fracture inlet pressure increases by 9% and the difference among 5 fractures is significantly reduced (Fig. 11); meanwhile, fluid distribution ratios of 5 fractures become similar to each other, where the difference is less than 5% (Fig.  12). However, the fracture widths of sub-central fractures are still obviously larger than the other 3 fractures because of the hinderance of their fracture lengths, and this is likely due to the limitation of model dimensions where the boundary effect cannot be ignored. Nevertheless, simulation results can already indicate that using XLE and reducing the perforation number can create more uniform hydraulic fractures in a stage when other parameters remain the same.

Conclusions
Using the cohesive zone model and the user subroutine UAMP, multi-cluster fracture propagation is simulated in ABAQUS that could capture the effects of fluid-solid coupling, stress shadow, and perforation friction change on fracturing fluid redistribution and competitive growth of multiple fractures within one fracturing stage. After validating the simulation method, impacts of different perforation designs, CLE or XLE, are further studied and compared. Main findings are concluded as follows.
(1) Dynamic fluid distribution among fractures within one fracturing stage can be simulated. Agreement between numerical and analytical results of fracture inlet pressures and fluid (2) Reducing the perforation number can raise the wellbore pressure that reduces the impacts of stress shadow effect and reservoir heterogeneity, and thus make fractures to propagate more uniform. Even for a fracturing case with five fractures in one stage, reducing the perforation number from 10 to 5 can reduce the difference of fluid distribution ratio from 60% to less than 5% without obviously increasing the wellhead operation pressure (less than 9%).
(3) The extended finite element method with the cohesive zone model can capture the change of propagation directions of fractures under the change of stress field, as well as the nonuniform propagation of fractures within one fracturing stage. However, due to the limitation of two-dimensional model and model dimensions, simulation results can be influenced by the boundary effect, and it is suggested to use three-dimensional simulation or larger model to eliminate such effect.
(4) XLE can maintain a uniform fracturing fluid flow into the multiple fractures in a stage, especially during the early time of fracturing. However, as time goes on, heterogeneity among fractures are gradually increased by the stress shadow effect and change of perforation friction. This limits the enhancement of XLE on uniform propagation of fractures. To obtain more uniform fractures or less fracture spacing is required, temporary plugging technique is recommended to used with the XLE design.

Declaration of Competing Interest
The authors declare that they have no conflict of interests.