Abstract
Intensive stage multicluster fracturing is the main technology of deep shale gas development, but this technology is easy to lead to the uneven propagation of multiple fractures. In order to deeply analyze the influence mechanism and law of proppant flow in the wellbore on multiple fractures and propagation, Firstly, based on the theoretical knowledge of computational fluid dynamics (CFD), a numerical model of solid-liquid two-phase flow of proppant migration in horizontal wellbore is established, which is calculated and solved by using the Euler-Euler multiphase-flow mixture model, and the effects of different influencing factors on the migration and distribution of proppant in wellbore are studied. Then, based on the fluid structure coupling theory, a planar three-dimensional (PL3D) numerical model of multicluster fracturing in horizontal wells is established. Taking the numerical simulation results of wellbore proppant migration as the initial condition, the influence law of proppant settlement on the fluid volume distribution among multiple clusters and the difference of multifracture geometry is clarified. The results show that the proppant content in the horizontal wellbore increases from heel end to toe end. The fracturing fluid viscosity and sand ratio are the main influencing parameters. The proppant settlement in the wellbore is also an important factor affecting the liquid and sand injection of multifracture. When considering the proppant settlement in the wellbore, the phenomenon of multifracture uneven propagation is more significant. The construction scheme of medium viscosity, large displacement fracturing fluid, large proppant size, and medium sand ratio in a certain range is conducive to promote the even propagation of multifracture. This study is helpful to guide the construction design of multicluster fracturing in horizontal wells, so as to improve the success rate of multicluster fracturing and greatly save the construction cost.
1. Introduction
So far, segmented multicluster perforation of horizontal wells is the main technical means of deep shale gas development. Multicluster fracturing technology in the section improves the construction efficiency, but according to near well temperature or acoustic monitoring on the filed, it is found that simultaneous initiation and even propagation of multifracture are not always assured, there are a large number of invalid perforations during construction, and there are obvious differences between fluid and sand volume of perforation at different positions [1–3].
Many scholars have done a lot of research on the migration law of proppant in the wellbore and the dynamic coupling between proppant and fracture propagation, whose research methods are mainly carried out from two aspects: experimental simulation and numerical simulation [4–8]. Wen et al. systematically studied the settlement and migration law of proppant in wellbore and fracture during shale gas slick water fracturing by using self-designed experimental instruments and established a critical velocity calculation model based on experimental data through multivariate nonlinear fitting. Their model considered the influencing factors such as proppant, fracturing fluid, horizontal well diameter, and sand ratio. It can accurately predict the critical sedimentation velocity and critical resuspension velocity [9, 10]. Kamga et al. studied the settlement law of proppant in horizontal wellbore during slick water fracturing based on experimental means. According to the experimental results, they found that fluid velocity and proppant properties (type and size) are the main factors affecting the uneven distribution of proppant in each perforation diameter [11, 12]. Bokane et al. studied the distribution of proppant in each cluster of fractures considering factors such as proppant particle size, density, fracturing fluid viscosity, and pumping displacement through large-scale laboratory experiments, and they found that high displacement, high viscosity and small particle size, and low-density proppant can improve the nonuniform sand injection effect of each cluster [13]. Yi et al. used the Eulerian multiphase flow model to fit Freddy Crespo’s laboratory test results, and verified the feasibility of Eulerian multiphase flow in simulating proppant migration in a wellbore [14]. Sophie et al. studied the distribution law of proppant among multiple perforation clusters during temporary plugging perforation based on the CFD/DEM method. The results show that even if the fluid is evenly distributed among all perforation clusters, the distribution of proppant may be very uneven. Proppant tends to accumulate near the toe perforation, and the proppant concentration increases with the outflow of fluid. According to the numerical simulation results, a method to quantify the proppant transportation effect on the wellbore scale is proposed. Hou et al. based on establishing the productivity prediction model of a shale gas well under the nonuniform distribution of proppant, the effects of the parameters of reservoir and fracturing fracture on the productivity of shale gas well are analyzed under the uniform distribution and the nonuniform distribution of proppant. The results show that the nonuniform distribution of proppant among different fracturing clusters makes the reduction of productivity more obvious [15]. In the numerical simulation of proppant migration and distribution in multiple fractures, predecessors have also studied in detail. Dontsov adopts the solid-liquid two-phase hydrodynamics method and introduces the constitutive equation of sand carrying fluid interaction established by Boyer to obtain a strict mathematical description of proppant transport. Dontsov’s model can describe the transition process of proppant sand carrying fluid from Poiseuille flow to Darcy flow. The model has been applied in KGD and P3D models [16, 17].
However, the current research on the law of proppant settlement in the wellbore and its influence on multifracture propagation is not clear, which needs further research and more reliable theoretical support. In view of this, by establishing a numerical model of solid-liquid two-phase flow of proppant migration in horizontal wellbore and a planar three-dimensional (PL3D) multifracture numerical model, and using the means of numerical simulation, this paper analyzes and studies the wellbore proppant migration law of multistage fracturing in deep shale horizontal wells and the influence of wellbore proppant accumulation on multifracture propagation.
2. The Establishment and Solution of Model
2.1. Wellbore Proppant Transport Model
The benchmark geometric models of multicluster perforation in the horizontal section and fracturing section of three-dimensional horizontal well are established by using ICEM CFD module of ANSYS software. As shown in Figures 1 and 2, the horizontal section of the geometric model is 1500 m in length, and the wellbore size is 5.5 inches. The length of fracturing section is 50 m, the number of perforations in a single section is 5 clusters per section, the cluster spacing is 10 m, the perforation density is 10 per meter, the perforation phase angle is 90°, the perforation diameter is 1 cm, and the perforation depth is 12 cm. The geometric model is meshed with a hexahedral structured mesh, and the mesh is refined for the flow domain where the pressure and flow rate change greatly, that is, the intersection of the wellbore and perforation and the wall of the entrance and exit.


After the grid model is established, it is imported into the FLUENT software for initial settings. The boundary conditions of the model are constant velocity inlet and constant pressure outlet, and the well wall is reflective wall. Taking the flow space of horizontal wellbore and perforation as the three-dimensional calculation domain of solid-liquid flow, a three-dimensional implicit transient solver based on pressure is selected, and the Eulerian-Eulerian multiphase flow mixture model is used for calculation and solution. The solution equations include the following: hydrodynamic equations (mass conservation equation, momentum conservation equation, and turbulence equation) and the force equation of proppant particles. The transport equations of solid mass, momentum, and particle temperature in particle dynamics are introduced to consider the interaction of particle collision and friction. Adopt the ~ ε turbulence model that characterizes the turbulence effect at low viscosity and high velocity. The finite volume method and discrete element coupling method are used to discretize the equation, and the coupling of pressure and velocity is solved by using the coupled simple algorithm. The Green-Gauss node method is used for the spatial discretization of the gradient, and the convective phase secondary upwind difference method is used to discretize other parameters. Meanwhile, in order to prevent the divergence of the equation, the simulation step size and underrelaxation factor are adjusted to make the equation converge. In this paper, the convergence residual is uniformly set to , and the simulation time is 10 min [14].
2.2. Numerical Model of Multifracture Propagation
The plane three-dimensional fracture propagation numerical model (PL3D) of multicluster fracturing in horizontal wells is established by Chen Ming [18–21]. The model is established based on the boundary element method, considering the effects of horizontal wellbore flow (laminar flow and turbulence), multicluster perforation friction, rock deformation (stress interference), fracture tip fracture, intrafracture flow, fluid leak off, and other factors. In order to study the multifracture morphology and propagation law when wellbore proppant settlement is considered, the following changes are made on the basis of the PL3D model and Dontsov model [16–21].
After proppant is taken into account, the fluid flow in the fracture is sand-carrying fluid. Proppant is regarded as a component of sand-carrying fluid; so, the mass conservation equation of sand-carrying fluid flow and proppant migration in fracture is
where is the fracture width, m; is the injection time, s; and are the flow in the fracture of sand-carrying fluid and proppant, respectively, m2/s; is the fluid leak off velocity, m/s; is the injection flow of cluster fracture, m3/s; and is the normalized proppant concentration at the entrance of cluster fracture, which is obtained from the numerical simulation of proppant migration in horizontal wellbore, dimensionless. x is the distance between the space point (a function related to the position coordinates x, y, z) and the origin, m; is the coordinate of of cluster fractures injection point [16, 17].
The calculation formula of flow in the fracture of sand carrying fluid and proppant is
where is the radius of proppant, m; is the acceleration of gravity, m/s2; are the density of proppant and fracturing fluid, kg/m3; is the unit vector in the -axis direction; is a dimensionless function describing the flow of sand carrying fluid; and are dimensionless functions describing proppant convection and settlement, respectively; and is the blocking function [16, 17].
At the same time, in order to quantitatively evaluate fluid volume distribution among multiclusters and the difference of multifracture morphology, the liquid inlet coefficient of variation was introduced based on the concept of distribution standard deviation. The larger the coefficient of variation, the greater the influence of this factor on the results, and the fracture propagation is more uneven [18–21]. where CV is the inlet liquid coefficient of variation, dimensionless; is the number of fractures, is the liquid volume of cluster n fracture, m3; and is the total liquid volume, m3.
When solving the model, the rock deformation is solved by the boundary element method, the fluid flow is solved by the finite volume method, and the proppant transport in the fracture is solved by the WENO (Weighted Essentially Non-Oscillatory) finite difference method. For the calculation equations involved, the model is solved by an explicit large step algorithm based on Runge-Kutta-Legendre (RKL), and the calculation efficiency and accuracy are greatly improved. It has good applicability for this study [18–21].
3. Simulation Scheme
Based on the construction parameters of intensive stage multicluster fracturing of typical deep shale gas horizontal well in a certain block in China, the numerical simulation research is carried out. The simulation scheme is divided into two steps: (1)Firstly, the wellbore proppant migration and settlement law under different construction conditions are analyzed by using the wellbore proppant migration model(2)According to the initial proppant concentration at the inlet of each shower hole obtained in step (1), in this study, we defined the initial proppant concentration as the volume fraction, that is, the total amount of proppant contained in a unit volume of fracturing. When analyzing the results, for the convenience of statistical analysis, the proppant concentration of each cluster is actually owned by the cluster. The sum of the perforation proppant concentrations is as follows. The numerical simulation of multifracture propagation considering wellbore proppant settlement and sand fracturing is carried out based on the PL3D model
The basic parameter settings of simulation are shown in Table 1.
4. Proppant Settlement Law at Different Positions of Wellbore
The simulation results of proppant settlement law at different positions of horizontal section wellbore are shown in Figure 3. Because proppant particles are affected by gravity, they mainly settle at the bottom of wellbore during migration, and the sand setting effect at the top and middle of wellbore is weak. As shown in Figure 4, the proppant volume fraction (concentration) at the bottom of the wellbore in the middle fracturing section increases from 0.1 to 0.6 along the heel end to toe end, and the proppant volume fraction at the middle and bottom of the wellbore remains below 0.1. There is no obvious change law; so, the follow-up study mainly focuses on the sand deposition law at the bottom of horizontal wellbore. Taking the bottom of the wellbore as the benchmark, the settlement law of proppant in the horizontal direction of the wellbore is analyzed. The results show that the proppant concentration gradually increases from heel end fracturing to toe end fracturing. With the increase of migration distance, the settlement rate gradually increases, and the settlement effect is more obvious in the middle and far end of the wellbore. As shown in Figure 4, the volume fraction of proppant at toe increases from 0.1 to 0.6 along the length of fracturing section, while the volume fraction of proppant at heel end and middle fracturing section changes only in the range of 0.1 to 0.2.


5. Effect of Wellbore Proppant Settlement on Multifracture Propagation
5.1. Comparative Analysis of Basic Cases
According to the simulation results of the benchmark case of wellbore proppant migration with horizontal section length of 1500 m, fracturing section length of 50 m, single-stage perforation of 5 clusters/section, fracturing fluid viscosity of 12 mPa·s, displacement of 12 m3 min-1, proppant particle size of 40/70 mesh, sand ratio of 10%, proppant density of 2310 kg m-3, using the PL3D model, the initial proppant concentration at the entrance of each cluster perforation (from the heel to the toe) is set to 0.422, 0.565, 0.641, 0.749, and 0.803 to simulate the migration and settlement of proppant in the horizontal wellbore. At the same time, in the control group, the average proppant concentration of 5 perforation clusters, 0.636, was taken as the uniform concentration at the entrance of the perforation cluster. The control group did not consider the influence of proppant migration and settlement in the horizontal wellbore. The fracture propagation morphology is shown in Figure 5. It can be seen from Figure 3 that due to the nonuniform settlement of proppant in the wellbore, the proppant concentration and the balance height of the sand bank were significantly different at the entrance of each cluster of perforations. At the same time, the existence of sand bank at the entrance of perforations prevents the fracturing fluid from entering the fracture from the perforation. The higher the balance height of the sand bank, the stronger the blocking effect on the fracturing fluid, the less fluid volume entering the fracture, and the more difficult the fracture is to propagate. Therefore, when considering the nonuniform settlement of proppant in the wellbore, the phenomenon of multifracture uneven propagation is more significant, and the restraint of intermediate cluster fractures is more stronger. As shown in Figure 6, under different conditions, each perforation cluster shows the advantage of heel end liquid inlet, and the fracture length of each cluster is positively correlated with the liquid inlet; that is, the greater the liquid inlet, the greater the fracture length. When considering the wellbore proppant settlement in the wellbore, the liquid inlet distribution proportion of each cluster of fractures is more uneven, and the liquid inlet difference coefficient is larger. When considering the wellbore settlement and not considering the wellbore settlement, the liquid inlet coefficient of variation of each cluster of fractures is 19.2% and 4.7%, respectively.

(a)

(b)

5.2. Influence of Different Construction Parameters
According to the simulation scheme in Section 2 above, with the help of the wellbore proppant migration model and PL3D model, the settlement law of proppant in wellbore under different construction conditions and its influence on multifracture propagation is analyzed.
5.2.1. Effect of Fracturing Section Length and Perforation Cluster Number
The length of the fracturing section and the number of perforation clusters (cluster spacing) are important parameters for segmented multicluster fracturing of horizontal wells. Different parameters will not only affect the construction cost and cycle but also affect the multifracture propagation due to the stress interference between perforation clusters. By modifying the numerical model, the migration law of proppant in wellbore under different section lengths and cluster numbers (cluster spacing) is studied. The results are as follows:
With a certain segment length, the proppant concentration increases from heel to toe. When the segment length is 40 m, the proppant concentration near the perforation cluster 5 at the heel end is about 0.35, while the volume content of perforation cluster 1 near the toe end increases to 0.47, the segment length increases, the proppant concentration at the perforation cluster at the same position increases relatively, and the proppant is more evenly distributed in each cluster (Figure 7(a)). When the section length increases from 40 m to 70 m, the advantage of liquid inlet of heel perforation cluster is gradually weakened, the liquid inlet coefficient of variation is reduced from 38% to 23% (Figure 7(b)), and the fracture propagation of each cluster is more even (Figure 8).

(a)

(b)

(a)

(b)

(c)

(d)
When the fracturing interval is a certain length, increasing the number of perforation clusters (reducing cluster spacing) will increase the nonuniform distribution of proppant in each cluster perforation, and the proppant settlement advantage is more obvious near the toe end perforation cluster. The proppant volume fraction near toe perforation cluster 3 is 0.57 in single segment 3 clusters, while the proppant content of toe end perforation cluster 8 increases to 1.22 in single segment 8 clusters. When there are 5 clusters in a single section, the content difference of perforating proppant in each cluster is the smallest, and the concentration difference is less than 0.02 (Figure 9(a)). With the increase of the number of perforation clusters, the advantage of liquid inlet at the heel end is more obvious, the inhibition of liquid inlet in the middle cluster is strengthened, the multifracture initiation is more difficult, and the phenomenon of multifracture uneven propagation is more serious (Figure 10). The liquid inlet difference of each cluster increases with the increase of the number of perforation clusters, the number of perforation clusters increases from 3 to 8, and the liquid inlet coefficient of variation increases from 3% to 55% (Figure 9(b)).

(a)

(b)

(a)

(b)

(c)

(d)
5.2.2. Effect of Fracturing Fluid Viscosity and Displacement
Because the fracturing fluid has a certain viscosity, the proppant will receive a force perpendicular to the movement direction during its movement. This force will make the proppant gather in the center of the flow line, so that the proppant can maintain the original movement direction, and it is difficult to change the movement direction into the perforation near the well wall. By simulating the concentration of proppant at each shower hole under different viscosity, it is found that with the increase of viscosity, the carryover of fracturing fluid to proppant particles is enhanced, and the proppant is more evenly distributed in different perforations. Under the condition of low viscosity 5 mPa·s, the proppant tends to be distributed to the toe end perforation cluster, the proppant concentration of heel end perforation cluster 5 is only 0.32, and the volume fraction of perforation cluster 1 near the toe end is 0.63. When the high viscosity is 100 mPa·s, the concentration of each cluster of proppants is in the range of 0.4 ~ 0.7, and the distribution is the most uniform (Figure 11(a)). The liquid inlet difference of each cluster increases with the increase of viscosity, and the liquid inlet coefficient of variation is 13% at 5 mPa·s and 39% at 100 mPa·s (Figure 11(b)). Fracture propagation is more uneven at high viscosity (Figure 12).

(a)

(b)

(a)

(b)

(c)

(d)
The larger the displacement, the greater the initial velocity of proppant, and the longer the equilibrium time required for proppant settlement under the same conditions, which is more conducive to the uniform distribution of proppant among perforation clusters. Under the condition of low displacement of 10 m3 min-1, the migration speed of proppant in the wellbore is small, and the settlement speed is fast. The volume fraction of perforation cluster 5 at the heel end is 0.35, and that of perforation cluster 1 at the toe end is 0.55. When the displacement increases to 16 m3 min-1, the volume fraction of each cluster of perforating proppant is in the range of 0.5 ~ 0.6, and the distribution is the most uniform (Figure 13(a)). The liquid inlet difference increases with the increase of displacement (Figure 13(b)). Under the condition of high displacement, the liquid inlet advantage at the heel end is more significant, and the inhibition of the middle cluster is stronger (Figure 14).

(a)

(b)

(a)

(b)

(c)

(d)
5.2.3. Effect of Proppant Particle Size and Initial Sand Ratio
The larger the proppant particle size, the greater the momentum and inertia it has, the less likely it is to deposit in the wellbore, and the less likely it is to change direction and enter the perforation. With the increase of particle size, the toe end settlement advantage of proppant is strengthened, and the nonuniformity of perforation distribution of each cluster is enhanced. At 20/40 mesh (500 μm), the concentration of proppant in heel end perforation cluster 5 is only 0.46, while the concentration of toe end perforation cluster 5 is 0.78.Under the condition of small particle size (100 mesh), the proppant is most evenly distributed among perforation clusters (Figure 15(a)). The liquid inlet coefficient of variation of each cluster increases with the increase of proppant particle size (Figure 15(b)). The smaller the proppant particle size, the more favorable for the uneven propagation of multiple fractures (Figure 16).

(a)

(b)

(a)

(b)

(c)

(d)
The larger the sand ratio, the stronger the interaction of solid-liquid two-phase molecules in unit volume fracturing fluid, the weaker the sand carrying capacity of fracturing fluid, and the shorter the effective migration distance of proppant particles, which is not conducive to the distant migration of proppant. When the sand ratio increases, the proppant gradually accumulates towards the toe end, and the nonuniformity of perforation distribution of each cluster increases. When the sand ratio is 30%, the volume fraction of proppant in heel end perforation cluster 5 is 0.46, and the volume fraction of toe end perforation cluster 1 reaches 1.22. Under the condition of low sand ratio (5%), the proppant is basically evenly distributed among each perforation cluster. At this time, the volume fraction of proppant of each cluster is in the range of 0.3 ~ 0.6. The concentration difference was the smallest and remained below 0.1 (Figure 17(a)). When the sand ratio is too large, it is easy to block the perforation at the toe end, so that the liquid inflow at this location is weakly close to zero, and the difference between the inflow of other cluster perforations is too large. The liquid inlet coefficient of variation reaches the minimum value of 13% when the sand ratio is 20% (Figure 17(b)). Increasing the sand ratio under certain conditions can reduce the liquid inlet difference coefficient of each cluster of fractures and promote the multiple fractures that propagate evenly (Figure 18).

(a)

(b)

(a)

(b)

(c)

(d)
6. Summary and Conclusions
Aiming at the problem of multicluster sand fracturing and multifracture propagation in deep shale horizontal wells, the numerical simulation research is carried out, and the following conclusions are obtained: (1)Proppant is mainly accumulated at the bottom of the wellbore, the settlement concentration increases from heel end to toe end in the horizontal fracturing section, and the settlement effect is more significant at the middle and far end of the wellbore(2)Fracturing fluid viscosity and proppant sand ratio are the main influencing parameters of proppant nonuniform settlement in wellbore. At the same time, under the same conditions, with the decrease of fracturing fluid viscosity, the increase of displacement, the increase of proppant particle size and sand ratio, and the concentration distribution of proppant in each cluster of perforation become more uneven(3)The proppant settlement in the wellbore is also an important factor affecting the fluid inflow of each perforation cluster and the propagation of multiple fractures(4)In the process of multifracture propagation, the heel end fluid inflow is dominant. When considering the nonuniform settlement of proppant in the wellbore, the fluid inflow difference of each cluster of multifracture is greater, and the phenomenon of uneven propagation is more significant(5)Considering the influence of nonuniform settlement and settlement of wellbore proppant, the construction scheme of selecting medium viscosity, large displacement fracturing fluid, large particle size, and medium sand ratio proppant in a certain range is conducive to promoting the balanced expansion of multiple fractures
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.