Numerical investigation on rules of fracture propagation during hydraulic fracturing in heterogeneous coal-rock mass

2.1. Governing equation for coal-rock mass deformation

The governing equations could be introduced by considering pore water pressure with Biot constitutive relation [14], and further adopting of strain-displacement relation and equilibration equation. For elastic porous medium with single-phase saturated under an isothermal state, the governing equation for coal-rock mass deformation is given by [32, 33]:

where $G$ and $v$ are the shear modulus and Poisson’s ratio of coal-rock mass, respectively; $u_j$ is the displacement of $j$th direction; $f_i$ denotes the volumetric body force of $i$th direction; $\alpha$ represents Biot coefficient; the footnote [$,i$] is the differentiation in regard to the coordinate $x_i$; and $p$ denotes pore pressure (negative for suction).

As can be seen from Eq. (1), it expresses the mechanical equilibrium of porous elastic medium where the effect of pore pressure on coal-rock mass deformation is accounted.

2.2. Governing equation for water seepage

The water seepage of saturated medium satisfies mass conservation equation, and further utilizing of Darcy seepage law for momentum equation [33, 34], the governing equation for water seepage is given by:

where $\varphi$ and $K$ are the porosity and bulk modulus of coal-rock mass, respectively; ${\beta }_w$, ${\mu }_w$ and ${\rho }_w$ are bulk modulus, dynamic viscosity coefficient and density of water, respectively; ${\epsilon }_v$ and $k$ are volumetric strain and permeability coefficient of coal-rock mass; $t$ is the time; $g$ denotes the gravitational acceleration, and $z$ represent the elevation.

2.3. Damage and permeability evolution equations

As mentioned, the elastic damage constitutive relationship would be employed to describe the damage evolution of coal-rock mass. As seen from Fig. 1, the relationship between strain and stress is assumed to be linear unless a damage threshold is satisfied. The coal-rock mass damage would initiate if its stress state reaches the Mohr-Coulomb criterion (abbreviated as C-C in short) or tensile stress criterion (abbreviated as T-C in short). In numerical algorithm, the T-C is always prior to C-C, namely, T-C is adopted first to ascertain whether tensile damage initiates in coal-rock mass or not, if the it doesn’t damage according to T-C, then the damage of coal-rock mass would be checked based on C-C. The mathematical functions of T-C and C-C can be given by [35]:

where $\phi$, $f_c$ and $f_t$ represent the internal frictional angle, compressive strength, and tensile strength of coal-rock mass, respectively; $F_1$ and $F_2$ denote damage initiation functions for T-C and C-C, respectively. It should be noted that Eqs. (3), (4) operate on effective stress space.

Fig. 1Elastic damage constitutive relations for coal-rock mass under uniaxial conditions

Based on elastic damage mechanics, as the evolution of damage, the coal-rock mass mechanical parameter, such as Young’s modulus $E$, would linearly degrade in accordance with following function [31]:

where $E_0$ denotes the initial Young’s modulus of coal-rock mass before damage initiates; while $D$ denotes the damage variable, which ranges from 0 to 1.

The constitutive relations of coal-rock mass under uniaxial tension and compressive conditions can be deduced from Fig. 1. If the coal-rock mass is under tension damage condition, its damage variable $D$ can be calculated by [33]:

where ${\epsilon }_{t0}$ represents the threshold strain of coal-rock mass when tensile damage initiates; $n$ is called the constitutive coefficient which equals to 2.0.

In analogy with the damage evolution under uniaxial tension condition, when the coal-rock mass satisfies the C-C under uniaxial compression condition, its damage variable $D$ can be calculated by [33]:

where ${\epsilon }_{c0}$ represents the threshold strain of coal-rock mass when shear damage initiates. With the degradation evolution of mechanical strength or stiffness of coal-rock mass, it would experience dilatancy, which would seriously change the hydraulic properties. The coal-rock mass porosity, $\varphi$, which is strongly related to its stress conditions, can be calculated as [34]:

where ${\varphi }_0$ represents the original porosity without stress; ${\varphi }_r$ represents the residual porosity in infinite stress state; ${\sigma }_e$ represents the effective stress calculated by ${\sigma }_e=\alpha p+{\sigma }_{ij}{\delta }_{ij}$/3, where ${\sigma }_{ij}$ ($i$, $j=$ 1, …, 3) represents the corresponding stress tensor; and ${\alpha }_{\phi }$ is called stress sensitivity coefficient with value being 5.0×10^-8 Pa^-1.

Furthermore, the permeability of coal-rock mass would also evolve during damage accumulation. Based on previous studies [34], the permeability is strongly dependent on the porosity of coal-rock mass, and then the permeability could be deduced according to the following exponential function as:

where $\beta$ is called index coefficient with value being 3.0; ${\alpha }_D$ is called damage-permeability coefficient with value being 5.0; $k_0$ denotes the permeability without stress; and $D$ could be calculated by Eqs. (6), (7).

2.4. Assignment of material heterogeneity

As the realistic coal-rock mass contains dispersed cementing materials and various mineral particles, which cause its mechanical properties to exhibit characteristics of heterogeneity. To characterize such heterogeneity of the coal-rock mass, assuming that it is composed of many mesoscopic elements, whose physical and mechanical parameters conform to a given Weibull distribution whose probability density function is [19]:

where $z$ represents the element’s mechanical parameter (e.g. elastic modulus), $\zeta$ represents the average value for all elements, and $m$ is called homogeneity index of element parameter $z$. According to characteristics of Weibull distribution, with the increase of $m$, the mechanical parameters of more elements will be closer to value of $\zeta$, and thus generating a more homogeneous coal-rock mass specimen. Previous results indicate that when $m$ is greater than 1000, the rock can be assumed as a homogeneous material [1].

2.5. Numerical algorithm

The aforementioned Eqs. (1)-(10) constitutes the mathematical model for hydraulic fracturing with seepage-damage coupling. The effect of water seepage on coal-rock mass deformation behavior is considered in Eq. (1), while the effect of coal-rock mass deformation on water seepage behavior is considered in Eq. (2). Moreover, the damage-induced evolution of porosity and permeability are considered in Eqs. (8), (9). Eq. (1) and Eq. (2) are both nonlinear, which are impossible to solve theoretically. In this respect, a numerical algorithm for the solution of mathematical model is proposed instead. The procedure of proposed numerical algorithm is as follows:

Step 1: After the geometry model is built, it is then discretized into numerous finite elements, whose mechanical parameters satisfy the Weibull distribution (see Eq. (10)) using a Monte Carlo simulation method coded by Matlab program. Moreover, to ensure the numerical model is in a quasi-static response, the prescribed boundary load on model is divided into several loading increments.

Step 2: For every loading increment, the Comsol program [36] is used to execute a fully coupled calculation according to governing Eqs. (1), (2), and numerical results of all finite elements (such as effective stress and displacement) are obtained for the following step.

Step 3: Substituting obtained effective stress of all finite elements into Eqs. (3), (4), and then the damage states of all elements can be determined. When an element is in damage state, damage variable $D$ at its corresponding node is calculated using Eq. (6) or Eq. (7), and then the elastic modulus $E$ and permeability $k$ at its node can be calculated using Eq. (5) and Eq. (9), thus forming a new numerical model with updated mechanical parameters. Furthermore, the second and third steps are cycled to judge the damage state of each element. The next loading increment is applied when a specific convergence criterion is satisfied: $\left|\right|\left(D_{j+1}-D_j\right)$/$D_{j+1}{\left|\right|}_{\infty }$ < 1×10^-5, where $j$ represents the solution step during numerical calculation, and $D$ denotes the damage variable matrix.

In a word, the above procedure of the numerical algorithm is coded into a Matlab program combined with Comsol program to solve mathematical model for hydraulic fracturing with seepage-damage coupling.

3.1. Influence of coal-rock mass heterogeneity on fracture propagation

To investigate how coal-rock mass heterogeneity influences the fracture propagation during hydraulic fracturing, two numerical models with homogeneity indexes of 1000 and 20 are established. Fig. 2 presents the geometry and corresponding boundary conditions for two numerical models. Related parameters adopted are presented in Table 1. The vertical stress ${\sigma }_v$ and horizontal stress ${\sigma }_h$ applied on boundaries are 15.45 MPa and 10.30 MPa, respectively. The initial water pressure in calculation domain is set to 0.0 MPa, while the pressurization rate on the wall of fracturing borehole is 0.5 MPa/step.

To realistically simulate the hydraulic fracturing process, the calculation scheme is divided into two steps. Firstly, the initial stress state before hydraulic fracturing is calculated by performing a mechanical-only solution. Secondly, a fully coupled seepage-damage solution is executed with the equilibrium state in the first step being the initial condition.

Fig. 2Geometry and corresponding boundary conditions for two numerical models

Fig. 3 illustrates two distribution diagrams for the first principal stress field and fracture propagation pattern with homogeneity indexes being 1000 and 20.

Fig. 3Distribution diagrams of stress field and fracture propagation pattern

a) First principal stress field, $m=$ 1000

b) Fracture propagation pattern, $m=$ 20

As seen from Fig. 3, the stress existed in heterogeneous model with $m$ being 20 exhibit roughness and discontinuity compared with homogeneous model with $m$ being 1000. In other words, coal-rock mass heterogeneity has a significant role in the modification of stress field; meanwhile, the path of fracture propagation in heterogeneous model exhibits a more zig-zag characteristics than that of the homogeneous model. Compared to numerical model with homogeneity index of 1000, there are more randomly distributed defects and grains in model with homogeneity index of 20. Therefore, the mechanical behavior of coal-rock mass at different locations of the numerical model with homogeneity index of 20 are more different. After the external load is applied, local stress concentrations would occur as a result of the different stress or deformation transfer abilities of coal-rock mass constitutes, thus causing the change of fracture propagation pattern and extension of local microcracks.

Moreover, numerical simulation results indicates that if the hydraulic pressure in fracturing borehole is increased to a certain pressure, two nearly symmetrical fractures (denoted by damaged elements) start to initiate. Such critical pressure is called fracture initiation pressure. The fracture initiation pressures for the homogeneous model ($m=$ 1000) and heterogeneous model ($m=$ 20) are 13.5 MPa and 12 MPa, respectively. Namely, with the decrease of homogeneity index, the fracture initiation pressure would also reduce, which conforms to related theoretical analysis and engineering practice [17, 37, 38]. Thus, the characteristics of coal-rock mass heterogeneity should be investigated before the operation of engineering practice of hydraulic fracturing.

3.2. Influence of confining pressure on fracture propagation

To study how confining pressure influences fracture propagation (including fracture propagation pattern, fracture initiation or breakdown pressure [19]), three heterogeneous numerical models with homogeneity index of 20 are used. The horizontal boundary stress ${\sigma }_h$ is kept at 10.30 MPa, while the ratios of vertical stress to horizontal stress (${\sigma }_v$/${\sigma }_h$) are set to 1.5, 1.25 and 1.0, respectively. The model geometry and its corresponding boundary conditions are shown in Fig. 2. Related parameters for hydraulic fracturing simulation are presented in Table 1. The fracture propagation patterns under varying confining pressures are illustrated in Fig. 4.

Fig. 4Fracture propagation patterns under varying confining pressures

a)${\sigma }_v$/${\sigma }_h=$ 1.5

b)${\sigma }_v$/${\sigma }_h=$ 1.25

c)${\sigma }_v$/${\sigma }_h=$ 1.0

As seen from Fig. 4, at a larger confining pressure ratio of ${\sigma }_v$/${\sigma }_h=$ 1.5, two nearly straight fractures propagate approximately orthogonal to direction of minimum principal stress (Fig. 4(a)). With the decrease of confining pressure ratio, the dominant fractures gradually extend, and they start to deviate from orientations of initial fractures; besides, some fractures would branch as hydraulic fracturing continues, and their propagation paths become devious and coarse. Especially when confining pressure ratio is 1.0 (namely, a hydrostatic pressure condition as shown in Fig. 5(c)), the randomly distributed defects and grains in coal-rock mass would significantly affect the path of fracture propagation; in other words, the micro fractures would initiate at the weakest region, causing the hydraulic fractures randomly propagate. In short, based on the results of numerical calculation, it could be concluded that fracture propagation is mainly controlled by confining pressure ratio, the hydraulic-induced fracture would always select the path with least resistance in coal-rock mass; and an irregular fracture path would occur due to the random locations of weak parts; perfectly planar fractures, which are parallel to any principal stress direction, are impossible in realistic engineering practice.

Fig. 5 presents simulation results of fracture breakdown pressure and initiation pressure under varying confining pressures. The breakdown pressures, $p_b$, are 20 MPa and 18.5 MPa for confining pressure ratios of 1.0 and 1.5, respectively, which are consistent with the trend of theoretical analysis calculated by $p_b=\text{3}{\sigma }_h$-${\sigma }_v$+$f_t$[17, 39]. With confining pressure ratio increasing, the fracture initiation or fracture breakdown pressure would decrease. The main reason is that the circumferential stress (positive for tension) at the top or bottom borehole boundary would increase with confining pressure ratio increasing, and then initial fractures initiate at these regions, which provide a path for water flow and further cause the extension or propagation of hydraulic fractures, and thus reducing the breakdown pressure for coal-rock mass.

Fig. 5Fracture breakdown and initiation pressure under varying confining pressures

3.3. Influence of beforehand hydraulic slotting on fracture propagation

Hydraulic fracturing or hydraulic slotting each has its own advantage of enhancing coalbed methane extraction. On basis of the two measures, scholars proposed the oriented fracturing technology with beforehand hydraulic slotting recently [40, 41]. The key of the oriented fracturing technology is to ascertain the influence rules of beforehand hydraulic slotting on fracture propagation. The influence of two major factors (including slotting length and principal stress difference) on fracture propagation are studied. The model geometry and its corresponding boundary conditions are presented in Fig. 6. Related parameters for hydraulic fracturing simulation are presented in Table 1. The horizontal stress ${\sigma }_h$ on boundary is 10.30 MPa. The calculation scheme is also divided into two steps as mentioned in section 3.1. The initial water pressure in computation domain is set to 0.0 MPa, while the pressurization rate in fracturing borehole is 0.5 MPa/step.

Fig. 6Model geometry and boundary conditions

3.3.1. Influence of slotting length on fracture propagation

In present numerical models, the angle between slotting and horizontal boundary stress is kept at 30°, the vertical stress ${\sigma }_v$ on boundary is 15.45 MPa, the slotting lengths are set to 25 mm, 50 mm and 125 mm for cases 1-3, respectively. The damage regions at the initial equilibrium step of case 1 are illustrated in Fig. 7(a), in which the tensile and shear damage are denoted by “○” and “△”, respectively. As seen from Fig. 7(a), regions with shear damage mainly occur at the slotting tips, while the regions with tensile damage occur at slotting tips, bottom of left slotting and upper of right slotting. With the evolution of damage, the mechanical strength of damage regions would reduce, which can be assumed as weak parts in coal-rock mass, thus affecting the fracture propagation pattern and fracture initiation pressure. For comparison and validation, the yield regions simulated by Phase² software [42] are shown in Fig. 7(b). As illustrated in Fig. 7, the damage regions are consistent with the yield regions. However, as the numerical model adopted in present study has considered the coal-rock mass heterogeneity, the stress contours present discontinuity and roughness compared with the homogeneous model simulated by Phase² software.

Fig. 7Damage and yield regions at the initial equilibrium step

a) Damage regions

b) Yield regions using Phase2 software

Fig. 8 presents the fracture propagation processes for different cases. As shown in Figs. 8(a), (b), the fracture propagation patterns of cases 1 and 2 are similar to the branch cracks propagation in compressive-shear stress state in references [43-45]. However, the differences are that the locations of fracture initiation are at a distance from slotting tips, and the directions of fracture initiation are approximately perpendicular to slotting directions.

As hydraulic fractures propagate, the propagation directions of fracture wings deviate to directions of vertical principal stresses. The fracture propagation pattern presents a “$S$” type, which is similar to the experimental results in references [43-45]. When the slotting length is 125 mm (Fig. 8(c)), due to the coal-rock mass heterogeneity and local stress concentrations at slotting location, more tensile damage regions would exist at the slotting tips and upper of right slotting, thus providing paths for water flowing, and then the fracture propagation pattern become more complex, where two primary branch cracks are generated at the slotting tips and upper of right slotting, causing the increase of the fracture breakdown pressure.

In brief, when the angle between hydraulic slotting and horizontal stress is 30°, and ratio of vertical stress to horizontal stress ${\sigma }_v$/${\sigma }_h=$ 1.5, hydraulic-induced fractures don’t directly propagate along the orientation of pre-existing hydraulic slotting but initiate at a distance from slotting tips, and they finally deviate to the orientation of maximum principal stress; the larger length of hydraulic slotting, the stronger tendency of fracture deviation, and the more difficult for oriented fracturing.

Fig. 8Fracture propagation processes for different cases

a) Case 1

b) Case 2

c) Case 3

3.3.2. Influence of principal stress difference on fracture propagation

In present numerical models, the angle between slotting and horizontal boundary stress is set to be 30°, and the slotting length is kept at 50 mm. To investigate influence of principal stress difference on fracture propagation, three cases are chosen, whose ratios of vertical stress to horizontal stress are 3, 1.5 and 1.25, respectively. The fracture propagation patterns with varying principal stress differences are illustrated in Fig. 9. The angle between the directions of hydraulic slotting and fracture propagation are called deviation angle.

As shown in Fig. 9, with the decrease of ratio of vertical stress to horizontal stress (${\sigma }_v$/${\sigma }_h$), the deviation angle between direction of hydraulic fracture and hydraulic slotting gradually reduces. When ${\sigma }_v$/${\sigma }_h$ is less than 1.25, the deviation angle is close to 0°. In other words, when the principal stress difference is lower than 2.58 MPa, the hydraulic fracture would propagate along the pre-existing hydraulic slotting. Theoretically, when the difference of principal stress is less than a critical value, the tensile damage regions mainly occur at the slotting tips other than other locations, which are beneficial for the fracture propagation along the pre-existing hydraulic slotting.

On basis of numerical results in sections 3.3.1 and 3.3.2, it can be inferred that the direction of fracture propagation with beforehand hydraulic slotting is mainly controlled by principal stress difference, the larger of principal stress difference, the more difficult of the oriented fracturing with beforehand hydraulic slotting. Thus, in practical engineering, a preliminary investigation of initial ground stress and reasonable design of hydraulic slotting parameters are required to achieve the goal of oriented fracturing with beforehand hydraulic slotting.

Fig. 9Fracture propagation patterns for varying principal stress differences

a)${\sigma }_v$/${\sigma }_h=$ 3

b)${\sigma }_v$/${\sigma }_h=$ 1.5

c)${\sigma }_v$/${\sigma }_h=$ 1.25

3.4. Influence of non-symmetric pressure gradient on fracture propagation

The experimental studies of Bruno [46, 47] indicated that the pore pressure had an influence on fracture propagation during hydraulic fracturing. The experimental method is as follows (Fig. 10). Firstly, three boreholes are drilled in a rock sample, and a constant hydraulic pressure is kept in borehole B2. Then, oil is gradually injected into borehole B1 to increase pressure to induce fractures, and the borehole B3 is not loaded by hydraulic pressure. The experimental result had illustrated that pore pressure had an influence on fracture propagation. However, the fracture initiation and propagation processes are not observed. Thus, the influence mechanism of non-symmetric pressure gradient on hydraulic fracturing is not understood fully and needs further study.

Fig. 10Hydraulic fracturing diagram of rock sample with three boreholes

Motivated by the experiment designed by Bruno, a numerical calculation model is built according to the experiment conducted as shown in Fig. 10. The three boreholes, each with a diameter of 4mm, are named B1, B2, B3. A small pre-existing crack is designed near the bottom of borehole B1. No initial stress is loaded on the model boundary during numerical simulation. The initial hydraulic pressure applied on borehole B1 is 1.0 MPa, and is increased at 0.5 MPa/step. The hydraulic pressure in borehole B2 is kept at 1.5 MPa or 0.8 MPa, and borehole B3 is free from hydraulic pressure. The related parameters used for hydraulic fracturing simulation are shown in Table 2 [46].

Figs. 11, 12 present the processes of fracture propagation driven by hydraulic pressure. As shown in Figs. 11(a)-12(a), the initial hydraulic fracture initiates from the pre-existing crack. When the length of hydraulic fracture reaches to a critical extent, the hydraulic fracture doesn’t extend downwards along pre-existing crack but deviates to borehole B2; meanwhile, at the top boundary of borehole B1, another hydraulic fracture which is symmetrical to the pre-existing crack is observed (Figs. 11(b)-12(b)). The fracture propagation patterns indicate that the fracture propagation path is not regular and fracture surface is coarse, which has a good similarity with the results observed in previous section 3.2.

Fig. 11Fracture propagation process (pressure in B2 is 0.8 MPa)

a) Fracture initiation stage

b) Fracture stable propagation stage

c) Final stage

d) Experimental result [46]

Moreover, as the distance between hydraulic-induced fracture and borehole B2 decreases, the influence of pore water pressure on fracture propagation becomes more obvious. As hydraulic pressure in borehole B2 increases, the hydraulic-induced fracture has a more obvious tendency of deviation to borehole B2 (Figs. 11(c)-12(c)).

Numerical results show that the hydraulic-induced fractures are mainly composed of elements under tensile damage state, thus the orientation of fracture propagation is nearly orthogonal with maximum tensile stress direction. In theory, when the hydraulic pressure in borehole B2 is unchanged, a maximum tensile stress would be generated in the perpendicular orientation of line connecting boreholes B1 and B2, resulting in that the hydraulic-induced fracture would propagate towards borehole B2. Furthermore, the breakdown pressures are 4.62MPa and 3.84MPa when the hydraulic pressures in borehole B2 are 0.8MPa and 1.5MPa, respectively. The greater the constant hydraulic pressure in borehole B2, the less the breakdown pressure in borehole B1. The aforementioned numerical results are consistent with the experimental results as illustrated in Figs. 11(d)-12(d).

In conclusion, a non-symmetric pressure distribution can be induced on a global scope due to the local pressure gradient between borehole B1 and B2. At the fracture initiation stage (Figs. 11(a)-12(a)), as the tip of pre-existing crack is relatively far apart from borehole B2, and then the non-symmetric pressure distribution has little effect on fracture propagation yet. However, at the fracture stable propagation stage (Figs. 11(b)-12(b)), when the hydraulic-induced fracture approaches to zones with non-symmetric pore pressure distribution, effect of pore pressure on propagation of fracture would gradually increase. The greater the constant hydraulic pressure in borehole B2, the more non-symmetric the pore pressure distribution. Thus, the above results have a beneficial guiding for oriented fracturing. In engineering practice, a reasonable arrangement of multiple boreholes can be arranged, and then hydraulic pressures are applied in several boreholes to generate a suitable non-symmetric pressure gradient, which could make the fracture propagate according to a specific direction.

Fig. 12Fracture propagation process (pressure in B2 is 1.5 MPa)

a) Fracture initiation stage

b) Fracture stable propagation stage

c) Final stage

d) Experimental result [46]