Regular Article
Meshless method simulation and experimental investigation of crack propagation of CBM hydraulic fracturing
^{1}
School of Mechanical & Electronic Information, China University of Geosciences (Wuhan), Wuhan 430074, China
^{2}
China Geological Equipment Group Co., Ltd., Beijing 100102, China
^{*} Corresponding author: wangyodan@163.com
Received:
15
March
2018
Accepted:
3
October
2018
For simulating CoalBed Methane (CBM) hydraulic fracturing using 3D 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 fluidstructure interaction problem. The geometric and mathematical models established above were described by 3D meshless Galerkin (EFG, ElementFree 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 MicroCT, 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.
© G. Wen et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 3D (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, FRANC3D, HYFRANC3D, or ABAQUS [1–4] are popular numerical simulation software based on the finiteelement 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 mechanicsbased 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 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 leastsquares 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 nearcracktip singularity and that crack propagation can be easily modeled since it requires hardly any remeshing [10, 11]. More simulating examples for arbitrary threedimensional 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 quasistatic 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 hydromechanical analysis for prediction of hydraulic fracture propagation in saturated porous media using EFG meshless method successfully simulates the leakoff 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 brittleelasticplastic 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 insitu 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 “springsplint” model illustrated in Figure 1. K_{1} and K_{2} are the equivalent spring coefficients of the coal rock between roof and floor, w is the crack width, and α is the opening angle.
Fig. 1. The approximate “springsplint” model for CBM hydraulic fracturing. 
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.
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 springsplint 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{2{K}_{\mathrm{IC}}^{2}}{E{f}_{\mathrm{t}}}.$$(1)
Similarly, since the elasto plastic deformation is dominated by the plane strain, the critical opening displacement is$${w}_{\mathrm{c}}=\frac{2\left(1{v}^{2}\right){K}_{\mathrm{IC}}^{2}}{E{f}_{\mathrm{t}}},$$(2)where K_{IC} is the fracture toughness, f_{t} is the tensile strength, E is the elastic modulus, and ν is the Poisson ratio of coal petrography.
Then with the COD criterion of W ≥ 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 constantheight rectangular section (KGD mode) can be employed to study fracture propagation based on fracture Mode I of coal rock illustrated in Figure 2. The socalled 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 crosssection, 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].
Fig. 2. The hydraulic fracturing shape of KGD model for coal rock. 
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 fluidstructure 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{\mathrm{d}m}{\mathrm{d}t}=\rho \frac{\mathrm{d}}{\mathrm{d}t}\underset{v}{\int}\mathrm{d}v+\rho \underset{\mathrm{\Omega}}{\oint}\left(u\cdot n\right)\mathrm{d}A=0,$$(3)where m is the mass of control volume, ρ is the mass density of control volume, v is the volume of control volume, Ω 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{\mathrm{d}M}{\mathrm{d}t}=\mathrm{\rho}\frac{\mathrm{d}}{\mathrm{d}t}\underset{v}{\int}u\mathrm{d}v+\rho \underset{\mathrm{\Omega}}{\oint}u\left(u\cdot n\right)\mathrm{d}A=\sum F,$$(4)where M is the momentum of control volume and ∑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:

The fluid volume should be incompressible;

The fluid belongs to Newtonian fluid;

The fluid loss should be ignored; and

The water flow should obey the rules of onedimensional (1D) laminar flow.
The velocity distribution of 1D 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{4{y}^{2}}{{w}_{x}}\right),$$(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.
Fig. 3. The control volume and fluid velocity distribution. 
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 constantheight 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,$$(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 higherorder differentials illustrated in Equation (7):$$\left(\frac{1}{3}\frac{1}{9}\frac{\rho}{\mu}{w}_{x}\frac{\mathrm{\partial}{w}_{x}}{\mathrm{\partial}t}\right)\frac{\mathrm{\partial}{p}_{x}}{\mathrm{\partial}x}\frac{{p}_{x}}{{w}_{x}}\frac{\mathrm{\partial}{\mathrm{w}}_{\mathrm{x}}}{\mathrm{\partial}\mathrm{x}}=0,$$(7)where w_{ x } is the fracture width of section x, μ 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{\mathrm{\partial}{q}_{x}}{\mathrm{\partial}x}+\frac{x}{{a}^{2}{x}^{2}}\frac{1}{\frac{1}{3}\frac{1}{9}\frac{\rho}{\mu}{w}_{\mathrm{o}}{w}_{0}^{\mathrm{\prime}}\left(1\frac{{x}^{2}}{{a}^{2}}\right)\frac{1}{9}\frac{\rho}{\mu}{w}_{0}^{2}{a\mathrm{\prime}}^{\frac{{x}^{2}}{{a}^{3}}}}{p}_{x}=0,$$(8)where w_{0} is the fracture width, ${w}_{0}^{\mathrm{\prime}}$ is the change rate of fracture width with time, and a′ is the change rate of fracture length with time.
For the brittleelasticplastic fracture of coal rock, the change rate of fracture length with time can be set as zero, i.e., α′ = 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}^{\mathrm{\prime}}}{\nu}\left(1\frac{{x}^{2}}{{a}^{2}}\right)}\right]}^{\frac{3}{2}},$$(9)where ν is the kinematic viscosity coefficient of fluid and 1/ν = ρ/μ.
2.4 The 2D processing method for 3D 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 3D to 2D process in which complex 3D fracture problem are transformed to easier 2D plane fracture problem.
Fig. 4. The slice along bedding direction of coal rock. 
Before the 3D to 2D process, the following operations should be considered:

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.

As the thickness of selected bedding slice along bedding direction is far less than that of coal reservoir, i.e., m << 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.
As shown in Figures 5 and 6, the model is very explicit after the 3D to 2D process on both mechanics and geometric sides.
Fig. 5. The 3D to 2D process for mechanics model. 
Fig. 6. The 3D to 2D process for geometric model. 
According to the principle of geometric transformation, the geometric transformation of 3D data points can be expressed as$$\left({x}^{\mathrm{\prime}}{y}^{\mathrm{\prime}}{z}^{\mathrm{\prime}}1\right)=\left(xyz1\right){\mathit{A}}_{\mathit{H}},$$(10) $$\mathbf{\Omega}\left(\mathbf{\Gamma}\right)\stackrel{{\mathit{A}}_{\mathit{H}}}{\Rightarrow}{\mathbf{\Omega}}_{\mathrm{xoy}}\left({\mathbf{\Gamma}}_{\mathrm{xoy}}\right).$$(11)
The solution domain is$${\mathrm{\Omega}}_{i}=\{\begin{array}{l}{x}_{i}={R}_{i}\mathrm{cos}{\alpha}_{i}\\ {y}_{i}={R}_{i}\mathrm{sin}{\alpha}_{i}\\ {z}_{i}={R}_{i}\mathrm{cot}{\beta}_{\mathrm{i}}\end{array}\left(i=0,1,2,\cdots ,n\right).$$(12)
The boundary condition is$${\mathrm{\Gamma}}_{\mathrm{xoy}}=\{\begin{array}{l}{R}_{i}=r\\ {R}_{i}=R\end{array}.$$(13)
The transformation matrix of vertical view is$${\mathit{A}}_{\mathit{H}}=\left[\begin{array}{llll}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& n& 1\end{array}\right].$$(14)
At the step i, the fracture propagation length will be$${l}_{i}=\sqrt{{\left({x}_{i}{x}_{0}\right)}^{2}+{\left({y}_{i}{y}_{0}\right)}^{2}+{\left({z}_{i}{z}_{0}\right)}^{2}}\left(i=0,1,2,\cdots ,n\right).$$(15)
As the fracture initiation always starts at the wall of wellbore, the initial solution domain is$${\mathrm{\Omega}}_{0}=\{\begin{array}{l}{x}_{0}=r\mathrm{cos}{\alpha}_{0}\\ {y}_{0}=r\mathrm{sin}{\alpha}_{0}\\ {z}_{0}=r\mathrm{cot}{\beta}_{0}\end{array}.$$(16)
In Equations (12), (13), (15) and (16), x_{ i } is xcoordinate of point at the solution domain, y_{ i } is ycoordinate of point, z_{ i } is zcoordinate 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, β_{ i } is the dip angle of sought point, α_{ i } is the azimuth angle of sought point, x_{0} is the xcoordinate of point of fracture initiation, y_{0} is the ycoordinate of point of fracture initiation, z_{0} is the zcoordinate of point of fracture initiation, β_{0} is the dip angle of point of fracture initiation, α_{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 2D fracture graph in bedding slice section will be transformed back to the 3D fracture propagation model by multiplying an inverse matrix $\overline{\mathit{A}}$ expressed in Equations (17) and (18) $${\mathrm{\Omega}}_{\mathrm{xoy}}\left({\mathrm{\Gamma}}_{\mathrm{xoy}}\right)\stackrel{\overline{A}}{\Rightarrow}\mathrm{\Omega}\left(\mathrm{\Gamma}\right),$$(17) $$\stackrel{\u0305}{\mathit{A}}=\overline{{\mathit{A}}_{\mathit{H}}}.$$(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.
3.1 The total program flow
The total program flowchart is illustrated in Figure 7. Visual Basic is used for the interactive window of input and output (shown in Fig. 13) while FORTRAN is for the fracture propagation program based on EFG method in the dashed box running in the background.
Fig. 7. The total program flowchart of numerical simulation. 
3.2 Program based on EFG method
The whole numerical simulation program of hydraulic fracturing propagation based on meshless method (including EFG method) is normally divided into three parts shown in Figure 8: the meshless method executive process, meshless method program module, and external program module. Both meshless method program module and external program module are developed in terms of the execution process of the program.
Fig. 8. The flowchart for overall structure of coal rock fracture propagation in hydraulic fracturing based on meshless methods. 
Meshless program module is used to describe the solving domain mathematically based on meshless method to establish meshless method basic equations of hydraulic fracture propagation for coal rock sample. It consists of general submodules: InPut for the external data input, GaussCoefficient for the setting of Gauss integral point coefficient, CellGaussPoints for the setting of all Gauss points in the background element and calculation of Jacob ratio, SuPPortDomain for the determination of single interpolation point support area, MLS_ShaPeFunc_2D for the calculation of shape function and derivative of single interpolation points, PointstiffessMatrix for the calculation of stiffness matrix of single Gauss integral point, EssentialBC for the implementation of penalty function 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 addin 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.
Fig. 9. The program flowchart for computing hydraulic pressure inside the fracture. 
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, XlnsertRNode XNODE, and XLOAD are involved during this process.
Fig. 10. The program flowchart for hydraulic fracture propagation based on Hillerborg model. 
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 × 10^{5} MPa, elastic modulus perpendicular to bedding direction is E = 4.3 × 10^{5} MPa, Poisson ratio is γ = 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 Scanner (MicroCT) 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.
Fig. 11. The results of electronic microscope scanning of existing fractures before hydraulic fracturing experiment. 
The existing fracture parameters of coal rock sample.
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.
Fig. 12. The geometric model of existing fractures in coal rock sample. 
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 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 × 10^{−6}m^{2}/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.
Fig. 13. The display interface and results for water injection pressure of 3 MPa in coal rock sample. 
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 (Nos. 1 and 2) and the new fracture (No. 3) are shown in Tables 2–4 respectively.
The No. 1 fracture propagation parameters.
The No. 2 fracture propagation parameters.
The No. 3 fracture propagation parameters.
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.
Fig. 14. The relationship of fracture parameter of coal rock and water injection pressure (triaxial pressure). 
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 concluded as follows.

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.

The expansion of hydraulic fractures is of inheritance since it is mainly based on the existing fractures and in the form of straight line without distortion normally.

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.

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

Acquiring the properties of coal rock sample using ESEM and MicroCT before water injection.

Hydraulic fracturing using customized JHLSI experiment instrument.

Scanning by ESEM and MicroCT again after water injection.

Comparing the scanned results before and after water injection; and

Analyzing the fracture propagation rules.
Fig. 15. The hydraulic fracturing experimental procedure of coal rock. 
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 MicroCT 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.
Fig. 16. The magnified results of electronic microscope scanning of the fractures after the 3 MPa water injection experiment. 
Figure 17 illustrates the internal characteristics of three fractures scanned by MicroCT after 3 MPa and 6 MPa water injection experiments.
Fig. 17. The results of MicroCT scanning of the fractures after the 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.
Fig. 18. The geometric model of the fractures in coal rock sample after water injection. 
The fracture parameters of coal rock sample after 3 MPa hydraulic fracturing.
The fracture parameters of coal rock sample after 6 MPa hydraulic fracturing.
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.
Fig. 19. The relationship of fracture parameters and water injection pressure (triaxial pressure) of coal rock. 
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.

The expansion of hydraulic fractures of the coal rock will mainly initiate from the existing fractures on the 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.

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.

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 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.

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 springsplint 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 MicroCT 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 springsplint 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:

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.

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.

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.

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.

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.

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
 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. [Google Scholar]
 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. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 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 USCanada Rock Mechanics Symposium, June 27–30, Salt Lake City, Utah, 8345713. [Google Scholar]
 Gu Y.T., Wang Q.X., Lam K.Y., Dai K.Y. (2007) A pseudoelastic local meshless method for analysis of material nonlinear problems in solids, Eng. Anal. Bound. Elem. 31, 771–782. [Google Scholar]
 Zhang Q.H. (2011) Theoretical analysis of numerical integration in Galerkin meshless methods, Bit Numer. Math. 51, 459–480. [CrossRef] [Google Scholar]
 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. [Google Scholar]
 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. [Google Scholar]
 Belytschko T., Lu Y.Y., Gu L. (1994) Elementfree Galerkin method, Int. J. Numer. Methods Fluids 37, 229–256. [Google Scholar]
 Belytschko T., Gu L., Lu Y.Y. (1994) Fracture and crack growth by elementfree Galerkinmethods, Model. Simul. Mater. Sci. Eng. 115, 277–286. [Google Scholar]
 Belytschko T., Lu Y.Y., Gu L. (1995) Crack propagation by elementfree Galerkin methods, Eng. Fract. Mech. 51, 295–315. [Google Scholar]
 Klysl P., Belytschko T. (1999) The elementfree Galerkin method for dynamic propagation of arbitrary 3D cracks, Int. J. Numer. Methods Eng. 44, 767–800. [Google Scholar]
 Kou X.D., Zhou W.Y. (2000) Using elementfree method to trace crack propagation, Chinese J. Rock Mech. Eng. 19, 18–23. [Google Scholar]
 Hu Y.J., Zhou W.Y., Lin P. (2003) Application of EFG method to threedimensional fracture mechanics, Rock Soil Mech. 24, 21–24. [Google Scholar]
 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. [Google Scholar]
 Shen M. (2006) Simulation of hydraulic fracturing of rock mass using elementfree method, Zhejiang University, Hangzhou. [Google Scholar]
 Oliae M.N., Pak A., Soga K. (2014) A coupled hydromechanical analysis for prediction of hydraulic fracture propagation in saturated porous media using EFG meshless method, Comput. Geotech. 55, 254–266. [Google Scholar]
 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. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Deng G.Z., Wang S.L., Huang B.X. (2004) Research on behavior character of crack development induced by hydraulic fracturing in coalrock mass, Chinese J. Rock Mech. Eng. 23, 3489–3493. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Ding S.D., Sun L.M. (1997) Fracture mechanic, China Machine Press, Beijing. [Google Scholar]
 Trunk B., Schober G., Helbling A.K., Wittmann F.H. (1999) Fracture mechanics parameters of autoclaved aerated concrete, Cem. Concr. Res. 29, 855–859. [Google Scholar]
 Ioannides A.M., Sengupta S. (2003) Crack propagation in Portland cement concrete beams – Implications for pavement design, 82nd Annual Meeting of the TransportationResearchBoard/Transportation RecordSeries 1853, 110–117. [CrossRef] [Google Scholar]
 Adachi J., Siebrits E., Peirce A., Desroches J. (2007) Computer simulation of hydraulic fractures, Int. J. Rock Mech. Min. Sci. 44, 739–757. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 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. [Google Scholar]
 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. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1. The approximate “springsplint” model for CBM hydraulic fracturing. 

In the text 
Fig. 2. The hydraulic fracturing shape of KGD model for coal rock. 

In the text 
Fig. 3. The control volume and fluid velocity distribution. 

In the text 
Fig. 4. The slice along bedding direction of coal rock. 

In the text 
Fig. 5. The 3D to 2D process for mechanics model. 

In the text 
Fig. 6. The 3D to 2D process for geometric model. 

In the text 
Fig. 7. The total program flowchart of numerical simulation. 

In the text 
Fig. 8. The flowchart for overall structure of coal rock fracture propagation in hydraulic fracturing based on meshless methods. 

In the text 
Fig. 9. The program flowchart for computing hydraulic pressure inside the fracture. 

In the text 
Fig. 10. The program flowchart for hydraulic fracture propagation based on Hillerborg model. 

In the text 
Fig. 11. The results of electronic microscope scanning of existing fractures before hydraulic fracturing experiment. 

In the text 
Fig. 12. The geometric model of existing fractures in coal rock sample. 

In the text 
Fig. 13. The display interface and results for water injection pressure of 3 MPa in coal rock sample. 

In the text 
Fig. 14. The relationship of fracture parameter of coal rock and water injection pressure (triaxial pressure). 

In the text 
Fig. 15. The hydraulic fracturing experimental procedure of coal rock. 

In the text 
Fig. 16. The magnified results of electronic microscope scanning of the fractures after the 3 MPa water injection experiment. 

In the text 
Fig. 17. The results of MicroCT scanning of the fractures after the 3 MPa and 6 MPa water injection experiments. 

In the text 
Fig. 18. The geometric model of the fractures in coal rock sample after water injection. 

In the text 
Fig. 19. The relationship of fracture parameters and water injection pressure (triaxial pressure) of coal rock. 

In the text 