Accepted: May 24, 2021 | Published Online: May 26, 2021

Citation: Li W, Xiaohui G, Bing S, Yunxing C (2021) Poromechanical Modeling of Hydraulic Fracture Propagation. Int J Earth Sci Geophys 7:046.

Copyright: © 2021 Li W, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

A poromechanical modeling of hydraulic fracture propagation is presented. It is established based on a constitutional model of fracture opening. Based on this constitutional model, a permeability tensor is established to model fracture propagation in different direction. Another is included into the computational equation. A full scheme of hydraulic fracturing simulation is presented in this paper.

Keywords

Hydraulic fracture, FEM, Poroelasticity, Damage, Fracture opening

Introduction

Hydraulic fracturing (HF) is an important technique in enhancing the permeability of petroleum and gas reservoirs. The mechanisms of fracture propagation are well understood by the analytical solutions (2D [1-5], P3D [6,7] and PL3D [8,9]) which are mainly dealing with the lubricant flow, elastic displacement of fracture walls and incomplete coupling between the fluid front and the fracture tip. However, the analytical solutions can only output temporal and spatial distribution of hydraulic fracture parameters such as fracture opening, length and fluid pressure on a predefined propagation path. Even some FEM methods [10-13] based on cohesive zone model also need to predefine the propagation path on lined nodes. These models ignore the temporal and spatial distribution of stress, damage, pore fluid pressure around the injection hole, which are significant in monitoring hydraulic fracture zone growth and assessment of permeability enhancement.

The poromechanical model can satisfy both the two needs for assessment of hydraulic fracture propagation. It is based on the coupling analysis between the porous flow and stress and damage utilizing the FEM, which is comprised of direct coupling and load transfer methods [14]. The former is also referred to as strongly coupled, where the final solution of the unknown multiple physical field variables are recovered by solving the simultaneous equations. However, many engineering problems do not satisfy the conditions of strong coupling, especially for some dynamical evolution problems such as damage-induced fracture - where it is difficult to ensure solution convergence. The load transfer method approaches a solution of the unknown field variables by successively solving the multi-physical field equations, during which one field variable is used as an input for the solution of another and repeated through the sequence of couplings until a tolerance for an equilibrium solution is reached. This load transfer, sequential, or leap-frog method represents only weak coupling, and since fluid-driven fractures are always evolving, hence this method is particularly appropriate for the nonlinearities in these problems. However, there are some key points in the poromechanical model to be dealt with, such as how to define the fracture opening using the continuum variables, how to deal with the strain energy loss resulting from hydraulic fracturing, how to control the direction of fracture propagation and how to apply the fluid load to simulate the continuum injection, etc. Previously, some hydraulic fracturing models [15-17] bypassed these issues.

In the present paper, we answer all these questions above. Basically, our poromechanical model is constituted of several components that (1) The fracturing opening is calculated based on damage localization employing a thickness of localization; (2) The strain energy loss resulting from fracturing is compensated by the fluid pressure invasion through poroelasticity coupling; (3) The direction of fracture propagation is controlled by the tensors of hydro-mechanical properties induced by hydraulic fracture opening; (4) The continuous injection with the fracture growth is achieved by a loading scheme of stepwise increased solution duration. As an example, the model is used representing hydraulic fracture propagation in a three-layer reservoir, the morphology of fracture zone and parameters such as length, width and fluid pressure are validated with the analytical solutions. And the model exhibits four stages of fracture propagation, which are fracture nucleation, kinetic propagation, steady propagation and propagation termination and represents the full coupling between fracture tip and fluid front.

Controlling Equations

We define the relations that enable the simulation of fluid-pressure effects on the propagation of a fluid-driven fracture. This involves both the transport of fluid and the mechanism of fracture expansion driven by that fluid.

Flow in porous media

In the reservoir formation, the porous flow satisfies Darcy' Law as

$$q\text{ = }-\text{K}\nabla P\text{ (1)}$$

In which, P is pore fluid pressure, q is the velocity of fluid flow, $\nabla$ is differential operator vector and K is the permeability, as,

$\mathrm{K}=\frac{1}{\mu }\left[\begin{array}{ccc}{K}_{xx}& 0& 0\\ 0& K_{yy}& 0\\ 0& 0& K_{zz}\end{array}\right]\text{ (2)}$

From mass conservation, the pore pressure is determined as

$\frac{d\xi }{dt}+{\nabla }^{T}q=\stackrel{⃛}{q}\text{ (3)}$

Where $\xi$ is the mass of fluid, $\stackrel{⃛}{q}$ is fluid generation rate of a unit solid volume and t is time.

According to the coupling between the compressibility of the solid and fluid [18,19], the differential change of fluid mass, $d\zeta$ can be represented as

$d\zeta \equiv \frac{1}{V_b}\frac{dm}{{\rho }_f}=(\phi C_{pp}dP+\phi C_{f}dP)-\phi C_{pc}d{P}_c\text{ (4)}$

in which ${\rho }_f$ is the mass density of the fluid, $dm$ is the differential increment of fluid mass, $V_b$ is the bulk volume and $\varphi$ is the porosity. Eq. 4 shows that the volume increment of fluid comprises three parts in which $\varphi C_{cc}dP$ and $\varphi C_{f}dP$ are the pore volume increment caused by the fluid pressure increment dP resulting from pore-elasticity and fluid compressibility, respectively, and $\varphi C_{pc}d{P}_c$ is the compressive volume increment of the pore volume caused by an increment of the confining stress $d{P}_c$. The coefficient $C_{pp}$, represents the internal expansibility in volume resulting from the increment of fluid pressure, $C_{pc}$, the internal contractility in volume resulting from the increment of confining stress, and $C_f$, the fluid compressibility resulting from increment of fluid pressure, they are defined as [20,21]

$$C_{pp}=\frac{1}{V_p}{\left(\frac{\partial V_p}{\partial P}\right)}_{Pc}$$, $C_{pc}=\frac{-1}{V_p}{\left(\frac{\partial V_p}{\partial P_c}\right)}_P$, $C_f=\frac{1}{{\rho }_f}{\left(\frac{\partial {\rho }_f}{\partial P}\right)}_{Pc}$     (5)

In which $V_p$ is the pore volume, $P_c={\sigma }_{ii}/3$ is the average stress and the sign convention is with compression defined positive. Further, we can define the storage coefficient of fluid by adding the first and third expressions in Eq.5, as $C=C_{pp}+C_f$. Rewriting Eq. 4 and substituting this into Eq. 3, yields

$\varphi C\frac{dP}{dt}\text{ = }{\nabla }^T(\mathrm{K}\nabla P)+\stackrel{⃛}{q}+\varphi C_{pc}\frac{d{P}_c}{dt}\text{ (6)}$

This states that in any element of the porous media, the increment of fluid volume comprises three parts, the net increment that flow in minus flow out, source generation and drainage resulting from external stress increment, which are corresponding to the three terms in the right hand of this equation.

Stress rebalance

The pore fluid pressure will enlarge the longitudinal strains, so the total strain is the superposition of confining stress induced strain and pore pressure induced longitudinal strains, as that

$\epsilon \text{ = }{\mathrm{D}}^{-1}\sigma +\beta \Delta P\text{ (7)}$

In which $\epsilon$ is the total strain vector, $\sigma$ is the (confining) stress vector, $\Delta P=P-P_0$ is an increment of pore fluid pressure, $P_0$ is the reference pore fluid pressure and $\beta$ is the linear expansion coefficient vector resulting from internal forces of fluid pressure increment. In the initial state $\beta$ is isotropic, $\beta ={\beta }_x={\beta }_y={\beta }_z$, which is defined as

$$\beta =\frac{1}{3}{C}_{bp}\text{ (8)}$$

While $C_{bp}$ is an bulk expansion coefficient, defined as [18,19]

$C_{bp}=\frac{1}{V_b}{\left(\frac{\partial V_b}{\partial P}\right)}_{P_c}\text{ (9)}$

By employing the concept of thermal elasticity, $\beta$ can also be calculated as

$\beta =\frac{\alpha }{3K_d}\text{ (10)}$

Where $K_d$ is the bulk modulus of the solid skeleton and $\alpha$ is the Biot coefficient, thus

$K_d=\frac{E}{3(1-2v)}=\lambda +\frac{2}{3}G,\alpha =1-K_d/K_s\text{ (11)}$

In which $K_s$ is the bulk modulus of matrix.

The stresses corresponding to the total strains in Eq. 7 are the effective stress ${\sigma }^{\prime }$, and can be written as

$${\sigma }^{\prime }\text{ = }\sigma +\text{α}\Delta P\text{ (12)}$$

Since the pore fluid pressure is not uniformly distributed, so the effective stress will cause stress redistribution. The tensor form of the stress differential equation can be written as

${\sigma }_{i,j}+\alpha P_{,j}+f_i=0\text{ (13)}$

Where $f_i$ is the body force per unit volume.

Fracture Opening and Anisotropy

As is mentioned previously, the key point using FEM to represent hydraulic fracture propagation is to establish the fracture opening in a continuum element, and based on this, to establish a series of second order Cartesian tensors of poromechanical properties such as damage, poro-elastic properties and permeability.

Fracture opening

Progressive fracturing in geo-materials transits a variety of scales and can be divided into three main stages: Evolution of distributed micro-damage, localization and subsequent macro-crack nucleation and macro-crack propagation [18,19,22]. For element of FEM, the entire fracture process can be represented by a damage variable and defined as

$D=\left\{\begin{array}{c}\begin{array}{cc}0,& {\epsilon }_I<{\epsilon }_{t0}\end{array}\\ \begin{array}{cc}1-\kappa \left(\frac{{\epsilon }_{tu}}{{\epsilon }_I}-1\right),& {\epsilon }_{t0}\le {\epsilon }_I<{\epsilon }_{tu}\end{array}\\ \begin{array}{cc}1,& {\epsilon }_{tu}\le {\epsilon }_I\end{array}\end{array}\right.\text{ (14)}$

In this, ${\epsilon }_{t0}$ is the threshold strain, representing the initiation of crack nucleation, ${\epsilon }_{tu}$ is the final strain when the fracture has transected the porous element, and $\kappa$ is a combined parameter, calculated as $\kappa ={\epsilon }_{t0}/({\epsilon }_{tu}-{\epsilon }_{t0})$, and ${\epsilon }_{{\rm I}}$ is the tensile strain controlling fracture opening. This tensile strain is in the same direction as the first effective principal stress ${{\sigma }^{\prime }}_I$ (Figure 1), while the effective stress following the cohesive law (Figure 2).

The magnitude of the fracture opening can be determined from strain and damage, by employing the thickness of the damage localization band, which is noted as $\delta$, typically 1-2 times the size of the cleat spacing. The element size is noted as L, thus it can be divided into two zones in the tensile direction: A concentrated damage zone of dimension $\delta$ and an undamaged zone of dimension $L-\delta$. Both zones are subject to the same stress, ${\sigma }^{\prime }$, thus

$\frac{w}{\delta }(1-D)E_0=\frac{w^t-w}{L-\delta }E={\sigma }^{\prime }\text{ (15)}$

Where w represents fracture opening, $w^t$ represents total elongation of the element and $E_0$ represents initial elastic modulus. From this equation, we have the relation

$w^t=\frac{(1-D)(L-\delta )+\delta }{\delta }w=\frac{(1-D)L+D\delta }{\delta }w\text{ (16)}$

Since the total strain can be expressed as

${\epsilon }_I=\frac{w^t}{L}\text{ (17)}$

Substituting Eq. 16 with 17 attains the hydraulic fracture opening, as

$$w={\epsilon }_I\frac{\delta }{D\frac{\delta }{L}+(1-D)}\text{ (18)}$$

According to Eq. 7, the relation between strain and stress can be written as

$${\epsilon }_I=\frac{{\sigma }_1}{E_0}+\beta \Delta P\text{ (19)}$$

Substituting Eq. 18 with Eq. 19, the constitutive relationship between fracture opening and fluid pressure is attained, thus

$w=\frac{\delta \left[\frac{{\sigma }_1}{E_0}+\beta \Delta P\right]}{D\frac{\delta }{L}+(1-D)}\text{ (20)}$

Where ${\sigma }_1$ is the minimum confining stress normal to fracture surfaces, is taken as negative.

Fracture induced anisotropy

The presence of a hydraulic fracture will bring about anisotropy to the hydro-mechanical properties such as damage, coefficients of poroelasticity, and permeability. These may be represented using a series of second order Cartesian tensors, which are defined in the global coordinates.

Anisotropy of damage: The damage tensor $D_{ij}$ in global coordinates can be converted from $D_{ij}^{\sigma }$, which is defined in the coordinates of the principal stresses $({\sigma }_1,{\sigma }_2,{\sigma }_3)$, as

$D_{pq}=M_{pi}{M}_{qj}{D}_{ij}^{\sigma }\text{ (21)}$

Where $M_{ij}$ is the conversion matrix, defined as

$\mathrm{M}\text{ = [}{M}_{ij}\text{] = [}{e^{\prime }}_i\cdot {\mathrm{e}}_j^{}\text{] (22)}$

and

$D_{11}^{\sigma }=D_n$, $D_{ij}^{\sigma }=0,(i,j\ne 1)\text{ (23)}$

In this $D_n$ represents the tensile damage in the direction of ${\sigma }_1$, which is derived from Eq. 14. Hence, the damage tensor in global coordinates are written as

$D_{pq}=M_{p1}{M}_{q1}{D}_n\text{ (24)}$

Anisotropy of poroelasticity coefficients: A hydraulic fracture also brings about anisotropy to the coefficients, $C_{pp}$ and $C_{bp}$. For these a set of damage related coefficients are introduced to achieve the anisotropy, as that

$C_{pp}=C_{pp}\psi$ and $C_{bp}=C_{bp}\sum \text{ (25)}$

In which $\psi$ and $\Sigma$ are the orthotropic coefficient vectors, with the components take the forms that

${\Psi }_i=\frac{1}{1-D_i}$ and ${\Sigma }_i=\frac{1}{1-D_i}$, $(i=1,2,3)$     (26)

Here $D_i$ represents damage in the global directions of x, y, z.

Anisotropy of permeability: The permeability tensor resulting from a series of joints sets or fractures is well known [23,24]. However, all those only consider fluid flow along the fracture, neglecting the normal flow, which represents the leak-off. We established a permeability tensor considering the both flows, as that

$$\mathrm{K}=M^T{K}^{\sigma }M\text{ (27)}$$

In which K and ${\mathrm{K}}^{\sigma }$ are the matrix of permeability coefficients in global coordinates and principal stress coordinates, respectively, and ${\mathrm{K}}^{\sigma }$ is written as

${\mathrm{K}}^{\sigma }=\frac{1}{\mu }\left[\begin{array}{ccc}(k_n-k_l)N_1{N}_1+k_l& (k_n-k_l)N_1{N}_2& (k_n-k_l)N_1{N}_3\\ (k_n-k_l)N_2{N}_1& (k_n-k_l)N_2{N}_2+k_l& (k_n-k_l)N_2{N}_3\\ (k_n-k_l)N_3{N}_1& (k_n-k_l)N_3{N}_2& (k_n-k_l)N_3{N}_3+k_l\end{array}\right]\text{ (28)}$

In which $N_1$, $N_2$ and $N_3$ are the three projections of the normal vector N of fracture surfaces on the three principal stress vectors; k_n and $k_l$ are the normal and tangent permeability, respectively. The normal permeability k_n can be recovered by modifying the in-situ permeability k with a modification $\xi$, as

$k_{\nu }=k\xi \text{ (29)}$

In which k is the in-situ permeability decreasing exponentially as the effective stress increases [25], calculated as

$k=k_0{e}^{-(\zeta {{\sigma }^{\prime }}_m/{\sigma }_0)}\text{ (30)}$

While $k_0$ is the intrinsic permeability tested in lab; ${{\sigma }^{\prime }}_m$ is the effective average stress, noted as compression positive; ${\sigma }_0$ is the reference stress, evaluated as the mean, maximum or minimum of magnitude vector of ${{\sigma }^{\prime }}_m$; $\zeta$ is a constant valued between 1.0 and 1.6.

The permeability along the fracture $k_l$ follows the cubic law [26,27] as

$k_l=\frac{w^2}{12}\text{ (31)}$

Here to all the components necessary for simulation of a hydraulic fracture propagation are established. In the next section, they are assembled by the weakly coupled FEM equations and run to work by the coupling analysis scheme.

Coupling Analysis

FEM format of coupling equations

The differential Eq. 6 and 13 can be transformed into FEM formats by utilizing Galerkin variational principle, so that in solid solution domain, the FEM format of Eq. 13 is written as

${\mathrm{K}}_{\text{s}}a=P_{\text{s}}+P_{\text{f}}\text{ (32)}$

Where ${\mathrm{K}}_{\text{s}}$ is the global stiffness matrix, a is the column matrix of the unknown nodal displacements, ${\mathrm{P}}_s$is the column matrix of solid load, ${\mathrm{P}}_f$ is column matrix of fluid pressure, calculated as

${\mathrm{P}}_f={\displaystyle \sum _e{\displaystyle {\int }_{{\Omega }_e}{\mathrm{B}}^T\mathrm{D}{\epsilon }_{\text{f}}}}d\Omega ={\displaystyle \sum _e{\displaystyle {\int }_{{\Omega }_e}{\mathrm{B}}^T\mathrm{D}\beta \Delta P}}d\Omega \text{ (33)}$

Where D is the elasticity stiffness matrix of the porous rock, B is element strain matrix, ${\epsilon }_f$ is column matrix of incremental strain, generated by the incremental pore pressure $\Delta P$.

In fluid solution domain, the FEM format of Eq. 6 can be written as

$\mathrm{C}\stackrel{•}{P}+{\mathrm{K}}_{\text{f}}\mathrm{P}=\mathrm{Q}\text{ (34)}$

Where P is the column matrix of the unknown pore pressure, $\stackrel{•}{P}\text{ = }d\mathrm{P}/dt$, ${\mathrm{K}}_{\text{f}}$ is the permeability matrix; C is the matrix of storage coefficient, Q is the column matrix of flow rate, which are respectively assembled as

$\mathrm{C}\text{ = }\varphi ({\mathrm{C}}_{pp}+{\mathrm{C}}_f)\text{ (35)}$

And

$\mathrm{Q}={\mathrm{Q}}_{\text{q}}^{}+{\mathrm{Q}}_{\text{g}}^{}+{\mathrm{Q}}_{{\text{p}}_{\text{c}}}^{}\text{ (36)}$

Where ${\mathrm{C}}_{pp}$ and ${\mathrm{C}}_f$ are matrixes of poroelasticity assembled using the coefficients defined in Eq.5; ${\mathrm{Q}}_q$, ${\mathrm{Q}}_g$, ${\mathrm{Q}}_{pc}$ are the column matrixes of flow rate, generation rate and confining stress change induced fluid content change, respectively. All the global matrixes are assembled as

${\mathrm{K}}_{\text{f}}\text{ = }{\displaystyle \sum _e{\mathrm{K}}^e}$, $C={\displaystyle \sum _e{C}^e}$, $\mathrm{Q}\text{ = }{\displaystyle \sum _e{\mathrm{Q}}^e}\text{ (37)}$

While

${\mathrm{K}}^e\text{ = }{\displaystyle {\int }_{{\Omega }^e}{\mathrm{B}}^T\mathrm{KB}}d\Omega$, ${\mathrm{C}}^e\text{ = }\varphi (C_{pp}+C_f){\displaystyle {\int }_{{\Omega }^e}{\mathrm{NN}}^{\text{T}}}d\Omega \text{ (38)}$

${\mathrm{Q}}_{\text{q}}^{\text{e}}\text{ = }{\displaystyle {\int }_{S_q}\mathrm{Nq}}dS$, $Q_g^e=\stackrel{⃛}{q}{\displaystyle {\int }_{{\Omega }^e}N}d\Omega$, $Q_{p_c}^e=\varphi C_{pc}\frac{d{P}_c}{dt}{\displaystyle {\int }_{{\Omega }^e}N}d\Omega \text{ (39)}$

In which N is shape function matrix. The weak coupling format of Eq. 32 and Eq. 34 is written as

$$\left[\begin{array}{cc}0& 0\\ 0& C\end{array}\right]\left\{\begin{array}{c}\stackrel{•}{a}\\ \stackrel{•}{\mathrm{P}}\end{array}\right\}+\left[\begin{array}{cc}{K}_s& 0\\ 0& C\end{array}\right]\left\{\begin{array}{c}a\\ P\end{array}\right\}\text{ = }\left\{\begin{array}{c}{P}_{\text{s}}\text{+}{P}_{\text{f}}\\ Q\end{array}\right\}\text{ (40)}$$

Analysis process

As is showed in the analysis flow chart in Figure 2, the coupling analysis can be divided into several calculation parts:

(1) Establish in-situ stress and permeability. This is very important, especially if no perforation made. As is known that the in-situ permeability depends on the initial in-situ stress distribution with injection hole bored, while the distribution of permeability around the injection hole determine the positions where the fracture propagation initiates.

(2) Porous flow analysis. For which, the transient analysis is carried out with fluid flux as surface load, permeability as input and pore pressure as the output.

(3) Stress adjustment analysis. For which, the pore fluid pressure input as the body force, the strains and effective stresses are attained by statics solution.

(4) Damage judgment. If no more new damaged elements generated, increase injection flow rate or injection time, and repeats the processes above; Else, calculate fracture opening, anisotropy of hydro-mechanical properties, renew the input parameters and repeat the processes above. Goes on these, the numeric hydraulic fracture propagates.

Fluid loading scheme

In order to identify fracture growth during continuous injection, a loading scheme of a stepwise increased solution duration is adopted. This requires the assumption that the fracture can completely be closed when hydraulic fluid is drained so that the parameters of hydro-mechanical properties can be elastically handled. As is showed in Figure 3a, $C_V$ represents the global elasticity rigidness of surrounding rock, it decreases with the fracture propagation going on, while the increasing fluid loading is carried out by the increasing injection volume V. Based on this, the loading scheme is that keeping the flow rate constant, the fluid loading duration is stepped into small time steps, and for each transient analysis of flow, the solution duration stepwise increases with the input parameters kept renewed. Further, in order to simulate a real injection process, a limit flow rate Q_m and a character time up to this pressure T₀ are set up, hence the flow rate load Q at each transient analysis is written as

$Q^{(n)}=\left\{\begin{array}{cc}\frac{Q_m}{T_0}{T}^{(n)},& T^{(n)}<T_0\\ Q_m,& T_0\le T^{(n)}\end{array}\right.\text{ (41)}$

Correspondingly, the flux loaded on the walls of the injection hole through coal bed, is written as

$q_0=Q_m/(\pi dh)\text{ (42)}$

where h represents coal bed thickness and d the hole diameter. In this paper $Q_m=10m^3/\min u.$

Example

FEM model

The example model (Figure 4) is taken from coal bed methane formation which is buried in 750 m depths, 200 m × 200 m in horizontal dimensions, 5 m thick for top and bottom beds, respectively, h = 10 m thick for coal bed, and d = 20 mm for injection hole diameter. The minor horizontal principal stress is ${\sigma }_h=7.8Mpa$, the major horizontal principal stress ${\sigma }_H=\lambda {\sigma }_h=15.6Mpa$, with the stress ratio is $\lambda =2.0$. The vertical stress is calculated as ${\sigma }_v=\rho gH$, in which $\rho$ is mass density of overlaying rocks, H is buried depth, g is gravitational acceleration. All boundary displacements are set to zero, since it is the requirement of implanting initial stresses, and on the other hand, it is more reasonable for real in-situ deformation state. Figure 4a shows the formation compositions and positions of cross-sections, Figure 4b and Figure 4c show the middle horizontal and vertical sections with injection hole crossed. The direction of horizontal major principal stress ${\sigma }_H$ is expressed using angle $\theta$, which rotates anti-clockwise from negative Z coordinate. Following the in-situ.

Property parameters

The formation property parameters include permeability, coefficients of poroelasticity and damage model parameters. These are usually assumed satisfying Weibull distribution with the probability density function as

$f(\Omega )=\frac{m}{\Omega }{\left(\frac{\Omega }{{\Omega }_0}\right)}^{m-1}\exp \left[-{\left(\frac{\Omega }{{\Omega }_0}\right)}^m\right]\text{ (43)}$

In which $\Omega$ represents element property parameters, ${\Omega }_0$ is the reference modulus, m is shape factor, representing the homogeneous degree of parameter distribution. The more m value is, the more uniform, in this paper, m = 10. The reference modulus ${\Omega }_0$ are valued as in Table 1, which probably amount to their average values. Hydraulic fluid property parameters are listed in Table 2.

Results and validations

Hydraulic fracture zones: The numeric results of hydraulic fracture zone are shown in Figure 5, Figure 6 and Figure 7, from which the basic features can be derived, that

(1) Hydraulic fracture always propagates in the direction of the major principal stress ${\sigma }_H$, with the fracture surfaces normal to the minimum principal stress ${\sigma }_h$ (Figure 5). This conforms to the routine knowledge and proves that the tensors established for the hydromechanical properties are right and can effectively control the direction of fracture propagation.

(2) The horizontal cross-section of fracture zone and fluid pressure contours are approximate to ellipse areas. This conforms to that computed by Liu [28] (Figure 8a) and those described in [1-13].

(3) The vertical section is approximate to an ellipse area, and slightly cut through the top and bottom beds in the vicinity of injection hole (Figure 6b). This conforms to that proposed by Peirce [29] for three-layered formation (Figure 8b).

Temporal variation of hydraulic fracture parameters: The parameters of hydraulic fracture include fracture length, opening and fluid pressure. The fracture lengths are measured from the resultant pictures in Figure 7. The fracture openings and fluid pressure at fracture mouth are measured by setting a monitor element at fracture mouth for every transient analysis. The comparisons of the temporal variation of these parameters via the analytical solutions are showed in Figure 9, Figure 10 and Figure 11. The analytical solutions proposed by Nordgren [2], Geertsman & Klerk [30] for KGD and PKN models in the case of high leak-off, are adopted, hence

$\left.\begin{array}{c}L(t)=\frac{q_0}{\pi h{C}_l}{t}^{1/2}\\ w(0,t)=4{\left[\frac{2(1-v)\mu q_0^2}{{\pi }^{3}Gh{C}_l}\right]}^{1/4}{t}^{1/8}\\ P(0,t)=1.135\eta {\left(\frac{G{q}_0\mu }{{(1-v)}^{3}L{(t)}^2}\right)}^{1/4}+S\end{array}\right\}\text{ (44)}$

Where G is coal bed shear constant, S is the principal stress normal to fracture surfaces, $S={\sigma }_h=7.8Mpa$, $\eta$ represents the effects of leak-off, $\eta =0.32\times {10}^{-15}$, $C_l$ is leak-off coefficient involving with permeability, porosity, in-situ stresses, fracture toughness and fluid viscosity, etc, is valued as $C_l=4.98\times {10}^{-4}$ [31], and $q_0$ is injection flux.

The comparisons show that numeric solutions are well conformed to the analytical solutions, while some slight differences indicate that the numeric solutions can exhibit more plentiful information about hydraulic fracture propagation. Base on these differences the process of hydraulic fracturing can be divided into four stages, which are conceptualized in Figure 12, that

Stage-I, fracture nucleation, during which, a macro embryo fracture takes its form in the close vicinity of injection hole. Since it is aggregated from distributive cracks, the gaps and bridging of them constitute fracture cohesive zones. In this stage the peak fluid pressure is used both in overcoming the traction of cohesive zone and support the confining normal stress;

Stage-II, kinetic propagation, during which, the sudden break of cohesive traction causes that the fluid pressure at fracture mouth has a big drop, the fracture opening goes fast up to a peak value and fracture length increases fast. These indicate a kinetic propagation.

Stage-III, steady propagation, during which, the fluid pressure at fracture mouth keeps constant while fracture length goes up fast and fracture opening goes to a constant slowly.

Stage-IV, propagation termination, during which, as fracture length increases, the injection flow rate cannot afford to going up of leak-off, so that gives rise to slow drop of fluid pressure, decrease of fracture opening and propagation termination in the end.

It is obvious that the analytical solutions ignore stage-I and stage-IV.

Spatial variation of parameters of hydraulic fracture: The analytical solutions of fracture opening and fluid pressure along fracture length are the so-called SCR asymptote by Adachi & Detournay, et al. [32-34], as that

$\left.\begin{array}{c}\Omega (\xi ,\tau )=A(\tau ){(1-\xi )}^{\alpha },0<\alpha <1\\ \prod (\xi ,\tau )=\frac{1}{4\gamma }A(\tau )\alpha \cot \pi \alpha {(1-\xi )}^{\alpha -1}\\ \alpha =2/3,A(\tau )=2^{1/3}{3}^{5/6}{\left(g_m{\gamma }^2\dot{\gamma }\right)}^{1/3}\end{array}\right\}\text{ (45)}$

In which $\Omega ,\prod ,\xi ,\gamma ,\tau ,g_m$ are the dimensionless forms of fracture opening, fluid pressure, position coordinate, fracture length, injection time and viscosity scaling, respectively. The comparisons between numerical solutions and the analytical solutions for fluid pressure and half width along the fracture length are showed in Figure 13 and Figure 14. The contrast between fluid pressure and half width along the fracture length is showed in Figure 15.

The comparison in Figure 13 shows that fluid pressure goes head of fracture tip and is almost evenly distributed along fracture hull. This is somewhat different from the analytical solutions. However, if we realize the premise of analytical solution, we would like to accept the numeric solution. Since mathematical difficulties, in analytical solutions the coupling between fracture tip and fluid front is generally assumed to be progressive, so that the fluid pressure distribution is always lagged behind fracture tip and that a pressure void always exists ahead of fluid front. So we can regard the analytical solution as the case of incomplete coupling between fracture tip and fluid front, which always occurs in the situations of large toughness, high viscosity and no leak-off; while in most cases, fluid pressure distribution goes ahead of fracture tip, complete coupling occurs, this is the case that numeric solution represents. The comparisons in Figure 14 and Figure 15 further indicate that there is cohesive zone ahead of the fracture tip.

Based on these, we can give out two types of coupling modes in fracture tip, the incomplete coupling model (Figure 16a) and complete coupling model (Figure 16b), which are corresponding to the analytical models and numeric solution in present paper, respectively. In the complete coupling model, D = 1 represents completely fractured zone, D = 0 represents elastically expansion zone of pore pressure, in between them, 1 > D > 0 represents cohesive fracture zone, the fracture opening is more like an aequilate 'bag', indicating that fluid pressure energy are mainly used in overcoming in-situ stress clamping and viscosity dissipation from fluid front invasion and leak-off.

Summary

In this paper we systematically clarified poromechanical modeling of hydraulic fracture propagation, which is proved to be correct and effectively. Especially, it is of the advantage that can be used in simulation of hydraulic fracturing in complex stress and inclined formation, which is the main aim we work to. Also some limitations should be noted here, that

(1) The fracture criterion and opening calculation method cannot be simply generalized to incomplete coupling model, in which case the toughness, viscosity and low leak-off are large,

(2) The thickness of localized fracture band is a material parameter, should be a specifically researched, and

(3) The patterned fracture zone shows that hydraulic fracture is not mere a plane but a band zone. This is verified by micro seismics monitored on the spot. Of course, the precision of fracture zone is element size related, it should not be larger than the size of RVE.

Acknowledgements

Many thanks to Special project for the 13th five-year plan, China, the project number is 2016ZX05067006-002.

---