Numerical investigation of fluid flow in fracture junctions with consideration to effect of fracture morphology

Because the fracturing fluid flow at fracture junctions directly affects the proppant distribution in fracture networks, it is very important to quantitatively investigate the fluid flow in fracture junctions. In this study, the fluid flow in a laboratory‐size fracture junction was investigated, and equations for calculating the fluid split ratio were established. The correlations between the pressure loss coefficients and the width ratio, approaching angle, and fluid split ratio were obtained by simulation. Using a custom‐made Visual Basic program, the fluid split ratios in an actual‐size fracture or complex fracture networks with several secondary or tertiary fractures were calculated. The results reveal that the fluid split ratio is significantly affected by the approaching angle and width ratio in laboratory‐size fractures. In contrast, in an actual‐size fracture junction, the fluid split ratio is mainly affected by the length of the fracture channel and the width ratio. The split ratios in laboratory‐size and real‐size fracture networks with several secondary or tertiary fractures are quite different, which indicates that the split ratio should be quantitatively estimated when investigating the proppant distribution by experiment or simulation. The findings of this study can be useful in the quantitative investigation of the proppant distribution in fracture networks.


| INTRODUCTION
The proppant distribution in fracture networks directly affects the fracturing effect. 1 Owing to the vital influence of fluid flow on the proppant distribution in fracture networks, 2 it is very important to investigate the fluid flow in complex fracture networks.
Typically, experiments and numerical simulations are used to investigate the proppant transport in fracture networks. When a sand-carrying fluid is injected, a fracture channel has been formed by the initial fluid. However, with the continuous injection of fracturing fluid, the hydraulic fracture propagates dynamically, which increases the complexity of investigating proppant transport in fracture networks. Thus, a fixed fracture size (width and length) is typically used to simplify the problem in experimental research or numerical simulations. Sahai et al 3 investigated the proppant flow from a primary fracture to a 90° secondary fracture and reported that the proppant flowed into the secondary fracture by turning the corner while traveling with the fluid or by the gravity effect after settling in a primary fracture. Han et al 4 investigated the proppant distribution in T-junction and crossing-junction fracture geometries using computational fluid dynamics (CFD) models. They analyzed the effects of the fluid viscosity, proppant diameter, proppant density, and injection rate of the fluid on the proppant migration. Alotaibi and Miskimins 5 conducted similar studies to Han et al. 4 The approaching angle affects the fluid split ratio and the volume fraction of the proppant, which enters the secondary fracture. Hence, this effect must be considered when investigating the proppant distribution in fracture networks, instead of simplifying complex fracture networks into an orthogonal fracture structure. Tong and Mohanty 6 and Li et al 7 considered the effects of the approaching angle on the proppant distribution in fracture networks. Although these studies can help in understanding the phenomenon of proppant migration in complex fracture networks, they did not quantitatively estimate the proppant transport at the fracture junction. In recent years, various studies have focused on this problem. Filippov et al 8 investigated the particle-laden fluid flow close to an individual perforation and demonstrated that the efficiency of particle transport to the fracture depends on two dimensionless parameters: the particle stokes number and the ratio of flow velocity in the wellbore to the fracture. Wen et al 9 used the fluid split ratio to approximate the volume fraction of the proppant, which enters into secondary fractures when the proppant transport in fracture networks is considered. Chang et al 2 established a simple mathematical model to calculate the proppant transport at an intersection. However, their model does not consider the effects of the approaching angle. From the above-mentioned studies, it is understood that the volume fraction of the proppant, which enters into secondary fractures, is affected by the fluid split ratio. In other words, the fluid split ratio at the fracture junctions is a key parameter for assessing the migration form, migration distance of the proppant in branch fractures, and whether the proppant can be carried into branch fractures. Thus, it is important to clarify the mechanisms of fluid flow in fracture junctions to quantitatively calculate the volume fraction of the proppant entering secondary fractures.
Earlier studies on the fluid split ratio have mostly focused on the fluid flow in pipe junctions. [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] However, few studies have calculated the fracturing fluid flow at the fracture junction. Additionally, laboratory-size fracture geometry has always been considered in experimental or numerical simulation studies on proppant transport in fracture networks, without verifying whether the size of fractures affects the fluid flow or proppant migration in fracture junctions. Thus, the objective of this study was to analyze the fluid flow at the fracture junction using a CFD method and analyze the effect of the fracture size on the fluid split ratio. The findings of this study will provide the basis for calculating the proppant transport in complex fracture networks.

LABORATORY-SIZE FRACTURE JUNCTION
A laboratory-size fracture is always used in the experimental method or numerical simulation studies on proppant transport in complex fracture networks. Thus, in this section, the factors influencing the fluid split ratio in a laboratory-size fracture junction were investigated using a CFD method, and the calculation model of the fluid split ratio was established.

| CFD model
The basic equations of fluid flow in a fracture junction are expressed as follows:

Continuity equation of the fluid:
where is fluid density (kg/m 3 ), V is fluid velocity (m/s), and t is the time (s).

Momentum conservation equation:
where P * is the partial pressure (Pa), and g is the gravitational acceleration (m/s 2 ).

Turbulence model
The k-ɛ model is typically used to simulate a fully developed turbulent flow away from the wall, while the k-ω model is more widely used in boundary layer problems. In this study, a shear-stress transport k-ω model (k-ω SST), 25,26 which combines the advantages of both models, was used owing to its excellent suitability to simulations of fluid flow in junctions. 24 (1) LIU et aL.

Solution method
The simulations were performed using the double-precision solver of the ANSYS FLUENT 18.1 CFD software. A pressure-based solver was adopted, and a coupled solution algorithm was applied for pressure-velocity coupling. 27 The second-order upwind scheme was used to discretize the momentum.

Boundary conditions
The velocity inlet boundary condition was adopted in the simulations. To ensure that the injected fluid had a fully developed turbulent flow, we simulated the fluid flow in a long straight fracture (E) and used the simulation results of the velocity and turbulence distribution at the outlet (E_outlet) as the injection parameters. Figure 2 shows the geometry of E (L e = 60d = 0.48 m, where d is the hydraulic diameter, Figure 3A shows the velocity distribution at the exit of E when the Re of the injected fluid was 32 000. The simulation results for the velocity at the centerline of E_out along the width and height directions are in good agreement with the experimental data, 23 as shown in Figure 3B. The error close to the wall is larger than that in the middle of the fracture, and the maximum error is less than 10%. The fluid velocity at the centerline of the computational domain E begins to stabilize from X = 0.2 m along the Xdirection, which indicates that the fluid transitioned into a fully developed turbulence, as shown in Figure 4. Table 1 presents the comparison between the pressure loss from X = 0.38 m to X = 0.48 m along the computational domain E calculated by Equation (3) and the simulation results, which also demonstrate the high accuracy of the simulation.
where p′ is the frictional pressure drop along the channel (Pa), V is the fluid velocity (m/s), is Darcy's friction coefficient (dimensionless), and l is the length of the calculated region along the X-direction (m).
The pressure outlet boundary condition was applied when investigating the split ratio of the fluid flow in the fracture junctions, and the outflow boundary condition was used when analyzing the pressure loss coefficients of the junction. A no-slip wall boundary condition was set.

Meshing strategy
Hexahedron elements were used to generate a structured grid, as shown in Figure 5. A boundary layer mesh was built close to the walls owing to the high boundary requirement of the k-ω SST turbulence model. The height of the first cell close to the wall was 5 × 10 −5 m to maintain the wall's YPLUS value at approximately one. It is well known that more accurate simulation results can be obtained with mesh encryption. However, because computer resources are limited, the rounding error caused by the computer's floating-point operations increases as the grid increases. Therefore, blindly increasing the meshes does not increase the simulation accuracy. The results were checked for grid independence in the course of a 90° rectangular junction with a width ratio of one. The total E_out E_in pressure drop was used to estimate the discretization uncertainty, which can be calculated using Equations (5) and (6). Table 2 presents the simulation results for the various cell numbers. As can be seen, the relative error of ΔPt is small when the cell number increases from 1 724 334 to 2 636 249. Thus, the third meshing strategy was applied to simulate fluid flow in a 90° fracture junction.
where Pt 1 , Pt 2 , and Pt 3 are the total pressure of the inlet channel, straight outlet channel, and branch channel, respectively (Pa); ΔPt is the total pressure drop (Pa); is the relative error of ΔPt (dimensionless); i refers to the case number.

| Influence factors of fluid split ratio
The factors affecting the fluid split ratio were investigated using a CFD method, and the required simulation parameters are listed in Table 3. The width ratio and fluid split ratio are defined by Equations (8) and (9), respectively: where subscripts 1, 2, and 3 denote the inlet fracture channel, straight outlet fracture channel, and secondary fracture, respectively; w * is the width ratio, the ratio of the width of the branch outlet channel to the width of the inlet channel (dimensionless); q * denotes to fluid split ratio, that is, the ratio of the fluid volume entering the secondary fracture to the volume of the injected fluid (dimensionless). 1. Effects of approaching angle on split ratio Figure 6 shows the impact of the approaching angle on the fluid split ratio. As can be seen, when the lateral branch angle is less than 90°, the split ratio decreases as the lateral branch angle increases; when the lateral branch angle is larger than 90°, the split ratio slightly increases with the branch angle. Figure 7 shows the impact of Re on the fluid split ratio with various lateral branch angles. As can be seen, when Re is less than 10 000, the split ratio decreases as Re increases, and is independent of Re when Re > 10 000.

Effects of Re on fluid split ratio
3. Effects of width ratio on fluid split ratio Figure 8 shows the influence of the width ratio on the split ratio. As can be seen, the split ratio increases with the width ratio, and the width ratio has a stronger influence on the split ratio when the lateral branch angle is small.  (10) and (11), respectively, according to the computational domain shown in Figure 1.  where P is the static pressure (Pa), V is the fluid velocity (m/s), L is the fracture length (m), K is the pressure loss coefficient (dimensionless), and is the energy shape factor, which is typically considered to be 1 in turbulent flow calculations (dimensionless). As expressed by Equations (10) and (11), when the size of the fracture geometry is small, the pressure loss coefficients are the key parameters affecting the fluid split ratio. To calculate the fluid split ratio in fracture junctions, it is important to construct formulas for the pressure loss coefficients influenced by various factors. In this study, the pressure loss coefficients were investigated based on the simulation results.

| Calculation model of fluid split ratio
(10)  Figure 9 shows the effects of different width ratios on the pressure loss coefficient K 12 . As can be seen, K 12 increases as the increases of the split ratio and is independent of the width ratio. Figure 10 shows the variation of K 12 with different approaching angles and a width ratio of one. As shown, K 12 is independent of approaching angle and is only affected by the split ratio. Thus, the pressure loss coefficient K 12 can be calculated according to the fluid split ratio, as expressed by Equation (12). Figures 11 and 12 show the relationship between K 13 and the approaching angle, when w * = 1. As can be seen, when the lateral branch angle is less than 90°, K 13 decreases as the fluid split ratio increases; when the lateral branch angle is more than 90°, K 13 increases with the split ratio. Additionally, K 13 increases with the lateral branch angle, and the linear relationship between K 13 and cos θ is shown in Figure 12. Figures 13 and 14 show the effects of the lateral branch angle and width ratio on K 13 . As can be seen, when the width ratio is smaller than 0.5, K 13 increases rapidly as the split ratio increases, whereas the effect of the width ratio on K 13 is small when the width ratio is more than 0.5. Additionally, K 13 increases with q * /w * and has a quadratic relationship with q * /w * when the lateral branch angle is 90°.

Pressure loss coefficient K 13
The simulation results reveal that K 13 is affected by q * / w * and cos , as shown in Figure 15. The relationship between q * / w * , cos , and K 13 can be fitted as expressed by Equation (13) using the CFTOOL toolbox of the MATLAB software. Figures 16 and 17 show the comparison between the results obtained for K 12 and K 13 using Equations (12) and (13) and the experimental data for the tube junction obtained from the literature, 11,13 respectively. As can be seen, the model is consistent with the experimental results. Figures 18 and 19 show the effect of the fracture scale on the fluid split ratio. As can be seen, the fluid split ratio increases with the length of the outlet channels. As the length of the outlet fracture channel increases, the fluid split ratio gradually approaches 0.5. When the length of the straight outlet

WITH SEVERAL JUNCTIONS
The custom-made Visual Basic program for calculating the fluid split ratio in fractures with more secondary or tertiary F I G U R E 1 5 Fit of q * / w * and cos with fractures is shown in Figure 22. The split ratios of laboratorysize and actual-size fractures with several secondary or tertiary fractures were analyzed. Figure 23 shows the geometry of fractures with more secondary fractures. Figure 24 shows the comparison between the split ratio of actual-size fractures and that of laboratory-size fractures with three secondary fractures. The size parameters of these fractures are listed in Table 4. As can be seen, when the fracture size is small, the approaching angle affects the fluid split ratio, and the fluid flowing across the outlet of the primary fracture accounts for the majority. In contrast, the split ratio of the secondary fracture in front of the primary fracture is more significant when the fracture size is large. Figure 25 shows the geometry of fractures with tertiary fractures. Figure 26 shows the comparison between the split ratio in laboratory-size fractures and that in actual-size fractures. The fracture contains three secondary fractures and two tertiary fractures in the first and third secondary fractures; the fracture parameters are listed in Table 5. As can be seen, there is a significant difference between the split ratio of the laboratory-size fractures and that of actual-size fractures. The fluid flowing across the secondary fractures and tertiary fractures located at the front of the primary fractures accounts for the majority when the fracture size is large, while the fluid flowing across the primary fracture and secondary fracture located at the front of the primary fracture accounts for the majority when the fracture size is small.

| DISCUSSION
From the above-mentioned investigations, it is understood that the fluid flow is affected by the approaching angle and width ratio in a small-scale fracture junction, but is only slightly affected by the approaching angle in a large-scale fracture junction. A study by Tong and Mohanty 6 demonstrated that, in a small-scale fracture junction, the volume  29,30 were not considered. Additionally, owing to the influence of the reservoir properties and in situ stress, the fluid flow during an actual slickwater hydraulic fracturing may be more complex. To better represent the proppant distribution in an actual fracture in the experimental or numerical simulation results obtained for proppant transport, it is worth discussing whether the outlet boundary condition of constant fluid flow or constant pressure should be adopted.