<figure><img src="image_1.png" /></figure>

HAL
open science

# Meshless method simulation and experimental investigation of crack propagation of CBM hydraulic fracturing

Guojun Wen, Haojie Liu, Hongbo Huang, Yudan Wang, Xinyu Shi

► To cite this version:

Guojun Wen, Haojie Liu, Hongbo Huang, Yudan Wang, Xinyu Shi. Meshless method simulation and experimental investigation of crack propagation of CBM hydraulic fracturing. Oil & Gas Science and Technology - Revue d'IFP Energies nouvelles, 2018, 73, pp.72. 10.2516/ogst/2018074. hal-01950488

HAL Id: hal-01950488

https://hal.science/hal-01950488

Submitted on 10 Dec 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Meshless method simulation and experimental investigation of crack propagation of CBM hydraulic fracturing

Guojun Wen<fnref n="1" />, Haojie Liu<fnref n="1" />, Hongbo Huang<fnref n="2" />, Yudan Wang<fnref n="1" />,*, and Xinyu Shi<fnref n="1" />

Received: 15 March 2018 / Accepted: 3 October 2018

**Abstract.** For simulating CoalBed Methane (CBM) hydraulic fracturing using 3-D meshless method, this paper analyzed the hydraulic fracturing mechanism and cracking form for coal rock and established the geometric and mathematical models of hydraulic fracturing propagation in coal rock in terms of the Hillerborg model on crack opening displacement theory. With the theoretical basis of hydromechanics, the formulas for calculating hydraulic pressure inside the fracture by numerical simulation were deduced from the analysis on this fluid-structure interaction problem. The geometric and mathematical models established above were described by 3-D meshless Galerkin (EFG, Element-Free Galerkin) method and compiled into the numerical simulation program using VB and FORTRAN programming language to simulate the fracture propagation for an actual coal rock sample with a drilling hole as an example. Then the physical simulation experiment of hydraulic fracturing propagation of coal seam was conducted on the same coal rock sample. Through the direct observation with naked eyes and detection by advanced instruments of ESEM and Micro-CT, the shape and parameters of cracks on the surface of and inside the coal rock sample were achieved, which indicated that experimental results are reasonably consistent with numerical simulation results.

## 1 Introduction

In recent decades CoalBed Methane (CBM) has become an important source of energy all over the world. Hydraulic fracturing is a popular and powerful technology mainly used in the CBM industry to stimulate reservoirs for improving the productivity. Although the productive benefit from hydraulically fractured wells is potentially large, it is hard to predict the fracture propagation and perform an optimal hydraulic fracturing process. The special structures and properties of coal seam (such as natural cleat, fracture growing, lower elastic modulus, higher Poisson ratio, frangibility, and strong heterogeneity), together with drilling disturbance and fluid flow make it difficult to figure out the complex and diverse 3-D (three dimensional) fracture propagation using mathematical analysis method without the help of computer.

Aimed to the principle and regulation of fracture propagation of oil and gas recovery for an optimal hydraulic fracturing process, a number of studies are carried out recently by computer numerical simulation or physical experimental simulation. Finite element methods are usually employed to the numerical investigation of fracture shape and propagation laws on hydraulic fracturing crack

of various kinds of rock mass. According to the relevant literatures, FRANC-3D, HYFRANC-3D, or ABAQUS [1-4] are popular numerical simulation software based on the finite-element methods for simulating the propagation of hydraulic fractures.

These studies are instrumental to understand the complex fracture growth to some extent. Results from Rahman show that two or more multiple transverse fractures may propagate but they involve different treatment pressure behaviors compared to a single fracture. Mofazzal Hossain and M.K. Rahman developed insights to explain how the complex fracture geometry grows under different stress and well conditions, and their effects on the injection pressure. Zuorong Chen developed the cohesive zone model to avoid the singularity in the crack tip stress field that is present in classic fracture mechanics. Yao Yao discussed the fundamental mechanical difference of line elastic fracture mechanics and cohesive fracture mechanics-based methods in according to the size of the fracture process zone and its effect on crack extension in ductile rock.

As a matter of fact, however, it was just an expedient way to simulate the fracture propagation of rock cracking, especially resulting from hydraulic fracture inside CBM wellbore, using finite element methods. As we well known, finite element methods were originally defined on meshes of data points. In such a mesh, each point has a fixed number

<footnote n="1">¹ School of Mechanical & Electronic Information, China University of Geosciences (Wuhan), Wuhan 430074, China</footnote>

<footnote n="2">² China Geological Equipment Group Co., Ltd., Beijing 100102, China</footnote>

<footnote>* Corresponding author: wangyodan@163.com</footnote>
of predefined neighbors, and this connectivity between neighbors can be used to define mathematical operators like the derivative for constructing the equations to simulate. But in simulations where the material being simulated can move around (as in computational fluid dynamics in hydraulic fracturing) or where large deformations of the material can occur (as in simulations of plastic or linear elasto materials of coal rock), the connectivity of the mesh can be difficult to maintain without introducing error into the simulation. If the mesh becomes tangled or degenerate during simulation, the operators defined on it may no longer give correct values. The mesh may be recreated during simulation (called remeshing), but this can also introduce error, since all the existing data points must be mapped onto a new and different set of data points.

As a much better alternative solution for simulating fracture propagation resulting from hydraulic fracturing, meshfree methods (also called element free methods), may remedy above mentioned shortfalls using finite element methods [5–8] since they are very suitable for: (1) simulations where creating a useful mesh from the geometry of a complex 3D object may be especially difficult or require human assistance, (2) simulations where nodes may be created or destroyed, such as in cracking simulations, (3) simulations where the problem geometry may move out of alignment with a fixed mesh, such as in bending simulations, and (4) simulations containing nonlinear material behavior, discontinuities or singularities.

Based on the use of moving least-squares interpolants with a Galerkin method, the Element Free Galerkin (EFG) method is particularly effective in line elastic fracture problems with the advantages of the circumvention of the need for element data, no any volumetric locking, the rate of convergence much faster than that of finite elements significantly, and a high resolution of localized steep gradients [9]. The implementation of the EFG method for problems of fracture and static crack growth based on an infinite plate with a central circular hole subjected to a unidirectional tensile load in the *x* direction shows that accurate stress intensity factors can be obtained without any enrichment of the displacement field by a near-crack-tip singularity and that crack propagation can be easily modeled since it requires hardly any remeshing [10, 11]. More simulating examples for arbitrary three-dimensional propagating cracks dynamically in elastic bodies also further verify that EFG method is particularly attractive in dynamic crack propagation under various tensile loading cases, such as quasi-static loading and step load, and even the combination of tension and torsion [12–15].

An expansive application on a crack within an infinite square plate using element free method attempted to simulate the hydraulic fracturing of rock mass under the external unidirectional compression loading without or with water pressure inside the rock mass [16]. Although this kind of condition seems like but totally not equal to the hydraulic fracturing in the wellbore for oil and gas recovery, the relevant results still verify the meshless method can commendably deal with the problems of rock fracturing under hydraulic pressure. In addition, a coupled hydro-mechanical analysis for prediction of hydraulic fracture propagation in

saturated porous media using EFG meshless method successfully simulates the leak-off of fluid from the fracture into the surrounding homogeneous and heterogeneous saturated soils [17].

From the analysis above, all of these significant researches indicate that all properties of EFG method seem to be customized for the requirements for simulating dynamic fracture propagation of hydraulic fracturing in coal seam for CBM recovery. However, these researches are hardly involved in directly applying EFG method to analyze fracture initiation and propagation of CBM hydraulic fracturing. Furthermore, up to now, many research scholars are still accustomed to numerical modeling of hydraulic fracturing in rock mass by finite element method [18, 19].

The main objective of this paper is to investigate the feasibility and validity of the EFG method on the crack propagation of CBM hydraulic fracturing by geometric and mathematic modeling after theoretic analysis on the cracks of an actual coal rock sample with a drilling hole. Then the physical simulation experiments are conducted on the same coal rock sample. Through the direct observation with naked eyes or detection by advanced instruments on the shape and parameters of cracks, the actual physical process of fracture propagation is scientifically described to verify and improve theoretical model and numerical model. The research results will be beneficial to the design of hydraulic fracturing in coal seam.

# 2 Geometric and mathematic modeling

## 2.1 Judging criteria for coal rock cracking

Although coal rock belongs to fragile material, the tip of existing crack inside the wellbore under hydraulic pressure will not be in the form of simple single brittle fracture, but rather in the form of the combination of brittle-elastic-plastic due to the crustal stress from coal reservoirs, especially buried deeply [20, 21]. This process can be considered as three procedures: (1) the fracture tip produces brittle fracture while it is subjected to the maximum fracture strength resulting from hydraulic pressure, (2) due to the constraint of in-situ stress, the fracture will propagate under the continuous action of hydraulic pressure till the pressure is unloaded or fracture propagation is complete, and (3) the fracture occurred will be partially closed after pressure unloading, thus forming the plastic fracture. This also can be vividly called a “spring-splint” model illustrated in Figure 1. K₁ and K₂ are the equivalent spring coefficients of the coal rock between roof and floor, *w* is the crack width, and *α* is the opening angle.

Within this model, the mechanic effect for a fracture propagation of coal rock can be stated as that the extended fractures are subjected to a closing force *σ* from a spring drive fractures. In this case, when the stability coefficient of tip ζ = *σ*/*w* exceeds the recovery coefficient *K* of spring, plastic fractures will be completed finally. Hence the fracture propagation laws based on plastic criteria can be used as the judging criteria in theory for hydraulic fracture propagation laws.
<figure><img src="image_3.png" /><figcaption>Fig. 1. The approximate “spring-splint” model for CBM hydraulic fracturing.</figcaption></figure>

The Hillerborg model based on Crack Opening Displacement (COD) theory has been proved that it is consistent with the mechanism of fracture propagation and quite suitable to the numerical simulation analysis for it [22–24]. Therefore the mechanic model based on the Hillerborg model is built to simulate the hydraulic fracture propagation of coal rock using meshless method. According to the spring-splint model, for a certain kind of coal rock, the critical opening displacement resulting from the brittle fracture or elasto plastic deformation can be derived from Equations (1) and (2) respectively.

Since the brittle fracture is dominated by the plane stress, the critical opening displacement is

$$w_c = \frac{2K_{IC}^2}{Ef_t}. \qquad (1)$$

Similarly, since the elasto plastic deformation is dominated by the plane strain, the critical opening displacement is

$$w_c = \frac{2(1 - v^2)K_{IC}^2}{Ef_t}, \qquad (2)$$

where $K_{IC}$ is the fracture toughness, $f_t$ is the tensile strength, $E$ is the elastic modulus, and $v$ is the Poisson ratio of coal petrography.

Then with the COD criterion of $w \ge w_c$, it can be judged that whether the fracture propagation directions follow the original fracture surfaces or not.

## 2.2 Modeling hydraulic fracturing shape of coal rock

The single fracture propagation model is the basis of analyzing and calculating complex fractures. There are three basic modes of crack extension, i.e., opening mode (Mode I), shearing mode (Mode II) and tearing mode (Mode III). Good agreement is found between the analytically and numerically estimated Mode I for the single fracture scenario [1]. That is the case the crack of coal rock by hydraulic fracturing is mainly Mode I.

Considering the differences on the natural properties of between the coal roof and coal floor resulting in high aspect ratio hydraulic fracture in coal seam, the constant-height rectangular section (KGD mode) can be employed to study fracture propagation based on fracture Mode I of coal rock

<figure><img src="image_4.png" /><figcaption>Fig. 2. The hydraulic fracturing shape of KGD model for coal rock.</figcaption></figure>

illustrated in Figure 2. The so-called KGD model (plane strain) was independently developed Geertsma and de Klerk at the same time the PKN model was formulated by Nordgren based on the PK model developed by Perkins and Kern. Variations of the KGD and PKN models were used routinely for treatment designs as recently as the 1990s, and are still used today for computer simulation of hydraulic fractures. The PKN model is applicable to long fractures of limited height and elliptical vertical cross-section, whereas the KGD model for width calculation is height independent, and is used for short fractures where plane strain assumptions are applicable to horizontal sections [25, 26].

## 2.3 Computing hydraulic pressure inside the fracture

For a specific coal seam, the fracture initiation and propagation law can be expressed by the functional relationship of the surrounding ground stress and the water pressure in coal seam since it belongs to typical fluid-structure interaction problem. As far as KGD hydraulic fracture model, the computing formula of pressure distribution of the water flow inside the fracture can be deduced using mass conservation and momentum conservation theorems.

### 2.3.1 Water flow motion equation inside the fracture

With regard to the water flow problem of high aspect ratio hydraulic fracture in coal seam, the control volume method is usually an efficient way to figure out the solution [27].

In terms of the mass conservation theorem, there is

$$\frac{dm}{dt} = \rho \frac{d}{dt} \int_V dv + \rho \oint_{\Omega} u(u \cdot n) dA = 0, \qquad (3)$$

where $m$ is the mass of control volume, $\rho$ is the mass density of control volume, $v$ is the volume of control volume, $\Omega$ is the surface of control volume, $u$ is the velocity vector of control volume, and $n$ is outer normal vector of control volume.

Also, in terms of the momentum conservation theorem, there is:

$$\frac{dM}{dt} = \rho \frac{d}{dt} \int_V u dv + \rho \oint_{\Omega} u(u \cdot n) dA = \sum F, \qquad (4)$$
where M is the momentum of control volume and $\Sigma F$ is the vector sum of forces on the control volume.

### 2.3.2 Modeling theoretic hydraulic pressure distribution inside hydraulic fracture

To simplify the computing model rationally, some recognized rules inside coal petrography fractures should be adopted on the premise of no excessive effect on computational accuracy as follows:

(1) The fluid volume should be incompressible;

(2) The fluid belongs to Newtonian fluid;

(3) The fluid loss should be ignored; and

(4) The water flow should obey the rules of one-dimensional (1-D) laminar flow.

The velocity distribution of 1-D laminar flow conforms to Poiseuille flow at time t in any section illustrated in Figure 3. There is:

$$u_{x,y} = u_0 \left( 1 - \frac{4y^2}{w_x} \right), \quad (5)$$

where $u_{x,y}$ is the flow velocity at any point inside the section x, $u_0$ is the maximum fluid velocity in the section x, and $w_x$ is the fracture width of the section x.

As the fracture propagation shape of coal rock is chosen to be the KGD model, the longitudinal profile of the fractures section should be the constant-height rectangular section while the transverse profile of the fractures section should be prolate ellipse expressed by expression (6) at time t is

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1, \quad (6)$$

where a and b are the long axis and minor axis of fracture ellipse section at time t respectively.

By taking an infinitesimal segment $dx$ for the control volume $x$ into the mass conservation equation and momentum conservation equation of fluid, the differential equation of hydraulic pressure inside the fracture can be obtained after ignoring higher-order differentials illustrated in Equation (7):

$$\left(\frac{1}{3} - \frac{1}{9}\frac{\rho}{\mu}w_x\frac{\partial w_x}{\partial t}\right)\frac{\partial p_x}{\partial x} - \frac{p_x}{w_x}\frac{\partial w_x}{\partial x} = 0, \quad (7)$$

where $w_x$ is the fracture width of section $x$, $\mu$ is the kinetic viscosity coefficient of fluid, $q_x$ and $p_x$ is the quantity of flow and the average head pressure of section $x$ respectively.

Substituting Equation (6) into Equation (7) gives

$$\frac{\partial q_x}{\partial x} + \frac{x}{a^2 - x^2} \frac{1}{\frac{1}{3} - \frac{1}{9}\frac{\rho}{\mu}w_o w'_0 \left(1 - \frac{x^2}{a^2}\right) - \frac{1}{9}\frac{\rho}{\mu}w_0^2 a'^2} p_x = 0, \quad (8)$$

where $w_0$ is the fracture width, $w'_0$ is the change rate of fracture width with time, and $a'$ is the change rate of fracture length with time.

<figure><img src="image_5.png" /><figcaption>Fig. 3. The control volume and fluid velocity distribution.</figcaption></figure>

For the brittle-elastic-plastic fracture of coal rock, the change rate of fracture length with time can be set as zero, i.e., $\alpha' = 0$. Then by introducing the boundary condition, i.e., $i_{x|x=0} = p_0$ ($p_0$ is the water injection pressure), Equation (8) can be simplified to Equation (9) that statements the distribution of hydraulic pressure inside the fracture of coal rock:

$$p_x = p_0 \left[ \frac{a^2 - x^2}{\frac{1}{3} + \frac{1}{9} \frac{w_0 w'_0}{v} \left(1 - \frac{x^2}{a^2}\right)} \right]^{\frac{3}{2}}, \quad (9)$$

where $v$ is the kinematic viscosity coefficient of fluid and $1/v = \rho/\mu$.

## 2.4 The 2-D processing method for 3-D fractures

For success solution of mechanics and mathematical models built above, coal rock is sliced along its bedding direction and then projected so that the coal seam with a certain thickness can be disposed as a plane approximately illustrated in Figure 4. This can be called 3-D to 2-D process in which complex 3-D fracture problem are transformed to easier 2-D plane fracture problem.

Before the 3-D to 2-D process, the following operations should be considered:

1. In order to prevent the nodes inside the solution domain overflow the solution domain due to the large deformation caused by crack propagation during hydraulic fracturing, the circumcircle of rectangular section in bedding direction is served as the boundary for background integral mesh. Since the circumcircle is theoretically the minimum regular boundary including the rectangular solution domain, any other envelope curves will be unreflected.

2. As the thickness of selected bedding slice along bedding direction is far less than that of coal reservoir, i.e., $m \ll c$, it can be regarded as an actual plane. Furthermore, it is assumed that the surface of bedding section is rather smooth without obvious uplift and fault, thus the vertical ground stress can be set as zero.
<figure><img src="image_6.png" /><figcaption>Fig. 4. The slice along bedding direction of coal rock.</figcaption></figure>

As shown in Figures 5 and 6, the model is very explicit after the 3-D to 2-D process on both mechanics and geometric sides.

According to the principle of geometric transformation, the geometric transformation of 3-D data points can be expressed as

$$ (x'y'z') = (xyz)A_H, \quad (10) $$

$$ \Omega(\Gamma) \stackrel{A_H}{\Rightarrow} \Omega_{xoy}(\Gamma_{xoy}). \quad (11) $$

The solution domain is

$$ \Omega_i = \begin{cases} x_i = R_i \cos \alpha_i \\ y_i = R_i \sin \alpha_i \quad (i = 0, 1, 2, \dots, n). \\ z_i = R_i \cot \beta_i \end{cases} \quad (12) $$

The boundary condition is

$$ \Gamma_{xoy} = \begin{cases} R_i = r \\ R_i = R \end{cases}. \quad (13) $$

The transformation matrix of vertical view is

$$ A_H = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & -n & 1 \end{bmatrix}. \quad (14) $$

At the step $i$, the fracture propagation length will be

$$ l_i = \sqrt{(x_i - x_0)^2 + (y_i - y_0)^2 + (z_i - z_0)^2} \quad (i = 0, 1, 2, \dots, n). \quad (15) $$

As the fracture initiation always starts at the wall of wellbore, the initial solution domain is

$$ \Omega_0 = \begin{cases} x_0 = r \cos \alpha_0 \\ y_0 = r \sin \alpha_0 \\ z_0 = r \cot \beta_0 \end{cases}. \quad (16) $$

In Equations (12), (13), (15) and (16), $x_i$ is x-coordinate of point at the solution domain, $y_i$ is y-coordinate of point, $z_i$ is

<figure><img src="image_7.png" /><figcaption>Fig. 5. The 3-D to 2-D process for mechanics model.</figcaption></figure>

z-coordinate of point, $R$ is the constraint boundary radius, $r$ is the radius of water injection hole, $R_i$ is the curvature radius at any node of the solution domain, $\beta_i$ is the dip angle of sought point, $\alpha_i$ is the azimuth angle of sought point, $x_0$ is the x-coordinate of point of fracture initiation, $y_0$ is the y-coordinate of point of fracture initiation, $z_0$ is the z-coordinate of point of fracture initiation, $\beta_0$ is the dip angle of point of fracture initiation, $\alpha_0$ is the azimuth angle of point of fracture initiation, $i$ is the step number for $i$ infinitesimal segment $dx$, and $l_i$ is the fracture length at $i$ step.

After fracture propagation is complete or the step is ended, the 2-D fracture graph in bedding slice section will be transformed back to the 3-D fracture propagation model by multiplying an inverse matrix $\bar{A}$ expressed in Equations (17) and (18)

$$ \Omega_{xoy}(\Gamma_{xoy}) \stackrel{\bar{A}}{\Rightarrow} \Omega(\Gamma), \quad (17) $$

$$ \bar{A} = \overline{A_H}. \quad (18) $$

## 3 Program design for numerical simulation

Since many studies have conducted on the theory of EFG method and its effectiveness in fracture problems, this section will not repeat them but only simply indicate the program flow of combined programming in FORTRAN and Visual Basic programming languages based on that.
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** The flowchart for overall structure of coal rock fracture propagation in hydraulic fracturing based on meshless methods.</figcaption></figure>

method with essential boundary conditions, NaturalBC for the implementation of natural boundary condition, SolveBand for the solution of asymmetric band matrix system equation, GetDisplacement for the calculation of actual

displacement, and Getstress for the solution of point stress and strain.

External program module, also called add-in program, is used to figure out the distribution of hydraulic pressure
inside the existing fracture and simulate the fracture propagation for coal rock sample. It consists of several customized submodels: 3DChange2D for the transformation of three dimensional fracture processing, XChangeNode for the modification of original fracture plane node information and generation of new fracture plane, XInsertRNode for the deletion of node information in background grid cell for the existing fracture plane and node insertion in fracture area, XNODE for the acquisition and loading of fracture surface node and modification of boundary condition, FMONODE for the calculation of fracture morphology and orientation parameters, MFLOAD for the calculation of fluid pressure on fracture surface, and XLOAD for the calculation of loading force in the fracture surface based on Hillerborg model.

### 3.2.1 Program for computing hydraulic pressure inside the fracture

The flowchart of hydraulic pressure inside the fracture illustrated in Figure 9 is established based on the mathematical model of Equation (9). Two submodules of FMONODE and MFLOAD are combined together to complete the loading of fluid pressure.

### 3.2.2 Program for simulating hydraulic fracture propagation

According to EFG process, hydraulic fracture propagation based on Hillerborg model can be programmed as the flowchart illustrated in Figure 10. Five submodules of 3DChange2D, XChangeNode, XInsertRNode, XNODE, and XLOAD are involved during this process.

## 4 Numerical simulation example

After the coding and compiling of all modules, the program can be used to simulate hydraulic fracturing.

### 4.1 Preparing coal rock sample for simulation

For comparison with physical experimental results at laboratory later, an actual coal rock sample of approximately 28 × 28 × 33 mm [3] with φ5 × 27 wellbore was designed to verify the developed mathematical model and EFG program as an example.

The physical properties of this coal rock sample are as follows: elastic modulus parallel to bedding direction is $E = 0.27 \times 10^5$ MPa, elastic modulus perpendicular to bedding direction is $E = 4.3 \times 10^5$ MPa, Poisson ratio is $\gamma = 0.2$, tensile strength parallel to bedding surface is $f_t = 0.3$ MPa, tensile strength perpendicular to bedding plane is $f_t = 1.2$ MPa, fracture toughness parallel to bedding direction is $K_{IC} = 0.8$ MPa m$^{1/2}$, and fracture toughness perpendicular to bedding direction is $K_{IC} = 2.4$ MPa m$^{1/2}$.

### 4.2 Distribution and characteristics of existing fractures of coal rock sample

The Environmental Scanning Electron Microscope (ESEM) of Quanta 200 and Micron Level Computerized Tomography

<figure><img src="image_9.png" /><figcaption>Fig. 9. The program flowchart for computing hydraulic pressure inside the fracture.</figcaption></figure>

Scanner (Micro-CT) of SkyScan1172 were used to detect the distribution and characteristics of existing fractures, including original fractures and artificial fractures resulting from drilling the hole, in the coal rock sample.

After the detection, two representative primary fractures were found in the coal rock sample illustrated in Figure 11. One is transverse fracture (named No. 1) at the depth of 10 mm in the wellbore, its dip angle is 0.43 rad, azimuth angle of 0, crack length is 12 mm, and crack width is 0.07 mm; the other is longitudinal fracture (named No. 2) at the depth of 10.5 mm in the wellbore, its dip angle is 1.28 rad, azimuth angle is 0, crack length is 13 mm, and crack width is 0.05 mm. These initial parameters are shown in Table 1.

### 4.3 Geometric description of existing fractures

In order to intuitionally calculate and analyze the distribution of existing fractures, their shape and parameters were described as geometric model illustrated in Figure 12.

### 4.4 Numerical simulation using developed program

Within the solving domain of bedding plane of projection, 11 × 10 + 10 + 10 nodes were placed. Each of them
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** The program flowchart for hydraulic fracture propagation based on Hillerborg model.</figcaption></figure>

is a center of establishing a rectangular influence domain where 4 × 4 Gauss points were placed for integral solution. The kinematic viscosity of water is $ν = 1.206 \times 10^{-6} \text{m}^2/\text{s}$ under the constant temperature of 20 °C.

After all setting, the numerical simulation was started in the developed program mentioned above. Figure 13 illustrates the simulating interface which is displaying the distribution and related parameters of hydraulic fractures

with triaxial confining pressure under the action of 3 MPa water injection.

## 4.5 Results of numerical simulation

Various water injection pressures of 0–6 MPa at interval of 1 MPa were adjusted to simulating hydraulic fracture propagation on coal rock sample under the action of increasing pressure. The simulated parameters of the existing fractures
<figure><img src="image_11.png" /><figcaption>The results of electronic microscope scanning of No.1 fracture(Transverse)
magnified 28 times and 56 times</figcaption></figure>

**Fig. 11.** The results of electronic microscope scanning of existing fractures before hydraulic fracturing experiment.

<table><caption>**Table 1.** The existing fracture parameters of coal rock sample.</caption><thead><tr><td>Fracture number</td><td>Fracture type</td><td>Fracture length (mm)</td><td>Fracture width (mm)</td><td>Dip angle of fracture tip (rad)</td><td>Azimuth of fracture tip (rad)</td></tr></thead><tbody><tr><td>1</td><td>Transverse</td><td>12</td><td>0.07</td><td>0.43</td><td>0</td></tr><tr><td>2</td><td>Longitudinal</td><td>13</td><td>0.07</td><td>1.28</td><td>0</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** The geometric model of existing fractures in coal rock sample.</figcaption></figure>
Fracture Propagation Numerical Simulation of Coal Petrography

<figure><img src="image_13.png" /><figcaption>Fig. 13. The display interface and results for water injection pressure of 3 MPa in coal rock sample.</figcaption></figure>

(Nos. 1 and 2) and the new fracture (No. 3) are shown in Tables 2–4 respectively.

By cubic polynomial fitting in Matlab, data of the
width, length, dip angle and azimuth angle changed with
various hydraulic pressure in above tables were fitted to
the curve graphs illustrated in Figure 14.

**4.6 Analysis on numerical simulation results**

Major findings from the analysis on the simulated data and
result curves of hydraulic fracturing in coal rock can be con-
cluded as follows.

(1) The expansion of hydraulic fractures is of step response on the water injection pressure since the existing fractures in coal rock sample wouldn't keep expanding in a certain water injection pressure, but rather stop at a certain stable state. Only when the water injection pressure increases to a certain higher situation, the expansion will continue to proceed again.

(2) The expansion of hydraulic fractures is of inheri-
tance since it is mainly based on the existing frac-
tures and in the form of straight line without
distortion normally.
<table><caption>**Table 2.** The No. 1 fracture propagation parameters.</caption><thead><tr><th>Step</th><th>Water injection pressure (MPa)</th><th>Confining pressure (MPa)</th><th>Fracture length (mm)</th><th>Fracture width (mm)</th><th>Dip angle of fracture tip (rad)</th><th>Azimuth of fracture tip (rad)</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>12.0</td><td>0.07</td><td>0.43</td><td>0</td></tr><tr><td>1</td><td>1</td><td>1</td><td>12.7</td><td>0.08</td><td>0.43</td><td>0</td></tr><tr><td>2</td><td>2</td><td>2</td><td>14.3</td><td>0.12</td><td>0.46</td><td>0</td></tr><tr><td>3</td><td>3</td><td>3</td><td>14.8</td><td>0.15</td><td>0.47</td><td>0</td></tr><tr><td>4</td><td>4</td><td>4</td><td>15.0</td><td>0.16</td><td>0.47</td><td>0</td></tr><tr><td>5</td><td>5</td><td>5</td><td>15.1</td><td>0.17</td><td>0.47</td><td>0</td></tr><tr><td>6</td><td>6</td><td>6</td><td>15.1</td><td>0.17</td><td>0.47</td><td>0</td></tr></tbody></table>

<table><caption>**Table 3.** The No. 2 fracture propagation parameters.</caption><thead><tr><th>Step</th><th>Water injection pressure (MPa)</th><th>Confining pressure (MPa)</th><th>Fracture length (mm)</th><th>Fracture width (mm)</th><th>Dip angle of fracture tip (rad)</th><th>Azimuth of fracture tip (rad)</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>13.0</td><td>0.05</td><td>1.28</td><td>0</td></tr><tr><td>1</td><td>1</td><td>1</td><td>13.0</td><td>0.05</td><td>1.28</td><td>0</td></tr><tr><td>2</td><td>2</td><td>2</td><td>13.2</td><td>0.05</td><td>1.31</td><td>0</td></tr><tr><td>3</td><td>3</td><td>3</td><td>13.4</td><td>0.05</td><td>1.32</td><td>0</td></tr><tr><td>4</td><td>4</td><td>4</td><td>14.2</td><td>0.06</td><td>1.46</td><td>0</td></tr><tr><td>5</td><td>5</td><td>5</td><td>14.7</td><td>0.06</td><td>1.52</td><td>0</td></tr><tr><td>6</td><td>6</td><td>6</td><td>15.2</td><td>0.08</td><td>1.55</td><td>0</td></tr></tbody></table>

<table><caption>**Table 4.** The No. 3 fracture propagation parameters.</caption><thead><tr><th>Step</th><th>Water injection pressure (MPa)</th><th>Confining pressure (MPa)</th><th>Fracture length (mm)</th><th>Fracture width (mm)</th><th>Dip angle of fracture tip (rad)</th><th>Azimuth of fracture tip (rad)</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>2</td><td>2</td><td>2</td><td>3.4</td><td>0.04</td><td>0.86</td><td>0</td></tr><tr><td>3</td><td>3</td><td>3</td><td>6.4</td><td>0.06</td><td>1.15</td><td>0</td></tr><tr><td>4</td><td>4</td><td>4</td><td>6.8</td><td>0.07</td><td>1.18</td><td>0</td></tr><tr><td>5</td><td>5</td><td>5</td><td>12.6</td><td>0.12</td><td>1.52</td><td>0</td></tr><tr><td>6</td><td>6</td><td>6</td><td>12.8</td><td>0.15</td><td>1.56</td><td>0</td></tr></tbody></table>

(3) For a given step hydraulic fracturing (this step is usually greater than that of the crack tip's fracture time of occurrence), the fracture propagation in the numerical simulation process is embedded in two stages. During the first stage, called fracture occurrence stage, a new long fracture will suddenly appear at the area without original fractures. During the second stage, called development stage, the fracture will propagate slowly in the water pressure loading process until the water pressure loading is ended or the fracture propagation is complete. It is clear that two stages were involved for the new fracture of No. 3 while only the second stage was for the existing fractures of No. 1 and No. 2.

## 5 Hydraulic fracturing experiments

For comparing with and verifying the numerical simulation results, the physical simulation experiments were conducted

on the coal rock sample used in numerical simulation above in the laboratory.

### 5.1 The experimental procedure

The basic experimental procedure is outlined as follows and illustrated in Figure 15.

(a) Preparing the coal rock sample for experiment with a drilling hole.

(b) Acquiring the properties of coal rock sample using ESEM and Micro-CT before water injection.

(c) Hydraulic fracturing using customized JHLS-I experiment instrument.

(d) Scanning by ESEM and Micro-CT again after water injection.

(e) Comparing the scanned results before and after water injection; and

(e) Analyzing the fracture propagation rules.
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** The relationship of fracture parameter of coal rock and water injection pressure ( triaxial pressure).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** The hydraulic fracturing experimental procedure of coal rock.</figcaption></figure>

The first two steps have been done for the numerical simulation. The water injection is a similar process as many other studies. The following parts will mainly focus on the results and analysis.

## 5.2 Coal rock sample scanning after hydraulic fracturing

After each water injection, ESEM and Micro-CT were used to acquire the properties for the hydraulic fractures.

Figure 16 illustrates the surface characteristics of three fractures scanned by ESEM after 3 MPa water injection

experiment. For ease of human observation, all of them are magnified 28 times, 56 times and 112 times.

Figure 17 illustrates the internal characteristics of three fractures scanned by Micro-CT after 3 MPa and 6 MPa water injection experiments.

## 5.3 Geometric description of hydraulic fractures

In the same way as in the numerical simulation, the shape and parameters of hydraulic fractures after water injection were described as geometric model illustrated in Figure 18 while the detailed parameters are listed in Tables 5 and 6.
<figure><img src="image_19.png" /><figcaption>Fig. 16. The magnified results of electronic microscope scanning of the fractures after the 3 MPa water injection experiment.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. The results of Micro-CT scanning of the fractures after the 3 MPa and 6 MPa water injection experiments.</figcaption></figure>

## 5.4 Interpolating process for data before and after hydraulic fracturing experiments

By cubic spline interpolating in Matlab [28] for the data before and after hydraulic fracturing experiments, the curve graphs of the shape and orientation parameters of three fractures with water injection pressure in above tables were illustrated in Figure 19.

## 5.5 Analysis on physical experiment results

Major findings from the analysis on above images, data and curves before and after the experiment can be concluded as follows.

(1) The expansion of hydraulic fractures of the coal rock will mainly initiate from the existing fractures on the
<figure><img src="image_21.png" /><figcaption>Fig. 18. The geometric model of the fractures in coal rock sample after water injection.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. The relationship of fracture parameters and water injection pressure ( triaxial pressure) of coal rock.</figcaption></figure>

wall of wellbore and then intersect with the internal porosity and fractures to form the fracture network. Less new fracture is able to be observed.

(2) For coal rock with only one single existing fracture, either transverse or longitudinal, hydraulic fracturing propagation will be primarily in the form of straight line along the existing fracture without the occurrence of larger distortion.

(3) For coal rock with both transverse and longitudinal fractures coexisting but no intersection, various phenomena resulting from increasing hydraulic pressure will be appeared: (1) Neither fracture propagated for lower water injection pressure (less than 1 MPa).

(2) While water injection pressure is increased (1–3 MPa), the transverse fracture will propagate significantly but longitudinal one won't. (3) While water
<table><caption>**Table 5.** The fracture parameters of coal rock sample after 3 MPa hydraulic fracturing.</caption><thead><tr><th>Fracture number</th><th>Fracture type</th><th>Fracture length (mm)</th><th>Fracture width (mm)</th><th>Dip angle of fracture tip (rad)</th><th>Azimuth of fracture tip (rad)</th></tr></thead><tbody><tr><td>1</td><td>Transverse</td><td>15</td><td>0.14</td><td>0.43</td><td>0.09</td></tr><tr><td>2</td><td>Longitudinal</td><td>13</td><td>0.07</td><td>1.28</td><td>0.09</td></tr><tr><td>3</td><td>New Transverse</td><td>7</td><td>0.12</td><td>1.16</td><td>0.61</td></tr></tbody></table>

<table><caption>**Table 6.** The fracture parameters of coal rock sample after 6 MPa hydraulic fracturing.</caption><thead><tr><th>Fracture number</th><th>Fracture type</th><th>Fracture length (mm)</th><th>Fracture width (mm)</th><th>Dip angle of Fracture Tip (rad)</th><th>Azimuth of Fracture Tip (rad)</th></tr></thead><tbody><tr><td>1</td><td>Transverse</td><td>15</td><td>0.16</td><td>0.43</td><td>0.09</td></tr><tr><td>2</td><td>Longitudinal</td><td>15</td><td>0.09</td><td>1.53</td><td>0.09</td></tr><tr><td>3</td><td>New transverse</td><td>13</td><td>0.15</td><td>1.54</td><td>0.88</td></tr></tbody></table>

injection pressure continues to be increased (3–6 MPa), both transverse and longitudinal fractures will propagate, and after they intersect with each other, the transverse propagation will be dominated and longitudinal one will cease basically.

(4) When water injection pressure is up to a certain value (about 6 MPa), new longitudinal fracture will be prone to occur for lower three axial confining pressure (about 1 MPa) while new transverse fracture for higher three axial confining pressure (about 3 MPa).

# 6 Conclusion

This paper conducted the numerical simulation and physical experiments for the fracture propagation of hydraulic fracturing in the coal seam. For the numerical simulation, a kind of new geometric model of spring-splint was built for the numerical simulation using EFG method and was applied on an actual coal rock sample with a drilling hole as an example. For the physical simulation, the same one coal rock sample was used for the experiment as well as the advanced detection instruments of ESEM and Micro-CT were used for the nuanced observation on the surface of and inside the coal rock sample.

Many findings are concluded from either numerical simulation or physical simulation experiments. Through comparison and analysis, the validity, feasibility and capability of spring-splint model and EFG method are illustrated by the reasonably good agreement between them. The results also show that since the development degree of transverse fractures (parallel to the bedding surface) is higher than that of longitudinal fractures (perpendicular to the bedding surface) in the coal seam, the extension of transverse fracture is dominated while the longitudinal crack extension is subsidiary during the process of hydraulic fracturing.

Based on the above results and analysis, the following major suggestions will be instrumental to some extent to

the design of hydraulic fracturing process in actual production of CBM:

1. The jet hole should be placed at the direction abounding in transverse fractures in the wellbore to the greatest extent, which will facilitate full extension of transverse fractures with lower water injection pressure, thus improving the efficiency of CBM recovery as well as reducing the cost of hydraulic fracturing.

2. Correspondingly, the jet hole should avoid the direction with longitudinal fractures in the wellbore as far as possible. This will not only be against the CBM recovery efficiency with higher cost of higher water injection pressure, but also induce some accidents such as collapse of the roof resulting from excessive expansion of longitudinal fracture.

3. If the direction of jet hole in the wellbore has both transverse and longitudinal fractures in the trend of separation, water injection pressure should be reasonably adjusted in a range to facilitate the extension of transverse fractures and inhibit the extension of longitudinal fractures.

4. If the direction of jet hole in the wellbore has both transverse and longitudinal fractures in the trend of intersection, water injection pressure should be increased to facilitate the sufficient extension of transverse fractures without influencing the stability of the roof of coal reservoir.

5. For thicker coal seam, hydraulic fracturing should firstly be conducted in the upper part of the reservoir to produce new transverse fracture. Then a certain water injection pressure will be loaded in the deeper part for producing longitudinal fracture after sealing the jet hole in the upper part. In this way, transverse and longitudinal fractures will intersect to form fracture network, thus benefiting the CBM recovery efficiency.

10. In the case of the existence of larger fault resulting from tectonic movement nearby the wall of wellbore,
no jet hole should be arranged at this place regardless of any kind of fractures to avoid the extension of hydraulic fracture along the fault, which can result in relative slip of fault plane to induce roof unstability of coal reservoir.

*Acknowledgments*. This research program was funded by the grant of National Natural Sciences Foundation of China (41672155, 61733016) and Excellent Youth Fund in Hubei Natural Sciences Foundation (2018CFA092).

## References

1. Rahman M.M., Hossain M.M., Crosby D.G., Rahman M.K., Rahman S.S. (2002) Analytical, numerical and experimental investigations of transverse fracture propagation from horizontal wells, *J. Pet. Sci. Eng.* **35**, 127–150.

2. Hossain M.M., Rahman M.K. (2008) Numerical simulation of complex fracture growth during tight reservoir stimulation by hydraulic fracturing, *J. Pet. Sci. Eng.* **60**, 86–104.

3. Chen Z.R., Bunger A.P., Zhang X., Jeffrey R.G. (2009) Cohesive zone finite element based modeling of hydraulic fractures, *Acta Mech. Solida Sin.* **22**, 443–452.

4. Yao Y., Gosavi S.V., Searles K.H., Ellison T.K. (2010) Cohesive fracture mechanics based analysis to model ductile rock fracture, *44th US Rock Mechanics Symposium and 5th US-Canada Rock Mechanics Symposium*, June 27–30, Salt Lake City, Utah, 83457–13.

5. Gu Y.T., Wang Q.X., Lam K.Y., Dai K.Y. (2007) A pseudo-elastic local meshless method for analysis of material nonlinear problems in solids, *Eng. Anal. Bound. Elem.* **31**, 771–782.

6. Zhang Q.H. (2011) Theoretical analysis of numerical integration in Galerkin meshless methods, *Bit Numer. Math.* **51**, 459–480.

7. Kennett D.J., Timme S., Angulo J., Badcock K.J. (2013) An implicit meshless method for application in computational fluid dynamics, *Int. J. Numer. Methods Fluids* **71**, 1007–1028.

8. Metsis P., Lantzounis N., Papadrakakis M. (2015) A new hierarchical partition of unity formulation of EFG meshless methods, *Comput. Methods Appl. Mech. Eng.* **283**, 782–805.

9. Belytschko T., Lu Y.Y., Gu L. (1994) Element-free Galerkin method, *Int. J. Numer. Methods Fluids* **37**, 229–256.

10. Belytschko T., Gu L., Lu Y.Y. (1994) Fracture and crack growth by element-free Galerkin methods, *Model. Simul. Mater. Sci. Eng.* **115**, 277–286.

11. Belytschko T., Lu Y.Y., Gu L. (1995) Crack propagation by element-free Galerkin methods, *Eng. Fract. Mech.* **51**, 295–315.

12. Klysl P., Belytschko T. (1999) The element-free Galerkin method for dynamic propagation of arbitrary 3-D cracks, *Int. J. Numer. Methods Eng.* **44**, 767–800.

13. Kou X.D., Zhou W.Y. (2000) Using element-free method to trace crack propagation, *Chinese J. Rock Mech. Eng.* **19**, 18–23.

14. Hu Y.J., Zhou W.Y., Lin P. (2003) Application of EFG method to three-dimensional fracture mechanics, *Rock Soil Mech.* **24**, 21–24.

15. Meng W.Y., Zhuo J.S. (2005) Displacement model for meshless method and its application to analysis of fracture problem, *Chinese J. Geotech. Eng.* **27**, 828–831.

16. Shen M. (2006) *Simulation of hydraulic fracturing of rock mass using element-free method*, Zhejiang University, Hangzhou.

17. Oliae M.N., Pak A., Soga K. (2014) A coupled hydro-mechanical analysis for prediction of hydraulic fracture propagation in saturated porous media using EFG mesh-less method, *Comput. Geotech.* **55**, 254–266.

18. Fries T.P., Schaetzer M., Weber N. (2014) XFEM simulation of hydraulic fracturing 3D in D with emphasis on stress intensity factors, *11th World Congress on Computational Mechanics (WCCM)/5th European Conference on Computational Mechanics (ECCM)/6th European Conference on Computational Fluid Dynamics (ECFD) II-IV*, 3282–3293.

19. Shi L.Y., Yu T.T., Tinh Q.B. (2015) Numerical modelling of hydraulic fracturing in rock mass by XFEM, *Soil Mech. Found. Eng.* **52**, 74–83.

20. Deng G.Z., Wang S.L., Huang B.X. (2004) Research on behavior character of crack development induced by hydraulic fracturing in coal-rock mass, *Chinese J. Rock Mech. Eng.* **23**, 3489–3493.

21. Lv Y.M., Li Z.P., Tang D.Z., Xu H., Chen X.Z. (2016) Permeability variation models for unsaturated coalbed methane reservoirs, *Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles* **71**, 32.

22. Ding S.D., Sun L.M. (1997) *Fracture mechanic*, China Machine Press, Beijing.

23. Trunk B., Schober G., Helbling A.K., Wittmann F.H. (1999) Fracture mechanics parameters of autoclaved aerated concrete, *Cem. Concr. Res.* **29**, 855–859.

24. Ioannides A.M., Sengupta S. (2003) Crack propagation in Portland cement concrete beams – Implications for pavement design, *82nd Annual Meeting of the Transportation-Research-Board/Transportation Record-Series* **1853**, 110–117.

25. Adachi J., Siebrits E., Peirce A., Desroches J. (2007) Computer simulation of hydraulic fractures, *Int. J. Rock Mech. Min. Sci.* **44**, 739–757.

26. Zhao Y.L., Zhang L.H., Feng G.Q., Zhang B.N., Kang B. (2016) Performance analysis of fractured wells with stimulated reservoir volume in coal seam reservoirs, *Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles* **71**, 8.

27. Li Z.L., Ren Q.W., Wang Y.H. (2005) Formula for water pressure distribution in rock or concrete fractures formed by hydraulic fracturing, *J. Hydraul. Eng.* **36**, 656–661.

28. Guo W., Yu R.Z., Zhang X.W., Hu Z.M. (2018) Physical and mathematical modeling of gas production in shale matrix, *Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles* **73**, 12.