Abstract

Bed load proppant transport is a significant phenomenon during slickwater hydraulic fracturing. However, the mechanism of bed load proppant transport is still unclear. In the present study, the proppant transport process during slickwater hydraulic fracturing was simulated with a coupled computational fluid dynamics- (CFD-) discrete element method (DEM) model, and the mechanism of the bed load proppant transport was analyzed. A model for calculating the mass flux of the bed load layer was proposed and verified with experimental results from the literature. The results show that bed load migration is an essential mechanism of proppant transport. When the shear force of fluid acting on the surface of the sand bank reaches the critical Shields number, the proppant in the upper layer of the sand bank begins to migrate in the form of bed load. The movement of the bed load layer increases the time for the sand bank to reach the equilibrium height. In addition, the mass flux of the bed load layer significantly affects the equilibrium height of the sand bank. The mass flux of the bed load layer decreases, and the equilibrium height increases as the proppant density, proppant diameter, or rolling friction coefficient and static friction coefficient of the proppant increase, but the mass flux of the bed load layer increases, and the equilibrium height decreases as the fluid viscosity increases.

1. Introduction

Slickwater hydraulic fracturing is a crucial technique for unconventional reservoir development. Due to the low viscosity of slickwater, the proppant settles out of suspension quickly to form a sand bank, and then bed load proppant transport becomes a significant mechanism of proppant transport [1]. Therefore, the study of bed load proppant transport is of great significance to analyze the proppant distribution in the fractures.

Experimental and numerical simulations are the key ways to study proppant distribution in fractures. Kern et al. [2] studied the proppant transport in fractures through experiments. They proposed that only when the sand dune reaches the equilibrium height can the proppant be further transported to the fracture in the form of fluidization. Woodworth and Miskimins [3] carried out the erosion slot flow experiments and performed a sensitivity analysis to show the sensitivities of various parameters, such as proppant density and fluid viscosity. Mack et al. [1] indicated that the static friction coefficient and coefficient of restitution of proppant influence the bed load proppant transport and further influence the proppant distribution in fractures. The advanced ceramic proppant can be transported more in-depth into the fracture due to its higher coefficient of restitution and lower friction coefficient compared to sand and ceramic proppant. Tsai et al. [4] analyzed the proppant transport in fractures with the Eulerian-Lagrangian model. They found that the proppant settling behavior is influenced by proppant density, diameter, and injection rate. Blyton et al. [5] simulated the sand-carrying fluid flow in fractures with the CFD-DEM model. The simulation results indicated that the relative velocity of proppant and fluid is affected by sand ratio and the ratio of proppant diameter to fracture width. Zhang et al. [6] analyzed the impacts of various factors such as injection velocity and proppant size and density on proppant transport by simulation with the CFD-DEM model. Baldini et al. [7] presented similar to those of Zhang et al. [6]. Li et al. [8] studied the influence of interaction parameters between proppants and walls on proppant transport with a CFD-DEM model. They found that the friction coefficient of the proppant affects the bed load proppant transport. A lower friction coefficient is beneficial for bed load proppant transport. In addition, many scholars have carried out experiments or numerical simulations of proppant transport in complex fracture networks. Sahai et al. [9] experimentally studied the proppant distribution in fracture networks. Li et al. [10] researched the laws of proppant distribution in complex fracture networks through many experiments and the Eulerian-Eulerian simulation method. They proposed that using a large injection rate, low particle diameter, and sand ratio can improve the effective proppant placement in fracture networks. Kong and McAndrew [11] studied the effects of fracture width, fluid viscosity, sand load, particle size, and particle density on proppant distribution in fracture networks with the Eulerian-Eulerian model. These studies mainly focus on the controversial issue of whether the proppant can enter secondary fractures.

From the literature review, one can observe that the proppant transport in the fractures has two mechanisms: suspension and bed load. Most research is about the influence of various factors on proppant distribution in fractures, while research on bed load proppant transport is limited. In this paper, we intend to conduct quantitative research on the transport of bed load proppants in fractures. This study is crucial for clarifying the proppant distribution in fractures.

2. Numerical Simulation of Bed Load Proppant Transport in Fractures

The Eulerian-Eulerian [1114], Eulerian-Lagrangian [4, 15, 16], and coupled CFD-DEM [57, 1720] models are commonly used for simulating proppant distribution in fractures. The Eulerian-Eulerian model treats particles as another fluid phase, and the interaction between the fluid and particles is calculated using the momentum exchange technique. Compared to the Eulerian-Eulerian model, the Eulerian-Lagrangian and CFD-DEM models can calculate the trajectory of each particle. However, the Eulerian-Lagrangian model is not suitable for simulating particle accumulation. Therefore, the coupled CFD-DEM model was chosen to simulate the bed load proppant transport during slickwater hydraulic fracturing, which is helpful for us to observe and analyze the bed load proppant transport.

In addition, during slickwater hydraulic fracturing, the fracturing fluid leak off will lead to the reduction of the fluid velocity. However, as the fracture model used in the simulation is small, the influence of the fracture fluid leak off on fluid velocity is little. Thus, the fracturing fluid leakage was ignored in the simulations.

2.1. Models

2.1.1. Particle Control Equations

By Newton’s second law, the particle movement is predominantly translational and rotational, and the expressions are as follows [8, 21]:
(1)midvidt=mig+i=1nFci+Fdi+Fb,Iidwidt=i=1nTij,
where vi refers to the translational velocity of particle i (m/s); wi refers to the angular velocity of particle i (rad/s); mi refers to the mass of particle i (kg); Ii refers to the moment of inertia of particle i (kg m2); Fci refers to the contact force acting on particle i (N); Tij refers to the torque acting on particle i by particle j or wall (N m); n refers to the number of total contacts for particle i (dless); Fdi refers to the drag force of the fluid acting on the particle (N); and Fb refers to the buoyancy of the particle (N).
In the simulations, the particle-moving pattern was calculated with the Hertz-Mindlin (no slip) model [22, 23].
(2)Fci=Fcni+Fcti+Fdni+Fdti,
where Fcni and Fcti are the contact forces in the normal and tangential directions (N), respectively, and Fdni and Fdti are the damping forces in the normal and tangential directions, respectively (N).
(3)Fcni=43ER1/2α3/2,1E=1v12E1+1v22E2,R=1R1+1R2,
where E refers to the effective Young’s modulus (N/m2); R refers to the effective particle radius (m); α refers to the displacement of the amount of normal overlap (m); E refers to Young’s modulus (N/m2); v refers to the Poisson ratio (dless); and R refers to the particle radius (m).
(4)Fdni=256λSnmvnrel,m=m1m2m1+m2,λ=lneln2e+π2,Sn=2ERα,α=R1+R2r1r2,vnrel=v1v2n,n=r1r2r1r2,
where m refers to the equivalent mass (kg); Sn refers to the normal stiffness (N/m); vnrel refers to the normal relative velocity component (m/s); m refers to the particle mass (kg); e refers to the coefficient of restitution (dless); v refers to the particle velocity before collision (m/s); n refers to the normal unit vector at the time of collision (dless); and r denotes the position vectors (m/s).
(5)Fcti=Stζ,St=8GRα,G=2v12G1+2v22G2,
where St is the tangential stiffness (N/m); ζ is the overlap along the tangent direction (m); G is the effective shear modulus (Pa); and G refers to the shear modulus of the particle (Pa).
(6)Fdti=256λSnmvtrel,vtrel=v1v2vnreln,
where vtrel denotes the relative velocity along the tangent direction (m/s).
(7)Tij=μrFcniRiwi,
where μr is the rolling friction coefficient (dless); Ri is the length from the centroid to contact point (m); and wi refers to the unit angular velocity vector at the contact point.
The CFD and the DEM modules are coupled by the forces of fluid applied on the proppant.
(8)Fdi=Vβεsvsvf,
where V refers to the volume of proppant (m3); εs denotes the volume fraction of proppant (dless); β is calculated with the Gidaspow drag model (kg/m3s) [5, 24], as shown in Equation (9); and the drag coefficient CD is considered by Equation (10).
(9)β=34CDεsρfvsvfdsεf1.65,αf>0.8,β=150εs2ufεfds2+1.75ρfεsdsvsvf,αf0.8,(10)CD=0.44,Res>1000,24εfRes1+0.15εfRes0.687,Res1000,
where μf is the fluid viscosity (Pa s); εs denotes the volume fraction of the solid phase (dless); ds is the proppant density (kg/m3); and Res is the particle Reynolds number, which can be calculated with Equation (11) (dless).
(11)Res=ρfvsdsμf.

2.1.2. Fluid Control Equations

In the CFD part, the fluid phase is solved directly by the Eulerian-Eulerian model [8, 11, 14].
(12)tεfρf+εfρfvf=0,tεfρfvf+εfρfvfvf=εfp+εfτf+εfρfg1Vcelli=1mFdi,
where εf refers to the volume fraction of the fluid (dless); ρf refers to the density of the fluid (kg/m3); vf refers to the average velocity of the fluid in a cell (m/s); t refers to the time (s); g refers to the gravitational acceleration (m/s2); τf refers to the stress-strain tensor (Pa); m refers to the number of particles in one fluid cell (dless); and Fdi refers to the momentum exchange between the fluid and solid phases (N).
The turbulence fluctuation was accounted for a standard kε turbulence model, as follows [14]:
(13)εfρfkt+εfρfvfk=εfμtσkk+Gk+χkεfρfε,εfρfεt+εfρfvfε=εfμtσεε+εfεkC1εGkC2ερfε+χε,
where k refers to the turbulent kinetic energy (m2/s2); ε refers to the turbulent dissipation rate (m2/s3); Gk represents the generation of turbulent kinetic energy (kg/(m3 s)); χk and χε is the turbulent exchange terms between the fluid phase and particles (kg/(m3 s)); μt refers to the turbulent eddy viscosity (Pa s); σk=1 and σε=1.3 denote the Prandtl numbers of turbulent kinetic energy and turbulent dissipation rate (dless), respectively; and C1ε=1.44 and C2ε=1.92 are empirical constants (dless).

2.2. Initial and Boundary Conditions

Figure 1 illustrates the computational domain of the fracture, where L=0.24m, W=0.006m, and H=0.032m. The sand-carrying fracturing fluid was injected into the fracture at the perforation and flowed out from the right outlet. In the simulations, the perforation was simplified to a small rectangle (0.004×0.006m), which was located at the center of the left boundary of the fracture.

In the CFD part, we used a structured mesh that was generated by hexahedral elements, as illustrated in Figure 2. The cell size was 0.002×0.002×0.002m, which ensured that the mesh size was larger than the particle diameter. The simulations applied a velocity boundary condition and a constant pressure boundary condition to the inlet and outlet, respectively. The timestep of the CFD solver was 8×104s.

In the DEM part, a particle factory was built at the position of the perforation. Randomly located particles were generated from the factory at a fixed rate. The generated rate of the particles was calculated based on the sand ratio (particle concentration) in the fracturing fluid. In addition, the velocity of the generated particles was the same as that of the injected fluid. The timestep of the DEM solver was 8×106s (29% of the Raleigh times), which is one hundredth of that in the CFD solver. The mesh size used in the DEM part is 2.5 times the particle radius.

The CFD and DEM modules are coupled through the drag force of fluid acting on particles, and the schematic of the CFD-DEM coupling calculation process is shown in Figure 3.

3. Results and Discussion

Figure 4 illustrates the result of proppant transport in the fracture simulated with the parameters in Table 1. The proppant settles out of suspension rapidly due to the effect of gravity when it is injected into the fracture, and a dune forms near the wellbore. The flow channel of the fluid above the sand bank decreases gradually as the height of the sand bank increases, which increases the shear force acting on the dune generated by the fluid. A bed load layer is generated when the shear force reaches the critical Shields number. However, the height of the dune increases gradually due to the continuous injection of proppant. When the shear force is enough to transport the injected proppants to move forward, the height of the dune stops increasing and reaches an equilibrium height. Due to the shear effect of fluid acting on the sand bank, there is an unsupported area at the front of the sand bank. The size of the unsupported area is related to the injection rate, the settlement velocity of the proppant, and the position of the perforation.

As seen from the process of proppant transport, the movement of the bed load layer only increases the time for the sand bank to reach the equilibrium height. The dune still reaches the equilibrium height near the wellbore rapidly, and the newly injected proppant moves forward with the bed load mechanism. The final morphology of the sand bank is less affected by the migration of the bed load layer. Figure 5 illustrates the mechanism of proppant migration in the fracture after the sand bank reaches the equilibrium height, in which H1, H2, and H3 denote the heights of the fluid flow area, the bed load layer, and the immobile sand bank, respectively.

4. Mass Flux of the Bed Load Layer

4.1. Model

The Shields number is a common parameter that governs bed load proppant transport, as shown in Equation (14) [1, 2527].
(14)S=τρsρfgds,
where τ is the shear stress of fluid acting on the surface of the dune (N); ρs refers to the proppant density (kg/m3); and S is the Shields number (m2).
When the Shields number is low, the sand bank is immobile, and when the Shields number reaches the critical level, the proppant up the sand bank begins to move in the form of a bed load under the action of the fluid drag force. As the shear stress can be calculated with Equation (15) in a fracture [1], the critical average velocity can be expressed as Equation (16), which denotes the average velocity of the fluid up the sand bank at which bed load begins.
(15)τ=ημfVf¯W,(16)Vfc¯=ScρsρfgdsWημf,
where Vf¯ refers to the average velocity of fluid up the sand bank (m/s); W refers to the fracture width (m); η refers to a coefficient, which is a constant (dless); Vfc¯ is the critical average velocity (m/s); and Sc is the critical Shields number, which varies between 0.03 and 0.06 [27].
Wang and Zheng [28] established the equation of mass flux of bed load layer as Equation (17), and as seen from it, the mass flux of the bed load is affected by the density of the fluid and proppant, the proppant diameter, the shear velocity of fluid acting on the bed, the static friction coefficient, and the rolling friction coefficient.
(17)Qcreep=ρf22ρsεs¯gμm2μmμdcdsg×U4,
where c1.5, μm refers to the static friction coefficient between proppants (dless); μd refers to the rolling friction coefficient of the proppant (dless); U refers to the shear velocity of fluid acting on the bed (m/s); and εs¯ denotes the average volume fraction of proppant in the bed load layer (dless).
The shear velocity can be calculated with Equation (18) [28, 29]. Then, the mass flux of the bed load layer can be calculated with Equation (19).
(18)U=τρf,(19)Qcreep=12ρsεs¯gμm2μmμdcdsg×ημfVf¯W2.

As the coefficient η is unknown, the simulation results of the equilibrium height with various parameters were used to regress the coefficient. Table 2 shows the simulation results of the equilibrium height with various parameters.

In the simulations, assume that all the injected proppants migrate in the form of bed load when the sand bank reaches the equilibrium height. The mass flux of the bed load layer can be calculated with
(20)Q=creepmsNsW,
where Ns is the number of injected proppants per unit time (1/s) and ms is the mass of proppant (kg).
The thickness of the bed load layer can be calculated as
(21)H2=msNsεs¯Wu¯creep,
where u¯creep is the average velocity of the bed load layer (m/s) and H2 is the thickness of the bed load layer (m).
The velocity distribution of the creep layer is complex [30, 31]. The Saint-Venant model was used to calculate the velocity of the bed load layer along the Z-direction [28] to simplify the calculation.
(22)ucreep=2UH2y+H2,
where ucreep is the velocity of the bed load layer (m/s) and y is the coordinate along the Y-direction (m). And the coordinate system is shown in Figure 5.
Then, the average velocity of the bed load layer can be expressed as
(23)u¯creep=H202U/H2y+H2dyH2=U=ημfVf¯Wρf.
Assume that the velocity of the fluid can be approximated as the proppant velocity in the bed load layer. Then, the average velocity of fluid up the bed load layer can be expressed as
(24)Vf¯=Qfinu¯creepWH21εs¯H1W,
where Qfin is the volume of injected fluid per unit time (m3/s).

According to the simulation results, the coefficient η can be calculated by combining Equations (19)–(24). Table 3 shows the results of the coefficient η based on the simulation results in Table 2 when the average volume fraction of proppant in the bed load layer is 0.5. As shown, the average result of η is 268.

During actual slickwater fracturing, the mass flux of the bed load layer can be calculated with Equation (25). The shear velocity, the thickness of the bed load layer, and the height of the fluid flow area can be calculated with Equations (26)–(28). Then, the equilibrium height of the sand bank can be calculated with Equation (29).
(25)Q=creepQpinρsεsW,(26)U=QpinρsεsW2ρsεs¯gμm2ρf2cdsgμmμd1/4,(27)H2=QpinεsWu¯creepεs¯=QpinεsWUεs¯,(28)H1=QfinQpinεs1εs¯/εs¯ημfW2ρfU2,(29)HEH=HQfinQpinεs1εs¯/εs¯ημfW2ρfU2,
where Qpin is the packing volume of injected proppant per unit time (m3/s); εs is the volume fraction of proppant dune, which is generally taken as 0.6 (dless); HEH is the equilibrium height of the sand bank (m); and H is the height of the fracture (m).

According to the simulation results shown in Table 2 and the equation of mass flux of the bed load layer, the mass flux of the bed load layer decreases, and the equilibrium height increases as the proppant density, proppant diameter, or the rolling friction coefficient and static friction coefficient of proppant increase, but the mass flux of the bed load layer increases, and the equilibrium height decreases as the fluid viscosity increases.

4.2. Model Validation

The experimental data of Wang et al. [29] were used to verify the built model. The size of the fracture used in their experiments was 2.44mL×0.305mH×0.00794mW, and the other parameters in their experiments were ρf=1000kg/m3, ρs=2650kg/m3, and dp=0.0006m. The static friction coefficient and rolling friction coefficient between the proppants (sand) in their experiments are not given. We took larger values than ceramic proppant [8], where μm=0.7 and μd=0.2. The results with various injection rates of fluid and proppant by the built model are shown in Table 4. Figures 6 and 7 show the comparison of H1 and H2 between the experimental data and the calculated results by the built model, respectively. As shown, the calculated results by the built model fit well with the experimental results when the sand ratio is larger than 0.5%.

5. Conclusions

In this study, bed load proppant transport was studied. From the research, the following conclusions can be obtained:

  • (1)

    Bed load proppant transport is an essential mechanism during slickwater hydraulic fracturing. The mass flux of the bed load layer mainly affects the equilibrium height of the sand bank

  • (2)

    The mass flux of the bed load layer increases with an increase in the fluid viscosity and with a decrease in the proppant density, proppant diameter, and friction coefficients of the proppant

  • (3)

    A model was built to calculate the mass flux of the bed load layer and the equilibrium height of the sand bank. The calculated results by the built model fit well with the experimental results when the sand ratio is larger than 0.5%

Data Availability

The results of the simulations or calculated by the built model (sand bank) used to support the findings of this study are included within the article. Previously reported (experimental results) data in Figures 5 and 6 were used to support this study and are available at doi:10.1016/S0301-9322(02)00152-0. This study is cited at relevant places within the text as reference [23].

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

We acknowledge the funding by the National Natural Science Foundation of China (Grant No. 51804175).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).