Hydraulic fracturing experiments with low-viscosity fluids, such as supercritical CO2, demonstrate the formation of complex fracture networks spread throughout the rocks. To study the influence of viscosity of the fracturing fluids on hydraulic fracture propagation, a hydromechanical-coupled cohesive zone model is proposed for the simulation of mechanical response of rock grains boundary separation. This simulation methodology considers the synergistic effects of unsteady flow in fracture and rock grain deformation induced by hydraulic pressure. The simulation results indicate a tendency of complex fracture propagation with more branches as the viscosity of fracturing fluids decrease, which is in accord with experimental results. The low-viscosity fluid can flow into the microfractures with extremely small aperture and create more shear failed fracture. This study confirms the possibility of effective well stimulations by hydraulic fracturing with low-viscosity fluids.

1. Introduction

Hydraulic fracturing has been one of the primary engineering tools for improving well productivity especially for unconventional reservoirs, such as shale or tight sand reservoirs [1]. CO2 has been considered to be a new king of fracturing fluid, which can be alternatives to the conventional aqueous fracturing fluid, due to its advantages of miscibility with hydrocarbons and enhanced fracturing properties [2]. Hydraulic fracturing experiments with supercritical CO2 have been conducted by some research institutions [38]. These experiment results provided valuable information to predict the initiation and propagation of artificial fracture in hydrofracturing with supercritical CO2. Compared to water, fracture induced by supercritical CO2 has a large number of branches (see Figure 1 [7]). Such complex fracture networks with many microfractures make the effective development of extremely low-permeability unconventional reservoirs possible. So far, the mechanism for the formation of complex fracture networks induced by hydrofracturing with low-viscosity supercritical CO2 is still not well studied.

In the past few decades, many numerical methods were proposed to investigate the hydraulic fracturing and fluid flow under micron-scale [912]. There are lots of numerical methods developed to investigate the hydraulic fracturing, such as boundary element method (BEM) [13, 14], discrete element method (DEM) [1517], and finite element method (FEM) [1820]. BEM is a very efficient method to simulate fracture propagation in rock formation as only the boundary of the domain is discretised, but it is restricted to homogeneous isotropic linear elastic rock and simple fracture geometry [21]. DEM treats the rock mass as an assembly of separate particles. The advantage of DEM is its ability to simulate relatively complex fracture geometry, and the disadvantages lie in the time-consuming calibration process and the restriction of triangular elements [22]. In the FEM community, cohesive zone models (CZM) and the extended finite element method (XFEM) can be cited. XFEM allows a fracture diverse to a new direction [23], but it is very difficult to using this method to solve multifracture problems. CZM is a domain fracture model in the FEM to study the mechanical responses at material interfaces such as grain boundary separation and crack growth [24, 25]. It removes the crack tip singularity of linear elastic fracture mechanics and can fit naturally into the FEM. Due to its high efficiency and easy implementation, CZM method is widely applied in hydraulic fracturing simulations to quantitatively analyse fracture behaviour in formations [2628]. Most of these methods only solve problems with very simply straight crack, and the path of fracturing propagation is given in advance. So far, few methods have capabilities to simulate the propagation of complex fractures with numerous branches. Thus, it is necessary to develop a numerical method to study the complex fracture network hydraulic fracturing process.

The objective of this study is to introduce a finite element simulation framework to study the complex hydraulic fracture propagation in reservoirs. The remainder of this paper is organized as follows. In Section 2, the numerical methods will be introduced, which includes a modified cohesive zone model and the mechanism of an unsteady fluid flow method in the fracture path. In Section 3, the numerical implementation will be presented. In Section 4, the hydrofracturing with different viscosities fluids was investigated by the hydromechanical-coupled simulation. Summary will be given in the last section.

2. Numerical Method

A fluid flow-dependent cohesive zone model is proposed and implemented in ABAQUS User Subroutine (UEL) [29].

2.1. Cohesive Interface Model

A large number of cohesive zone models and various mechanisms have been developed to study the separation between two joining interfaces [30]. The constitutive equation relating the traction vector (Tn, Tt) to relative displacement (Δn, Δt) is the most conveniently characterized by a scalar interplanar potential ΦΔn,Δt by setting
The value of Φ represents the work done per unit area in separating the interface. Subscripts n,t denote the normal and tangential directions, respectively. In this paper, the function of the separation potential Φ is developed by Xu and Needleman [31]:
where δn and δt are characteristic lengths. Under normal loading, the interface has a separation work фn, and the normal traction stress reaches a maximum valueσmax=ϕn/δnexp1 at an interface separation Δn=δn. Under purely shear loading, the interface has work of separation фt, and the tangential traction reaches a maximum value τmax=ϕt2/δtexp1 atΔt=δt/2. For Mode I loading condition in plane strain, the fracture energy is GIc=ϕn, and fracture toughness KIc=EGIc/1v2.
The cohesive interface simulations are often limited by the occurrence of an elastic snap-back instability in the Newton-Raphson iterations. Such numerical convergence problem can be solved by introducing a small viscosity in the cohesive constitutive equations [32].
where ζn and ζt are viscosity parameters that govern viscous energy dissipation under normal and tangential loading.

The traction-displacement relations are plotted in Figure 2.

2.2. Fluid Flow in the Fractures

In the real fracturing process, the fluid flow is in an unsteady state in the flow path. It is necessary to develop an unsteady state solution method suitable for the coupling between fluid and solid.

For laminar flow, the basic equations for flow from fracture element j to i is
where w is the equivalent aperture of the fracture element, L is the length of the flow channel element, and μ is the viscosity of injection fluid; i and j are element numbers.
For the leak-off flow from element i, the equation is defined as
where cl is the fluid leak-off coefficients and pr denotes the pore pressures in the adjacent rock.

In our method, there will be at most four other fracture elements connecting on one fracture element due to the existence of fracture junctions (see Figure 3).

The fluid volume changes in one fracture element equal to the sum of the inflow, leaf off flow from the fracture surface to rock, and fracture volume change. Thus, the expression of fluid volume changes in fracture element i, in a t step, is given by
Based on the volume changes of fluid, the pressure changes in the fracture element i (connecting to at least one fracture element) are given by
where Kf is the bulk modulus of injection fluid.
The time step size (Δt) is chosen according to the following equation:
where n is the number of the fracture elements in the model. The time step value of Equation (9) ensures that in most cases, only a small fraction of a volume of a fracture element is filled with inflow fluid in a time step.

2.3. Numerical Implementation

The rock-like material can be simulated as a dense assembly of deformable grains that interact at their contact points, and the grain boundaries can represent flaws in rock and allow for numerical replication of fracture propagation along grain boundaries. A two-dimensional (2D) Voronoi-grain model with cohesive elements that represents the reservoir rock is created, and the model generation procedures are shown in Figure 4. The formed Voronoi polygonal grains are meshed with quadrangular shapes. The advantage of the Voronoi random polygonal grain model is that it shows a wide variety of particle shapes and volumes, which guarantees that the macroscopic behaviour of rock is not biased by particle distribution [33]. Grain boundaries are highlighted by red solid lines consisting of four node 2D user-defined elements (UEL) where the proposed cohesive zone model and fluid flow are incorporated. To illustrate how user-defined cohesive elements are coupled with neighbouring grains, a detailed sketch of the grain boundary is shown in Figure 5. The initial width of cohesive elements is zero, so that the top and bottom cohesive faces share the same coordinates. In addition, a sketch of junction point at grain boundaries is presented in Figure 6. Note the width of these cohesive elements are exaggerated in Figures 5 and 6.

The failed cohesive elements (UEL) are the flow network path (or fracture) for fluid flow. To solve the fluid pressure, a fracture judging algorithm is proposed to initialize and update the flow network in the coupled hydromechanical process. At the beginning of simulation, all cohesive elements are not considered as the flow network except the injection element. The judging algorithm will read the positions of them and determine their connecting status (as shown in Figure 2). During the simulation, flow network is constantly updated at each increment. Only when the two conditions are satisfied can a cohesive element be added to the flow network: (1) this cohesive element is failed; (2) at least one of its connecting cohesive elements is failed.

In each time increment, the UEL subroutine is called to return the assembled stiffness matrix. Given Equation (7) and Equation (8), in order to collect necessary values for the gradient term, the COMMON BLOCK technique is adopted for data storage and acquiring. Values of pressure and width in both tn1 and tn are stored in four COMMON BLOCKs. These data are identified with UEL element number. All UEL cohesive elements are initialized based on boundary condition in current COMMON BLOCK sets, P_CURRENT and W_ CURRENT. At the next increment, the data stored in current COMMON BLOCK sets are transferred to old COMMON BLOCK sets, P_OLD and W_ OLD. Then, the pressure data in each element is updated by solving Equation (7) based on the corresponding data in old COMMON BLOCK sets, and the width is updated by accumulating cohesive element nodes displacement. The UEL element can be called multiple times during every time increments. To distinguish new increment from old increment with the second call, a signal flag variable is adopted to differentiating increments.

Figure 7 shows the computational scheme of the coupled procedure between the flow solver and the solid solver. The approach is fully coupled in the sense that the mechanical solving in solid module is affected by the fluid pressures, and the flow calculations in flow solver are affected by the change of fracture width and length due to the rock deformation and fracturing process. The flow solver obtains the shape of the flow path from the node displacements. The calculated fluid pressure is then taken as a load acting on the fracture up and bottom surfaces. With the surface load, the solid solver can update the node displacement and judge whether the failure of cohesive elements happens. By repeating the above procedures between flow and solid solvers, a hydramechanical coupling analysis by Abaqus is achieved.

3. Results and Discussion

To investigate influences of the low viscosity of supercritical CO2 on fracture geometry, a 25m×25m anisotropic model is established. There are two kinds of boundary conditions in this model. First, the edges of the model are applied in situ stresses, maximum and minimum horizontal principal stresses, σH and σh. Second, the hydraulic flow is continuing pumped into a fracture element (point) in the centre of the model. The boundary condition is shown in Figure 8.

The parameters of rock and fluid come from the experiment date in Zhang et al.’s work [7]. Two types of fracturing fluids is adopted, water at room temperature and supercritical CO2. The viscosity of supercritical CO2 is around 5% of that of water. Detailed mechanical and hydraulic parameters in simulation are listed in Table 1.

Hydraulic fracture propagation in rocks with water at room temperature and supercritical CO2 is shown in Figure 9. Noted that the width of fractures is amplificated. The fracture geometry agrees well with the results of experiments [7]. The injection of supercritical CO2 leads to more complex fracture pattern (see Figure 9(b)), and only long main fractures with simple geometry are observed (see Figure 9(a)). The pressure distribution along fracture paths is also presented in Figure 9. In hydraulic fracturing with water, the pressure falls rapidly from injection point to crack tips. In comparison, high pressure spread widely along all the microfracture.

The crack failure type is also investigated. As shown in Figure 10, the tensile fail cracks are dominated in hydraulic fracturing with water, but shear fail cracks increase significantly in hydraulic fracturing with supercritical CO2. The numbers of tensile and shear failed cracks are presented in Figure 11. The total number of failed cracks in hydraulic fracturing with supercritical CO2 is much more than hydraulic fracturing with water (246 to 101), which corresponds to the complex geometry and lots of microfractures. In addition, the number of shear failed cracks in hydraulic fracturing with supercritical CO2 is 87, which account for 35% of total failed cracks.

With super low viscosity, the supercritical CO2 can flow fast along the narrow fracture paths and rise the pressure in the microfractures and exert surface load on the fracture surfaces. The widely spread high pressure can press the rock grains easily and then induce shear failure of cracks. That is the main reason for forming complex branched fractures induced by supercritical CO2 fracturing.

4. Conclusions

This study proposes a finite element framework for simulating hydraulic fracture with a complex crack network. An unsteady flow solution is embedded in cohesive zone model and implemented in ABAQUS UEL subroutine. Then, the hydromechanical-coupled method is applied in a Voronoi grain-based model, and the fracture pattern in hydraulic fracturing with water and supercritical CO2 is investigated. The simulation results reproduce the same results of previous experimental results [7]. With supercritical CO2, high pressure can spread widely along microfractures due to the super low viscosity. The high pressure in microfractures may extrude rock grains and induce lots of shear failed cracks. The shear cracks make fractures more branched and promote to form complex fracture networks. Our research provides a better understanding of fracture propagation in hydraulic fracturing with low-viscosity fluid. In the present research, the seepage effect is not considered and will be studied in future work.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflict of interest regarding this publication.


This work is funded by the National Science and Technology Major Project of China (2017ZX05009005).

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.