Introduction

Shale gas plays a significant role in unconventional resources, while shale gas reservoirs have an extremely low permeability and porosity. Hydraulic fracturing is an effective method to improve unconventional resource production (Ni et al. 2019; Xie et al. 2019). Previous studies have demonstrated that shale reservoirs are apt to form complex fracture networks, because these hydraulic fractures interact with pre-existing natural fractures. The hydraulic fracture path, which occurs as mixed-mode extension, is not straight due to the collective effect of shear and normal stresses on fracture propagation. Several factors highly influence hydraulic fracture initiation and propagation, such as Young’s modulus, Poisson’s ratio, in situ stresses and fracturing fluid rate. It is important to understand the fracture mechanism to optimize its design and improve production.

In practical hydraulic fracturing, the rock is not only subjected to tensile stress but also to shear stress, so mixed-mode loading is common in engineering applications. Mixed-mode failure always occurs in brittle or quasi-brittle materials. Therefore, it is important to investigate the structural integrity of fracture components in rock samples under mixed-mode loading. Many researchers have examined mixed fractures in different materials and methods. Four methods are recommended to obtain the mode I fracture toughness by the International Society for Rock Mechanics (ISRM), which include chevron-bend (CB), short-rod (SR), cracked chevron-notched Brazilian disc (CCNBD) and semi-circular bend (SCB) specimens. However, these methods are only widely used to determine the pure mode I fracture toughness. The semi-circular bend (SCB) test has attracted much attention for testing different materials, such as rocks (Zuo et al. 2019), asphalt mixtures (Ameri et al. 2012) and polymethyl methacrylate (PMMA) (Ayatollahi et al. 2006). To propagate mixed-mode fractures, the method of changing the initial fracture direction has been applied in semi-circular bend (SCB) tests. The SCB test to analyse mixed-mode fractures has the advantage of being easily applied to shale samples, and the results are reliable.

In numerous experimental and theoretical analysis studies, the influence of natural fractures on hydraulic fracture initiation, propagation and interaction has been investigated (Warpinski and Teufel 1987; Renshaw and Pollard 1995; Zhou et al. 2008; Olson et al. 2012). Blanton (1982) used Devonian shale and hydrostone to demonstrate how hydraulic fractures intersect with natural fractures under a high stress and large incidence angles. Renshaw and Pollard (1995) proposed a direct standard method to estimate whether hydraulic fractures would pass through nearby orthogonal natural fractures. Gu et al. (2011) extended the criteria of Renshaw and Pollard to include hydraulic fractures crossing frictional interfaces at nonorthogonal angles, as well as used laboratory experiments to validate the criteria. The effect of natural fractures on hydraulic fractures was evaluated through energy-based analysis (Yao et al. 2018). Many numerical methods, such as the displacement discontinuity method (DDM) (Zhang and Jeffrey 2007), the finite-element method (FEM) (Guo et al. 2015, 2017), the extended finite-element method (XFEM) (Mohammadnejad and Khoei 2013; Shi et al. 2016, 2017; Suo et al. 2020), and the discrete-element method (DEM) (Wasantha and Konietzky 2016), have been widely used to investigate hydraulic fracture interactions with natural fractures. The extended finite-element method (XFEM) not only inherits the advantages of the finite-element method but also overcomes its disadvantages. XFEM does not need to re-mesh fractures, and a refined fracture tip node is used to solve propagation and interaction problems. Thus, the extended finite-element method (XFEM) effectively addresses hydraulic fracture initiation, propagation and interaction.

In this study, a two-dimensional (2D) hydraulic fracturing mixed-mode model is developed  based on the extended finite-element method (XFEM) and used to simulate hydraulic fracture initiation and propagation and damage evolution. Fluid leak-off and diffusion into shale samples are considered. Mixed-mode hydraulic fracture initiation, propagation and interaction with different types of natural fractures are modelled. In addition to comparing different factors influencing hydraulic fracturing propagation, the optimal values of these parameters in hydraulic fracturing are also determined. The interaction between two different natural fractures and hydraulic fractures is examined.

Mixed-mode fracture initiation and propagation

XFEM

The extended finite-element method (XFEM) is a numerical method that has been developed in recent years to solve the problems of discontinuous mechanics. XFEM was proposed by Ted Belytschko and Mullen (1978) and was derived from the traditional finite-element method. The extended finite-element method (XFEM) can overcome the shortcomings of the finite-element method. First, there is no need to re-mesh fractures during the fracture expansion process. Second, the crack tip does not need to be refined. Third, there is no need to transfer data, because the same grid is used (Yan et al. 2019). The basic equation is as follows:

$$u = \mathop \sum \limits_{I = 1}^{N} N_{I} \left( x \right)\left[ {u_{I} + H\left( x \right)a_{I} + \mathop \sum \limits_{\alpha = 1}^{4} F_{\alpha } \left( x \right)b_{I}^{\alpha } } \right],$$
(1)

where NI(x) is the conventional shape function, uI is the displacement vector of the continuous part, H(x) is the jump function, bIα is the crack tip node expansion degree of freedom, Fα(x) is the crack tip progressive displacement function, and I refers to the node set for all grid nodes. The conventional shape function is available for all nodes, such as traditional finite-element nodes. The jump function is used for those nodes whose shape function support is cut by the crack interior. The crack tip shape function is only applicable for the crack tip.

Fluid–solid coupling

At the beginning of hydraulic fracturing, there is no hydraulic fracture. With increasing injected fracturing fluid volume, crack propagation occurs along the maximum horizontal stress. It is important to understand the theory of the hydraulic fracturing process. Using the finite-element method, the pressure field p is as follows:

$${\text{p}}\left( {\text{s}} \right) = \mathop \sum \limits_{{I \in s_{n} }} N_{I}^{p} \left( s \right)p_{I} = N^{p} \left( s \right)P,$$
(2)

where Sn is the fluid node set, \(N_{I}^{p} \left( s \right)\) is fluid node I corresponding to pressure shape function p\(,{ }N^{p} \left( s \right)\) is the shape function matrix, and P is the node pressure degree of freedom vector. Additionally, the crack width w has the following discrete form:

$${\text{w}} = \mathop \sum \limits_{{I \in S_{w} }} N_{I}^{w} u_{I} = N^{w} U,$$
(3)

where Sw is the collection of solid elements containing fluid junctions, \(N_{I}^{w}\) indicates the transform of the node degree of freedom vector into a shape function of the crack width, and \(N^{w}\) refers to the shape function matrix. Based on the weak forms of the stress equilibrium and lubrication equations, the corresponding discretized forms can be written as

$${\text{KU}} - {\text{QP}} - F^{xt} = 0$$
(4)
$$Q^{T} {\text{U}} + {\text{HP}} + {\text{S}} = 0,$$
(5)

where K is the overall stiffness matrix, Q is the coupling matrix to efficiently convert the fluid node pressure vector P into a solid node, F refers to the external load vector, H is the fluid flow matrix, and S is the source term. The fluid–solid coupling equations of hydraulic fracturing can be solved with the Newton–Raphson method.

Crack initiation and propagation criteria

Several theories have been proposed for mix-mode fractures. Three criteria, namely, maximum energy release (Gmax), minimum strain energy density (Smin) and maximum tangential stress (MTS), are the most commonly used criteria for mixed-mode fractures. The maximum circumferential stress (MTS) criterion is applied to the propagation of hydraulic fractures. When the stress intensity factor value at the crack tip is greater than the toughness of the rock, the fracture will expand and extend along the direction of the maximum circumferential stress. The expression is as follows:

$$K_{I} \ge K_{IC} ,$$
(6)

where KI is the crack tip fracture toughness, and KIC is the mixed-mode rock toughness. All the criteria mentioned above depend on the stress field that exists prior to initiation of fracture propagation. The stress field near the crack tip in a homogeneous and isotropic linear elastic material in polar coordinates (Fig. 1) is defined as

$$ \begin{aligned}\sigma_{rr} &= \frac{1}{{\sqrt {2\pi r} }}cos\frac{\theta }{2}\left[ {K_{I} \left( {1 + sin^{2} \frac{\theta }{2}} \right) + \frac{3}{2}K_{II} \left( {sin\theta - 2tan\frac{\theta }{2}} \right)} \right] \\ &\quad+ Tcos^{2} \theta + O\left( {r^{\frac{1}{2}} } \right)\end{aligned} $$
(7)
$$\sigma_{\theta \theta } = \frac{1}{{\sqrt {2\pi r} }}cos\frac{\theta }{2}\left( {K_{I} cos^{2} \frac{\theta }{2} - \frac{3}{2}K_{II} sin\theta } \right) + Tcos^{2} \theta + O\left( {r^{\frac{1}{2}} } \right)$$
(8)
$$\sigma_{r\theta } = \frac{1}{{\sqrt {2\pi r} }}cos\frac{\theta }{2}\left[ {K_{I} sin\theta + K_{II} \left( {3cos\theta - 1} \right)} \right] - Tcos^{2} \theta + O\left( {r^{\frac{1}{2}} } \right),$$
(9)

where \(\sigma_{rr,}\)\(\sigma_{\theta \theta }\) and \(\sigma_{r\theta }\) are the polar stresses. The higher-order terms \(O\left( {r^{\frac{1}{2}} } \right)\) can be considered negligible near the crack tip. KI, KII and T play a part in the mode I/II stress intensity factor and T-stress. Erdogan and Sih (1974) proposed a brittle material standard for the traditional maximum tangential stress (MTS): (1) crack propagation occurs radially along the direction perpendicular to the maximum tension (θ0) within the plane crack tip initiation (for example, the tangential stress (\(\sigma_{\theta \theta }\)) reaches a maximum) (2) when the maximum tangential stress (\(\sigma_{\theta \theta max}\)) is considered to be the critical material constant value \({ }\sigma_{\theta \theta crit}\). Based on this standard, with \(\sigma_{r\theta }\) = 0, the higher-order terms in r are ignored and only the singular stress terms are applied. The non-dimensional form of the T-stress (T*) for the specimen is as follows:

$$T^{*} \left( {\frac{a}{R},\frac{S}{R},{\upbeta }} \right) = {\text{T}}\frac{2RT}{P},$$
(10)

where P is the loading pressure. T* increases with increasing β shift from mode I fractures to mode II fractures. Especially for larger fracture angles or higher load conditions of mode II fractures, the magnitude of T* is more significant (Ayatollahi 2006). The critical distance rc can be regarded as a parameter to describe the extent of the damage area formed around the crack tip. According to the maximum tangential stress (MTS) criteria and the study of Ayatollahi, the KI/KIC and KII/KIC ratios during mix-mode loading are as follows:

$$\frac{{K_{I} }}{{K_{IC} }} = [{\cos}\frac{{\theta_{0} }}{2}\left( {cos^{2} \frac{{\theta_{0} }}{2} - \frac{3}{2}\frac{{Y_{II} }}{{Y_{I} }}{\sin}\theta_{0} } \right) + \sqrt {\frac{{2r_{c} }}{a}} \frac{{T^{*} }}{{Y_{I} }}sin^{2} \theta_{0} ]^{ - 1}$$
(11)
$$\frac{{K_{II} }}{{K_{IC} }} = [{\cos}\frac{{\theta_{0} }}{2}\left( {\frac{{Y_{I} }}{{Y_{II} }}cos^{2} \frac{{\theta_{0} }}{2} - \frac{3}{2}{\sin}\theta_{0} } \right) + \sqrt {\frac{{2r_{c} }}{a}} \frac{{T^{*} }}{{Y_{II} }}sin^{2} \theta_{0} ]^{ - 1} .$$
(12)
Fig. 1
figure 1

Stress field in polar coordinate system

Full size image

Damage evolution under mixed-mode failure is defined based on the Benzeggagh–Kenane fracture criterion (Benzeggagh and Kenane 1996), and when the critical fracture energies during deformation along the primary and secondary shear directions are similar, i.e., Gsc = Gtc, the criteria are given by

$$G_{n}^{C} + \left( {G_{s}^{C} - G_{s}^{C} } \right)\left( {\frac{{G_{S} }}{{G_{T} }}} \right)^{\eta } = G^{C} ,$$
(13)

where the mixed-mode fracture energy is defined as follows:

$$G^{C} = G_{n} + G_{s} + G_{t} ,$$
(14)

where Gn, Gs, and Gt refer to the work performed by traction and its conjugate relative displacement in the normal, primary, and secondary shear directions, respectively. Gnc, Gsc and Gtc refer to the critical fracture energies required to cause failure in the normal, primary, and secondary shear directions, respectively.

Mixed-mode cohesive zone model

The Fig. 1 shows stress components around the crack tip in the stress plane. Figure 2 shows a graphic of the mixed mode in the cohesive zone model (Zhang et al. 2015). The two unshaded  triangles represent pure normal and pure shear deformation responses. The other shaded triangle between the two pure modes is the mixed-mode response. Point A is the fracture initiation point, and the fracture gradually develops down to point B, which represents the complete failure for the mixed-mode response. All intermediate vertical planes (including vertical axes) represent the damage response under mixed-mode conditions with different modes of mixing. The position of the mixed-mode response between the pure normal and shear responses is controlled by the mode ratio as follows:

$$m_{1} = \frac{{G_{n} }}{{G_{T} }}$$
(15)
$$m_{2} = \frac{{G_{s} }}{{G_{T} }}$$
(16)
$$m_{3} = \frac{{G_{t} }}{{G_{T} }}.$$
(17)
Fig. 2
figure 2

Schematic of the bilinear, mixed-mode cohesive zone model

Full size image

Numerical model setup

Mixed-mode verification using the semi-circular bend (SCB) test

There exist several experimental techniques to calculate the fracture toughness under mixed-mode propagation. The semi-circular bend (SCB) test has the advantage of easy specimen preparation and simple loading condition. The semi-circular bend (SCB) test has an initial crack (a) of different lengths and approach angles. Figure 3a shows the different initial crack angles to investigate the fracture toughness. Vertical load P is exerted on the specimen, and two rollers support the sample with a distance of 2S. When we vary the initial crack angle (α), different values of the fracture toughness of modes I and II are obtained. A pure mode I fracture occurs when the crack inclination angle is α = 0°. With increasing initial crack angle α, mode I decreases, while mode II increases. When the crack inclination angle α reaches a specific value, the fracture becomes a pure mode II fracture. The crack parameters KI and KII are functions of a/R and α, respectively, in the semi-circular bend test. The modes I and II stress intensity factors of the SCB specimen can be described as follows:

$$K_{If} = \frac{P}{2RB}\sqrt {\pi \alpha } Y_{I}$$
(18)
$$K_{IIf} = \frac{P}{2RB}\sqrt {\pi \alpha } Y_{II} ,$$
(19)
Fig. 3
figure 3

a Semi-circular bend (SCB) specimen under mixed-mode loading. b Finite-element mesh for mixed-mode fracture propagation

Full size image

where P is the loading pressure of the sample; B is the sample thickness; and KIf and KIIf are the stress intensity factors of modes I and II, respectively. Normally, YI and YII are functions of a/R, S/R and α. For the mixed-mode numerical simulations, a/R is 0.3, S/R = 0.72 and α = 30°.

The typical mesh pattern to simulate mixed-mode fracture propagation is shown in Fig. 3b, whereby α is 30°, and the force is applied at the top of the sample. The initial fracture length divided by the sample radius R equals 0.3. The different areas have mesh grids with different densities to promote calculation convergence. The meshing of an incrementally extending fracture was performed simply in ABAQUS by defining the location of the new crack tip. The basic numerical simulation parameters are listed in Table 1.

Table 1 Parameters for the numerical simulation
Full size table

Figure 4a shows (Zuo et al. 2019) the three-point bending test result on semi-circular bend sandstone specimens, which indicated that the concentrated force has a guiding effect on crack propagation, and we use the same parameters to perform the numerical simulations. Figure 4b shows the mixed-mode fracture trajectory with α = 30°. There are three stages in the fracture process: crack initiation, crack propagation and complete damage. The crack propagation path is similar to the experimental results of Zuo et al. (2019). The parameters KI and KII are applied in the mixed-mode hydraulic fracturing model.

Fig. 4
figure 4

Experimental and simulation results of the crack path in an SCB specimen with α = 30°

Full size image

Numerical model establishment

There are two types of wells in hydraulic fracturing: vertical and horizontal wells. Horizontal wells are widely used in shale gas reservoirs, because they increase the reservoir exposure areas and improve well production. Figure 5a shows a typical horizontal well geometry pattern. The hydraulic fracturing fluid from the wellbore is injected into the shale reservoir to generate and propagate fractures. Figure 5b shows the 2D model mesh and geometry pattern of a single hydraulic fracture. A 2D model is built to simulate hydraulic fracture propagation using the extended finite-element method (XFEM) based on the cohesive zone model in ABAQUS. The hydraulic fracture propagation region is divided into two parts, and the mesh sizes differ. The mesh size of the hydraulic fracture propagation part is more refined than that of the other parts, which promotes convergence in the software.

Fig. 5
figure 5

Two-dimensional mechanical model of horizontal shale gas well pre-crack

Full size image

Table 2 summarizes the input simulation parameters. In our model, the shale gas reservoir length is 50 m, and the thickness is 50 m. The initial hydraulic fractures are emitted from the centre of the left edge (the wellbore). The borehole wall is a fixed boundary, while the other boundaries are loading boundaries. The entire numerical simulation process lasts 50 s.

Table 2 Simulation parameters
Full size table

Numerical model verification

The hydraulic fracture length is considered to increase in one direction in the KDG model (Geertsma 1969). The KDG model is a simplified hydraulic fracturing model, and it does not contain mesh along cracks and is less flexible; however, XFEM enables calculation of each node distance to obtain the whole crack displacement. This example is used to validate the efficiency of the proposed method. The basic parameters are the same as in Table 2. The change in hydraulic fracture half-length with the fracturing time is shown in Fig. 6 between the KDG and XFEM models. The difference between the calculated result and the classical model is very small, which supports the correctness of the proposed model.

Fig. 6
figure 6

Comparison of the XFEM and KDG model results

Full size image

Results and discussion

In this section, the effect of different parameters on the fracture propagation are investigated. According to previous results of fracturing numerical modelling, the value range of the Young's modulus is generally 12–34GPa (Zhao et al. 2018) and the value range of Poisson's ratio is generally 0.15–0.25 (Dahi-Taleghani and Olson 2011; Li et al. 2020). For this reason, in this work the value of Young's modulus is taken in the range 20–30 GPa and the value of Poisson's ratio is taken in the range 0.16 –0.23. The in-situ stress differences are generally in the range 0–10 MPa and the fracturing fluid injection rate is usually defined in the range 0.8–6 m3/min (Guo et al. 2011; Wang et al. 2019; Li et al. 2019). Therefore, in this study the in-situ stress differences are evaluated in the range 0–10 MPa and the fracture fluid injection rate is evaluated in the range 0.9 –1.2 m3/min.

The influence of Young’s modulus and Poisson’s ratio on the fracture length

The basic rock mechanics parameters of Young’s modulus and Poisson’s ratio are significant hydraulic fracturing parameters. Young’s modulus of the rock is the ratio of stress and axial strain, and Poisson’s ratio is the radial expansion divided by the vertical contraction. The impact of Poisson’s ratio on the fracture half-length with time is shown in Fig. 7. In this model, the reservoir Poisson’s ratio changes from 0.16 to 0.23, which indicates that the fracture length is negatively correlated with Poisson's ratio. The simulation results based on different elastic modulus values are shown in Fig. 8. The reservoir Young's modulus is changed to 20 GPa, 25 GPa and 30 GPa. As Young’s modulus increases, the half-length of the hydraulic fracture gradually increases. In the initial stage, the fracture half-length did not change substantially, while it slowly increased in the later stages. The fracture half-length changes more when Young’s modulus is increased from 25 to 30 GPa compared to when it is increased from 20 to 25 GPa. Compared with Poisson's ratio, Young's modulus has a greater effect on the crack growth length, and Poisson’s ratio is negatively related to Young's modulus. In this study, Young’s modulus and Poisson’s ratio are key parameters controlling fracture geometry.

Fig. 7
figure 7

Influence of Poisson’s ratio on the fracture half-length with time

Full size image
Fig. 8
figure 8

Influence of Young’s modulus on the fracture half-length with time

Full size image

The influence of in situ stress differences on the fracture length

The impact of the in situ stress on hydraulic fracturing propagation is examined in this section. There are three principals in situ stresses: the minimum horizontal stress (ϭh), the maximum principal stress (ϭH) and the overburden stress (ϭv). In-situ stress differences lead to an uneven reservoir stress distribution and compaction of natural fractures, resulting in more discontinuities. For enhanced fracture growth, the minimum horizontal stress is set to σh = 20 MPa and the maximum horizontal stress is σH = 25 MPa (Suo et al. 2018a, b). The fracture half-length gradually increases with increasing in situ stress. Figure 9 shows the in situ stress changes over time. At the beginning of hydraulic fracturing, the different maximum principal stresses have little influence on the fracture length. Over time, the difference in length gradually increases. In addition, the figure illustrates that for different stress difference (the maximum principal stress minus the minimum principal stress) values, the horizontal stress differences have no significant impact on the fracture shape, and the fracture length remains uniform.

Fig. 9
figure 9

Influence of the horizontal stress difference on the fracture half-length with time

Full size image

The influence of the fracturing fluid rate on the fracture length

In engineering applications in the field, the proper fracturing fluid injection rate will yield more production. The numerical simulation results of the fracture length at different fracturing fluid rates are shown in Fig. 10. The injection flow rate ranges from 0.015 to 0.02 m3/s in the model. Over time, the pressure inside the fracture increases promptly, paving the way for the fracture extending along the length direction. The results show that the injection flow rate has little effect on the fracture length.

Fig. 10
figure 10

Influence of the injection fluid rate on the fracture half-length with time

Full size image

Hydraulic fractures intersecting with frictional natural fractures

There are many natural fractures in shale reservoirs. Activating pre-existing natural fractures during fracture treatment results in a complex fracture network that can improve production. Open natural fractures intersected by hydraulic fractures will acquire a new orientation, and thus, the numerous fracture interactions result in a complex network fracture. In this section, the impact of two types of natural fractures intersecting with hydraulic fractures at different angles is simulated using the extended finite-element method (XFEM). The geometry of the hydraulic fracturing model is shown in Fig. 11. The reservoir is 50 m×50 m, and the initial hydraulic fracture is set in the middle of the left edge, perpendicular to the edge. The initial length of the hydraulic fracture is 1 m, and the fracturing fluid Q is also injected at this point. The natural fracture length is 6 m. The input parameters are listed in Table 3.

Fig. 11
figure 11

Half of the symmetrical hydraulic fracturing model

Full size image
Table 3 Input XFEM coupled model parameters
Full size table

Three possibilities can be occurred when the hydraulic fracture interacts with natural fracture: arresting, crossing, or offsetting. Several criteria have been proposed by researchers, while in our study the extended Renshaw and Pollard criteria are used to reveal the interactions between hydraulic and natural fractures; when the hydraulic fracture pressure is lower than the rock tensile stress, fractures cannot propagate. When the maximum principle stress is equal to the rock tensile stress, a new fracture is initiated on the other side of the natural fracture (Wang et al. 2018). At the same time, no-slip conditions should be fulfilled along the frictional natural fracture fault surface. In spite of these conditions, the hydraulic fracture will cross or branch into the deviated frictional natural fracture. Equations (20) and (21) define the tensile and frictional fracture criteria respectively.

$$\sigma_{1} = T_{0}$$
(20)
$$\tau_{\beta } < S_{0} - \mu_{f} \sigma_{\beta y} ,$$
(21)

where β is the intersection angle between the hydraulic and natural fractures;

\(\tau_{\beta }\) is the combined shear stress along the natural fracture (NF) surface under the action of remote and local crack tip stresses;

\(\sigma_{\beta y}\) is the combined normal stress;

\(T_{0}\) is the rock tensile strength;

\(S_{0}\) is the cohesion force of the friction natural fracture; and.

\(\mu_{f}\) is the frictional coefficient of the natural fracture surface.

The von Mises stress distributions of a single hydraulic fracture intersecting at different angles with a frictional natural fracture are shown in Fig. 12. In Fig. 12a, there are two stress concentration regions at the tip of the natural fracture, and a stress balance is maintained in the whole model. As shown in Fig. 12b, c, the von Mises stress value of the upper part is higher than that of the lower part in the case where the fracture is easily diverted upward. The stress at an intersection angle of β = 60° is higher than that at β = 45°. With decreasing interaction angle, the pressure inside the crack is also reduced, which indicates that the fracture propagation resistance is lowers at a smaller intersection angle.

Fig. 12
figure 12

Von Mises stress distributions at different angles of hydraulic fracture intersect with frictional natural fracture

Full size image

Hydraulic fractures intersecting with cemented natural fractures

When a hydraulic fracture interacts with a cemented natural fracture, there are only two possible scenarios, namely, the hydraulic fracture crosses the natural fracture or the hydraulic fracture propagates along the cemented natural fracture. Dahi-Taleghani and Olson (2011) proposed that oblique crossover transfer can only happen along one direction. Because of the propagation and interaction of fractures, mixed patterns always occur. In this paper, a standard based on the energy release rate method considering KI and KII is applied to obtain the cross/prohibition behaviour, which has been verified by researchers (Dahi-Taleghani and Olson 2011). A fracture propagates when the maximum energy release rate is higher or equal to the critical Gc. The energy release rate in a particular direction θ can be written as

$$G_{0} = \frac{{K_{I\theta }^{2} + K_{II\theta }^{2} }}{{E^{*} }}$$
(22)
$$K_{I\theta } = \frac{1}{2}cos\left( {\frac{\theta }{2}} \right)\left[ {K_{I} \left( {1 + cos\theta } \right) - 3K_{II} sin\theta } \right]$$
(23)
$$K_{II\theta } = \frac{1}{2}cos\left( {\frac{\theta }{2}} \right)\left[ {K_{II} sin\theta + K_{II} \left( {3cos\theta - 1} \right)} \right]$$
(24)
$$E^{*} = \left\{ {\begin{array}{*{20}c} {E\quad \left( {\text{plane stress}} \right)} \\ {\frac{E}{{1 - \vartheta^{2} }} \quad \left( {\text{plane strain}} \right)} \\ \end{array} } \right.,$$
(25)

where E is Young’s modulus, θ is defined in the local polar coordinate system at the fracture tip, ν is Poisson’s ratio, KI is the mode I fracture toughness and KII is the mode II fracture toughness.

The hydraulic fracture interacts with the cemented natural fracture according to the same numerical model. Equation (22) is used to calculate the different angles of energy release, and then the different values of the critical Gc are compared. The von Mises stress distributions at different angles are shown in Fig. 13. When the interaction angle is 90°, the hydraulic fracture crosses the cemented natural fracture directly. When the interaction angle gradually decreases, the hydraulic fracture is arrested by the cemented natural fracture, and a new branch is then generated. The upper cemented natural fracture remains closed at the same time and there is no effect on the intersection.

Fig. 13
figure 13

Von Mises stress distributions at different angles of hydraulic fracture intersecting with cemented natural fracture

Full size image

Conclusions

  1. 1.

    A two-dimensional mixed-mode hydraulic fracture propagation model is established using the extended finite-element method (XFEM). The model considers fluid-soil coupling, fluid flow within the fracture, fluid leak-off and rock deformation. KI and KII are investigated in a semi-circular bend test and used in the mixed-mode hydraulic fracture propagation model.

  2. 2.

    Parameter sensitivity analysis of hydraulic fracture initiation and propagation is conducted. The maximum horizontal stress difference has little effect on the fracture length, and the minimum principal stress is the main factor limiting fracture growth. The lower the minimum principal stress is, the longer the fracture extension length is. A high elastic modulus and low Poisson’s ratio contribute more to the crack length. The higher the injection rate is, the longer the fracture growth length is.

  3. 3.

    The interactions between hydraulic fractures and two types of natural fractures at different angles are investigated. A smaller coefficient of the frictional natural fracture enables the hydraulic fracture to cross the natural fracture. A small interaction angle and low natural cementing fracture strength allow hydraulic fractures to propagate along cemented natural fractures.