# DEM-CFD Modeling of Proppant Pillar Deformation and Stability during the Fracturing Fluid Flowback

Authors: Haiyan Zhu; Jiadong Shen; Fengshou Zhang; Bo Huang; Liaoyuan Zhang; Wei Huang; John D. McLennan
Journal: Geofluids
 (2018)

---

## Abstract

In this study, proppant pillar deformation and stability during the fracturing fluid flowback of channel fracturing was simulated with DEM-CFD- (discrete element method-computational fluid dynamics-) coupling method. Fibers were modeled by implementing the bonded particle model for contacts between particles. In the hydraulic fracture-closing period, the height of the proppant pillar decreases gradually and the diameter increases as the closing stress increases. In the fracturing fluid flowback period, proppant particles could be driven away from the pillar by the fluid flow and cause the instability of the proppant pillar. The proppant flowback could occur easily with large proppant pillar height or a large fluid pressure gradient. Both the pillar height and the pillar diameter to spacing ratio are key parameters for the design of channel fracturing. Increasing the fiber-bonding strength could enhance the stability of the proppant pillar.


---

## Body

## 1. Introduction

Channel fracturing is a relatively new stimulation technique first proposed by Gillard et al. in 2010 [1]. It is mainly composed of three technical parts: the pulse pumping technique, the multicluster perforating process, and the injection of fracturing fluid mixed with fibers. The major difference between channel fracturing and conventional fracturing lies in the pattern of proppant placement. In channel fracturing, the fibers are expected to keep proppant pulses cohesive and prevent them from spreading when traveling through the fracture slots. As a result, the proppant after channel fracturing could assemble as clusters in the fracture. When the fractures close during fracturing fluid flowback or the production process, those proppant clusters are analogous to pillars which resist the closing process. The void space between proppant pillars form the fluid channels which greatly improve the fracture conductivity during the production [2]. An investigation based on the results of more than 1000 times channel fracturing found that more than 99.9% of the channel-fracturing jobs fully completed the proppant placement and by average channel fracturing could save 43% of the proppant compared with the conventional technique implemented in adjacent wells [3, 4]. The oil production with channel fracturing could also be greatly enhanced. For instance, the channel-fracturing operation of tight oil and gas reservoirs in Ordos Basin, China, has produced 2.4 times as much oil as conventional fracturing, and the gas well production is 4-5 times as high as that of conventional fracturing [5].After the hydraulic-fracturing operation, the fracturing fluid flows back to the wellbore which results in additional drag force to the proppant particles. Consequently, the proppant pillar might lose its stability due to the flow of proppant, resulting in sanding and significant decrease of fracture conductivity. Therefore, how to avoid the flowback of proppant has attracted great attention. Channel-fracturing technology adds fibers to the fracturing fluid by which the rheological properties of the proppant are changed. The fracturing fluid mixed with low concentration short fibers can be affected by the formation temperature and exhibit some adhesion [6, 7]. Figure 1 shows the network structure of fibers and proppant mixture. The bonding strength of fibers could effectively glue the proppant particle together, and thus, the fiber-proppant cluster is able to better support the hydraulic fractures, prevent the proppant flowback, and also reduce the risk of sanding. In the channel-fracturing tests of Eagle Ford and other shale fields, sand production was improved and fibers showed better proppant flowback resistance [5, 8, 9].Figure 1
Network structure of fibers and proppant mixture [6, 7].The network structure of fibers and proppant mixture and the adhesion of fibers make the mechanical mechanism of the proppant pillar greatly different from the one for the traditional proppant placement scheme. The existing literatures mainly focus on the mechanical process of the flowback of fiber-free proppant. Vreeburg et al. [10] performed a laboratory study on two types of back-production to help clarify the effect of curing temperature, water production rate, proppant size, and stress cycling on the integrity of resin-coated proppant (RCP) packs. The number of applied stress cycles and the initial RCP pack strength appear to be the dominant factors that govern proppant back-production. Goel and Shah [11] investigated proppant flowback phenomena using a large-scale fracturing simulator. Experimental results show that flowback initiates at lower cleanup rates when the closure stress increases or when the fracture width increases and flowback initiates at higher cleanup rates when the sand size increases. Smith et al. [12] investigated the effects of various factors on the proppant flowback, including the differential-fracture closure, leak-off of fracturing fluid, and fracturing fluid rheology, and claimed that those effects should be combined in the proppant transport analysis. Aidagulov et al. [13] presented a quantitative model to predict proppant flowback based on treating both the proppant pack and the reservoir as poro-elasto-plastic media. But with adoption of such a continuum model, interactions between flowback proppant particles would not be observed directly. Daneshy [14] suggested that the main cause of the proppant flowback was the presence of shear fractures which led to the formation of many randomly distributed tight proppant packs. The three requirements for proppant flowback were motion initiation, motion maintenance, and infinite conductivity along the return path, and gravity played a very important role on proppant flowback. Hu et al. [15] established a proppant mechanical model to predict critical flowback velocity. They concluded that adopting a lower velocity before the fracture closes and a higher velocity after the fracture closes helps to discharge fracturing fluid in time and effectively prevent proppant flowback. Qi and Jiao [16] presented a model to predict proppant flowback in a fractured gas well. However, with the hypothesis of a one-dimensional flow and ignoring the influence of gravity, this model cannot demonstrate real interactions between liquid and proppants. McLennan et al. [17] described the experimental results of a radially convergent flow and a linear flow into a wellbore through an initially packed fracture. Results of a radial flow through a transverse fracture suggested the creation of flow channels, while results of a linear flow through a longitudinal fracture suggested more proppant removal. The above studies have pointed out that effects such as fluid viscosity, fluid pressure gradient, and closing stress are crucial to the proppant flowback for the conventional proppant placement scheme. They laid an important foundation for the research of proppant flowback with channel fracturing.Other studies explored the evolution of fracture width and conductivity in channel fracturing. The deformation characteristics of the proppant pillar are described using numerical and analytical methods though some of them are not quite reasonable. Gomaa et al. [18] described experimental and numerical models to generate pillar-propped fractures based on fingering phenomena observed in fluid injection. Experimental and numerical results confirmed that increasing the injection rate reduced the main channel width and increased branching. Neto and Kotousov [19] developed a simplified semianalytical method for calculating the residual opening of fractures partially filled with proppant. They used the one-dimensional Terzaghi’s consolidation model to describe proppant response. Both pore water and particles are assumed to be incompressible and changes of void volume are linked to deformation. Zheng et al. [20] derived an analytical model to calculate permeability and conductivity in channel fracturing. The proppant pillar is treated as a cylindrical indenter which is rigid. Guo et al. [21] established an analytical model to describe fracture aperture change and conductivity for cuboid shaped pillars. Yan et al. [22] developed an analytical model to represent the physical deformation of channel-fracturing fractures, and the Darcy-Brinkman equation was applied to simulate the flow in pillars and fluid channels. Both Guo et al. and Yan et al. calculate the axial deformation of the proppant pillar using Hooke’s law.Later, Deng et al. [23] adopted the discrete element method to study interactions between shale and proppant. Effects of factors such as shale modulus and proppant size on fracture aperture were numerically modelled and learned. Fan et al. [24] developed a DEM-LB workflow to understand the interaction between reservoir depletion, proppant compaction, and single-/multiphase flows in a hydraulic fracture. Furthermore, Zhang et al. [25] developed an integrated DEM-CFD modeling workflow to model proppant embedment and fracture conductivity. Fracture conductivity after fracture closing was evaluated by modeling fluid flow through the proppant pack by use of DEM coupled with CFD. Results showed that the fracture conductivity increased with the increase of proppant concentration or proppant size and decreased with the increase of fracture-closing stress or degree of shale hydration. Shale-hydration effect was confirmed to be the main reason for the large amount of proppant embedment. With increasingly frequent application of DEM-CFD, it was believed that DEM-CFD was the most suitable computational method for modeling two-way solid–fluid interactions (e.g., Tomac and Gutierrez [26]). This paper utilizes the DEM-CFD-coupling method, by which the cementation between the fibers and the proppant particles is mimicked by the tensile and shear bonding strength between discrete element particles, to set up the hydromechanically coupled model for the interactions of proppant-fiber-fracturing fluid. The model aims to reveal the mechanism of proppant pillar deformation and stability during fracturing fluid flowback by combing the analysis at both micro- and macroscales. The effects of pressure gradient, fluid viscosity, pillar height, pillar diameter to spacing ratio, and bonding strength of fibers are investigated, and the amount of flowback proppant and the spreading area of proppant are studied. This work could provide a potential guideline and theoretical background for the design of channel fracturing and also the optimization of field fracturing operation.


## 2. Theoretical Background and Numerical Model Setup

### 2.1. DEM-CFD-Coupled Model for Proppant Pillar Deformation and Stability during the Fracturing Fluid Flowback of Channel Fracturing

#### 2.1.1. Proppant Particle-Fiber Interaction

Romero and Feraud [27] first studied the stability of a proppant pack reinforced with fibers and presented a fiber-reinforced failure criterion in their study. They set up a laboratory model to demonstrate the effect of fiber properties on proppant pack stability. Two modes of failure were distinguished in their experiments: collapse of the arch at the perforation and total failure of the proppant pack creating a channel within the fracture. The fracturing fluid used in the channel fracturing contains fibers that have wrapping and restraining effects on proppant, which can enhance the sand-carrying capacity of the fracturing fluid. When proppant and fibers are injected into the fractures with fracturing fluid, the bonding and friction between the fibers and proppant make the proppant cluster as a whole, which achieves the balance under the action of closing stress and lateral stress. During the flowback, the proppant pillar is subjected to shear deformation by the fluid force; however, the existence of fibers helps to prevent the collapse of the proppant pillar.In this study, we apply the bonded particle model in DEM to model the interaction between proppant particles and fibers [28, 29]. DEM is an ideal tool to model the mechanical behaviors of a granular assembly such as a proppant pack. The macroscale mechanical response emerges from the interaction of particles or grains through their contact behaviors at the particle scale, either frictional or with cementation. Since the shape of fibers is irregular, it is technically difficult to model fibers together with proppant particles. Thus, we indirectly model the fibers by implementing the bonded particle model which adds bonding strength between proppant particles but without explicitly representing the fibers. The theory of DEM and bonded particle model are well documented in the literatures and hence are not presented here for brevity.


#### 2.1.2. DEM-CFD Coupling for Proppant Particle-Fracturing Fluid Interaction

By considering force balance for particles within a cell, the driving force applied to a single particle is given by (1). The driving force consists of two kinds of forces caused by the flow. One is the viscous force induced by fluid viscous properties and characterized by shear force on the particle surface. The other is the force acting on the particles by pressure gradient and characterized by a normal force on the particle surface [30].
(1)fdi=−fint1−ϕ+∇pπ6dpi3,where fdi is the driving force applied to particle i, fint is the interaction force per unit volume between particles and fluid, dpi is the diameter of particle i, ϕ is the porosity of this cell, and ∇p is the pressure gradient.In the first term of (1), fint represents the interaction force per unit volume between particles and fluid. The minus sign means that the force applied to the fluid is opposite in direction with fint. The sign of fint depends on the relative moving direction of particles in flow. If the relative moving direction is opposite between particles and the flow, fint is negative and the driving force is positive, which means particles are pushed forward by the flow. fint equals 0 when particles are static relative to the flow. In the second term of (1), the minus sign means that pressure decreases with positive flow when ∇p is negative. The force acting on the particles by the pressure gradient pushes particles forward if ∇p is negative.When fracturing fluid begins to flow back, the flow direction is opposite to the relative moving direction of proppants and pressure decreases with positive flow. Thus, both viscous force and the force acting on the particles by the pressure gradient work as positive forces to particles.Interaction force per unit volume between particles and fluid in thej direction is given by
(2)fintj=βintjv¯j−uj,where uj is fluid velocity in the j (j=x,y,z) direction, v¯j is average velocity of particles in the j direction, and βintj is the fluid-particle friction coefficient. βintj is defined by different expressions according to different porosity values. When porosity is lower than 0.8, βintj is given by
(3)βintj=1501−ϕ2d¯p2ϕμf+1.751−ϕd¯pρfv¯j−uj,where d¯p is the average diameter of particles, μf is the kinematic viscosity of fluid, and ρf is the density of fluid. When porosity is higher than 0.8, βintj is given by
(4)βintj=34ρfCD1−ϕϕ−1.7d¯pv¯j−uj,where CD is the drag coefficient (dimensionless).Pressure gradient∇pj also has two expressions accordingly. When porosity is lower than 0.8, ∇pj is given by Ergun’s equation:
(5)∇pj=1501−ϕ2d¯p2ϕ2μf+1.751−ϕd¯pϕρfv¯j−ujv¯j−uj.When porosity is higher than 0.8,∇pj is given by the Wen-Yu equation:
(6)∇pj=34ρfCD1−ϕϕ−2.7d¯pv¯j−ujv¯j−uj.


### 2.2. Numerical Model Setup

Figure2 shows the numerical model setup. The model consists of one proppant pillar which is sandwiched by two shale plates. The rock plate and proppant pillar are made of blue particles with a face-centered cubic (FCC) structure and yellow particles with random distribution, respectively. Taking the 6 mm high proppant pillar as an example, the simulation procedure can be summarized in the following steps.
(1)
Generate a cubic-shape rock sample with the FCC structure. The initial model size is 20 mm × 20 mm × 30 mm. The microscale properties of rock are assigned to all contacts.
(2)
A small confining stress of 0.5 MPa is applied to the sample by a servocontrol procedure.
(3)
A 6 mm thick layer in the middle of the sample is deleted to generate a fracture slot perpendicular to thez direction.
(4)
A cylindrical proppant pillar with a height of 6 mm and diameter of 10 mm is generated in the middle of the fracture slot.
(5)
The two rock plates are gradually loaded with the compression stress of 41.4 MPa while the proppant pillar is deformed.
(6)
After the loading, the fluid grid is set up for the proppant layer and fluid flow calculation is carried out.Figure 2
(a) Model perspective view and (b) coupled DEM-CFD analysis for proppant pillar deformation and stability. The blue and yellow particles represent the rock matrix and proppant pillar, respectively, and the green meshes represent the CFD fluid grids.
(a)
(b)For the calculation of fluid-particle coupling, 10 fluid cells are set up along the two horizontal directions and 1 fluid cell is set up along the vertical direction of the proppant pillar height. In order to apply the fluid boundary condition, one additional layer of fluid cells is needed at each of the six faces of the calculation domain. It means that in each of the three directions (x, y, z), there are two additional boundary layers. The final cell pattern is 10+2×10+2×1+2=12×12×3 given that the initial cell pattern is 10×10×1. Therefore, the gird mesh has 12×12×3=432 fluid cells in total. With the two additional boundary layers, the final model size in both the x and y directions is 10+2×2mm=24mm (Figure 3). Cells outside the red line in Figure 4 are boundary cells. Only cells within this “boundary” ran the fluid-particle interaction computation.Figure 3
Illustration of model size: (a) grid mesh without particles; (b) size of boundary cells and flow area.
(a)
(b)Figure 4
Coupled DEM-CFD model: (a) vertical view; (b) front view.
(a)
(b)In fluid flow calculation, a negative pressure gradient is applied along thex direction by applying positive fluid pressure on the left boundary at the inlet and applying zero fluid pressure on the right boundary at the outlet. The flowback pattern of the proppant pillar at different times can then be simulated.


### 2.3. Calibration of the Numerical Model with Experimental Results

#### 2.3.1. Microscale Properties of Rock

The microscale properties of rock are firstly calibrated by comparing the measured rock macroscale properties in numerical modeling with the given experimental results (uniaxial compression strength 215 MPa, elastic modulus 30 GPa, and Poisson’s ratio 0.28). The summary of calibrated microscale properties of rock is listed in Table1.Table 1
Summary of microscale properties of rock.
Microparameters
Linear contacts
Apparent modulusEc (GPa)
5.7
Normal stiffnesskn (N/m)
2.6ks
Shear stiffnessks (N/m)
2DEc
Parallel bonds
Apparent modulusEc (GPa)
5.7
Normal stiffnesskn (GPa/m)
2.6ks
Shear stiffnessks (GPa/m)
EcD
Normal bond strength (MPa)
38
Shear bond strength (MPa)
38
Radius multiplier
1
Friction
0.5
Density (kg/m3)
2650
Particle radius (D/2) (mm)
0.2


#### 2.3.2. Microscale Properties of the Proppant Pillar

In order to calibrate the microscale properties of the proppant pillar, lab experiments of proppant pillar compression are carried out first with the following procedure.(1)
Put proppant particles and fibers (weight ratio 1 : 0.004) into a beaker.
(2)
Use a glass rod to mix the proppant and fibers (Figure5(a)).
(3)
Add glue into the proppant-fiber mixture and stir the sample until it is solidified.
(4)
Fill the sample into a metal mold (hollow cylinder with inner diameter of 10 mm and height of 10 mm) and make proppant pillars (Figure5(b)).
(5)
Heat the proppant pillars in an oven for 1 hour with the temperature of 60 degrees, then cool them down under ambient condition.
(6)
Fill the array of proppant pillars in the testing chamber designed based on the API standard for proppant conductivity testing; load the pillars at the given normal loading stress and keep the loading for 1 minute until reaching the steady state (Figure5(c)).
(7)
Calculate the change of pillar height based on the recorded variation of vertical displacement.
(8)
Disassemble the experimental setup and record the profiles of proppant pillars after the testing (Figure6(a) and 6(b)).
(9)
Repeat step (1) to (8) for a different closing stress.Figure 5
(a) Proppant-fiber mixture, (b) array of proppant pillars, and (c) testing chamber with proppant pillars placed.
(a)
(b)
(c)Figure 6
Comparison of proppant pillar profiles from experiments and from numerical modeling: (a) lateral view before the testing in the experiment, (b) top view after the testing in the experiment, (c) lateral view before the testing in the numerical modeling, and (d) top view after the testing in the numerical modeling.
(a)
(b)
(c)
(d)To calibrate the microscale properties of proppant pillars, a numerical sample of the proppant pillar was first made and then compressed with constant stress from the two ends (Figure6(c) and 6(d)). Numerical tests were carried out with different compression stresses, and the pillar height and pillar diameter after the tests were compared with the experimental results for the model calibration, as shown in Figures 7 and 8. The final calibrated microscale properties of proppant particles are listed in Table 2.Figure 7
Comparison of numerical modeling results and experimental results for pillar height at different closing stress.Figure 8
Comparison of numerical modeling results and experimental results for pillar diameter at different closing stresses.Table 2
Summary of microscale properties of proppant.
Proppant size combination
40/70
30/50
20/40
Microparameters
Linear contacts
Apparent modulusEc (MPa)
100
Normal stiffnesskn (N/m)
2DEc
Shear stiffnessks (N/m)
kn
Parallel bonds
Apparent modulusEc (MPa)
100
Normal stiffnesskn (MPa/m)
Ec/D
Shear stiffnessks (MPa/m)
kn
Normal bond strength (Pa)
3000
Shear bond strength (Pa)
3000
Radius multiplier
1
Friction
0.9
Density (kg/m3)
2650
Particle radius (D/2) (mm)
0.21–0.42
0.3–0.6
0.42–0.84


## 2.1. DEM-CFD-Coupled Model for Proppant Pillar Deformation and Stability during the Fracturing Fluid Flowback of Channel Fracturing

### 2.1.1. Proppant Particle-Fiber Interaction

Romero and Feraud [27] first studied the stability of a proppant pack reinforced with fibers and presented a fiber-reinforced failure criterion in their study. They set up a laboratory model to demonstrate the effect of fiber properties on proppant pack stability. Two modes of failure were distinguished in their experiments: collapse of the arch at the perforation and total failure of the proppant pack creating a channel within the fracture. The fracturing fluid used in the channel fracturing contains fibers that have wrapping and restraining effects on proppant, which can enhance the sand-carrying capacity of the fracturing fluid. When proppant and fibers are injected into the fractures with fracturing fluid, the bonding and friction between the fibers and proppant make the proppant cluster as a whole, which achieves the balance under the action of closing stress and lateral stress. During the flowback, the proppant pillar is subjected to shear deformation by the fluid force; however, the existence of fibers helps to prevent the collapse of the proppant pillar.In this study, we apply the bonded particle model in DEM to model the interaction between proppant particles and fibers [28, 29]. DEM is an ideal tool to model the mechanical behaviors of a granular assembly such as a proppant pack. The macroscale mechanical response emerges from the interaction of particles or grains through their contact behaviors at the particle scale, either frictional or with cementation. Since the shape of fibers is irregular, it is technically difficult to model fibers together with proppant particles. Thus, we indirectly model the fibers by implementing the bonded particle model which adds bonding strength between proppant particles but without explicitly representing the fibers. The theory of DEM and bonded particle model are well documented in the literatures and hence are not presented here for brevity.


### 2.1.2. DEM-CFD Coupling for Proppant Particle-Fracturing Fluid Interaction

By considering force balance for particles within a cell, the driving force applied to a single particle is given by (1). The driving force consists of two kinds of forces caused by the flow. One is the viscous force induced by fluid viscous properties and characterized by shear force on the particle surface. The other is the force acting on the particles by pressure gradient and characterized by a normal force on the particle surface [30].
(1)fdi=−fint1−ϕ+∇pπ6dpi3,where fdi is the driving force applied to particle i, fint is the interaction force per unit volume between particles and fluid, dpi is the diameter of particle i, ϕ is the porosity of this cell, and ∇p is the pressure gradient.In the first term of (1), fint represents the interaction force per unit volume between particles and fluid. The minus sign means that the force applied to the fluid is opposite in direction with fint. The sign of fint depends on the relative moving direction of particles in flow. If the relative moving direction is opposite between particles and the flow, fint is negative and the driving force is positive, which means particles are pushed forward by the flow. fint equals 0 when particles are static relative to the flow. In the second term of (1), the minus sign means that pressure decreases with positive flow when ∇p is negative. The force acting on the particles by the pressure gradient pushes particles forward if ∇p is negative.When fracturing fluid begins to flow back, the flow direction is opposite to the relative moving direction of proppants and pressure decreases with positive flow. Thus, both viscous force and the force acting on the particles by the pressure gradient work as positive forces to particles.Interaction force per unit volume between particles and fluid in thej direction is given by
(2)fintj=βintjv¯j−uj,where uj is fluid velocity in the j (j=x,y,z) direction, v¯j is average velocity of particles in the j direction, and βintj is the fluid-particle friction coefficient. βintj is defined by different expressions according to different porosity values. When porosity is lower than 0.8, βintj is given by
(3)βintj=1501−ϕ2d¯p2ϕμf+1.751−ϕd¯pρfv¯j−uj,where d¯p is the average diameter of particles, μf is the kinematic viscosity of fluid, and ρf is the density of fluid. When porosity is higher than 0.8, βintj is given by
(4)βintj=34ρfCD1−ϕϕ−1.7d¯pv¯j−uj,where CD is the drag coefficient (dimensionless).Pressure gradient∇pj also has two expressions accordingly. When porosity is lower than 0.8, ∇pj is given by Ergun’s equation:
(5)∇pj=1501−ϕ2d¯p2ϕ2μf+1.751−ϕd¯pϕρfv¯j−ujv¯j−uj.When porosity is higher than 0.8,∇pj is given by the Wen-Yu equation:
(6)∇pj=34ρfCD1−ϕϕ−2.7d¯pv¯j−ujv¯j−uj.


## 2.1.1. Proppant Particle-Fiber Interaction

Romero and Feraud [27] first studied the stability of a proppant pack reinforced with fibers and presented a fiber-reinforced failure criterion in their study. They set up a laboratory model to demonstrate the effect of fiber properties on proppant pack stability. Two modes of failure were distinguished in their experiments: collapse of the arch at the perforation and total failure of the proppant pack creating a channel within the fracture. The fracturing fluid used in the channel fracturing contains fibers that have wrapping and restraining effects on proppant, which can enhance the sand-carrying capacity of the fracturing fluid. When proppant and fibers are injected into the fractures with fracturing fluid, the bonding and friction between the fibers and proppant make the proppant cluster as a whole, which achieves the balance under the action of closing stress and lateral stress. During the flowback, the proppant pillar is subjected to shear deformation by the fluid force; however, the existence of fibers helps to prevent the collapse of the proppant pillar.In this study, we apply the bonded particle model in DEM to model the interaction between proppant particles and fibers [28, 29]. DEM is an ideal tool to model the mechanical behaviors of a granular assembly such as a proppant pack. The macroscale mechanical response emerges from the interaction of particles or grains through their contact behaviors at the particle scale, either frictional or with cementation. Since the shape of fibers is irregular, it is technically difficult to model fibers together with proppant particles. Thus, we indirectly model the fibers by implementing the bonded particle model which adds bonding strength between proppant particles but without explicitly representing the fibers. The theory of DEM and bonded particle model are well documented in the literatures and hence are not presented here for brevity.


## 2.1.2. DEM-CFD Coupling for Proppant Particle-Fracturing Fluid Interaction

By considering force balance for particles within a cell, the driving force applied to a single particle is given by (1). The driving force consists of two kinds of forces caused by the flow. One is the viscous force induced by fluid viscous properties and characterized by shear force on the particle surface. The other is the force acting on the particles by pressure gradient and characterized by a normal force on the particle surface [30].
(1)fdi=−fint1−ϕ+∇pπ6dpi3,where fdi is the driving force applied to particle i, fint is the interaction force per unit volume between particles and fluid, dpi is the diameter of particle i, ϕ is the porosity of this cell, and ∇p is the pressure gradient.In the first term of (1), fint represents the interaction force per unit volume between particles and fluid. The minus sign means that the force applied to the fluid is opposite in direction with fint. The sign of fint depends on the relative moving direction of particles in flow. If the relative moving direction is opposite between particles and the flow, fint is negative and the driving force is positive, which means particles are pushed forward by the flow. fint equals 0 when particles are static relative to the flow. In the second term of (1), the minus sign means that pressure decreases with positive flow when ∇p is negative. The force acting on the particles by the pressure gradient pushes particles forward if ∇p is negative.When fracturing fluid begins to flow back, the flow direction is opposite to the relative moving direction of proppants and pressure decreases with positive flow. Thus, both viscous force and the force acting on the particles by the pressure gradient work as positive forces to particles.Interaction force per unit volume between particles and fluid in thej direction is given by
(2)fintj=βintjv¯j−uj,where uj is fluid velocity in the j (j=x,y,z) direction, v¯j is average velocity of particles in the j direction, and βintj is the fluid-particle friction coefficient. βintj is defined by different expressions according to different porosity values. When porosity is lower than 0.8, βintj is given by
(3)βintj=1501−ϕ2d¯p2ϕμf+1.751−ϕd¯pρfv¯j−uj,where d¯p is the average diameter of particles, μf is the kinematic viscosity of fluid, and ρf is the density of fluid. When porosity is higher than 0.8, βintj is given by
(4)βintj=34ρfCD1−ϕϕ−1.7d¯pv¯j−uj,where CD is the drag coefficient (dimensionless).Pressure gradient∇pj also has two expressions accordingly. When porosity is lower than 0.8, ∇pj is given by Ergun’s equation:
(5)∇pj=1501−ϕ2d¯p2ϕ2μf+1.751−ϕd¯pϕρfv¯j−ujv¯j−uj.When porosity is higher than 0.8,∇pj is given by the Wen-Yu equation:
(6)∇pj=34ρfCD1−ϕϕ−2.7d¯pv¯j−ujv¯j−uj.


## 2.2. Numerical Model Setup

Figure2 shows the numerical model setup. The model consists of one proppant pillar which is sandwiched by two shale plates. The rock plate and proppant pillar are made of blue particles with a face-centered cubic (FCC) structure and yellow particles with random distribution, respectively. Taking the 6 mm high proppant pillar as an example, the simulation procedure can be summarized in the following steps.
(1)
Generate a cubic-shape rock sample with the FCC structure. The initial model size is 20 mm × 20 mm × 30 mm. The microscale properties of rock are assigned to all contacts.
(2)
A small confining stress of 0.5 MPa is applied to the sample by a servocontrol procedure.
(3)
A 6 mm thick layer in the middle of the sample is deleted to generate a fracture slot perpendicular to thez direction.
(4)
A cylindrical proppant pillar with a height of 6 mm and diameter of 10 mm is generated in the middle of the fracture slot.
(5)
The two rock plates are gradually loaded with the compression stress of 41.4 MPa while the proppant pillar is deformed.
(6)
After the loading, the fluid grid is set up for the proppant layer and fluid flow calculation is carried out.Figure 2
(a) Model perspective view and (b) coupled DEM-CFD analysis for proppant pillar deformation and stability. The blue and yellow particles represent the rock matrix and proppant pillar, respectively, and the green meshes represent the CFD fluid grids.
(a)
(b)For the calculation of fluid-particle coupling, 10 fluid cells are set up along the two horizontal directions and 1 fluid cell is set up along the vertical direction of the proppant pillar height. In order to apply the fluid boundary condition, one additional layer of fluid cells is needed at each of the six faces of the calculation domain. It means that in each of the three directions (x, y, z), there are two additional boundary layers. The final cell pattern is 10+2×10+2×1+2=12×12×3 given that the initial cell pattern is 10×10×1. Therefore, the gird mesh has 12×12×3=432 fluid cells in total. With the two additional boundary layers, the final model size in both the x and y directions is 10+2×2mm=24mm (Figure 3). Cells outside the red line in Figure 4 are boundary cells. Only cells within this “boundary” ran the fluid-particle interaction computation.Figure 3
Illustration of model size: (a) grid mesh without particles; (b) size of boundary cells and flow area.
(a)
(b)Figure 4
Coupled DEM-CFD model: (a) vertical view; (b) front view.
(a)
(b)In fluid flow calculation, a negative pressure gradient is applied along thex direction by applying positive fluid pressure on the left boundary at the inlet and applying zero fluid pressure on the right boundary at the outlet. The flowback pattern of the proppant pillar at different times can then be simulated.


## 2.3. Calibration of the Numerical Model with Experimental Results

### 2.3.1. Microscale Properties of Rock

The microscale properties of rock are firstly calibrated by comparing the measured rock macroscale properties in numerical modeling with the given experimental results (uniaxial compression strength 215 MPa, elastic modulus 30 GPa, and Poisson’s ratio 0.28). The summary of calibrated microscale properties of rock is listed in Table1.Table 1
Summary of microscale properties of rock.
Microparameters
Linear contacts
Apparent modulusEc (GPa)
5.7
Normal stiffnesskn (N/m)
2.6ks
Shear stiffnessks (N/m)
2DEc
Parallel bonds
Apparent modulusEc (GPa)
5.7
Normal stiffnesskn (GPa/m)
2.6ks
Shear stiffnessks (GPa/m)
EcD
Normal bond strength (MPa)
38
Shear bond strength (MPa)
38
Radius multiplier
1
Friction
0.5
Density (kg/m3)
2650
Particle radius (D/2) (mm)
0.2


### 2.3.2. Microscale Properties of the Proppant Pillar

In order to calibrate the microscale properties of the proppant pillar, lab experiments of proppant pillar compression are carried out first with the following procedure.(1)
Put proppant particles and fibers (weight ratio 1 : 0.004) into a beaker.
(2)
Use a glass rod to mix the proppant and fibers (Figure5(a)).
(3)
Add glue into the proppant-fiber mixture and stir the sample until it is solidified.
(4)
Fill the sample into a metal mold (hollow cylinder with inner diameter of 10 mm and height of 10 mm) and make proppant pillars (Figure5(b)).
(5)
Heat the proppant pillars in an oven for 1 hour with the temperature of 60 degrees, then cool them down under ambient condition.
(6)
Fill the array of proppant pillars in the testing chamber designed based on the API standard for proppant conductivity testing; load the pillars at the given normal loading stress and keep the loading for 1 minute until reaching the steady state (Figure5(c)).
(7)
Calculate the change of pillar height based on the recorded variation of vertical displacement.
(8)
Disassemble the experimental setup and record the profiles of proppant pillars after the testing (Figure6(a) and 6(b)).
(9)
Repeat step (1) to (8) for a different closing stress.Figure 5
(a) Proppant-fiber mixture, (b) array of proppant pillars, and (c) testing chamber with proppant pillars placed.
(a)
(b)
(c)Figure 6
Comparison of proppant pillar profiles from experiments and from numerical modeling: (a) lateral view before the testing in the experiment, (b) top view after the testing in the experiment, (c) lateral view before the testing in the numerical modeling, and (d) top view after the testing in the numerical modeling.
(a)
(b)
(c)
(d)To calibrate the microscale properties of proppant pillars, a numerical sample of the proppant pillar was first made and then compressed with constant stress from the two ends (Figure6(c) and 6(d)). Numerical tests were carried out with different compression stresses, and the pillar height and pillar diameter after the tests were compared with the experimental results for the model calibration, as shown in Figures 7 and 8. The final calibrated microscale properties of proppant particles are listed in Table 2.Figure 7
Comparison of numerical modeling results and experimental results for pillar height at different closing stress.Figure 8
Comparison of numerical modeling results and experimental results for pillar diameter at different closing stresses.Table 2
Summary of microscale properties of proppant.
Proppant size combination
40/70
30/50
20/40
Microparameters
Linear contacts
Apparent modulusEc (MPa)
100
Normal stiffnesskn (N/m)
2DEc
Shear stiffnessks (N/m)
kn
Parallel bonds
Apparent modulusEc (MPa)
100
Normal stiffnesskn (MPa/m)
Ec/D
Shear stiffnessks (MPa/m)
kn
Normal bond strength (Pa)
3000
Shear bond strength (Pa)
3000
Radius multiplier
1
Friction
0.9
Density (kg/m3)
2650
Particle radius (D/2) (mm)
0.21–0.42
0.3–0.6
0.42–0.84


## 2.3.1. Microscale Properties of Rock

The microscale properties of rock are firstly calibrated by comparing the measured rock macroscale properties in numerical modeling with the given experimental results (uniaxial compression strength 215 MPa, elastic modulus 30 GPa, and Poisson’s ratio 0.28). The summary of calibrated microscale properties of rock is listed in Table1.Table 1
Summary of microscale properties of rock.
Microparameters
Linear contacts
Apparent modulusEc (GPa)
5.7
Normal stiffnesskn (N/m)
2.6ks
Shear stiffnessks (N/m)
2DEc
Parallel bonds
Apparent modulusEc (GPa)
5.7
Normal stiffnesskn (GPa/m)
2.6ks
Shear stiffnessks (GPa/m)
EcD
Normal bond strength (MPa)
38
Shear bond strength (MPa)
38
Radius multiplier
1
Friction
0.5
Density (kg/m3)
2650
Particle radius (D/2) (mm)
0.2


## 2.3.2. Microscale Properties of the Proppant Pillar

In order to calibrate the microscale properties of the proppant pillar, lab experiments of proppant pillar compression are carried out first with the following procedure.(1)
Put proppant particles and fibers (weight ratio 1 : 0.004) into a beaker.
(2)
Use a glass rod to mix the proppant and fibers (Figure5(a)).
(3)
Add glue into the proppant-fiber mixture and stir the sample until it is solidified.
(4)
Fill the sample into a metal mold (hollow cylinder with inner diameter of 10 mm and height of 10 mm) and make proppant pillars (Figure5(b)).
(5)
Heat the proppant pillars in an oven for 1 hour with the temperature of 60 degrees, then cool them down under ambient condition.
(6)
Fill the array of proppant pillars in the testing chamber designed based on the API standard for proppant conductivity testing; load the pillars at the given normal loading stress and keep the loading for 1 minute until reaching the steady state (Figure5(c)).
(7)
Calculate the change of pillar height based on the recorded variation of vertical displacement.
(8)
Disassemble the experimental setup and record the profiles of proppant pillars after the testing (Figure6(a) and 6(b)).
(9)
Repeat step (1) to (8) for a different closing stress.Figure 5
(a) Proppant-fiber mixture, (b) array of proppant pillars, and (c) testing chamber with proppant pillars placed.
(a)
(b)
(c)Figure 6
Comparison of proppant pillar profiles from experiments and from numerical modeling: (a) lateral view before the testing in the experiment, (b) top view after the testing in the experiment, (c) lateral view before the testing in the numerical modeling, and (d) top view after the testing in the numerical modeling.
(a)
(b)
(c)
(d)To calibrate the microscale properties of proppant pillars, a numerical sample of the proppant pillar was first made and then compressed with constant stress from the two ends (Figure6(c) and 6(d)). Numerical tests were carried out with different compression stresses, and the pillar height and pillar diameter after the tests were compared with the experimental results for the model calibration, as shown in Figures 7 and 8. The final calibrated microscale properties of proppant particles are listed in Table 2.Figure 7
Comparison of numerical modeling results and experimental results for pillar height at different closing stress.Figure 8
Comparison of numerical modeling results and experimental results for pillar diameter at different closing stresses.Table 2
Summary of microscale properties of proppant.
Proppant size combination
40/70
30/50
20/40
Microparameters
Linear contacts
Apparent modulusEc (MPa)
100
Normal stiffnesskn (N/m)
2DEc
Shear stiffnessks (N/m)
kn
Parallel bonds
Apparent modulusEc (MPa)
100
Normal stiffnesskn (MPa/m)
Ec/D
Shear stiffnessks (MPa/m)
kn
Normal bond strength (Pa)
3000
Shear bond strength (Pa)
3000
Radius multiplier
1
Friction
0.9
Density (kg/m3)
2650
Particle radius (D/2) (mm)
0.21–0.42
0.3–0.6
0.42–0.84


## 3. Proppant Pillar Deformation and Damage Laws during Fracture Closing and Fracturing Fluid Flowback

### 3.1. Fracture-Closing Period

The calibrated microscale properties (Tables1 and 2 are used as input parameters in the model shown in Figure 2. Numerical modeling of fracture closing and fracturing fluid flowback is then calculated. Figures 9 and 10 plot the profiles of a proppant pillar at two different closing stresses, 20.7 MPa and 41.4 MPa, respectively. The height of the proppant pillar is 6 mm, and the diameter is 10 mm. With the increase of the closing stress, the height of the proppant pillar decreases gradually and the diameter increases. The pillar is compressed into a pancake shape. The compression from the shale plate changes the stress state in the proppant pillar and crushes down the pillar. Some particles near the edge of the pillar falls down and rolls away from the pillar.Figure 9
Profiles of the proppant pillar at 20.7 MPa closing stress: (a) lateral view and (b) top view.
(a)
(b)Figure 10
Profiles of the proppant pillar at 41.4 MPa closing stress: (a) lateral view and (b) top view.
(a)
(b)


### 3.2. Fracturing Fluid Flowback Period

The fluid flowback calculation with coupled DEM-CFD analysis was carried out after the fracture closing. Figure11 shows the profiles of the proppant pillar at different flow times. The initial height and diameter of the proppant pillar are 6 mm and 10 mm, respectively. The fluid viscosity is 0.001 Pa·s and the pressure difference is 64 kPa. The black arrows in the plot represent the fluid velocity vector at each fluid cell, and the velocity magnitude is indicated by the length of the arrow. The particles for rock plates are not shown in the plot. With the increase of flow time, the proppant pillar starts to spread and particles around the proppant pillar diffuse outward and form a single layer. After 0.0105 s, the proppant particles gradually escaped from the right side of the flow field and start to flow back. The number of flowback particles in this paper is defined as the number of proppant particles escaping from the flow field.Figure 11
Profiles of the proppant pillar at different flow times: (a) right after fracture closing, (b) flow time of 0.00015 s, (c) flow time of 0.0005 s, and (d) flow time of 0.0105 s. The red square in (b) represents the location of a monitoring fluid cell.
(a)
(b)
(c)
(d)We select the first row of fluid grids at the entry (shown in Figure11(b)) to monitor the change of flow velocity. Fluid velocities reach stable values at about 0.04 s flow time (Figure 12). Before 0.005 s, fluid velocities show fluctuation which indicates that the flow is subjected to unstable change of proppant ahead due to the spreading of particles. Between 0.025 s and 0.035 s, fluid velocities show more drastic fluctuation caused by fast flowback of a number of proppants. The results above indicate that a larger pillar diameter is prone to cause proppant particle spreading as well as unstable fluid flow. On the other hand, a smaller pillar diameter will lead to wider fluid channel and favor the stable fluid flow; however, the pillar-supporting capacity might not be met to resist the fracture-closing stress. Therefore, a proper size of proppant pillar is needed which can satisfy both the stability requirement and conductivity enhancement.Figure 12
Fluid velocity of the 10 fluid cells at the entry.


## 3.1. Fracture-Closing Period

The calibrated microscale properties (Tables1 and 2 are used as input parameters in the model shown in Figure 2. Numerical modeling of fracture closing and fracturing fluid flowback is then calculated. Figures 9 and 10 plot the profiles of a proppant pillar at two different closing stresses, 20.7 MPa and 41.4 MPa, respectively. The height of the proppant pillar is 6 mm, and the diameter is 10 mm. With the increase of the closing stress, the height of the proppant pillar decreases gradually and the diameter increases. The pillar is compressed into a pancake shape. The compression from the shale plate changes the stress state in the proppant pillar and crushes down the pillar. Some particles near the edge of the pillar falls down and rolls away from the pillar.Figure 9
Profiles of the proppant pillar at 20.7 MPa closing stress: (a) lateral view and (b) top view.
(a)
(b)Figure 10
Profiles of the proppant pillar at 41.4 MPa closing stress: (a) lateral view and (b) top view.
(a)
(b)


## 3.2. Fracturing Fluid Flowback Period

The fluid flowback calculation with coupled DEM-CFD analysis was carried out after the fracture closing. Figure11 shows the profiles of the proppant pillar at different flow times. The initial height and diameter of the proppant pillar are 6 mm and 10 mm, respectively. The fluid viscosity is 0.001 Pa·s and the pressure difference is 64 kPa. The black arrows in the plot represent the fluid velocity vector at each fluid cell, and the velocity magnitude is indicated by the length of the arrow. The particles for rock plates are not shown in the plot. With the increase of flow time, the proppant pillar starts to spread and particles around the proppant pillar diffuse outward and form a single layer. After 0.0105 s, the proppant particles gradually escaped from the right side of the flow field and start to flow back. The number of flowback particles in this paper is defined as the number of proppant particles escaping from the flow field.Figure 11
Profiles of the proppant pillar at different flow times: (a) right after fracture closing, (b) flow time of 0.00015 s, (c) flow time of 0.0005 s, and (d) flow time of 0.0105 s. The red square in (b) represents the location of a monitoring fluid cell.
(a)
(b)
(c)
(d)We select the first row of fluid grids at the entry (shown in Figure11(b)) to monitor the change of flow velocity. Fluid velocities reach stable values at about 0.04 s flow time (Figure 12). Before 0.005 s, fluid velocities show fluctuation which indicates that the flow is subjected to unstable change of proppant ahead due to the spreading of particles. Between 0.025 s and 0.035 s, fluid velocities show more drastic fluctuation caused by fast flowback of a number of proppants. The results above indicate that a larger pillar diameter is prone to cause proppant particle spreading as well as unstable fluid flow. On the other hand, a smaller pillar diameter will lead to wider fluid channel and favor the stable fluid flow; however, the pillar-supporting capacity might not be met to resist the fracture-closing stress. Therefore, a proper size of proppant pillar is needed which can satisfy both the stability requirement and conductivity enhancement.Figure 12
Fluid velocity of the 10 fluid cells at the entry.


## 4. Discussions

### 4.1. Effect of Pressure Gradient

To investigate the effect of pressure gradient, a set of numerical simulations were run by gradually increasing the pressure difference between the two ends from 64 kPa to 512 kPa, while the height and diameter of the proppant pillar are kept as 6 mm and 10 mm, respectively, and the fluid viscosity is 0.001 Pa·s. Figure13 plots the evolution of the number of flowback proppant particles and proppant spreading area for different pressure gradients. The flowback process can be divided into two stages. At the first stage, the number of flowback proppant particles has an early increase but the growth rate decreases, while the difference of the flowback proppant particles among the four different fluid pressure gradients is minor. At the second stage, the difference of the number of flowback proppant particles gets larger and larger. Note that the difference of spreading area for different fluid gradients is minor. The reason is that the calculation of spreading area is based on the farthest particle in the domain. Figure 14 plots the evolution of the fluid velocity at the monitoring fluid cell (red cell in Figure 11(b)) for different pressure gradients. As the calculation starts, the fluid velocity gradually increases and reaches the plateau value which is proportional to the pressure difference. The difference of fluid velocity for the first stage of the flowback is small which accounts for the small difference of flowback proppant particles as shown in Figure 13.Figure 13
Evolution of the number of flowback proppant particles and proppant spreading area for different pressure gradients.Figure 14
Evolution of the fluid velocity at the monitoring fluid cell for different pressure gradients.Figure15 plots the evolution of the proppant pillar profile for different pressure gradients. In order to make better comparison for pillar profiles, proppants which run outside the flow field are not shown and similar figures are treated in the same way. For different pressure gradients, the profiles of proppant pillars are similar at the first stage. At the second stage, more proppants flow back in the high-pressure gradient case than in the low-pressure gradient case.Figure 15
Evolution of the proppant pillar profile for different pressure gradients: (a) pressure difference 64 kPa, (b) pressure difference 128 kPa, (c) pressure difference 256 kPa, and (d) pressure difference 512 kPa.
(a)
(b)
(c)
(d)


### 4.2. Effect of Fluid Viscosity

To study the effect of fracturing fluid viscosity on the flowback of proppant particles, the proppant pillar with 6 mm height and 10 mm diameter was run with different fluid viscosities, while the pressure difference between the inlet and outlet of the flow field was fixed as 64 kPa. Figure16 shows the evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities. It can be seen that when the pressure boundary condition is used, the amount of flowback particles and the spreading area decreases with the increase of fluid viscosity. According to (5), with the constant pressure gradient, increase of the fluid viscosity causes the decrease of the velocity difference vj¯−uj. However, from (3), the size of the βintj is difficult to determine (one part of (3) increases and the other decreases). Therefore, the force exerted on the particles by the fluid is difficult to determine. Figure 17 plots the evolution of fluid velocity for the monitoring fluid cell (blue cell in Figure 11(b)) for different fluid viscosities. It shows that when the viscosity of the fluid is large, the overall flow rate is small, which indirectly explains the phenomenon that the number of particles and the spreading area are smaller when the viscosity is larger. Figure 18 plots the evolution of the proppant pillar profile for different fluid viscosities. For lower viscosity cases, proppant pillars break down into separated packs, while for higher viscosity cases, proppant pillars are able to retain integrity during the flowback period.Figure 16
Evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosity.Figure 17
Evolution of the fluid velocity for the monitoring fluid cell for different fluid viscosities.Figure 18
Evolution of the proppant pillar profile for different fluid viscosities: (a) 0.001 Pa·s, (b) 0.01 Pa·s, and (c) 0.02 Pa·s.
(a)
(b)
(c)It should be noted that the trends in Figure16 will be reversed if velocity boundary is applied at the inlet. Figure 19 plots the evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities, with a fixed fluid velocity of 3.64 m/s applied at the inlet. The results show that larger fluid viscosity causes more flowback proppant particles and larger spreading area, opposite with the results in Figure 16. At the beginning of calculation, forces applied by the fluid on particles varied due to different particle sizes. Particle movement was further complicated by collision, and thus, average particle velocity might change dramatically. This kinematical disorder in turn caused fluctuation at the beginning of the blue line by particle-fluid interaction term in (2).Figure 19
Evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities, with velocity boundary applied at the inlet.Figure20 plots the evolution of the proppant pillar profile for different fluid viscosities, with velocity boundary applied at the inlet. Since the flowback process is different for the pressure boundary case and velocity boundary case, pillar profiles at different flow times are used. For the velocity boundary case, proppant pillars are more possible to break down when fluid viscosity is higher, which shows an exactly reverse law to the pressure boundary case.Figure 20
Evolution of the proppant pillar profile for different fluid viscosities, with velocity boundary applied at the inlet: (a) 0.001 Pa·s, (b) 0.01 Pa·s, and (c) 0.02 Pa·s.
(a)
(b)
(c)


### 4.3. Effect of Proppant Pillar Height

With the given proppant pillar diameter, an initially tall proppant pillar is expected to provide larger closing aperture and more conductive channels for the fluid flow after fracture closing than the initially short proppant pillar. However, the tall proppant pillar might also be subjected to the large amount of flowback proppant particles since during the fracture-closing process, particles at the edge of the proppant pillar are prone to be stripped off. Figure21 shows the proppant pillar profile after the flowback for two different proppant pillar heights, 6 mm and 8 mm. The case with the 8 mm pillar height has larger spreading area and more flowback proppant particles than the case with the 6 mm pillar height.Figure 21
Proppant pillar profile after the flowback: (a) proppant height 6 mm, before the flowback; (b) proppant height 6 mm, after the flowback of 0.00075 s; (c) proppant height 8 mm, before the flowback; and (d) proppant height 8 mm, after the flowback of 0.00075 s. The proppant diameter is 10 mm, the pressure difference is 64 kPa, and the viscosity is 0.001 Pa·s.
(a)
(b)
(c)
(d)The lab experiments by Gillard et al. [1] and Nguyen et al. [31] show that the conductivity of proppant pillars can be up to two orders of magnitude higher than the one with the traditional proppant placement scheme. The proppant pillar profile after the conductivity measurement (Figures 22 and 23) shows that the proppant particles keep bonded together during the flow testing process and the flowback of proppant particles cannot be clearly identified. Our modeling results however show that the flowback could occur with a large proppant pillar height or large fluid pressure gradient. As a result, the conductivities measured in the numerical models could be largely dropped due to the occupancy of particles in the fluid channels. The field tests with channel fracturing shows the increase of oil/gas production [3–5] compared with the traditional proppant placement scheme, but the improvement is much less remarkable than the experimental results by Gillard et al. [1] and Nguyen et al. [31].Figure 22
Proppant pillar profile after the conductivity measurement (Gillard et al. [1]).Figure 23
Proppant pillar profile after the conductivity measurement; four types of proppant were used (Nguyen et al. [31]).
(a)
(b)Due to the flowback of proppant particles, the proppant particles move into the channels between pillars, which results in the single layer or a few layers of proppant particles between pillars (Figure21). Part of a proppant pillar could still remain stable after flowback, and the initially open fluid channels might be filled with loosely packed proppant particles. The work by Fredd et al. [32] shows that the single-layer proppant in the fracture could enhance the fracture conductivity. Therefore, the fluid channels between pillars, though filled with a single layer or a few layers of proppant, may still effectively increase the fluid conductivity and make channel fracturing preferable than conventional-fracturing methods.


### 4.4. Effect of Pillar Diameter to Spacing Ratio

In channel fracturing, the fracturing fluid pulses, containing either the fracturing fluid without proppant or proppant-laden fluid, alternately pumped into the fractures, forming loose and discontinuous proppant clusters. Figure24 shows the two-dimensional schematic of high-conductivity fluid channels. The volumes of the fracturing fluid without proppant and the proppant-laden fluid can be simplified as two fluid columns in the plot while the volume ratio of the two fluids is related to the ratio of the pillar diameter and spacing by the following equation:
(7)D/L=d/x.Figure 24
Two-dimensional schematic of high-conductivity fluid channels.Among them,d is the diameter of the proppant pillar, x is the pillar spacing, L is the length of the fracturing fluid without proppant, and D is the length of proppant-laden fluid.The pillar diameter to spacing ratio (d/x) is therefore a key parameter for the design of channel-fracturing treatment. If the pillar diameter to spacing ratio is small, the pillar column might not be stable which causes small fracture-closing apertures and the fracture conductivity is reduced. On the contrary, large d/x leads to a narrow fluid channel and increasing amount of flowback particles. To optimize the value of d/x and determine the appropriate volume ratio of two fracturing fluids, simulation series by fixing the pillar spacing but varying the pillar diameter is carried out.Figure25 plots the percentage of flowback proppant particles for different d/x by using the pillar height of 6 mm and diameter of 10 mm. The fluid viscosity is 0.001 Pa·s, and the pressure difference is 64 kPa. There is nearly no flowback particles when the d/x ratio is smaller than 0.43. The percentage of flowback proppant particles reaches 57% when the d/x=0.70. If the d/x ratio is more than 0.70, the percentage of flowback proppant particles is close to 100%. The results suggest that the optimized d/x ratio is between 0.43 and 0.70, with the assigned input parameters in this case. However, it is worth mentioning that the analysis above is based on the assumption that the pillars are placed in arrays, while the real distribution of proppant pillars in field operation is expected to be less regular, given the fact that there could be multiple circumstances which cannot all be considered in the numerical models.Figure 25
Percentage of flowback proppant particles for different pillar diameter to spacing ratios.Figure26 plots the evolution of the proppant pillar profile for different pillar diameter to spacing ratios. With increasing diameter, proppant pillars tend to break down at earlier flow times and fill the flow field. There are almost no flowback proppants when d/x is smaller than 0.43 (Figure 26(a) and 26(b)), which can be deemed as well-controlled flowback conditions. Proppant flowback becomes a severe problem as d/x grows higher.Figure 26
Evolution of proppant pillar profile for different pillar diameter to spacing ratios: (a) 0.25, (b) 0.43, (c) 67, (d) 0.70, (e) 0.80, (f) 0.90, (g) 1.00, (h) 1.50, (i) 2.33, and (j) 4.00.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


### 4.5. Effect of Fiber-Bonding Strength

The effect of fiber-bonding strength on the stability of the proppant pillar was also studied. Figure27 shows the evolution of the number of flowback proppant particles and proppant spreading area for different fiber-bonding strengths, with the proppant pillar height of 6 mm and the diameter of 10 mm. The fluid viscosity is 0.001 Pa·s, and the pressure difference is 64 kPa. The greater the bond strength, the stronger the cohesive effect of the simulated fibers on the proppant particles. With the enhancement of cohesion, the number of flowback particles from the proppant pillar is reduced and the spreading area is also decreased. Figure 28 plots the evolution of the proppant pillar profile for different fiber-bonding strengths. Comparing with bonding strengths of 3E + 3 and 3E + 6 Pa, the proppant pillar remains quite stable during the flowback period when the bonding strength is 3E + 7 Pa. It can be seen that reasonably adjusting the physicochemical properties of the fibers increases the cohesion between the proppant particles and helps to enhance the stability of the proppant pillar in the high-permeability fracture.Figure 27
Evolution of the number of flowback proppant particles and proppant spreading area for different fiber-bonding strength.Figure 28
Evolution of proppant pillar profile for different fiber-bonding strengths: (a) 3E + 3 Pa, (b) 3E + 6 Pa, and (c) 3E + 7 Pa.
(a)
(b)
(c)


## 4.1. Effect of Pressure Gradient

To investigate the effect of pressure gradient, a set of numerical simulations were run by gradually increasing the pressure difference between the two ends from 64 kPa to 512 kPa, while the height and diameter of the proppant pillar are kept as 6 mm and 10 mm, respectively, and the fluid viscosity is 0.001 Pa·s. Figure13 plots the evolution of the number of flowback proppant particles and proppant spreading area for different pressure gradients. The flowback process can be divided into two stages. At the first stage, the number of flowback proppant particles has an early increase but the growth rate decreases, while the difference of the flowback proppant particles among the four different fluid pressure gradients is minor. At the second stage, the difference of the number of flowback proppant particles gets larger and larger. Note that the difference of spreading area for different fluid gradients is minor. The reason is that the calculation of spreading area is based on the farthest particle in the domain. Figure 14 plots the evolution of the fluid velocity at the monitoring fluid cell (red cell in Figure 11(b)) for different pressure gradients. As the calculation starts, the fluid velocity gradually increases and reaches the plateau value which is proportional to the pressure difference. The difference of fluid velocity for the first stage of the flowback is small which accounts for the small difference of flowback proppant particles as shown in Figure 13.Figure 13
Evolution of the number of flowback proppant particles and proppant spreading area for different pressure gradients.Figure 14
Evolution of the fluid velocity at the monitoring fluid cell for different pressure gradients.Figure15 plots the evolution of the proppant pillar profile for different pressure gradients. In order to make better comparison for pillar profiles, proppants which run outside the flow field are not shown and similar figures are treated in the same way. For different pressure gradients, the profiles of proppant pillars are similar at the first stage. At the second stage, more proppants flow back in the high-pressure gradient case than in the low-pressure gradient case.Figure 15
Evolution of the proppant pillar profile for different pressure gradients: (a) pressure difference 64 kPa, (b) pressure difference 128 kPa, (c) pressure difference 256 kPa, and (d) pressure difference 512 kPa.
(a)
(b)
(c)
(d)


## 4.2. Effect of Fluid Viscosity

To study the effect of fracturing fluid viscosity on the flowback of proppant particles, the proppant pillar with 6 mm height and 10 mm diameter was run with different fluid viscosities, while the pressure difference between the inlet and outlet of the flow field was fixed as 64 kPa. Figure16 shows the evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities. It can be seen that when the pressure boundary condition is used, the amount of flowback particles and the spreading area decreases with the increase of fluid viscosity. According to (5), with the constant pressure gradient, increase of the fluid viscosity causes the decrease of the velocity difference vj¯−uj. However, from (3), the size of the βintj is difficult to determine (one part of (3) increases and the other decreases). Therefore, the force exerted on the particles by the fluid is difficult to determine. Figure 17 plots the evolution of fluid velocity for the monitoring fluid cell (blue cell in Figure 11(b)) for different fluid viscosities. It shows that when the viscosity of the fluid is large, the overall flow rate is small, which indirectly explains the phenomenon that the number of particles and the spreading area are smaller when the viscosity is larger. Figure 18 plots the evolution of the proppant pillar profile for different fluid viscosities. For lower viscosity cases, proppant pillars break down into separated packs, while for higher viscosity cases, proppant pillars are able to retain integrity during the flowback period.Figure 16
Evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosity.Figure 17
Evolution of the fluid velocity for the monitoring fluid cell for different fluid viscosities.Figure 18
Evolution of the proppant pillar profile for different fluid viscosities: (a) 0.001 Pa·s, (b) 0.01 Pa·s, and (c) 0.02 Pa·s.
(a)
(b)
(c)It should be noted that the trends in Figure16 will be reversed if velocity boundary is applied at the inlet. Figure 19 plots the evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities, with a fixed fluid velocity of 3.64 m/s applied at the inlet. The results show that larger fluid viscosity causes more flowback proppant particles and larger spreading area, opposite with the results in Figure 16. At the beginning of calculation, forces applied by the fluid on particles varied due to different particle sizes. Particle movement was further complicated by collision, and thus, average particle velocity might change dramatically. This kinematical disorder in turn caused fluctuation at the beginning of the blue line by particle-fluid interaction term in (2).Figure 19
Evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities, with velocity boundary applied at the inlet.Figure20 plots the evolution of the proppant pillar profile for different fluid viscosities, with velocity boundary applied at the inlet. Since the flowback process is different for the pressure boundary case and velocity boundary case, pillar profiles at different flow times are used. For the velocity boundary case, proppant pillars are more possible to break down when fluid viscosity is higher, which shows an exactly reverse law to the pressure boundary case.Figure 20
Evolution of the proppant pillar profile for different fluid viscosities, with velocity boundary applied at the inlet: (a) 0.001 Pa·s, (b) 0.01 Pa·s, and (c) 0.02 Pa·s.
(a)
(b)
(c)


## 4.3. Effect of Proppant Pillar Height

With the given proppant pillar diameter, an initially tall proppant pillar is expected to provide larger closing aperture and more conductive channels for the fluid flow after fracture closing than the initially short proppant pillar. However, the tall proppant pillar might also be subjected to the large amount of flowback proppant particles since during the fracture-closing process, particles at the edge of the proppant pillar are prone to be stripped off. Figure21 shows the proppant pillar profile after the flowback for two different proppant pillar heights, 6 mm and 8 mm. The case with the 8 mm pillar height has larger spreading area and more flowback proppant particles than the case with the 6 mm pillar height.Figure 21
Proppant pillar profile after the flowback: (a) proppant height 6 mm, before the flowback; (b) proppant height 6 mm, after the flowback of 0.00075 s; (c) proppant height 8 mm, before the flowback; and (d) proppant height 8 mm, after the flowback of 0.00075 s. The proppant diameter is 10 mm, the pressure difference is 64 kPa, and the viscosity is 0.001 Pa·s.
(a)
(b)
(c)
(d)The lab experiments by Gillard et al. [1] and Nguyen et al. [31] show that the conductivity of proppant pillars can be up to two orders of magnitude higher than the one with the traditional proppant placement scheme. The proppant pillar profile after the conductivity measurement (Figures 22 and 23) shows that the proppant particles keep bonded together during the flow testing process and the flowback of proppant particles cannot be clearly identified. Our modeling results however show that the flowback could occur with a large proppant pillar height or large fluid pressure gradient. As a result, the conductivities measured in the numerical models could be largely dropped due to the occupancy of particles in the fluid channels. The field tests with channel fracturing shows the increase of oil/gas production [3–5] compared with the traditional proppant placement scheme, but the improvement is much less remarkable than the experimental results by Gillard et al. [1] and Nguyen et al. [31].Figure 22
Proppant pillar profile after the conductivity measurement (Gillard et al. [1]).Figure 23
Proppant pillar profile after the conductivity measurement; four types of proppant were used (Nguyen et al. [31]).
(a)
(b)Due to the flowback of proppant particles, the proppant particles move into the channels between pillars, which results in the single layer or a few layers of proppant particles between pillars (Figure21). Part of a proppant pillar could still remain stable after flowback, and the initially open fluid channels might be filled with loosely packed proppant particles. The work by Fredd et al. [32] shows that the single-layer proppant in the fracture could enhance the fracture conductivity. Therefore, the fluid channels between pillars, though filled with a single layer or a few layers of proppant, may still effectively increase the fluid conductivity and make channel fracturing preferable than conventional-fracturing methods.


## 4.4. Effect of Pillar Diameter to Spacing Ratio

In channel fracturing, the fracturing fluid pulses, containing either the fracturing fluid without proppant or proppant-laden fluid, alternately pumped into the fractures, forming loose and discontinuous proppant clusters. Figure24 shows the two-dimensional schematic of high-conductivity fluid channels. The volumes of the fracturing fluid without proppant and the proppant-laden fluid can be simplified as two fluid columns in the plot while the volume ratio of the two fluids is related to the ratio of the pillar diameter and spacing by the following equation:
(7)D/L=d/x.Figure 24
Two-dimensional schematic of high-conductivity fluid channels.Among them,d is the diameter of the proppant pillar, x is the pillar spacing, L is the length of the fracturing fluid without proppant, and D is the length of proppant-laden fluid.The pillar diameter to spacing ratio (d/x) is therefore a key parameter for the design of channel-fracturing treatment. If the pillar diameter to spacing ratio is small, the pillar column might not be stable which causes small fracture-closing apertures and the fracture conductivity is reduced. On the contrary, large d/x leads to a narrow fluid channel and increasing amount of flowback particles. To optimize the value of d/x and determine the appropriate volume ratio of two fracturing fluids, simulation series by fixing the pillar spacing but varying the pillar diameter is carried out.Figure25 plots the percentage of flowback proppant particles for different d/x by using the pillar height of 6 mm and diameter of 10 mm. The fluid viscosity is 0.001 Pa·s, and the pressure difference is 64 kPa. There is nearly no flowback particles when the d/x ratio is smaller than 0.43. The percentage of flowback proppant particles reaches 57% when the d/x=0.70. If the d/x ratio is more than 0.70, the percentage of flowback proppant particles is close to 100%. The results suggest that the optimized d/x ratio is between 0.43 and 0.70, with the assigned input parameters in this case. However, it is worth mentioning that the analysis above is based on the assumption that the pillars are placed in arrays, while the real distribution of proppant pillars in field operation is expected to be less regular, given the fact that there could be multiple circumstances which cannot all be considered in the numerical models.Figure 25
Percentage of flowback proppant particles for different pillar diameter to spacing ratios.Figure26 plots the evolution of the proppant pillar profile for different pillar diameter to spacing ratios. With increasing diameter, proppant pillars tend to break down at earlier flow times and fill the flow field. There are almost no flowback proppants when d/x is smaller than 0.43 (Figure 26(a) and 26(b)), which can be deemed as well-controlled flowback conditions. Proppant flowback becomes a severe problem as d/x grows higher.Figure 26
Evolution of proppant pillar profile for different pillar diameter to spacing ratios: (a) 0.25, (b) 0.43, (c) 67, (d) 0.70, (e) 0.80, (f) 0.90, (g) 1.00, (h) 1.50, (i) 2.33, and (j) 4.00.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


## 4.5. Effect of Fiber-Bonding Strength

The effect of fiber-bonding strength on the stability of the proppant pillar was also studied. Figure27 shows the evolution of the number of flowback proppant particles and proppant spreading area for different fiber-bonding strengths, with the proppant pillar height of 6 mm and the diameter of 10 mm. The fluid viscosity is 0.001 Pa·s, and the pressure difference is 64 kPa. The greater the bond strength, the stronger the cohesive effect of the simulated fibers on the proppant particles. With the enhancement of cohesion, the number of flowback particles from the proppant pillar is reduced and the spreading area is also decreased. Figure 28 plots the evolution of the proppant pillar profile for different fiber-bonding strengths. Comparing with bonding strengths of 3E + 3 and 3E + 6 Pa, the proppant pillar remains quite stable during the flowback period when the bonding strength is 3E + 7 Pa. It can be seen that reasonably adjusting the physicochemical properties of the fibers increases the cohesion between the proppant particles and helps to enhance the stability of the proppant pillar in the high-permeability fracture.Figure 27
Evolution of the number of flowback proppant particles and proppant spreading area for different fiber-bonding strength.Figure 28
Evolution of proppant pillar profile for different fiber-bonding strengths: (a) 3E + 3 Pa, (b) 3E + 6 Pa, and (c) 3E + 7 Pa.
(a)
(b)
(c)


## 5. Conclusions

In this study, proppant pillar deformation and stability during the flowback of channel fracturing was simulated with the DEM-CFD-coupling method. The effect of fibers was considered by implementing the bonded particle model for contacts between particles. The major conclusions of this study can be summarized as follows.In the fracturing-closing period, the height of the proppant pillar decreases gradually and the diameter increases as the closing stress increases. In the flowback period, proppant particles could be driven away from the pillar by the fluid flow and cause the instability of proppant pillar.The flowback could occur easily with a large proppant pillar height or large fluid pressure gradient. Increasing the initial proppant pillar height gives more conductive channels for the fluid flow after fracture closing, but on the other hand, a tall proppant pillar might also be subjected to the large amount of flowback proppant particles.Besides the proppant height, the pillar diameter to spacing ratio is another key parameter for the design of channel fracturing. A small pillar diameter to spacing ratio causes the instability of the proppant pillar and small fracture-closing aperture; on the contrary, a large pillar diameter to spacing ratio leads to narrow fluid channels and increasing amounts of flowback particles.The profile of proppant placement inside the fracture in the field should be somehow between the ideal competent proppant pillar arrangement and the uniform layer of the proppant based on the traditional proppant placement scheme. Adjusting the physicochemical properties of the fibers increases the cohesion between the proppant particles and helps to enhance the stability of the proppant pillar in the high-permeability fracture.



---

*Source: 3535817-2018-08-02.xml*