Fracture propagation laws of staged hydraulic fracture in fractured geothermal reservoir based on phase field model

Genbo Peng

doi:10.1007/s40789-023-00636-y

Fracture propagation laws of staged hydraulic fracture in fractured geothermal reservoir based on phase field model

Download PDF

Research Article

Open Access

Published: 15 September 2023

0 Accesses

International Journal of Coal Science & Technology Volume 10, article number 52, (2023)

Genbo Peng ¹

1.

School of Environment, Liaoning University, Shenyang, China

Abstract

Hydraulic fracturing is widely used in geothermal resource exploitation, and many natural fractures exist in hot dry rock reservoirs due to in-situ stress and faults. However, the influence of natural fractures on hydraulic fracture propagation is not considered in the current study. In this paper, based on the phase field model, a thermo-hydro-mechanical coupled hydraulic fracture propagation model was established to reveal the influence of injection time, fracturing method, injection flow rate, and natural fracture distribution on the fracture propagation mechanism. The results show that fracture complexity increases with an increase in injection time. The stress disturbance causes the fracture initiation pressure of the second cluster significantly higher than that of the first and third clusters. The zipper-type fracturing method can reduce the degree of stress disturbance and increase fracture complexity by 7.2% compared to simultaneous hydraulic fracturing. Both low and high injection flow rate lead to a decrease in fracture propagation time, which is not conducive to an increase in fracture complexity. An increase in the natural fracture angle leads to hydraulic fracture crossing natural fracture, but has a lesser effect on fracture complexity. In this paper, we analyzed the influence of different factors on initiation pressure and fracture complexity, providing valuable guidance for the exploitation of geothermal resources.

1.Introduction

With the rapid development of economic globalization, human demand for energy is increasing. The long-term use of conventional energy sources such as coal, oil and natural gas in large quantities has led to the emergence of environmental problems such as the greenhouse effect and land acidification, and it is urgent to find clean and renewable alternative energy sources. Geothermal energy, as a clean and renewable energy source, has good potential for development, and its development and exploitation scale are increasing in countries all over the world (Kepinska 2003; Trumpy et al. 2016; Vonsée et al. 2019). Geothermal resources are mainly divided into two major categories: hydrothermal type geothermal resources and hot dry rock type geothermal resources. Hydrothermal type geothermal resources are buried in 200–3000 m, mainly stored in high permeability pores or fractured media. Hot dry rock type geothermal resources are buried in 3000–10,000 m with a temperature higher than 180 °C. Currently, hydrothermal type geothermal resources that can be effectively developed account for only about 10% of the explored geothermal resources, and more geothermal energy is stored in the hot dry rock type geothermal resources (Xu et al. 2018). The research of hot dry rock started in the 1970s, and after decades of research and development, it has shown great application value and good development prospects (Bruhn 2002; Górecki et al. 2003; Ziabakhsh-Ganji et al. 2018; Leveni et al. 2019; Anderson and Rezaie 2019). Morton Smith from Los Alamos National Lab was the first to propose the idea and concept for the development of hot dry rock, and its development method is the Enhanced Geothermal System (EGS), the use of hydraulic fracturing and other methods to improve permeability in hot dry rock reservoirs, thereby increasing the volume of heat exchange between the fluid and the high-temperature reservoir. Hot dry rock reservoirs have natural fractures and high temperatures, how to adequately connect the natural fractures is the key to the efficient construction of EGS.

In recent years, scholars have studied the propagation mechanism of hydraulic fractures under natural fractures through experiments. It is generally believed that four methods after the hydraulic fracture meets the natural fracture, such as crossing, opening, arrest and offset (Zhang et al. 2021a). However, the defects of the current experimental research include: (1) the mechanical parameters of the experimental samples are different, which have impact on the experimental results, and it is difficult to get the ideal experimental law. (2) The experimental loading conditions cannot fully conform to the real formation requirements, resulting in the deviation of experimental results from the real formation. (3) The size difference between the experimental and the formation is huge, and the size effect is obvious, and the complex working conditions such as multiple wells and multiple fractures cannot be simulated due to the size limitation.

To better study the fracture propagation mechanism, numerical simulation methods are widely used. For fracture propagation problems, discrete models and continuous are widely used. The discrete models including the boundary element method (BEM) (Rezaei et al. 2019; Cong et al. 2021, 2022a), the discrete element method (DEM) (Shimizu et al. 2011; Chamanzad et al. 2017; Cong et al. 2022b; Li et al. 2023) and the displacement discontinuity method (DDM) (Jiao et al. 2015). The discrete model can simulate the fracture initiation and propagation, but the matrix is assumed to be discontinuous, which is insufficient for the simulation of continuous media, and the fracture propagation results are affected by the grid size. The continuous model includes damage evolution law (Shojaei et al. 2014; Li et al. 2017, 2019) and phase field methods (Miehe et al. 2010b; Mikelić et al. 2015; Miehe and Mauthe 2016). The damage evolution law that examines the stress condition of rock using failure criteria in relation to fracture initiation and propagation. The phase field model describes the fracture initiation and propagation by establishing a variational-based energy minimization framework, this is a better way than damage evolution law (Li et al. 2018) by solving the challenges of fracture tip discontinuity, or grid size sensitivity for fracture propagation, and hence the phase field model has been widely used for fracturing simulation in porous media. Bourdin et al (2012) used a phase field model for the first time to simulate hydraulic fracturing by assuming that the material of the model is a homogeneous, impermeable, and nonporous reservoir. Wheeler et al (2014) applied phase field model on porous media by the augmented Lagrangian approach and filled the gap for hydraulic fracturing simulation using the phase-field models. Xia et al (2017) and Zhou et al (2019) considered the impact of material heterogeneity on fracture propagation in 2D and 3D phase field model and realized hydraulic fracturing simulation.

Based on the phase field model, this paper establishes a fracture propagation model for hydraulic fracturing of hot dry rocks of natural fracture development. The influence of fracturing method, injection flow rate and natural fracture distribution on fracture propagation and morphology are analyzed respectively. The model adopts the simulation technique that realizes the study of fracture propagation and meets natural fractures in the hydraulic fracturing process of EGS and explores the fracture propagation mechanism of hydraulic fracturing. The research results can provide a theoretical basis for the optimal design of EGS.

2.Mathematical model

2.1 Energy functional

According to Zhou et al. (2019) the total potential energy ψ of porous media rocks consists of kinetic energy ψ_k, elastic energy ψ_ε, fracture energy, external energy, fluid pressure, traction f and body force b. The total potential energy ψ is given by:

$$\begin{aligned} \psi \left( {u,\varGamma } \right) = & \;\int_{\varOmega }^{{}} {\psi _{k} } {\text{d}}\varOmega - \int_{\varOmega }^{{}} {\psi _{\varepsilon } } {\text{d}}\varOmega + \int_{\varOmega }^{{}} {\alpha p \cdot \left( {\nabla \cdot u} \right)} {\text{d}}\varOmega \\ & - \int_{\varGamma }^{{}} {G_{\text{c}} } {\text{d}}S + \int_{\varOmega }^{{}} {b \cdot u} {\text{d}}\varOmega + \int_{{\partial \varOmega _{{hi}} }}^{{}} {f \cdot u} {\text{d}}S \\ \end{aligned}$$

(1)

To distinguish fracture and rock matrix, defining ϕ(x, t) ∈ [0,1], for fracture domain ϕ = 1, for matrix domain ϕ = 0. Length scale parameter l₀ is used to control ϕ change from 0 to 1, and the outer boundary is ∂Ω, Ω is an arbitrary bounded computational domain. as shown in Fig. 1.

The kinetic energy ψ_k is:

$$\psi _{k} = \psi _{k} \left( {\dot{u}} \right) = \frac{1}{2}\rho \dot{u}_{i} \dot{u}_{i}$$

(2)

The elastic energy ψ_ε is:

$$\psi_{\varepsilon } \left( \varepsilon \right) = \frac{1}{2}\lambda \varepsilon_{ii} \varepsilon_{jj} + \mu \varepsilon_{ii} \varepsilon_{jj}$$

(3)

where λ and μ are Lame constants.

2.2 Phase field approximation

To divide the fractures and the rock matrix, ϕ(x, t) ∈ [0,1] is defined to divide the rock into fracture region (ϕ = 1) and matrix region (ϕ = 0), and the regularized fracture length l₀ is used to control the transition of ϕ from 0 to 1.

The fracture surface density per unit volume of rock is (Miehe et al. 2010a):

$$\gamma \left( {\phi ,\nabla \phi } \right) = \frac{{\phi^{2} }}{{2l_{0} }} + \frac{{l_{0} }}{2}\frac{\partial \phi }{{\partial x_{i} }}\frac{\partial \phi }{{\partial x_{i} }}$$

(4)

Fracture energy ψ_k is:

$$\int_{\varGamma }^{{}} {G_{\text{c}} } {\text{d}}S = \int_{\varOmega }^{{}} {G_{\text{c}} } \left[ {\frac{{\phi ^{2} }}{{2l_{0} }} + \frac{{l_{0} }}{2}\frac{{\partial \phi }}{{\partial x_{i} }}\frac{{\partial \phi }}{{\partial x_{i} }}} \right]{\text{d}}\varOmega$$

(5)

where, G_c is the critical energy release rate, N/m.

The elastic energy ψ_ε is decomposed into two components, tension and compression, with a strain tensor of the form ε = ε₊ + ε_-, where ε₊ and ε_- are the tensile and compressive tensors, respectively.

To avoid singularities in the elastic energy density, the model parameter k is defined and 0 < k << 1. Assuming that the phase field affects the stretching part of the elastic energy density, the elastic energy density is (Borden et al. 2012):

$$\psi_{\varepsilon } \left( \varepsilon \right) = \left[ {\left( {1 - k} \right)\left( {1 - \phi } \right)^{2} + k} \right]\psi_{\varepsilon }^{ + } \left( \varepsilon \right) + \psi_{\varepsilon }^{ - } \left( \varepsilon \right)$$

(6)

2.3 Governing equations for evolution of the phase field

By using Eqs. (2), (5) and (6), the expression for the energy can be obtained as:

$$\begin{aligned} L = \psi \left( {u,\varGamma } \right) = & \frac{1}{2}\int_{\varOmega }^{{}} {\rho \dot{u}_{i} \dot{u}_{i} } {\text{d}}\varOmega \\ & - \int_{\varOmega }^{{}} {\left\{ {\left[ {\left( {1 - k} \right)\left( {1 - \phi } \right)^{2} + k} \right]\psi _{\varepsilon }^{ + } \left( \varepsilon \right) + \psi _{\varepsilon }^{ - } \left( \varepsilon \right)} \right\}} {\text{d}}\varOmega \\ & + \int_{\varOmega }^{{}} {\alpha p \cdot \left( {\nabla \cdot u} \right)} {\text{d}}\varOmega - \int_{\varGamma }^{{}} {G_{c} } \left[ {\frac{{\phi ^{2} }}{{2l_{0} }} + \frac{{l_{0} }}{2}\frac{{\partial \phi }}{{\partial x_{i} }}\frac{{\partial \phi }}{{\partial x_{i} }}} \right]{\text{d}}\varOmega \\ & + \int_{\varOmega }^{{}} {b_{i} u_{i} } {\text{d}}\varOmega + \int_{{\partial \varOmega _{{hi}} }}^{{}} {f_{i} u_{i} } {\text{d}}S \\ \end{aligned}$$

(7)

Calculating the first-order variation of the energy functional L and setting its value to 0 to obtain:

$$\left\{ \begin{gathered} \frac{{\partial \sigma_{ij} }}{{\partial x_{i} }} + b_{i} = \rho \ddot{u}_{i} \hfill \\ \left[ {\frac{{2l_{0} \left( {1 - k} \right)\psi_{\varepsilon }^{ + } }}{{G_{\text{c}} }} + 1} \right]\phi - l_{0}^{2} \frac{{\partial^{2} \phi }}{{\partial x_{i}^{2} }} = \frac{{2l_{0} \left( {1 - k} \right)\psi_{\varepsilon }^{ + } }}{{G_{\text{c}} }} \hfill \\ \end{gathered} \right.$$

(8)

where: $\ddot{u}_{i} = \frac{{\partial^{2} u}}{{\partial t^{2} }}$, σ_ij are the stress components, which are calculated as:

$$\sigma_{ij} = \left[ {\left( {1 - k} \right)\left( {1 - \phi } \right)^{2} + k} \right]\frac{{\partial \psi_{\varepsilon }^{ + } }}{{\partial \varepsilon_{ij} }} + \frac{{\partial \psi_{\varepsilon }^{ - } }}{{\partial \varepsilon_{ij} }}$$

(9)

To prevent fractures already generated by the phase field model from reverting to an unopened state, the strain history field H needs to be introduced (Borden et al. 2012):

$$H\left( {x,t} \right) = \mathop {\max }\limits_{{s \in \left[ {0,t} \right]}} \psi_{\varepsilon }^{ + } \left( {\varepsilon \left( {x,s} \right)} \right)$$

(10)

substituting ψ_ε⁺ by H (x, t) in Eq. (10), the strong forms are rewritten as (Borden et al. 2012):

$$\left\{ \begin{gathered} \frac{{\partial \sigma_{ij}^{{}} }}{{\partial x_{i} }} + b_{i} = \rho \ddot{u}_{i} \hfill \\ \left[ {\frac{{2l_{0} \left( {1 - k} \right)H}}{{G_{\text{c}} }} + 1} \right]\phi - l_{0}^{2} \frac{{\partial^{2} \phi }}{{\partial x_{{}}^{2} }} = \frac{{2l_{0} \left( {1 - k} \right)H}}{{G_{\text{c}} }} \hfill \\ \end{gathered} \right.$$

(11)

2.4 Governing equations for fluid flow

Define χ_R and χ_F as two linear functions, where the matrix region χ_R = 1 and χ_F = 0; the fracture region χ_R = 0 and χ_F = 1; the transition region expression is (Lee et al. 2016) (Fig. 2):

$$\chi_{R} = \frac{{c_{2} - \phi }}{{c_{2} - c_{1} }}$$

(12)

$$\chi_{F} = \frac{{\phi - c_{1} }}{{c_{2} - c_{1} }}$$

(13)

where: c₁ and c₂ are the boundary of the transition region, the fracture region when ϕ ≤ c₁, the matrix region when ϕ ≥ c₂, and the transition region when c₁ < ϕ < c₂. The fluid and solid properties of the transition domain are obtained by interpolating the matrix and fracture domains using the functions χ_R and χ_F.

In this paper, we assume that the fluid flow within a porous medium conforms to Darcy's law, and the fluid pressure control equation is (Yu et al. 2019):

$$\rho S\frac{\partial p}{{\partial t}} - \nabla \cdot \left( {\rho v + \rho g} \right) = q_{\text{m}} - \rho \alpha \chi_{R} \frac{{\partial \varepsilon_{\text{vol}} }}{\partial t}$$

(14)

where, q_m is fluid source term, kg/(m³/ s); ε_vol = ∇·u is volumetric strain of Ω, dimensionless; g is gravitational acceleration, m/s²; p is pressure, Pa; ρ = ρ_Rχ_R + ρ_Fχ_F, ρ_R and ρ_F represent rock and fluid density respectively, kg/m³; α = α_Rχ_R + χ_F, α_R is rock Biot coefficient.

The Darcy velocity is expressed as (Xu et al. 2019):

$$v = - \frac{K}{\mu }\nabla p$$

(15)

where K = K_Rχ_R + K_Fχ_F, K_R and K_F represent rock and fracture permeability respectively, m²; μ is the fluid viscosity, Pa s.

The storage coefficient S is calculated by the formula (Zhou et al. 2019):

$$S = \varphi_{R} c + \frac{{\left( {\alpha - \varphi_{R} } \right)\left( {1 - \alpha } \right)}}{{K_{\text{vol}} }}$$

(16)

where S is the storage coefficient, 1/Pa; compressibility coefficients c = c_Rχ_R + c_Fχ_F, c_R and c_F denote rock and fluid compressibility coefficients, respectively; φ_R is the rock porosity; and K_vol is the bulk modulus.

Therefore, the fluid pressure control equation is rewritten as:

$$\rho S\frac{\partial p}{{\partial t}} - \nabla \cdot \frac{\rho K}{\mu }\left( {\nabla p + \rho g} \right) = q_{\text{m}} - \rho \alpha \chi_{R} \frac{{\partial \varepsilon_{\text{vol}} }}{\partial t}$$

(17)

2.5 Thermal stress calculation equations

The heat exchange between fluid and rock within the fracture generates thermal stress that effect fracture initiation and propagation. The thermal stress due to temperature changes are calculated (Zhang et al. 2021b).

$$\sigma_{ij} = \lambda \varepsilon_{\text{vol}} \delta_{ij} + 2\mu \varepsilon_{ij} + 3\alpha_{T} K_{\text{vol}} \left( {T - T_{0} } \right)\delta_{ij}$$

(18)

where δ_ij is the Kronecker symbol, uncaused; α_T is the coefficient of linear thermal expansion, 1/K; T₀ is the initial temperature of the formation, K.

The equation for heat transfer between fluid and rock within a porous medium is (Li. 2018):

$$c_{p,eff} \frac{\partial T}{{\partial t}} + \rho c_{p} \frac{{\partial (Tu_{i} )}}{{\partial x_{i} }} + \frac{\partial }{{\partial x_{i} }}\left( { - \lambda_{eff} \frac{\partial T}{{\partial x_{i} }}} \right) = Q$$

(19)

where, c_p,eff is equivalent isobaric specific heat capacity, J/(m³ K); λ_eff is equivalent thermal conductivity, W/(m K); Q is heat source, W/m³; c_p,eff and λ_eff are given by:

$$c_{\text{p,eff}} = \rho_{R} c_{p,R} (1 - \varphi_{R} ) + \rho_{F} c_{p,F} \varphi_{R}$$

(20)

$$\lambda_{\text{eff}} = \lambda_{R} (1 - \varphi_{R} ) + \lambda_{F} \varphi_{R}$$

(21)

where, c_p,R and c_p,F represent rock and fluid isobaric specific heat capacity, respectively, J/m³ K; λ_R and λ_F represent rock and fluid thermal conductivity, W/(m K).

2.6 Validation

To verify the accuracy of the hydraulic fracturing model, we compare the simulation result with the experimental result carried out by Zhang et al (2017). The established physical model is shown in Fig. 3a, the model size is 0.2 m × 0.2 m, the stress of 10 MPa and 8 MPa are loaded at the top and right boundaries of the model, the left and bottom boundaries are supported by rolls, the injection fluid is water, the injection flow rate is 30 mL/s, the initial pressure is 0.1 MPa, and the specimen is heated at 333.15 K for 24 h before the experiment. Therefore, the fracture fluid and specimen temperature are assumed to be 333.15 K. The radius of the circular hole is 0.0075 m. The specific simulation parameters are shown in Table 1.

Table 1 Parameter setting of phase field model fracturing simulation verification (Feng et al. 2021)

Parameter	Value	Parameter	Value
c₁	0.4	c₂	0.6
Young’s modulus E	33 GPa	Porosity φ_R	0.01
Poisson’s ratio N	0.35	Biot coefficient α	0.01
Critical energy release rate G_c	130 N/m	Permeability K_R	0.001 mD

The hydraulic fracturing experimental fracture morphology of Zhang et al (2017) and the fracture morphology of numerical simulation in this paper are shown in Fig. 3b and c. Due to the stress difference is 2 MPa, the simulated results and the experimental results fractures all propagation along the direction of the maximum principal stress, and the fracture morphology and propagation paths are the same. By comparing with the experimental results, it is verified that the model in this paper has good accuracy.

3.Results and discussion

In this paper, a two-dimensional model is used to simulate the hydraulic fracture propagation problem during geothermal mining. The boundary loads of 2 MPa and 1 MPa are applied in the x-direction and y-direction, respectively, and other boundaries are fixed. The fracturing fluid is water, and the model has two wells with one injection and one production, and the fracturing method is zipper-type hydraulic fracturing of staged fracturing, each well is divided into two stages, each stage has three clusters, and the cluster spacing is 12 m. The natural fractures in the reservoir are randomly distributed, and the natural fracture length ranges from 8 to 12 m. The model size is 120 m × 120 m. The grid number of the model is 59254. The specific model parameters are shown in Table 2 (Figs. 4, 5).

Table 2 Fracturing fluid and formation parameter data (Li et al. 2022)

Parameter	Value	Parameter	Value
Young’s modulus E	20 GPa	Poisson's ratio N	0.3
Fracture length	5 m	Length and width	120 m
Stage spacing	12 m	Cluster spacing	12 m
Critical energy release rate of interface	183 N/m	Critical energy release rate of bulk	500 N/m
σ_y	1 MPa	σ_x	2 MPa
Permeability K_R	1 × 10⁻¹⁸m²	Fluid viscosity	0.001 Pa s
Biot coefficient α	0.05	Porosity φ_R	0.05
Fluid density ρ	1000 kg/m³

3.1 Hydraulic fracture propagation morphology study

This section analyzes the fracture morphology during hydraulic fracturing. The model is fractured by zipper-type hydraulic fracture, and the fracturing sequence is the first stage of the first well, the first stage of the second well, the second stage of the first well, and the second stage of the second well, respectively. The injection flow rate is 8 kg/(m³ s), and the fracturing time of each stage is 4 s. Figure 6 shows the fracture propagation at different times.

By comparing the fracture propagation at different times, with the increase of time, the first and second stages hydraulic fractures propagate and meet with the natural fractures, and then the hydraulic fractures open the natural fractures and propagate along the natural fracture direction. When the fracture propagates to the tip of the natural fracture, the hydraulic fracture deflects in the direction of the maximum principal stress due to the stress difference. When the third and fourth stages are fractured, the inter-stage and inter-cluster stress interference exist, resulting in difficulty in propagate the hydraulic fracture along the natural fracture direction when it meets the natural fracture, at which the hydraulic fracture cross the natural fracture. To compare the fracture complexity, the ratio of the number of fracture grids (ϕ > 0.95) to the total number of grids in the model is defined as the fracture area ratio according to Li et al (2022), and the fracture complexity can be reflected according to the fracture area ratio. Figures 7 and 8 reflect the fracture initiation pressure per cluster per stage and the fracture complexity, respectively.

From the fracture initiation pressure curve of Fig. 7, the fracture initiation pressure of second cluster is significantly higher than that of first and third clusters, which is caused by the inter-cluster interference during fracture propagation. With the increase of fracture stages, there is a small increase of fracture initiation pressure. Comparing the fracture complexity (Fig. 8), the fracture complexity gradually increases with the increase of fracturing time, but the growth of fracture complexity gradually decreases. Combined with the fracture morphology cloud diagram (Fig. 6), the fracture propagation is mainly concentrated in the first and second stages, and the fractures in the third and fourth stages are propagated much less than the first and second stages. In the second stage propagation, the fracture complexity increases the most significantly by 0.058, but in the third and fourth stages, the fracture complexity only increases 0.014 and 0.009, respectively. Analyzing the reasons for this, we can see that: (1) in the second stage of propagation, the pressure in the first stage still causes the fracture to continue to propagate; (2) the hydraulic fracture activates many natural fractures during the second stage propagation, and the above two reasons work together to increase the fracture complexity in the second stage. From third stage to fourth stage, from 8 to 16 s, the fracture complexity increases by only 0.023, which is due to the higher stress disturbance leads to an increase in the fracture initiation pressure every cluster. The increase in the fracture initiation time, and a decrease in the fracture propagation time, which eventually leads to a decrease in the growth of the fracture complexity.

3.2 Study on the fracture propagation of different fracturing methods

This section focuses on the effect of zipper-type hydraulic fracturing and simultaneous hydraulic fracturing on the fracture propagation morphology. The injection flow rate is 8 kg/(m³ s), and the injection time is 4 s per stage. The schematic diagrams of the two fracturing methods are shown in Figs. 4 and 9, and the fracturing fracture morphology is shown in Fig. 10.

The fracture propagation morphology of the two fracturing methods are basically the same, both mainly along the direction of maximum principal stress, and there is connection between stages and clusters, the difference is that the fracture propagation of zipper-type hydraulic fracturing is more adequate, the fracture length is longer, and the fracture propagation of the second and third stages of zipper-type hydraulic fracturing is higher than that of simultaneous hydraulic fracturing.

To compare the differences between the two fracturing methods, we investigate the differences between the fracturing methods and fracture complexity, as shown in Figs. 11 and 12. To visually compare the differences in fracture initiation pressure, we take the average of the three clusters of initiation pressure as the fracture initiation pressure of the stage and compare the fracture complexity at the end of the second fracturing stage and the end of the fourth fracturing stage of zipper-type hydraulic fracturing with simultaneous hydraulic fracturing. By comparing the fracture initiation pressure (Fig. 11), the simultaneous hydraulic fracturing has a higher fracture initiation pressure than zipper-type hydraulic fracturing, which in turn affects the fracture complexity less than zipper-type hydraulic fracturing. The third and fourth stages of simultaneous hydraulic fracturing have significantly higher pressure than the first and second stages due to inter-stage stress interference.

3.3 Study on fracture propagation morphology by injection flow rate

The injection flow rate is increased from 4 to 10 kg/(m³ s), and the fracturing method is zipper-type hydraulic fracturing. To ensure the same injection volume under different injection flows, each injection time is set to t = 8, 5.3, 4, and 3.2 s, respectively. The fracture morphology after hydraulic fracturing is shown in Fig. 13.

Similarly, we choose the average value of three clusters of fracture initiation pressure as the fracture initiation pressure of this stage. The comparison of fracture initiation pressure and fracture complexity under different injection flow rate conditions are shown in Figs. 14 and 15, respectively. The comparison of fracture morphology at different injection flow rates shows (Fig. 14) that the difference in fracture morphology and orientation is small. Comparing with the injection flow rate of 4 kg/(m³ s), it is easier to open the natural fractures and connect with each other when the injection flow rate is 6 kg/(m³ s) and 8 kg/(m³ s). When the injection flow rate is 10 kg/(m³ s), the first cluster of the second stage cross the natural fractures and connects with the first cluster of first stage, therefore, the increased injection flow rate of hydraulic fractures easily leads to the propagation of hydraulic fractures cross the natural fractures.

The first and second stages fracture initiation pressure increases with higher injection flow rate due to the higher injection flow rate in the fracture will result in higher frictional, therefore, higher pressure at the injection point required to reach the initiation pressure at the fracture tip and higher injection flow rates require higher pressures. However, the inter-stage interference between the third and fourth stages leads to higher fracture initiation pressure at injection flow rates of 4 kg/(m³ s) and 6 kg/(m³ s) initiation pressure than at injection flow rate of 8 kg/(m³ s). It takes shorter time to reach the initiation pressure at injection flow rate of 8 kg/(m³ s), which in turn leads to the highest fracture complexity (Fig. 15). More fracture complexity can be obtained with a medium injection flow rate, which is also consistent with Tan's experimental results (Tan et al. 2017). Although the injection time is longer at an injection flow rate of 4 kg/(m³ s), the longer time to reach the fracture initiation pressure results in a shorter hydraulic fracture propagation time, therefore, slower fracture propagation leading to a lower fracture complexity. As the injection flow rate increases, the hydraulic fracture propagation time is longer, and the fracture complexity is higher. Due to the shorter simulation time at injection flow rate of 10 kg/(m³ s), the hydraulic fracture propagation time is shorter, resulting in insufficient fracture propagation and lower fracture complexity.

3.4 Study of natural fracture distribution on fracture propagation

This section mainly analyzes the effect of natural fracture distribution pattern on fracture propagation. The natural fracture angles are 0° and 90°, respectively. The schematic diagram of natural fracture distribution pattern is shown in Figs. 16 and 17.

By comparing the natural fracture random distribution, natural fracture angles are 0° and natural fracture angles are 90° fracture patterns (Fig. 18a and c), it can be seen that: the natural fracture distribution pattern has a great influence on the fracture propagation. When the natural fracture angles are 0°, a large number of natural fractures are opened by the hydraulic fractures and propagated along the natural fracture direction because the natural fracture direction is the same as the maximum principal stress direction, while when the natural fracture angles are 90°, the hydraulic fractures cross the natural fracture obviously, and almost all the fractures in the second and fourth stages cross the natural fractures and propagate along the maximum principal stress direction. And when the natural fracture is random distributed, the randomness of hydraulic fracture propagation is significantly higher than that when the angles are 0° and 90° (Fig. 18).

The fracture initiation pressure in the second and third stages shows natural fracture angle 0° > random distribution > natural fracture angle 90°, which is caused by the increased stress disturbance due to the horizontal propagation of the fracture. When the fourth stage initiation fracturing, the natural fracture angle 0° < random distribution < natural fracture angle 90°, and it can also be seen from the fracture morphology that the fracture propagation is higher when the natural fractures are 0° than when the natural fractures are 90° (Fig. 19).

The comparison of the fracture complexity shows (Fig. 20) that the fracture complexity is higher than the random distribution when the natural fractures are 0° and 90°. This is because both 0° and 90° open the natural fractures to a higher degree than random fractures, resulting in a slightly higher final fracture complexity than the random distribution. The fracture complexity is basically the same when the natural fractures are 0° and 90°. This is because when the natural fractures are 0°, the fracture mainly propagates horizontally, and the fracture length is longer, but the fracture width is smaller. However, when the natural fractures are 90°, the first and third stages of fractures propagate along the natural fracture direction due to the influence of the natural fracture, and the fractures continue to propagate along the natural fracture direction with difficulty due to the stress difference, and the fractures are deflected to the direction of the maximum principal stress, resulting in smaller fracture lengths but increased fracture widths. The fracture complexity is basically the same when the natural fractures are 0° and 90°.

4.Conclusions

In this paper, based on the phase field model, a hydraulic fracture propagation model is constructed to study the influences of fracturing time, fracturing method, injection flow rate and natural fracture distribution on fracture morphology, fracture initiation pressure and fracture complexity, and the following conclusions are obtained:

(1)
With the increase of fracturing time, the fracture initiation pressure of each cluster increases, but due to the inter-cluster interference problem, the fracture initiation pressure of the second cluster of each stage is higher than that of the first and third clusters; the fracture complexity increases with the increase of time, but the growth rate gradually decreases.
(2)
Compared with simultaneous fracturing, the use of zipper-type fracturing can reduce the fracture initiation pressure and obtain a higher fracture complexity. The average fracture initiation pressure is reduced by 2.9% and the fracture complexity is increased by 7.2% for the zipper-type fracturing.
(3)
The fracture complexity increases first and then decreases, reaching the maximum at the injection flow rate of 8 kg/(m³ s). Both low and high injection flow rates result in short hydraulic fracture propagation time and lead to low fracture complexity.
(4)
When the natural fractures are horizontal to the direction of the maximum principal stress, the hydraulic fractures are more likely to open the natural fractures and propagate along the natural fractures, and when the natural fractures are perpendicular to the direction of the maximum principal stress, the hydraulic fractures are easy to propagate cross the natural fractures. The natural fracture distribution has less influence on the fracture complexity.

References

[1]	Anderson A, Rezaie B (2019) Geothermal technology: trends and potential role in a sustainable future. Appl Energy 248:18–34
[2]	Borden M, Verhoosel C, Scot M, Hughes T, Landis C (2012) A phase-field description of dynamic brittle fracture. Comput Methods Appl Mech Eng 217:77–95
[3]	Bourdin B, Chukwudozie C, Yoshioka K (2012) A variational approach to the numerical simulation of hydraulic fracturing. In: SPE annual technical conference and exhibition. OnePetro
[4]	Bruhn M (2002) Hybrid geothermal–fossil electricity generation from low enthalpy geothermal resources: geothermal feedwater preheating in conventional power plants. Energy 27:329–346
[5]	Chamanzad M, Ramezanzadeh A, Tokhmechi B, Norouzi H (2017) Comparison of different hydraulic fracture growth models based on a carbonate reservoir in Iran. J Chem Petrol Eng 51:95–104
[6]	Cong Z, Li Y, Liu Y, Xiao Y (2021) A new method for calculating the direction of fracture propagation by stress numerical search based on the displacement discontinuity method. Comput Geotech 140:104482
[7]	Cong Z, Li Y, Pan Y, Liu B, Shi Y, Wei J, Li W (2022a) Study on CO₂ foam fracturing model and fracture propagation simulation. Energy 238:121778
[8]	Cong Z, Li Y, Tang J, Martyushev A, Yang F (2022b) Numerical simulation of hydraulic fracture height layer-through propagation based on three-dimensional lattice method. Eng Fract Mech 264:108331
[9]	Feng Y, Haugen K, Firoozabadi A (2021) Phase-field simulation of hydraulic fracturing by CO₂, water and nitrogen in 2D and comparison with laboratory data. JGR Solid Earth 126:e2021JB022509
[10]	Górecki W, Kozdra T, Kuźniak T, Myśko A, Strzetelski W (2003) Geothermal-energy resources in the Polish Lowlands and the possibility of their industrial utilization. Appl Energy 74:53–64
[11]	Jiao Y, Zhang H, Zhang X, Zhang L, Li B, Jiang H (2015) A two-dimensional coupled hydromechanical discontinue model for simulating rock hydraulic fracturing. Int J Numer Anal Meth Geomech 39:457–481
[12]	Kepinska B (2003) Current state and prospects of geothermal-energy implementation in Poland. Appl Energy 74:43–51
[13]	Lee S, Wheeler MF, Wick T (2016) Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model. Comput Methods Appl Mech Eng 305:111–132
[14]	Leveni M, Manfrida G, Cozzolino R, Mendecka B (2019) Energy and exergy analysis of cold and power production from the geothermal reservoir of Torre Alfina. Energy 180:807–818
[15]	Li Y, Jia D, Rui Z, Peng J, Fu C, Zhang J (2017) Evaluation method of rock brittleness based on statistical constitutive relations for rock damage. J Petrol Sci Eng 153:123–132
[16]	Li P, Zhu Q, Gu S, Ni T (2018) A phase field method to simulate crack nucleation and crack propagation in rock-like materials. Eng Mech 35:41–48
[17]	Li Y, Long M, Zuo L, Li W, Zhao W (2019) Brittleness evaluation of coal based on statistical damage and energy evolution theory. J Petrol Sci Eng 172:753–763
[18]	Li X, Hofmann H, Yoshioka K, Luo Y, Liang Y (2022) Phase-field modelling of interactions between hydraulic fractures and natural fractures. Rock Mech Rock Eng 55:6227–6247
[19]	Li Y, Hubuqin WuJ, Zhang J, Yang H, Zeng B, Xiao Y, Liu J (2023) Optimization method of oriented perforation parameters improving uneven fractures initiation for horizontal well fracturing. Fuel 349:128754
[20]	Li X (2018) Research on flow-heat transfer and rock damage during CO₂ fracturing based on multiphysics coupling. China University of Petroleum (Beijing)
[21]	Miehe C, Mauthe S (2016) Phase field modeling of fracture in multi-physics problems. Part III. Crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media. Comput Methods Appl Mech Eng 304:619–655
[22]	Miehe C, Hofacker M, Welschinger F (2010a) A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Eng 199(45–48):2765–2778
[23]	Miehe C, Welschinger F, Hofacker M (2010b) Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. Int J Numer Meth Engng 83:1273–1311
[24]	Mikelić A, Wheeler MF, Wick T (2015) Phase-field modeling of a fluid-driven fracture in a poroelastic medium. Comput Geosci 19:1171–1195
[25]	Rezaei A, Siddiqui F, Bornia G, Soliman M (2019) Applications of the fast multipole fully coupled poroelastic displacement discontinuity method to hydraulic fracturing problems. J Comput Phys 399:108955
[26]	Shimizu H, Murata S, Ishida T (2011) The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution. Int J Rock Mech Min Sci 48:712–727
[27]	Shojaei A, Dahi Taleghani A, Li G (2014) A continuum damage failure model for hydraulic fracturing of porous rocks. Int J Plast 59:199–212
[28]	Tan P, Jin Y, Han K, Hou B, Chen M, Guo X, Gao J (2017) Analysis of hydraulic fracture initiation and vertical propagation behavior in laminated shale formation. Fuel 206:482–493
[29]	Trumpy E, Botteghi S, Caiozzi F, Donato A, Gola G, Montanari D, Pluymaekers MPD, Santilano A, van Wees JD, Manzella A (2016) Geothermal potential assessment for a low carbon strategy: A new systematic approach applied in southern Italy. Energy 103:167–181
[30]	Vonsée B, Crijns-Graus W, Liu W (2019) Energy technology dependence: a value chain analysis of geothermal power in the EU. Energy 178:419–435
[31]	Wheeler MF, Wick T, Wollner W (2014) An augmented-Lagrangian method for the phase-field approach for pressurized fractures. Comput Methods Appl Mech Eng 271:69–85
[32]	Xia L, Yvonnet J, Ghabezloo S (2017) Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media. Eng Fract Mech 186:158–180
[33]	Xu T, Hu Z, Li S et al (2018) Enhanced geothermal system: International progresses and research status of China. Acta Geol Sin 92(9):1936–1947
[34]	Xu Z, Wu K, Song X, Li G, Zhu Z, Sun B (2019) A unified model to predict flowing pressure and temperature distributions in horizontal wellbores for different energized fracturing fluids. SPE J 24:834–856
[35]	Yu H, Zhu Y, Jin X, Liu H, Wu H (2019) Multiscale simulations of shale gas transport in micro/nano-porous shale matrix considering pore structure influence. J Nat Gas Sci Eng 64:28–40
[36]	Zhang X, Lu Y, Tang J, Zhou Z, Liao Y (2017) Experimental study on fracture initiation and propagation in shale using supercritical carbon dioxide fracturing. Fuel 190:370–378
[37]	Zhang J, Li Y, Pan Y, Wang X, Yan M, Shi X, Zhou X, Li H (2021a) Experiments and analysis on the influence of multiple closed cemented natural fractures on hydraulic fracture propagation in a tight sandstone reservoir. Eng Geol 281:105981
[38]	Zhang W, Wang C, Guo T, He J, Zhang L, Chen S, Qu Z (2021b) Study on the cracking mechanism of hydraulic and supercritical CO₂ fracturing in hot dry rock under thermal stress. Energy 221:119886
[39]	Zhou S, Zhuang X, Rabczuk T (2019) Phase-field modeling of fluid-driven dynamic cracking in porous media. Comput Methods Appl Mech Eng 350:169–198
[40]	Ziabakhsh-Ganji Z, Nick HM, Donselaar ME, Bruhn DF (2018) Synergy potential for oil and geothermal energy exploitation. Appl Energy 212:1433–1447

About this article

Cite this article

Peng, G. Fracture propagation laws of staged hydraulic fracture in fractured geothermal reservoir based on phase field model.Int J Coal Sci Technol 10, 52 (2023).

https://doi.org/10.1007/s40789-023-00636-y

Received

13 June 2023
Revised

28 July 2023
Accepted

10 August 2023
DOI

https://doi.org/10.1007/s40789-023-00636-y
Share this article

Copy to clipboard

About issue

Editors-in-Chief

Managing Editor

Associate Editors

Publishing model

Fracture propagation laws of staged hydraulic fracture in fractured geothermal reservoir based on phase field model

Abstract

1.Introduction

2.Mathematical model

2.1 Energy functional

2.2 Phase field approximation

2.3 Governing equations for evolution of the phase field

2.4 Governing equations for fluid flow

2.5 Thermal stress calculation equations

2.6 Validation

3.Results and discussion

3.1 Hydraulic fracture propagation morphology study

3.2 Study on the fracture propagation of different fracturing methods

3.3 Study on fracture propagation morphology by injection flow rate

3.4 Study of natural fracture distribution on fracture propagation

4.Conclusions

References

About this article

Cite this article

Share this article

Keywords

For Authors

Explore