Shock and Vibration

Volume 2018, Article ID 2748408, 14 pages

https://doi.org/10.1155/2018/2748408

## Adaptive Finite Element-Discrete Element Analysis for Microseismic Modelling of Hydraulic Fracture Propagation of Perforation in Horizontal Well considering Pre-Existing Fractures

^{1}State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Beijing 100083, China^{2}School of Mechanical & Civil Engineering, China University of Mining and Technology, Beijing 100083, China^{3}State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, China

Correspondence should be addressed to Yang Ju; nc.ude.btmuc@yuj

Received 4 November 2017; Accepted 1 March 2018; Published 17 April 2018

Academic Editor: Longjun Dong

Copyright © 2018 Yongliang Wang et al. 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.

#### Abstract

Hydrofracturing technology of perforated horizontal well has been widely used to stimulate the tight hydrocarbon reservoirs for gas production. To predict the hydraulic fracture propagation, the microseismicity can be used to infer hydraulic fractures state; by the effective numerical methods, microseismic events can be addressed from changes of the computed stresses. In numerical models, due to the challenges in accurately representing the complex structure of naturally fractured reservoir, the interaction between hydraulic and pre-existing fractures has not yet been considered and handled satisfactorily. To overcome these challenges, the adaptive finite element-discrete element method is used to refine mesh, effectively identify the fractures propagation, and investigate microseismic modelling. Numerical models are composed of hydraulic fractures, pre-existing fractures, and microscale pores, and the seepage analysis based on the Darcy’s law is used to determine fluid flow; then moment tensors in microseismicity are computed based on the computed stresses. Unfractured and naturally fractured models are compared to assess the influences of pre-existing fractures on hydrofracturing. The damaged and contact slip events were detected by the magnitudes, -values, Hudson source type plots, and focal spheres.

#### 1. Introduction

Hydrofracturing technology has been extensively used to stimulate low-permeable unconventional oil and gas reservoirs to increase production rates [1]. The fracture growth is hidden from the designers and microseismicity is often used to infer hydraulic fractures state. Microseismic output can also be computed from changes of the stresses in the numerical model and compared against field microseismicity, in which the fracturing process needs to be properly addressed in numerical simulation for detecting the mechanism of hydrofracturing [2]. A rising tide of evidence from practical treatments in the field suggests that hydrofracturing cracks may initiate and grow in a complicated way that is highly affected by the complexity of naturally fractured reservoir formations [3–5]. If the influencing mechanisms and factors are not fully considered, clearly understood, and resolved, it will not be possible to predict, control, and optimize the hydraulic fracturing efficiency, which would ultimately negatively impact oil and gas production.

The microseismic event can be triggered by different types of failure modes, termed Mode I, Mode II, and Mode III. The damaged event caused by tensile failure, Mode I, represents sources with a volume change. The contact slip event caused by shear slip failure, Mode II and Mode III, represents sliding and tearing, both of which represent slip most likely appearing on the pre-existing fracture surfaces. Microseismic events arise due to displacements in the earth’s crust, resulting in faint earth tremors. The information at the source of the displacements can be calculated and stored in a moment tensor, which provides information on the magnitude and orientation of an event. Moment tensors can be determined in several ways, including stress drop, kinetic energy, or internal forces. Once the moment tensor is determined, the information may be presented in a variety of forms, such as Hudson source type plots and focal mechanisms as beach ball plots [6]. The concept behind integrating geomechanics and microseismic prediction for the finite element method (FEM) [7] is based on the approach proposed by Hazzard and Young [8]. The integrated geomechanics and microseismic prediction was applied to a North Sea field, where subsidence prediction was used to calibrate the coupled flow-geomechanical model of the field [9]. This paper will detect monitoring fractures in the reservoir layer with the potential increased risk of microseismicity. The approach adopted in this paper is outlined by Angus et al. [7]. Seismic events are monitored in space and time, where a microseismic event is predicted to occur within an element when the effective stress satisfies the Mohr-Coulomb yield envelope. Based on the differential stress tensor Angus calculates a pseudo-scalar seismic moment (i.e., stress drop) which can be used to infer microseismicity. In order to describe a seismic event in terms of moment tensors, the physical processes are connected to the system of equivalent forces. In an unfaulted medium, the equivalent force system would produce the same displacements in the far field as at the source [10]. Observations of the displacement field can provide information on the equivalent force system. Failure can be thought of as a sudden change in the stress-strain relation. Prior to an event, the stress field satisfies the equations of equilibrium. At failure time the stress field will change, causing dynamic motions which release elastic waves. Many researchers express the equivalent forces in terms of a stress “glut” from [11]. However, this paper uses the true stress, which is not known, to express a moment tensor through a stress “drop” [12], that is, the temporal difference in stress which occurs before and after an event. The method calculates the moment tensor through the stress drop and then uses this information to output the relevant results to obtain the event number, the location coordinates, “beach ball” plot, and so forth.

In situ monitoring and laboratory experiments have been used to identify the morphologies of hydrofracturing cracks [13–17]. Microseismic and nondestructive testing methods have been used to investigate the hydrofracturing crack morphologies [18–22]. The laboratory experiments show that the dip angle affects crossing behaviours between hydraulic and pre-existing fractures to form the complex fracture networks [23]. Some other techniques, that is, microseismic measurements, mine-back experiments, and surface tiltmeter, also have been effectively applied [24–28]. The significant effects of pre-existing fractures on hydraulic fractures were investigated [29–31]. However, the actual hydrofracturing behaviours and evolutionary process are difficult to be detected, so numerical methods have been developed as alternatives.

Recently, some numerical methods and models considering the effects of hydromechanical (HM) coupling and pre-existing fractures have been developed to simulate the hydraulic fracturing. For analysing the effects of pre-existing microcracks on the supercritical CO_{2} fracture network [32], the authors of this paper have established the three-dimensional (3D) bonded particle models (BPMs), based on the discrete element method (DEM) [33]. For investigating the hydrofracturing cracks in heterogeneous media, the authors of this paper developed the continuum-based discrete element method (CDEM) [34]. The conventional FEM has limitations to deal with the fracture propagation; furthermore, the mesh refinements have to be adopted [35]. The coupled finite element- (FE-) discrete element (DE) method makes the cracks flexibly propagate between the elements; the adaptive analysis and remeshing strategy could conquer the problems in solving fracture propagation by conventional FEM [36]. Using the FE-DE method and local remeshing technique, the software package ELFEN TGR [37] has been developed and will be applied in this paper.

In order to represent the actual fractured characteristics of a rock, this study uses the discrete fracture network (DFN) model for modelling naturally fractured reservoir. The unfractured and fractured models will be compared with identical conditions to study the effects of pre-existing fractures on hydrofractures and microseismicity. The remainder of this paper is organized as follows. In Section 2, the numerical simulation using the adaptive FE-DE method and the governing equations are described in detail; the adaptive remeshing technique of fractures propagation is presented. Section 3 presents the microseismic modelling techniques for evaluation and visualization of moment tensors. Section 4 presents the numerical fracturing models in perforated horizontal wells, DFN model for naturally fractured reservoir. Section 5 outlines the numerical results, the analysis of the crack growth, and distribution in naturally fractured reservoir, and the damaged and contact slip events were detected by magnitudes, -values, Hudson source type plots, and focal spheres. Concluding remarks are summarised in Section 6.

#### 2. Adaptive Finite Element-Discrete Element Method for Hydraulic Fracturing

##### 2.1. Governing Equations

In this study, the adaptive FE-DE algorithm combining the FEM and DEM, complementary HM coupling, and fracture network mechanism are introduced by transfer between physical quantities of the solid stress, fluid pressure, and fracture network. The three main sets of governing equations are for structure field, seepage field, and network field.

The main governing equations are derived assuming that the equilibrium of stresses with an appropriate constitutive model is able to mimic both tensile and shear failure, and the mass conservation of fluid flow inside the fracture region with a flow constitutive response is able to recover parallel plate flow theory. The update of the mechanical stresses satisfies the momentum balance equation, with the assumption that fluid acceleration relative to the solid and the convective terms can be ignored. The mechanical governing equation is given by the following [35]:where is the spatial differential operator; is the effective stress tensor; is the Biot coefficient; is the identity tensor; is the pore fluid pressure in the rock formation; is the wet bulk density; and is the gravity vector.

The liquid seepage equation combines mass conservation along with Darcy’s law and is given as follows [37]:where is the intrinsic permeability of the porous media; is the viscosity of the pore liquid; is the pore liquid pressure; is the density of the pore liquid; is the porosity of the porous media; is the bulk stiffness of the pore liquid; is the bulk stiffness of the solid grains; is the volumetric strain of the porous media.

The fracture fluid flow governing equation of the network field is given by the following [35]:where is the intrinsic permeability of the fractured region; is the viscosity of the fracturing fluid; is the fracturing fluid pressure; is the density of the fracture fluid; is the storage coefficient which is effectively a measure of the compressibility of the fractured region when a fluid is present; and is the aperture strain rate. Assuming parallel plate theory, the intrinsic permeability of a fractured region is given by the following:where is the fracture aperture. The storage term is given by the following:where is the fracture normal stiffness and is the bulk modulus of the fracturing fluid.

The mass transport equation is based on the convection equation and given aswhere is the volumetric concentration; is the Darcy flux for the fluid phase in the fracture; is the effective width of the fracture; is the fluid flux through the fracture walls.

The three main equations governing a production stage are derived assuming that the equilibrium of stresses is established based on an appropriate constitutive model for mechanical analysis, the mass conservation is based on the Darcy’s law for gas seepage analysis, and fracture region considers gas flow characteristics for gas network analysis. Note that the mass transport field analyses are inactive during the gas production stage. The update of the mechanical stresses satisfies the momentum balance equation with the assumption that the fluid acceleration relative to the solid and the convective terms can be neglected. The mechanical equation is given as (1).

The gas seepage equation combines mass conservation along with Darcy’s law and is given as follows [37]:where is the viscosity of the pore gas; is the pore gas pressure; is the density of the pore gas; is the porosity of the porous media; is the bulk stiffness of the pore gas; is the mass of adsorbed gas per unit volume.

Similar to the gas seepage equation, the gas network equation combines mass conservation along with Darcy’s law and is given by the following:where is the gas compressibility; is the gas compressibility factor.

##### 2.2. Adaptive Remeshing for Fractures Propagation

The geometry insertion technique is a means of introducing new geometry lines in two-dimensional (2D) models or geometry areas in 3D models into the FE discretisation, which do not necessarily follow the edges of the FE mesh. The newly introduced geometry lines or areas simulate fracture growth during hydraulic stimulation. In addition, a local remeshing algorithm is implemented, which only operates around the fracture tip region, thus avoiding the need for a computationally exhaustive global remeshing.

At the material level each element response follows continuum damage theory [37]. The schematic of the damage process zone and its numerical treatment introduce the opening of a crack mouth with a tensile damage zone ahead of the fracture tip. The elements at the tip are fully damaged with the remaining elements below this threshold and still capable of supporting load. The predicted fracture surface is not forced to follow the element edges and can follow the stress state as dictated by the simulation. Due to the assumption of failure, the failure direction is defined as orthogonal to the maximum tensile principal stress. When the failure path exceeds a user-specified length, all points along the path are used to form a geometric entity, and this finally leads to the fracture surface. According to the above loop for fracture propagation, the hydraulic fractures effectively extend based on the remeshing.

The important modelling aspect is the capability to deal with meshing the new geometry in a computationally attractive manner. A local meshing methodology is employed which makes the need for a traditionally expensive global remeshing redundant, which is also likely to introduce dispersion of key material variables such as stress, strain, and damage indicators. The local meshing zone is defined via a patch region which is adjacent to the fracture tip, and it is always very small compared to the problem size, so in relative terms the computational cost is low. It is very important that the remeshing is only performed locally at a fracture tip. The mesh in all other regions of the problem domain remains unchanged. Clearly, fracturing is a very dynamic process in that the model is constantly changing, so this procedure of following the fracture tip via a patch region and subsequently only remeshing locally is continually being updated during hydraulic stimulation. This is extremely important for this class of problem as it has been observed that for some tight gas reservoir stimulated fracture lengths reach values of many hundreds of metres [37]. Here the key techniques are the error estimation and adaptive remeshing for the above-mentioned fracture propagation and local remeshing. The error estimation method is based on the superconvergent patch recovery (SPR), in which the stress recovery of element using the higher-precision stress results in Gaussian integration points as the natural superconvergent points in triangular and quadrilateral element [38, 39]. There is a parent-child relationship between the network and fracture surface nodes which allows, for example, the fluid pressure from the network field to be transferred as an external load to a corresponding node on the fracture surface of the structure field. During a local remeshing around a fracture tip this mesh topology must be maintained, and this is achieved by initially stitching the mesh back to a nondiscrete body, performing the local remeshing on the resulting bonded domain and then displacing the elements back to their original spatial location after the remeshing. This is achieved by monitoring the displacements of the fracture surfaces such that they can be returned to their position after the remeshing.

#### 3. Microseismic Modelling Techniques

##### 3.1. Evaluation of Moment Tensors

The differential stress from before and after events is calculated over a single time step, and the eigendecomposition of the drop in the effective stress tensor is carried out as follows:where is the eigenvector and is the eigenvalue. The eigenvalues are used to scale the eigenvector, resulting in the definition of the moment tensor as follows:The moment tensor is rotated into its principle direction using the unit eigenvectors; then the moment tensors could be computed [37].

##### 3.2. Visualization of Moment Tensors

The Hudson plot is concerned with the relative sizes of the three principle moments and representing the probability density of these sizes [40]. The aim is to separate the force system from the information related to orientation in the moment tensor. This is done to ensure that the model is capable of interpreting the single or combined physical processes, for example, shear, tensile, explosion, or implosion, each of which would have their own force system and principle moments. A parameterisation method is used to separate the physical processes [40], and two parameters are used:(1) is used to describe the constant volume in the source and provides information on shear failure.(2) is used for the changing volume.

The parameters can be plotted for each source type on a 2D diagram, known as the Hudson source type plot as shown in Figure 1. For a specific region of the plane, an equal area distribution of the source is produced. When there is no prior knowledge of a specific source type, the probability of the source type which corresponds to a range of and values is directly proportional to the region area. The Hudson source type plot can thus be used to determine the type of mechanism occurring for every event. This can be plotted against known distances from the source and a frequency plot for each mechanism. The decomposition of a moment tensor can be done in various ways and therefore is not unique. However, a method is used to determine the probability density of the source types which compose moment tensor [40]. The probability of each source type is then plotted in a Hudson source type plot. Together with the orientations of the principal moments, it is then possible to obtain a full picture of the event.