Using Cohesive Zone Model to Simulate the Hydraulic Fracture Interaction with Natural Fracture in Poro-Viscoelastic Formation
Hydraulic fracturing is a widely used production stimulation technology for conventional and unconventional reservoirs. The cohesive element is used to explain the tip fracture process. In this paper, the cohesive zone model was used to simulate hydraulic fracture initiation and propagation at the same time rock deformation and fluid exchange. A numerical model for fracture propagation in poro-viscoelastic formation is considered. In this numerical model, we incorporate the pore-pressure effect by coupling fluid diffusion with shale matrix viscoelasticity.
Yu Suo1, *, Zhixi Chen1, Hao Yan1, 2, Daobing Wang3 and Yun Zhang
1Minerals and Energy Resources Engineering, University of New South Wales, Sydney 2052, Australia. 2State Key Laboratory of Coal Resources and Safe Mining, China University of Mining & Technology, Xuzhou 221116, China. 3School of Mechanical Engineering, Beijing Key Laboratory of Pipeline Critical Technology and Equipment for Deepwater Oil & Gas Development, Beijing Institute of Petrochemical Technology, Beijing 102617, China
Received: 18 January 2019 / Accepted: 26 March 2019 / Published: 1 April 2019
The numerical procedure for hydraulically driven fracture propagation uses a poro-viscoelasticity theory to describe the fluid diffusion and matrix creep in the solid skeleton, in conjunction with pore-pressure cohesive zone model and ABAQUS was used as a platform for the numerical simulation. The simulation results are compared with the available solutions in the literature. The higher the approaching angle, the higher the differential stress, tensile stress difference, injection rate, and injection fluid viscosity, and it will be easier for hydraulic fracture crossing natural fracture. These results could provide theoretical guidance for predicting the generation of fracture network and gain a better understanding of deformational behavior of shale when fracturing.
Unconventional reservoirs, such as shale, coal seam and tight-gas sand, are highly reliant on hydraulic fracturing to increase the production. Although it is difficult to exploit, these low-permeability reservoirs contain considerable hydrocarbon resources such as Marcellus Shale , Barnett Shale  and Eagle Ford Shale . The purpose of hydraulic fracturing is to open natural fracture and increase the number of flow paths, to connect the natural fractures in the reservoir with the hydraulic fractures, opening the existing natural fractures, and increasing the path volume in the fractures.
Natural fracture plays an important role in hydraulic fracturing such as increasing shale gas flow path, enhancing productivity, and reducing propped fracture length . Natural fractures have an important influence in hydraulic fracturing. It cannot only affect the propagation of fracturing in reservoirs, but it also affects the efficiency of hydraulic fracturing treatment. The benefit of natural fractures is that it can expand the flow path volume and increase the production of hydrocarbons. On the negative side, it may cause additional resources to leak. Therefore, understanding of how the hydraulic fractures interact with natural fractures is necessary and is beneficial to our future simulations and improvements.
Overall, continuous research on hydraulic fracturing is necessary and worthwhile. Hydraulic fracture intersection with natural fracture was investigated through both laboratory and mine-back experimental approaches [4,5,6,7]. When hydraulic fractures intersect with natural fractures, there are three different developments. Hydraulic fractures cross over natural fractures without changing the volume or shape of natural fractures; hydraulic fractures divert into the natural fracture and expand the volume of the flow path; hydraulic fractures enter natural fractures and form a complex fracture network structure.
Biot first proposed viscoelastic theory to study the linear viscoelasticity and anisotropy of porous media . The influence of the permeability on the deformation and the settlement problem of the loaded column was discussed. A micromechanical method for explaining the theory of pore viscoelasticity of Biot is proposed by Abousleiman et al. . The problem of drilling under plane-strain deformation is discussed and the effects of combined poro- and visco-elasticity are investigated.
Simakin and Ghassemi  put forward another poro-viscoelastic model by considering the relaxation of the deviatoric stress and the symmetric effective stress. These constitutive models can be used to simulate the coupling process between fluid flow and creep deformation of matrix rocks, but they cannot be directly applied to simulate hydraulic fracturing, especially the interaction of hydraulic fracturing with natural fractures. Hydraulic fracturing is a complex process and it contains rock deformation, fracturing fluid flow, and fracture imitation and propagation.
Linear Elastic Fracture Mechanics (LEFM), which is known as a method to describe fracture mechanics in brittle rock, made outstanding contributions. However, for soft rock, LEFM has no advantage. Under formation temperature and pressure, the brittle rock in shale exhibits ductility after hydraulic fracturing treatments. Therefore, it may not be the best choice to use elastic theory to analyze the shale failure.
The cohesive zone model (CZM) is a model in fracture mechanics in which fracture formation is considered to be a gradual phenomenon. The separation of the surface involved in the crack occurs on the extended crack tip or cohesive region and resists traction by bonding. By using CZM, the crack initiation and propagation in many kinds of materials could be predicted and it makes a detailed description of fracture process zone become possible . CZM assumed the fracture process zone may be separated by traction-separation law (TSL), which is determined by the material of the fracture treatment zone.
The presence of cohesive traction eliminates the singularity of the crack tip, and this singularity significantly affects the convergence of the solution. The standard model for describing the fracture tip processing region is the assumed bond portion extending between the fracture surfaces until they break at an appropriate level of traction [12,13]. CZM has great advantages in predicting crack propagation orientation in different kinds of materials (such as metals, concrete, and rocks ). Furthermore, it has been used for simulating hydraulic fracture interactions with natural fracture [15,16,17].
In this paper, using this method to simulate the hydraulic fracturing from initiation and propagation, the rock deformation and fluid exchange between porous infiltrated media and fractures are being coupled. The cohesive element is used to explain the tip fracture process. The numerical procedure for hydraulically driven fracture propagation uses a poro-viscoelasticity theory to describe the fluid diffusion and matrix creep in the solid skeleton. Moreover, pore-pressure CZM and ABAQUS was used as a platform for the numerical simulation.
When the rock material is subjected to an external load greater than the long-term strength, and the external stress conditions are constant, the rock deformation gradually increases with time. The rock incurs creep when the internal strain reaches the accelerated creep threshold, and the strain will eventually reach the strength of the fracture and finally break. It is generally considered that the rock has a long-term effect on the external constant load creep phenomenon, so the viscoelastic model is suitable.
The How law is used to model the stress and strain for an elastic material. Nevertheless, viscoelastic materials generally show elasticity and viscous deformation over time. The rock deformation and fluid diffusion are involved in the physical process and influence each other. The poro-viscoelasticity theory has been developed and used for the fluid-rock coupling effect . In this paper, the Maxwell model will be used to illustrate the viscoelasticity shown in Equation (1):
Poro-viscoelasticity is the synthesis of poro-elasticity and viscoelasticity. Biot  pioneered the approach by using the correspondence principle, which allows for straightforward conversion from elastic to viscoelastic behavior. The constitutive equations for a poro-elastic rock are regarded as:
where K is rock bulk Young’s modulus, G is shear Young’s modulus, ɑ is the Biot’s effective stress coefficient which measures the ability of the pore pressure to resist compressional stresses, M is the Biot’s Modulus. The corresponding total stress components added to pore pressure of the fluid (p) value is the components of effective stress tensor.
The difference between the constitutive equation of poro-viscoelasticity and that of poro-viscoelasticity is the viscous term (see right-hand side of Equation (1)). The deviatoric and symmetric components of stress act on the time developed by the viscous behavior, by adding these two elements using the finite method. A fully coupled hydro-mechanical process, the increment of total stress at the current stage is showed as:
where σ(vis)ij is the pseudo-stress tensor which denotes the viscous influence of the strain history. Besides the above equation, a fluid diffusion equation should be used,
where k is the permeability. Equations (5) and (6) show a combination for poro-viscoelasticity. They can be used to simulate fluid flow under plane-strain conditions as well as the interaction between viscoelastic deformation of rock mass and fluid flow during mechanical processes.
In the ABAQUS platform, it assumes that the time domain viscoelasticity is determined by a Prony series expansion. In the Prony series the parameters are defined for each term. On the other hand, ABAQUS can calculate the terms in the Prony series using time-dependent creep test data, time-dependent relaxation test data, or frequency-dependent cyclic test data. The parameters used in Prony series are acquired from nonlinear regression based on power-law equations acquired from experimental data .
Compare Poro-Elastic Model to Poro-Viscoelastic Model
In numerical simulation, a 2D plane-strain crack propagation problem is considered. The formation is regarded as homogeneously poro-viscoelastic. Only half of the model is meshed due to its symmetry. The initial perforated crack is located at the center, while the predefined crack path is located on the x-axis. The model is subjected to far-field in situ stress, as shown in Figure 1.
Figure 1. 2D crack in plane-strain condition.
To investigate the effect of poro-viscoelastic formation on fracture propagation, results from poro-viscoelastic formation are compared with poro-elastic formation. The input parameters are the same except for cohesive strength and viscosity. The parameters used for numerical simulations are presented in Table 1. Figure 2 presents fracture propagation geometry in poro-elastic and poro-viscoelastic formations, respectively.
Poro-viscoelastic model fracture width is larger than poro-elastic, while the half fracture height is smaller than poro-elastic. Since the hydraulic energy is partly dissipated during the subsequent ongoing viscous deformation, more time is needed to build up the required input to initiate fracture growth. In the meantime, the fracture is widened due to decaying modulus of the bulk. The model shows that the poro-viscoelastic is differently poro-elastic so this principle for shale needs to be considered in our simulation.
Figure 2. Fracture geometry with time.
Table 1. All the parameters used in numerical simulation.
Cohesive Zone Model
The interaction between hydraulic fracture and natural fracture from initiation to the propagation was modeled in the CZM by introducing a layer of elements (possible path) with zero thickness. Fully modeling the hydraulic fracturing process requires solving a coupled system of governing equations [20,21]. A well is drilled and then a hydraulic fracture with injection fluid enters the reservoir. The fracture process zone (unbroken cohesive zone) is defined within the upper and lower surfaces in Figure 3. The fracturing fluid is filled with fracture where no traction from shale fracture exists.
The bilinear model is used in the CZM. Three stages can be observed in the process of opening fracture. (1) At first, the opening increases at the tip. Following with this change, the traction of cohesive element also increases. Traction changes from 0 to Tmax, traction in this point is the largest tensile strength, the opening at this point is called critical opening (δn), and interface bonding degradation starts at this point. (2) As the opening increases, traction decreases with the opening until 0; during this time, the cohesive element is continually damaged. The opening when traction reaches 0 is called failure opening (δt). (3) Exceeding this point, the cohesive element was totally damaged. Because of the fluid flow, the opening will continually increase.
Figure 3. The hydraulic fracture tips.
Cohesive using stiffness degradation to describe the damage evolution process of the element:
The D is the damage index to the elastic linear material
where δm is the maximum effective relative displacement, δ0 is normal loading displacement and δf is shear loading displacement.
There is a fluid flow in the viscous surface gap. By Newtonian rheology, we assume that the fluid is incompressible. The lubrication equations control the tangential flow in the tip, which is formulated from Poiseuille’s law:
The ordinary flow is explained:
where pt is the top surface pore pressures, pb is bottom surfaces of the fracture pore pressure, and Ct and Cb define the corresponding fluid leak-off coefficients.
Interaction Model Setup
ABAQUS is finite-element software for structural and field analysis in the mechanical, civil, and electronics industries, which mainly include CAE, STANDARD, and EXPLICIT . In this program, a hydraulic fracture is injected and interacts with a natural facture in rock, and is analyzed and simulated by it. Existing cracks appear to be closed, and there is no initial separation in the entire path. The Newtonian fluid, which is incompressible, continues to enter the fracture at a constant injection flow rate (Q0) and causes the initiation and expansion of fractures created by hydraulic fracturing.
There is no fluid leak-off through the impermeable surfaces of the fracture, so that the simulation only models the flow in the fracture radius direction. Six-node coupled displacement and pore-pressure axisymmetric cohesive elements are used to model the fracture, and the rock is modeled with 4-node axisymmetric elastic elements. The cohesive elements at the injection point are defined as initially open to allow entry of the fluid, which makes the process of initial flow and fracture growth come true. Since the model we built is for tensile fractures, the cohesive element is continually undergoing damage until 0, and fails in purely normal deformation conditions without considering shear deformation.
Four categories of parameters were set up, including Young’s modulus in rock properties; tensile strength; maximal horizontal stress and minimum horizontal stress in in situ stress; and viscosity and injection rate in pumping parameters. Some of these values are not fixed; a preliminary simulation can be performed by setting initial values, and then the results can be changed by adjusting the values.