A novel hydraulic fracturing model for the fluid-driven fracture propagation in poroelastic media containing the natural cave

Jin, Jie; Wang, Xiaoqiang; Liu, Xiaohua; Xu, Yanmei; Lu, Detang

doi:10.1063/5.0160672

Hydraulic fracturing is an efficient technology to extract hydrocarbon within natural caves. However, these caves can markedly affect the fracture propagation behavior. This paper proposes a novel hydraulic fracturing model to simulate the fracture propagation in poroelastic media containing the natural cave, utilizing the strengths of the phase-field method. By coupling the Reynolds flow with cubic law in fracture domain, free flow in cave domain, and low-permeability Darcy flow in reservoir domain, the fracture-cave-reservoir flow governing equations are established. The Biot poroelasticity theory and fracture width are the links of hydro-mechanical coupling. The smooth phase-field is introduced to diffuse not only the sharp fracture but also the sharp cave edge. The fully coupling model is solved by a staggered scheme, which independently solves the pressure field and displacement field in inner cycle, and then independently solves the phase field in outer cycle. The proposed model is verified by comparing with the Khristianovic–Geertsma–de Klerk (KGD) model and Cheng's hydraulic fracturing model. Then, the interaction between hydraulic fracture and natural cave is investigated through several two-dimensional and three-dimensional cases. The result shows that the cave effect can make the hydraulic fracture deflect and raise its propagation velocity. Increasing the fracture-cave distance, injection rate, and in situ stress difference can all decline the cave effect. The displayed cases also substantiate the capability and efficiency of the proposed model.

I. INTRODUCTION

Saturated geomaterials, such as carbonate rocks, often contain the natural caves, which results in multiscale pore structures. The geometry of pore structures can influence the mechanical properties of the geomaterials and also cause complex fracture growth path.^1–3 For example, fracture-caved carbonate reservoirs are a common kind of carbonate reservoirs.^4,5 They are composed of fractures, caves, and tight matrix, which leads to strong heterogeneity.^6,7 Oil/gas is primarily stored within the natural caves. Hydraulic fracturing is an efficient technology to stimulate productivity of fracture-caved carbonate reservoirs, because it can connect the wellbore and natural caves by fluid-driven fracture.^8,9

Hydraulic fracture (HF) propagation is an extremely complicated process, including branching, deflection, merging, and so on. The existence of natural cave can cause strong stress disturbance, which will play a noticeable effect on the HF propagation behavior.^10–12 However, the inter-action mechanism between HF and natural cave is still vague. Moreover, according to reports in recent years, the hydraulic fracturing may create unfavorable transport channel accidentally. This circumstance can increase the possibility of groundwater being contaminated.^13,14 Therefore, accurately predicting the HF propagation behavior in the existence of natural cave is essential and challenging, which contributes to the economic benefits and environmental protection. The primary objective of this paper is to develop a reliable numerical model to address this challenge.

Even now, various numerical techniques are developed to simulate the fracture initiation and propagation. In general, these techniques can be divided into discrete and continuous approaches. In the discrete approaches, the displacement discontinuity is directly or indirectly introduced into the computation. The extended finite element method (XFEM) is a widely used discrete approach.^15–18 It adds the enrichment functions to depict the displacement discontinuity of fracture. In addition, the popular discrete approaches include cohesive element method,¹⁹ discrete element method (DEM),^20,21 displacement discontinuity method (DDM),^22,23 and so on. In the continuous approaches, the intact region and fractured region are viewed as a whole, and the crack is smeared at a certain area.^24,25 In recent years, the phase-field method (PFM) has become the most common continuous approach. PFM originated in the late 1990s. Francfort and Marigo²⁶ developed a variational model for the quasi-static cracks evolution. After that, Bourdin et al.^27,28 established the early phase-field model based on the variational framework. Then, Miehe et al.²⁹ established a solid–fluid coupling hydraulic fracturing model in porous media by PFM. So far, PFM has been utilized by many academic groups to establish hydraulic fracturing model.^30–34 The highlight of PFM is that it applies a continuous scalar field to deal with the sharp fracture. Hence, the complex fracture evolution is governed by the scalar field. This treatment can eliminate the laborious tracing procedure of the fracture topological structure, which helps reduce the complexity of numerical implementation. PFM is constructed on the thermodynamic framework, and it is easily expanded to the multi-physical field problem.^35–37 The strengths of PFM include the following two aspects: first, it can naturally detect and simulate the complex fracture evolution behavior, such as deflection, branching, and coalescence. The additional treatment for stress intensity factor is not needed. Second, PFM is conducted on the fixed topology, which can avoid remeshing.

In the past few years, some academic groups have studies HF propagation behavior with cave. Youn³⁸ established a hydro-mechanical coupling model to study HF propagation in the existence of cave by XFEM. He used the level-set method to characterize the fracture and cave. The result indicated that HF will deflect and spread along the cave when it approaches the cave. After that, Cheng et al.³⁹ improved Youn's model by adding a ramp function in blending nodes. They studied the influence of cave location and stress difference for the interaction mechanism. However, in Youn's and Cheng's models, the fluid flow in the cave is not considered. Luo et al.⁴⁰ proposed that the cave is not empty and there is fluid or cement in it. They established a free flow-seepage-mechanical coupling model by XFEM. The simulation result revealed that the fluid pressure in cave can affect the HF propagation path. Liu et al.⁴¹ studied HF propagation with natural cave by PFM. They changed Young's modulus of cave to distinguish between cave and matrix. However, the matrix and cave are both impermeable in their model. Azarov et al.⁴² investigated the HF growth in a poroelastic media containing a cave. The result showed the stress arising near the cave can change the HF growth path. Yet, their model did not couple the fluid in the cave. Liu et al.⁴³ simulated the HF evolution with natural fractures and natural cave by the element partition method. The cave is represented by an ellipse equation. Overall, the fluid flow in the cave is not considered in most studies. Therefore, it is necessary to couple the cave fluid in the hydraulic fracturing model. In addition, the HF propagation behavior is very complex near the cave, which results in great difficult to treat the stress intensity factor. If expand to the three-dimensional issue, it becomes more complicated, and there are a few studies now. Hence, we employ the PFM to establish our hydraulic fracturing model, due to its strengths in dealing with the complex fracture propagation behavior.

In this study, we establish a novel hydraulic fracturing model for the fluid-driven fracture propagation in poroelastic media containing the natural cave. First, the fracture-cave-reservoir flow equations are derived. In fracture domain, the fluid flow obeys Reynolds flow with cubic law; in cave domain, it conforms to free flow, while in reservoir domain, it follows low-permeability Darcy flow. Then, the PFM is used to diffuse the sharp cave edge and fracture. We reshape the total energy functional by adding the work done by fracturing fluid pressure and considering the pore pressure effect in elastic strain energy term. The coupled equations are discretized by the finite element method and solved by a staggered scheme. After that, the interaction mechanism between HF and natural cave is investigated in detail. The strong capability of the proposed model in dealing with 2D, 3D, and multiple cave issues is presented.

This paper is outlined as follows. The mathematical models and governing equations are presented in Sec. II. Then, the discrete process and numerical implementation scheme are shown in Sec. III. Next, the proposed model is verified in Sec. IV. After that, the interaction mechanism is studied by some numerical examples in Sec. V. The performance of the proposed model in multiple cave issues is displayed in Sec. VI. Finally, the concluding remarks of this work are drawn in Sec. VII.

II. MATHEMATICAL MODELS

A poroelastic domain

Ω

is shown in the Fig. 1. It contains three subdomains, namely, reservoir domain

Ω_{r}

⁠, fracture domain

Γ

⁠, and cave domain

Ω_{c}

⁠. The poroelastic domain

Ω

is in

R^{n}

(⁠

n

= 2 or 3). The reservoir domain

Ω_{r}

is a undamaged domain and

Ω_{r} \subset R^{n}

⁠. The fracture domain

Γ

is discontinuous, which has two surfaces

Γ^{+}

and

Γ^{-}

⁠. In mathematical modeling, the fracture is assumed to coincide and

n_{Γ^{+}} = - n_{Γ^{-}} = n_{Γ}

⁠. The fracture domain

Γ

is deemed as an internal boundary of

Ω

⁠. The geometry of natural cave is usually irregular. To simplify the mathematical modeling, the geometry of cave is assumed as circle in two-dimensional space and sphere in three-dimensional space. In addition, there are two outer boundaries, which are the Dirichlet boundary

\partial_{D} Ω

and the Neumann boundary

\partial_{N} Ω

⁠. It is worth noting that the boundary conditions on

\partial_{D} Ω

and

\partial_{N} Ω

are different in fluid flow model and mechanical model. For the fluid flow model in Sec. II A, the boundary conditions are

p_{r} = p_{n} on \partial_{D} Ω,

(1)

V_{r} \cdot n = q_{n} on \partial_{N} Ω,

(2)

where

p_{r}

is the fluid pressure in the reservoir,

p_{n}

is the prescribed pressure on boundary,

V_{r}

is the fluid flow velocity in the reservoir,

n

is outside normal vector of the poroelastic domain, and

q_{n}

is the prescribed fluid flux on boundary.

FIG. 1.

View large Download slide

The schematic of poroelastic domain.

For the mechanical model in Sec. II D, a specified displacement $u^{*}$ is imposed on the Dirichlet boundary $\partial_{D} Ω$ ⁠. Meanwhile, a specified traction force $τ$ is imposed on the Neumann boundary $\partial_{N} Ω$ ⁠.

Before obtaining the governing equations, some assumptions are made as follows:

The poroelastic domain is linearly elastic and isotropic.
There is only single-phase fluid in the poroelastic domain.
The fluid is viscous and slightly compressible.

A. The reservoir flow model

First, we establish the flow model in the reservoir. The fluid infiltration and transportation in the reservoir follow Darcy's law.^44–46 Hence, the fluid flow velocity is

V_{r} = - \frac{k_{r}}{μ} \nabla p_{r},

(3)

where

V_{r}

is the fluid flow velocity in the reservoir,

k_{r}

is the reservoir permeability,

μ

is the fluid viscosity, and

p_{r}

is the fluid pressure in the reservoir.

Biot⁴⁷ first proposed the interaction between solid deformation and fluid flow. Then, some scholars^25,30 developed his theory. According to their work, the Biot effective stress and volumetric strain are introduced to the mass conservation equation. For the slightly compressible fluid in the reservoir, the corresponding equation is as follows:

\frac{1}{M} \frac{\partial p_{r}}{\partial t} + α \frac{\partial ε_{vol}}{\partial t} + \nabla \cdot V_{r} = q_{r},

(4)

where

M

is Biot modulus,

t

is the time,

α

is Biot coefficient, and

ε_{vol}

is volumetric strain.

q_{r}

is the fluid source or sink term in the reservoir, and its unit is volumetric flow velocity per unit volume.

The Biot modulus and volumetric strain are given as follows:

\frac{1}{M} = \frac{α - ϕ}{K_{s}} + \frac{ϕ}{K_{f}},

(5)

ε_{vol} = ε_{x} + ε_{y} + ε_{z} = \nabla \cdot u,

(6)

where

ϕ

is the porosity,

K_{s}

is the rock bulk modulus,

K_{f}

is the fluid bulk modulus, and

u

is the displacement vector.

Then, upon substituting Eq. into Eq. , the mass conservation equation becomes

\frac{1}{M} \frac{\partial p_{r}}{\partial t} - \nabla \cdot \frac{k_{r}}{μ} \nabla p_{r} = q_{r} - α \frac{\partial ε_{vol}}{\partial t} .

(7)

In this paper, we use the finite element method (FEM) to discretize the governing equations. Hence, Eq. (7) should be changed into weak formulation.

We define a test function

ψ

⁠. Multiplying Eq. by the test function

ψ

and integrating over the reservoir domain

Ω_{r}

⁠, the weak formulation is acquired

\int_{Ω_{r}} ψ \frac{1}{M} \frac{\partial p_{r}}{\partial t} d V - \int_{Ω_{r}} ψ \nabla \cdot \frac{k_{r}}{μ} \nabla p_{r} d V = \int_{Ω_{r}} ψ q_{r} d V - \int_{Ω_{r}} ψ α \frac{\partial ε_{vol}}{\partial t} d V .

(8)

The second term of Eq. is reformulated through integration by parts.

\begin{matrix} \int_{Ω_{r}} ψ \nabla \cdot \frac{k_{r}}{μ} \nabla p_{r} d V = - \int_{Ω_{r}} \frac{k_{r}}{μ} \nabla p_{r} \cdot \nabla ψ d V + \int_{\partial_{N} Ω} ψ \frac{k_{r}}{μ} \nabla p_{r} \cdot n d S - \int_{Γ^{+}} ψ \frac{k_{r}}{μ} \nabla p_{r} \cdot n_{Γ^{+}} d S - \int_{Γ^{-}} ψ \frac{k_{r}}{μ} \nabla p_{r} \cdot n_{Γ^{-}} d S - \int_{\partial Ω_{c}} ψ \frac{k_{r}}{μ} \nabla p_{r} \cdot n_{c} d S \end{matrix} .

(9)

The second term on the right of the equation indicates the influence of the Neumann boundary. The third and fourth terms on the right of the equation indicate the fluid flux from fracture to reservoir. According to the mass conservation, these two terms can be considered as the fluid loss in fracture domain during hydraulic fracturing. Hence, in this paper, this fluid loss is assumed as the leak-off from fracture to reservoir. Define

q_{l}

as the leak-off flux.

q_{l} = - 〚 \frac{k_{r}}{μ} \nabla p_{r} 〛 \cdot n_{Γ} .

(10)

The notation $〚〛$ denotes the “discontinuity” of the vector $\nabla p$ at the fracture face.

The fifth term on the right of the equation indicates the fluid exchange between reservoir and cave. Define

q_{c r}

as the exchange flux between reservoir and cave.

q_{c r} = - \frac{k_{r}}{μ} \nabla p_{r} \cdot n_{c} .

(11)

Equation becomes

\int_{Ω_{r}} ψ \nabla \cdot \frac{k_{r}}{μ} \nabla p_{r} d V = - \int_{Ω_{r}} \frac{k_{r}}{μ} \nabla p_{r} \cdot \nabla ψ d V - \int_{\partial_{N} Ω} ψ q_{n} d S + \int_{Γ} ψ q_{l} d S + \int_{\partial Ω_{c}} ψ q_{c r} d S,

(12)

where

q_{n}

is the prescribed fluid flux on boundary,

q_{l}

is the leak-off flux from fracture to reservoir, and

q_{c r}

is the fluid exchange flux between reservoir and cave.

Upon substituting Eq. into Eq. , the governing equation becomes

\begin{matrix} \int_{Ω_{r}} ψ \frac{1}{M} \frac{\partial p_{r}}{\partial t} d V + \int_{Ω_{r}} \frac{k_{r}}{μ} \nabla p_{r} \cdot \nabla ψ d V = \int_{Ω_{r}} ψ q_{r} d V - \int_{Ω_{r}} ψ α \frac{\partial ε_{vol}}{\partial t} d V - \int_{\partial_{N} Ω} ψ q_{n} d S \\ + \int_{Γ} ψ q_{l} d S + \int_{\partial Ω_{c}} ψ q_{c r} d S \end{matrix} .

(13)

B. The fracture flow model

The fluid flow velocity in fracture is vastly greater than the velocity in the reservoir.⁴⁸ Due to the fracture width is much smaller than the fracture length, the fracture can be thought of as a plane object. Hence, we assume that the fluid flow in the fracture domain conforms to the Reynolds equation with cubic law. The cubic law^49,50 is as follows:

w V_{f} = - \frac{w^{3}}{12 μ} \nabla_{Γ} p_{f},

(14)

where

w

is the fracture width and

V_{f}

is the fluid flow velocity.

The Reynolds equation is as follows:

\frac{\partial w}{\partial t} + \nabla_{Γ} \cdot (w V_{f}) + q_{l} = q_{f} + q_{c f},

(15)

where

q_{f}

is the fluid source or sink term in the fracture.

Then, Eq. can be changed by the definition of surface divergence

\frac{\partial w}{\partial t} + (n_{Γ} \times \nabla) \cdot (n_{Γ} \times w V_{f}) + q_{l} = q_{f} + q_{c f} .

(16)

Multiplying Eq. by a test function

ψ

and integrating over the fracture domain

Γ

⁠, the weak formulation can be written as

\int_{Γ} ψ \frac{\partial w}{\partial t} d S + \int_{Γ} ψ (n_{Γ} \times \nabla) \cdot (n_{Γ} \times w V_{f}) d S + \int_{Γ} ψ q_{l} d S - \int_{Γ} ψ q_{c f} d S - \int_{Γ} ψ q_{f} d S = 0 .

(17)

The second term of Eq. is complicated, and discretizing it by FEM is very difficult. Hence, vector algebra identity is applied to simplify it.⁵¹

\int_{Γ} ψ (n_{Γ} \times \nabla) \cdot (n_{Γ} \times w V_{f}) d S = - \int_{Γ} (w V_{f} \cdot \nabla_{Γ} ψ) d S - \int_{\partial Γ_{c}} ψ w V_{f} \cdot n_{c} d s .

(18)

The second term on the right side of above formula only appears when fracture is connected with cave. Where $Γ_{c}$ is the boundary between fracture and cave, ds is the line integral.

Substituting Eq. into Eq. , then Eq. becomes

\int_{Γ} ψ (n_{Γ} \times \nabla) \cdot (n_{Γ} \times w V_{f}) d S = \int_{Γ} \frac{w^{3}}{12 μ} (\nabla_{Γ} p_{f} \cdot \nabla_{Γ} ψ) d S - \int_{\partial Γ_{c}} ψ w V_{f} \cdot n_{c} d s .

(19)

Finally, the governing equation in the fracture domain becomes

\int_{Γ} ψ \frac{\partial w}{\partial t} d S + \int_{Γ} \frac{w^{3}}{12 μ} (\nabla_{Γ} p_{f} \cdot \nabla_{Γ} ψ) d S + \int_{Γ} ψ q_{l} d S - \int_{Γ} ψ q_{f} d S - \int_{\partial Γ_{c}} ψ w V_{f} \cdot n_{c} d s = 0 .

(20)

C. The cave flow model

Caves form through a long-term geological movement, such as sedimentation, tectogenesis, and dissolution.^52,53 Scholars classify caves into two categories: one is unfilled cave, and the other is filled cave.⁴⁵ In this paper, we only consider the unfilled cave. Hence, there are no any loose substances inside the cave, and the fluid flow conforms to free flow. The Stokes equation is a very good appropriate choice to describe the fluid flow in unfilled cave. Many scholars combined Darcy–Stokes equations to establish flow model in porous domain and cave domain.^54–58 Their simulation results show the pressure spreads very quickly in cave, and there is almost no pressure drop inside cave. Hence, some scholars proposed the equipotential body model in cave domain. This model assumes the pressure inside the cave is equal everywhere. Many literatures have proven the correctness of this model.^59–65 In this paper, we choose the equipotential body model to describe the fluid flow in cave domain.

The flow governing equation in cave domain is

V c_{f} \frac{d p_{c}}{d t} = Q_{c},

(21)

where V is the cave volume,

c_{f}

is the fluid compressibility coefficient, and cf. = 1/K_f. Q_c is the total volume flow rate into or out of cave.

Q_c is related to the flow velocity at the boundary of reservoir cave and fracture cave.

Q_{c} = - \int_{\partial Ω_{c}} V_{r} \cdot n_{c} d S - \int_{\partial Γ_{c}} w V_{f} \cdot n_{c} d s,

(22)

where

\partial Γ_{c}

is the boundary between fracture and cave. w is the fracture width, and ds is the line integral.

In summary, the flow governing equations in three domains have been established. Adding up Eqs. and , the governing equation for non-free flow domain (fracture and reservoir domain) is obtained

\begin{matrix} \int_{Ω_{r}} (ψ \frac{1}{M} \frac{\partial p}{\partial t} + \frac{k_{r}}{μ} \nabla p \cdot \nabla ψ) d V + \int_{Γ} \frac{w^{3}}{12 μ} (\nabla_{Γ} p \cdot \nabla_{Γ} ψ) d S \\ = - \int_{Γ} ψ \frac{\partial w}{\partial t} d S - \int_{Ω_{r}} ψ α \frac{\partial ε_{vol}}{\partial t} d V - \int_{\partial_{N} Ω} ψ q_{n} d S + \int_{\partial Ω_{c}} ψ q_{c r} d S \\ + \int_{\partial Γ_{c}} ψ w V_{f} \cdot n_{c} d s + \int_{Ω_{r}} ψ q_{r} d V + \int_{Γ} ψ q_{f} d S \end{matrix} .

(23)

Due to the hydraulic communication, the fluid pressure is continuous. Therefore, a single pressure variable $p$ can represent $p_{r}$ and $p_{f}$ ⁠. Meanwhile, the leak-off flux from fracture to reservoir is eliminated. So, there is no leak-off term in Eq. (23).

In Sec. II D, we will establish the mechanical model.

D. Energy functional

The research on brittle material deformation and fracture can be traced back to Griffith.⁶⁶ According to Griffith's criterion, fractures begin to extend when the energy release rate

G

equates to the fracture toughness

G_{c}

⁠. Based on Griffith's criterion and the variational fracture approach,²⁷ we establish the total energy functional

F

of a single-phase solid. It is composed of the following parts: elastic strain energy, fracture surface energy, and work function of external force

P_{ext}

⁠.

F = \int_{Ω \ Γ} Ψ_{ε} d V + \int_{Γ} G_{c} d S - P_{ext},

(24)

where

Ψ_{ε}

is elastic strain energy density.

Ψ_{ε} = \frac{1}{2} D ε \cdot ε .

(25)

In the above formula, $D$ is the elastic matrix, and $ε$ is the strain.

However, when the poroelastic material is filled with fluid, the effect of pore pressure should be taken into account. Following the Biot poroelasticity theory,⁴⁷ the effective stress is introduced:

σ_{eff} = σ - α p I

⁠. Then, the strain energy density can be modified to

Ψ_{m ε} = \frac{1}{2} D (ε - \frac{α}{n κ} p I) \cdot (ε - \frac{α}{n κ} p I),

(26)

where

Ψ_{m ε}

is the modified strain energy density,

I

is the identity tensor,

κ

is the bulk modulus of material, and

n

represents the space dimension,

n = 2 or 3

⁠.

In the mechanical model, a specified traction force

τ

is imposed on the Neumann boundary. In addition, a body force

b^{*}

is imposed on the whole poroelastic domain. Hence,

P_{ext}

is expressed as follows:

P_{ext} = \int_{\partial_{N} Ω} τ \cdot u d S + \int_{Ω} b^{*} \cdot u d V .

(27)

During hydraulic fracturing, fracturing fluid is injected to increase the fracture volume. The energy that is contributed by fracturing fluid pressure must be accounted for the total energy functional. The work done by the fracturing fluid pressure is expressed as follows:

P_{inject} = \int_{Γ} p 〚 u 〛 \cdot n_{Γ} d S .

(28)

In summary, the total energy functional

F

becomes

\begin{matrix} F = \int_{Ω \ Γ} \frac{1}{2} D (ε - \frac{α}{n κ} p I) \cdot (ε - \frac{α}{n κ} p I) d V + \int_{Γ} G_{c} d S - \int_{\partial_{N} Ω} τ \cdot u d S \\ - \int_{Ω} b^{*} \cdot u d V + \int_{Γ} p 〚 u 〛 \cdot n_{Γ} d S \end{matrix} .

(29)

After establishing the total energy functional $F$ ⁠, we will use the phase-field method (PFM) to regularize it.

E. Regularization by phase-field method

In PFM, a scalar field

v

is introduced to diffuse the sharp fracture topology, which can eliminate the difficulties in the explicit tracing procedure of fracture surface. Specifically speaking, v = 1 represents the intact region and v = 0 represents the completely fractured region. When 0 < v < 1, it denotes the intermediate region linking the intact region and completely fractured region. Hence, displacement field is continuous in the whole region. Bourdin et al.²⁷ proposed an exponential function to approximate the phase-field scalar

v .

The function is as follows:

v = {\begin{matrix} 0 if d_{T} (x) \leq 0.5 d_{f} \\ 1 - e^{- \frac{d_{T} - 0.5 d_{f}}{2 ϵ}} otherwise, \end{matrix}

(30)

where

d_{f}

is the thickness of the initial fracture,

d_{T}

is the distance from any point

x

to the fracture, and

ϵ

is a regularization parameter, which decides the size of the intermediate region.

The interior of the cave is completely damaged. Yet, the edge between cave and reservoir is sharp. The PFM is still suitable for characterizing the cave. Following the above method, the function for phase-field scalar

v

in cave domain is proposed.

v = {\begin{matrix} 0 if d_{T} (x) \leq 0.5 d_{c} \\ 1 - e^{- \frac{d_{T} - 0.5 d_{c}}{2 ϵ}} otherwise, \end{matrix}

(31)

where

d_{c}

is the cave diameter, and

d_{T}

is the distance from any point

x

to the cave center.

The phase-field representation of fracture and cave is shown in Fig. 2. The sharp fracture and cave edge are diffused by phase field.

FIG. 2.

$The phase-field representation of fracture and cave.$

View large Download slide

The phase-field representation of fracture and cave.

Upon above settings, the fracture surface energy can be approximated in the following form:

\int_{Γ} G_{c} d S ≃ \frac{1}{2} G_{c} \int_{Ω} (\frac{{(1 - v)}^{2}}{ϵ} + ϵ {| \nabla v |}^{2}) d V .

(32)

Then, PFM is utilized to regularize the modified strain energy density

Ψ_{m ε}

⁠. Based on the modeling hypothesis that the damage occurs on the subpore scale,^24,34 we select to let the phase-field variable to affect the stress term and not the fluid pressure term.

Ψ_{m ε} = \frac{1}{2} D (v ε - \frac{α}{n κ} p I) \cdot (v ε - \frac{α}{n κ} p I) .

(33)

Through tensor operation,

Ψ_{m ε}

can be further converted into

Ψ_{m ε} = \frac{1}{2} v^{2} σ_{i j} ε_{i j} - vap \cdot ε_{i i} = v^{2} Ψ_{ε} - v α p \cdot (\nabla \cdot u) .

(34)

Next, the work done by the fracturing fluid pressure needs to be regularized.

Define a discontinuous vector

Φ

that is admitted traces

Φ^{+}

and

Φ^{-}

on each side of

Γ

⁠.

〚 Φ 〛 \cdot n_{Γ}

represents the difference value between the two sides of fracture surface. The integral of the difference value in fracture domain can be transformed to the integral in the whole poroelastic domain by the following approximation:^67,68

\int_{Γ} 〚 Φ 〛 \cdot n_{Γ} d S ≃ lim_{ϵ \to 0} \int_{Ω} Φ \cdot n_{Γ} | \nabla v | d V .

(35)

\nabla v

is a vector that points toward the orientation where

v

changes fastest. The region with the fastest change of

v

is the transition region. So, the fracture normal

n_{Γ}

can be approximated as

n_{Γ} \approx \frac{\nabla v}{| \nabla v |} .

(36)

Substituting Eq. into Eq. , the integral transformation becomes

\int_{Γ} 〚 Φ 〛 \cdot n_{Γ} d S ≃ lim_{ϵ \to 0} \int_{Ω} Φ \cdot \nabla v d V .

(37)

The displacement vector

u

is also discontinuous on each side of

Γ

⁠. So let

Φ = p u

⁠, the work done by the fracturing fluid pressure is regularized

\int_{Γ} p 〚 u 〛 \cdot n_{Γ} d S = lim_{ϵ \to 0} \int_{Ω} p u \cdot \nabla v d V .

(38)

After the above phase-field regularization, the total energy functional

F

can be expressed as

\begin{matrix} F = \int_{Ω} v^{2} Ψ_{ε} d V - \int_{Ω} v α p \cdot (\nabla \cdot u) d V + \frac{1}{2} G_{c} \int_{Ω} (\frac{{(1 - v)}^{2}}{ϵ} + ϵ {| \nabla v |}^{2}) d V \\ - \int_{\partial_{N} Ω} τ \cdot u d S - \int_{Ω} b^{*} \cdot u d V + \int_{Ω} p u \cdot \nabla v d V \end{matrix} .

(39)

F. The decomposition of elastic strain energy

As can be seen in Eq. ,

F

is symmetric in compression and tension. However, the crack evolution under compression is not allowed.^9,69 Therefore, the elastic strain energy needs to be decomposed. In this paper, we follow the decomposition method proposed by Amor et al.⁷⁰ In this method, the strain is decomposed into spherical part

ε_{S}

and deviatoric part

ε_{D}

⁠,

ε = ε_{S} + ε_{D},

(40)

{\begin{matrix} ε_{S} = \frac{1}{n} t r (ε) I \\ ε_{D} = ε - \frac{1}{n} t r (ε) I, \end{matrix}

(41)

where

n

is the space dimension, and

I

is the identity tensor.

With this decomposition, the strain energy density becomes

Ψ_{ε} = \frac{κ}{2} t r {(ε)}^{2} + μ ε_{D} \cdot ε_{D},

(42)

where

κ

is the bulk modulus of material and

κ = λ + 2 μ / n

⁠,

λ

⁠, and

μ

are Lamé constants.

The crack is expected to open in the tensile regions and prohibited from interpenetration in compressed regions. Therefore, Amor et al.⁷⁰ decomposed the trace of strain tensor into positive and negative part: $t r^{+} (ε) = max (t r (ε), 0)$ ⁠, $t r^{-} (ε) = max (- t r (ε), 0)$ ⁠.

On the basis, the strain energy density is distinguished into tension, compression, and shear part.

Ψ_{ε} = \frac{κ}{2} t r^{+} {(ε)}^{2} + \frac{κ}{2} t r^{-} {(ε)}^{2} + μ ε_{D} \cdot ε_{D} .

(43)

Only the strain energy related to tension and shear part drives the crack evolution. Therefore, the phase-field variable only affects the two parts. As a result, the total energy functional

F

is written as

\begin{matrix} F = \int_{Ω} [(v^{2} + η) (\frac{κ}{2} t r^{+} {(ε)}^{2} + μ ε_{D} \cdot ε_{D}) + \frac{κ}{2} t r^{-} {(ε)}^{2}] d V - \int_{Ω} v α p \cdot (\nabla \cdot u) d V + \frac{1}{2} G_{c} \int_{Ω} (\frac{{(1 - v)}^{2}}{ϵ} + ϵ {| \nabla v |}^{2}) d V - \int_{\partial_{N} Ω} τ \cdot u d S - \int_{Ω} b^{*} \cdot u d V + \int_{Ω} p u \cdot \nabla v d V \end{matrix},

(44)

where

0 < η ≪ 1

⁠, and it is an additional parameter to ensure numerical stability when v tends to 0.

G. Minimization of energy functional

Crack generation and propagation occur when the total energy functional

F

reaches minimum value.^25,71 In this paper, we employ the variational approach to minimize

F

⁠. Specifically, the first-order variation of

F

is set to zero.

δ F = δ_{u} F + δ_{v} F = 0 .

(45)

In the above formula

\begin{matrix} δ_{u} F = \int_{Ω} [(v^{2} + η) (κ t r^{+} (ε) δ_{i j} + 2 μ ε_{D}) + κ t r^{-} (ε) δ_{i j}] {(δ u_{i})}_{, j} d V \\ - \int_{Ω} v α p δ_{i j} {(δ u_{i})}_{, j} d V - \int_{\partial_{N} Ω} τ_{i} δ u_{i} d S - \int_{Ω} b_{i}^{*} δ u_{i} d V + \int_{Ω} p δ u_{i} v_{, i} d V \end{matrix},

(46)

\begin{matrix} δ_{v} F = \int_{Ω} 2 (\frac{κ}{2} t r^{+} {(ε)}^{2} + μ ε_{D} \cdot ε_{D}) v δ vdV - \int_{Ω} α p \cdot (\nabla \cdot u) δ vdV \\ + G_{c} \int_{Ω} (\frac{v δ v}{ϵ} - \frac{δ v}{ϵ} + ϵ v_{, i} {(δ v)}_{, i}) d V + \int_{Ω} p u_{i} {(δ v)}_{, i} d V \end{matrix},

(47)

where the operator

*_{, i}

represents the derivative to Cartesian space, such as

{(δ v)}_{, i} = \partial (δ v) / \partial x_{i}

⁠.

For any variable displacement

δ u_{i}

⁠, it meets Eq. . In addition, for any variable phase-field

δ v

⁠, it also meets Eq. . Therefore, two governing equations are obtained.

\begin{matrix} \int_{Ω} [(v^{2} + η) (κ t r^{+} (ε) δ_{i j} + 2 μ ε_{D}) + κ t r^{-} (ε) δ_{i j}] {(δ u_{i})}_{, j} d V \\ - \int_{Ω} v α p δ_{i j} {(δ u_{i})}_{, j} d V - \int_{\partial_{N} Ω} τ_{i} δ u_{i} d S - \int_{Ω} b_{i}^{*} δ u_{i} d V + \int_{Ω} p δ u_{i} v_{, i} d V = 0 \end{matrix} .

(48)

\begin{matrix} \int_{Ω} 2 (\frac{κ}{2} t r^{+} {(ε)}^{2} + μ ε_{D} \cdot ε_{D}) v δ vdV - \int_{Ω} α p \cdot (\nabla \cdot u) δ vdV \\ + G_{c} \int_{Ω} (\frac{v δ v}{ϵ} - \frac{δ v}{ϵ} + ϵ v_{, i} {(δ v)}_{, i}) d V + \int_{Ω} p u_{i} {(δ v)}_{, i} d V = 0 \end{matrix} .

(49)

III. NUMERICAL IMPLEMENTATION

This section will delineate the numerical implement of the proposed multi-field coupling model. First, the FEM is used to discretize the spatial domain and the finite difference method is used to discretize the time domain. After that, a staggered scheme is designed to solve the discretized equations. The details are as follows.

A. Governing equation discretization

A structured eight-node hexahedral element is used to discretize the spatial domain in our work. Hence, the numerical algorithm can simulate the fracture propagation in three-dimensional space. If we assume the plane strain conditions, the numerical algorithm is also suitable for simulating the fracture propagation in two-dimensional space. The code of this numerical algorithm is completed by C language.

Using the eight-node element, the displacement

u

⁠, phase-field

v

⁠, and fluid pressure

p

are converted to

u = N_{u} \hat{d}, v = N_{v} \hat{v}, p = N_{p} \hat{p},

(50)

where

N_{u}, N_{v}, and N_{p}

are the interpolation function matrices for

u, v

⁠, and

p .

\hat{d}, \hat{v}

⁠, and

\hat{p}

are the vectors composed of node values

u_{i}, v_{i}

⁠, and

p_{i}

⁠. Their definitions are as follows:

N_{u} = [\begin{matrix} N_{1} & 0 & 0 & \dots & N_{n} & 0 & 0 \\ 0 & N_{1} & 0 & \dots & 0 & N_{n} & 0 \\ 0 & 0 & N_{1} & \dots & 0 & 0 & N_{n} \end{matrix}],

N_{v} = N_{p} = [\begin{matrix} N_{1} & N_{2} & \dots & N_{n} \end{matrix}],

\hat{d} = {[\begin{matrix} u_{1} & u_{2} & \dots & u_{n} \end{matrix}]}^{T},

u_{i} = {[\begin{matrix} u_{i x} & u_{i y} & u_{i z} \end{matrix}]}^{T},

\hat{v} = {[\begin{matrix} v_{1} & v_{2} & \dots & v_{n} \end{matrix}]}^{T},

\hat{p} = {[\begin{matrix} p_{1} & p_{2} & \dots & p_{n} \end{matrix}]}^{T},

where

n

represents the count of node, n = 8, and

N_{i}

represents the interpolation function of node i.

In this paper, we use the Voigt notations. The strain, gradient of phase field, and gradient of pressure are converted to

ε = B_{u} \hat{d}, \nabla v = B_{v} \hat{v}, \nabla p = B_{p} \hat{p},

(51)

where

B_{u}

is the strain matrix, and

B_{v}

and

B_{p}

are the gradients matrices. They are defined as

B_{u} = [\begin{matrix} N_{1, x} & 0 & 0 & \dots & N_{n, x} & 0 & 0 \\ 0 & N_{1, y} & 0 & \dots & 0 & N_{n, y} & 0 \\ 0 & 0 & N_{1, z} & \dots & 0 & 0 & N_{n, z} \\ N_{1, y} & N_{1, x} & 0 & \dots & N_{n, y} & N_{n, x} & 0 \\ 0 & N_{1, z} & N_{1, y} & \dots & 0 & N_{n, z} & N_{n, y} \\ N_{1, z} & 0 & N_{1, x} & \dots & N_{n, z} & 0 & N_{n, x} \end{matrix}],

B_{v} = B_{p} = [\begin{matrix} N_{1, x} & N_{2, x} & \dots & N_{n, x} \\ N_{1, y} & N_{2, y} & \dots & N_{n, y} \\ N_{1, z} & N_{2, z} & \dots & N_{n, z} \end{matrix}] .

Similarly, the variable displacement

δ u

and variable phase-field

δ v

are converted to

δ u = N_{u} δ \hat{d}, δ v = N_{v} δ \hat{v} .

(52)

By substituting Eqs. into Eqs. and , the discrete form of

δ_{u} F

and

δ_{v} F

is obtained.

R^{u} = K_{u} \hat{d} - f_{u} = 0,

(53)

R^{v} = K_{v} \hat{v} - f_{v} = 0,

(54)

where

R^{u}

and

R^{v}

are the residuals,

K_{u}

and

K_{v}

are the stiffness matrices, and

f_{u}

and

f_{v}

are the coupling vectors. Their forms are as follows:

K_{u} = \int_{Ω} B_{u}^{T} [(v^{2} + η) (κ D_{t r} B_{u} + 2 μ B_{u} - \frac{2 μ}{3} D_{t r} B_{u}) + κ D_{t r} B_{u}] d V,

f_{u} = \int_{Ω} v α p B_{u}^{T} I d V + \int_{\partial_{N} Ω} N_{u}^{T} τ d S + \int_{Ω} N_{u}^{T} b^{*} d V - \int_{Ω} p N_{u}^{T} \cdot \nabla vdV,

K_{v} = \int_{Ω} [2 N_{v}^{T} (\frac{κ}{2} t r^{+} {(ε)}^{2} + μ ε_{D} \cdot ε_{D}) N_{v} + \frac{G_{c}}{ϵ} N_{v}^{T} N_{v} + G_{c} ϵ B_{v}^{T} B_{v}] d V,

f_{v} = \int_{Ω} N_{v}^{T} α p \cdot (\nabla \cdot u) d V + \int_{Ω} \frac{G_{c}}{ϵ} N_{v}^{T} d V - \int_{Ω} B_{v}^{T} p u d V,

where

D_{t r}

is the matrix to realize the trace of strain tensor.

D_{t r} = [\begin{matrix} 1 & 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}],

I = {[\begin{matrix} 1 & 1 & 1 & 0 & 0 & 0 \end{matrix}]}^{T} .

Next, we will discretize the fluid flow model. First, the source or sink term needs to be regularized. During hydraulic fracturing, the injection point is in the fracture domain. To simulate the hydraulic fracturing more realistically, the source or sink term

q_{f}

is generally composed of the actual injection rate and the Dirac delta function. In PFM, the fracture is represented by a smooth scalar field. Therefore,

q_{f}

also needs to be regularized, which is shown as follows:

q_{f} = {\begin{matrix} \sum_{i = 1}^{n} Q_{f, i} \frac{e^{- | x - x_{i} |}}{2 π ϵ^{2}} two dimension \\ \sum_{i = 1}^{n} Q_{f, i} \frac{e^{- | x - x_{i} |}}{4 π ϵ^{3}} three dimension, \end{matrix}

(55)

where

x_{i}

is the location of injection point, and

Q_{f, i}

is the actual injection rate, with the unit of

m^{2} / s

in two dimension and

m^{3} / s

in three dimension.

There are two flow equations in our model. We use a staggered scheme to solve them. First, we solve Eq. to get the pressure p. Then, we calculate the fluid velocity at the boundary grids of cave. After that, we solve Eq. to get the cave pressure p_c. Then, we update the fluid velocity at the boundary grids of cave. Subsequently, we proceed to next iteration step and solve Eq. again. When

p^{l + 1} \to p^{l}

and

p_{c}^{l + 1} \to p_{c}^{l}

⁠, the iteration ends.

\begin{matrix} \int_{Ω_{r}} (ψ \frac{1}{M} \frac{\partial p^{l + 1}}{\partial t} + \frac{k_{r}}{μ} \nabla p^{l + 1} \cdot \nabla ψ) d V + \int_{Γ} \frac{w^{3}}{12 μ} (\nabla_{Γ} p^{l + 1} \cdot \nabla_{Γ} ψ) d S \\ = - \int_{Γ} ψ \frac{\partial w}{\partial t} d S - \int_{Ω_{r}} ψ α \frac{\partial ε_{vol}}{\partial t} d V - \int_{\partial_{N} Ω} ψ q_{n} d S + \int_{\partial Ω_{c}} ψ q_{c r}^{l} d S \\ + \int_{\partial Γ_{c}} ψ w V_{f}^{l} \cdot n_{c} d s + \int_{Ω_{r}} ψ q_{r} d V + \int_{Γ} ψ q_{f} d S \end{matrix},

(56)

V c_{f} \frac{d p_{c}^{l + 1}}{d t} = - \int_{\partial Ω_{c}} q_{c r}^{l + 1} d S - \int_{\partial Γ_{c}} w V_{f}^{l + 1} \cdot n_{c} d s .

(57)

Equation (56) involves the fracture width $w$ and the volumetric strain $ε_{vol} .$ The fracture width is the difference of displacement on both sides of the fracture face: $w = 〚 u 〛 \cdot n_{Γ}$ ⁠. Since the fracture width and the volumetric strain are the function of displacement, Eq. (56) contains the pressure field and the displacement field.

The fixed stress split method⁷² is employed to make, Eq. (56), decoupling from deformation. It is an unconditional stable algorithm, which is appropriate for the multi-field coupling problem. Specifically, the fracture width and the volumetric stress are assessed using the value of previous iteration step when solves Eq. (56).

The relationship between volumetric strain and volumetric stress is

σ_{vol}^{j + 1} = σ_{vol}^{j} + κ (ε_{vol}^{j + 1} - ε_{vol}^{j}) - α (p^{j + 1} - p^{j}) .

(58)

Then, Eq. becomes

\begin{matrix} \int_{Ω_{r}} ((\frac{ψ}{M} + ψ \frac{α^{2}}{κ}) \frac{\partial p^{l + 1, j + 1}}{\partial t} + \frac{k_{r}}{μ} \nabla p^{l + 1, j + 1} \cdot \nabla ψ) d V + \int_{Γ} \frac{{(w^{j})}^{3}}{12 μ} (\nabla_{Γ} p^{l + 1, j + 1} \cdot \nabla_{Γ} ψ) d S \\ = - \int_{Γ} ψ \frac{\partial w^{j}}{\partial t} d S - \int_{Ω_{r}} ψ α \frac{\partial ε_{vol}^{j}}{\partial t} d V + \int_{Ω_{r}} ψ \frac{α^{2}}{κ} \frac{\partial p^{j}}{\partial t} d V - \int_{\partial_{N} Ω} ψ q_{n} d S \\ + \int_{Ω_{r}} ψ q_{r} d V + \int_{Γ} ψ q_{f} d S + \int_{\partial Ω_{c}} ψ q_{c r}^{l} d S + \int_{\partial Γ_{c}} ψ w^{j} V_{f}^{l} \cdot n_{c} d s, \end{matrix}

(59)

where

j

represents the iteration step of fixed stress split method.

Section III B provides a detailed explanation of the implementation scheme.

Because the location of the fracture is variable during hydraulic fracturing, it is difficult to directly conduct numerical integration on the fracture domain. Hence, a regularized method is used to transform the integral in fracture domain $Γ$ ⁠. Correspondingly, the integral in reservoir domain $Ω_{r}$ is also transformed.

Finally, the flow governing equation in reservoir and fracture domain can be discretized as

R^{p} = [M_{p} + Δ t S_{p}] {[\hat{p}]}_{n + 1}^{l + 1, j + 1} - [M_{p}] {[\hat{p}]}_{n} - Δ t {[f_{p}]}_{n + 1} = 0,

(60)

where

R^{p}

is the residuals,

M_{p}

is the compressibility matrix,

S_{p}

is the transportation matrix, and

f_{p}

is the coupling vector. Their forms are as follows:

M_{p} = \int_{Ω_{d}} (\frac{N_{p}^{T} N_{p}}{M} + v^{2} N_{p}^{T} N_{p} \frac{α^{2}}{κ}) d V,

S_{p} = \int_{Ω_{d}} [\frac{k_{r}}{μ} B_{p}^{T} B_{p} + \frac{{(w^{j})}^{3}}{12 μ} | \nabla v | B_{Γ p}^{T} B_{Γ p}] d V,

\begin{matrix} f_{p} = \int_{Ω_{d}} N_{p}^{T} (- \frac{\partial w^{j}}{\partial t} | \nabla v | - α \frac{\partial ε_{vol}^{j}}{\partial t} v^{2} + \frac{α^{2}}{κ} \frac{\partial p^{j}}{\partial t} v^{2} + \frac{(1 - v^{2})}{M} \frac{\partial p^{j}}{\partial t} + q_{r} v^{2} + q_{f} | \nabla v |) d V \\ - \int_{\partial_{N} Ω_{d}} N_{p}^{T} q_{n} d S + \int_{\partial Ω_{c}} N_{p}^{T} q_{c r}^{l} d S + \int_{\partial Ω_{c f}} N_{p}^{T} q_{c f}^{l, j} d S \end{matrix} .

The flow governing equation in cave domain can be discretized as

p_{n + 1}^{l + 1} = p_{n} - \frac{Δ t}{V c_{f}} (\int_{\partial Ω_{c}} q_{c r}^{l + 1} d S + \int_{\partial Ω_{c f}} q_{c f}^{l + 1, j} d S) .

(61)

Please see Appendix for the discretizing process.

B. Numerical implementation scheme

In Sec. III A, we have obtained three discrete governing equations: Eqs. (53), (54), and (60) The process of solving the three equations is called module_u, module_v, and module_p, respectively. The three modules are interrelated closely. Module_u solves the displacement field $u$ ⁠. Based on the solution of the displacement field, relevant parameters such as the stress, strain, and fracture width can be calculated. Module_u supplies the fracture width and volumetric strain to module_p. In addition, it supplies strain energy to module_v. Module_p is designed to solve the pressure field $p$ ⁠. It supplies fluid pressure to module_u and module_v. Module_v is designed to solve the phase field $v$ ⁠. It provides indicator function for module_p and phase-field variable for module_u. The Newton–Raphson method is applied in the three modules to solve the respective equations.

There are two ways to solve the proposed multi-field coupling model, which are monolithic approach and staggered scheme. The monolithic approach is capable of simulating the real physical process. Yet, it needs more precise mathematical techniques to obtain the convergence, because of the non-convexity of energy functional.^73,74 The staggered scheme is robust, but it needs more costs of computing.^75,76 In this paper, we choose the staggered scheme to solve the multi-field coupling model.

The detail of the staggered scheme is shown in Fig. 3. First, we construct the initial $v$ field for the cave domain and the perforated fracture domain. Next, the boundary conditions and initial conditions for the multi-field coupling model are established. After that, the solving cycle of the staggered scheme begins. The staggered scheme involves three cycles: external cycle, middle cycle, and inner cycle. In the middle and inner cycles, the time is separated and the physical field iterates at each time step. Further explained, the displacement field $u$ and pressure field $p$ are solved by keeping the phase-field $v$ fixed in the inner cycle. Then, the phase-field $v$ is solved by keeping $u$ and $p$ fixed in the middle cycle.

FIG. 3.

View large Download slide

The staggered scheme for the multi-field coupling model.

In the inner cycle, the fixed stress split method is applied to solve $u$ and $p$ ⁠. First, module_p is executed to obtain the fluid pressure $p$ through retaining the volumetric stress fixed. Then, module_u is executed to obtain the displacement $u$ by the updated $p$ from module_p. After that, it iterates between module_p and module_u until the error $ε_{inn}$ is less than the predetermined tolerance $ε_{t}$ ⁠. In the middle cycle, the phase-field $v$ is solved using the convergent solutions of $u$ and $p$ from the inner cycle. The irreversible condition is set here to ensure that the fracture does not recover to the unbroken state. If the error $ε_{mid}$ is less than the tolerance $ε_{t}$ ⁠, the computation in new time step begins. In addition, the final convergent solutions of $u, p$ ⁠, and $v$ in current time step are used as the initial values in next time step.

IV. MODEL VALIDATION

This section offers three numerical cases to verify the proposed model. The proposed model couples three systems including fracture, cave, and reservoir. However, as far as we know, there is no analytical solution of the hydraulic fracturing model when natural cave exists. Hence, the basic validation idea is as follows: first, we compare our model with the Khristianovic–Geertsma–de Klerk (KGD) model, which is a classic case for verifying the hydraulic fracturing. Through this comparison, the correctness can be verified when cave is not considered in our model. Next, in order to verify the correctness when cave is considered, we compare our model with the model of Cheng et al.³⁹ through two numerical cases.

A. Comparison with KGD model

The KGD model was first established in the study of Khristianovic, Geertsma, and de Klerk.^77,78 In the KGD model, a single fracture propagates under the constant fluid driven. It is a plane-strain issue and the fracture height is invariant. The reservoir is assumed impermeable, so the fluid leakage from fracture is not existent. In addition, the reservoir is linear elastic, isotropic, and homogeneous. There are many analytical solutions and semi-analytical solutions in the KGD model. In this paper, we choose the analytical solutions of Spence and Sharp⁷⁹ to verify our model. The characteristic of Spence's model is injecting viscous fluid into an impermeable reservoir. Hence, their model is in viscosity-controlled regime.

The geometry and boundary conditions for simulating KGD issue are shown in Fig. 4. In the KGD model, a preexisting fracture propagates in an infinite domain. In our modeling, the computation domain is set as $200 \times 200 m^{2},$ and its center is placed an initial fracture, with the half-length equaling 2.5 m. The four boundaries are restricted for the displacement field. The above setting is to simulate the infinite domain approximately. An injection point is set at the center of fracture, with the rate $Q = 3 \times 10^{- 4} m^{2} / s .$ The total injection time is 12 s. To simulate the impermeable reservoir, the reservoir permeability is set very small with the value of $2 \times 10^{- 22} m^{2} .$ The initial pore pressure is $1 \times 10^{- 6} MPa .$ Other important parameters are as follows: critical energy release rate $G_{c} = 100 N / m,$ fluid viscosity $μ = 1 \times 10^{- 3} Pa \cdot s,$ Young's modulus $E = 15$ GPa, Poisson's ratio $v = 0.2,$ and time interval $Δ t = 0.05 s .$

FIG. 4.

View large Download slide

Geometry and boundary conditions for simulating KGD issue.

We focus on the evolution of fluid pressure, fracture half-length, and fracture width vs time. Figure 5 shows the comparison results. We can find the fluid pressure at injection point rises nonlinearly prior to the fracture initiation in Fig. 5(a). Then, it gradually decreases with the fracture propagation. The analytical solution and numerical solution match well in terms of fluid pressure after the fracture initiation. The numerical solutions of fracture half-length and fracture width are also matched with the analytical solutions, which are shown in Figs. 5(b) and 5(c). Although these two matches are not as good as the fluid pressure match, these results are acceptable.

FIG. 5.

$The comparison results between analytical solution and numerical solution: (a) fluid pressure at injection point vs time; (b) fracture half-length vs time; and (c) fracture width at injection point vs time.$

View large Download slide

The comparison results between analytical solution and numerical solution: (a) fluid pressure at injection point vs time; (b) fracture half-length vs time; and (c) fracture width at injection point vs time.

There are some reasons that cause the matching difference. First, the phase-field method makes several assumptions for regularizing the hydraulic fracturing model. Second, the computational domain cannot be infinite. When the fracture length becomes larger, the fixed boundary will influence the numerical result.

To sum up, the comparison results between our model and the KGD model are favorable. They verify the correctness of our model when cave is not considered.

B. Comparison with Cheng's model

In this verification example, the natural cave is considered in the reservoir. We compare the proposed model with the model of Cheng et al..³⁹ They established a numerical model to simulate the fracture propagation when natural cave exists by XFEM. In their study, four numerical cases were set to study the impact of cave location. Here, we choose to compare with the first and second numerical cases to verify our model.

Figure 6(a) shows the geometric structure of the first numerical case in Cheng's model. The reservoir size is $50 \times 50 m^{2} .$ The right and bottom boundaries of reservoir are fixed. The maximum and minimum horizontal principal stresses are applied to the left boundary and top boundary, respectively. The reservoir center places a natural cave, with the radius of 3 m. In their model, the cave is empty and there is no fluid in the cave. An initial fracture is set at the left boundary, whose length is 1.8 m. The main input parameters in our model and Cheng's model are shown in Table I.

FIG. 6.

$The first numerical case. (a) The geometric structure of Cheng' model. (b) The fracture propagation path in our model and Cheng's model.$

View large Download slide

The first numerical case. (a) The geometric structure of Cheng' model. (b) The fracture propagation path in our model and Cheng's model.

TABLE I.

Main parameters in our model and Cheng's model.

Parameters	Our model	Cheng' s model
Young's modulus E	35 GPa	35 GPa
Poisson's ratio v	0.23	0.23
Reservoir permeability k_r	2 mD	2 mD
Fluid injection rate Q	0.05 m²/s	0.05 m²/s
Fluid viscosity μ	10 $mPa \cdot s$	10 $mPa \cdot s$
Horizontal stress difference $σ_{H} - σ_{h}$	5 MPa	5 MPa
Biot coefficient $α$	1	/
Porosity ϕ	0.04	0.04
Critical energy release rate G_c	200 N/m	/

Parameters	Our model	Cheng' s model
Young's modulus E	35 GPa	35 GPa
Poisson's ratio v	0.23	0.23
Reservoir permeability k_r	2 mD	2 mD
Fluid injection rate Q	0.05 m²/s	0.05 m²/s
Fluid viscosity μ	10 $mPa \cdot s$	10 $mPa \cdot s$
Horizontal stress difference $σ_{H} - σ_{h}$	5 MPa	5 MPa
Biot coefficient $α$	1	/
Porosity ϕ	0.04	0.04
Critical energy release rate G_c	200 N/m	/

It is worth noting that Cheng et al. did not establish the fluid flow model in cave domain. Hence, there is no pressure inside cave, while in our model, the injection fluid can enter the cave and the cave pressure will change vs time. To better match their model, our idea is keeping the cave pressure fixed. Under this setting, the cave pressure remains constant and is always equal to the initial pressure. The similarity is that there is no pressure change inside cave domain of the two models. In this way, the influence of cave pressure change on solid deformation can be neglected.

Figure 6(b) presents the numerical results of our model and Cheng's model. The fracture propagation path in Cheng's model is depicted by the white dotted line, while the fracture propagation path in our model is depicted by the phase-field evolution. We can see that the fracture propagates horizontally to the right side and finally is connected with the cave. The result of the two models matches well.

The second numerical case is presented in Fig. 7. The only variation is that the cave location is changed. The cave moves up, and the cave angle, which can be seen in Fig. 7(a), is $10 °$ ⁠. The input parameters remain unchanged. Figure 7(b) presents the numerical results. The fracture propagation path in Cheng's model is depicted by the white dotted line, while the fracture propagation path in our model is depicted by the phase-field evolution. We can find that the fracture first propagates to the right side and then deflects up. Finally, the fracture is connected with cave. In a word, the matching result of the two models is good. Yet, there is a little difference of the propagation path between the two models. This is because our model is built upon PFM, while their model is upon XFEM. There may be a little calculation difference between the two methods.

FIG. 7.

$The second numerical case. (a) The geometric structure of Cheng' model. (b) The fracture propagation path in our model and Cheng's model.$

View large Download slide

The second numerical case. (a) The geometric structure of Cheng' model. (b) The fracture propagation path in our model and Cheng's model.

To sum up, the correctness of our model is fully verified through the comparison with KGD model and Cheng's model.

C. Sensitivity analysis of the regularization parameter

The regularization parameter $ϵ$ can influence the numerical simulation. Making a reasonable choice for it is crucial in our model. Therefore, we will discuss the reasonable selection of $ϵ$ here.

In our model, the total energy functional [Eq. (29)] is regularized to Eq. (39) by introducing the diffusive interface to replace the sharp interface. Equation (39) converges to Eq. (29) as the regularization parameter $ϵ \to 0$ ⁠. Hence, the error between regularized model and actual model decreases as $ϵ$ decreases. The ideal thought is to choose the regularization parameter $ϵ$ as small as possible. However, previous studies have shown that $ϵ$ is related to the grid element size (h_w) in the FEM.^27,28 When $ϵ \leq h_{w}$ ⁠, the phase-field method will construct a sharp fracture profile. In this situation, the fracture surface energy is overrated. Hence, we need to make $ϵ > h_{w}$ in the numerical simulation. However, $ϵ$ cannot be chosen too large because large $ϵ$ will produce softening effects, which can result in the undervaluation of the strain energy. In the meanwhile, fracture and cave are overly smeared to the intact region, which can interfere with the fracture–cave interaction.

To explore the reasonable choice of

ϵ

in our model, we conduct a sensitivity analysis of

ϵ

⁠. The geometric structure is shown in Fig. 8(a). The computation domain is 4

\times

4 m, and an initial fracture is set at its center. The coordinates of point 0 are (0, 0). The position of cave's center is (3, 2), while its radius is 0.18 m. The critical energy release rate is G_c = 200 N/m, Young's modulus is E = 17 GPa, Poisson's ratio is 0.2, and the injection rate is 5

\times

10⁻⁵ m²/s. A structured mesh is used in our model, and the grid element size h_w is 0.005 m. We set four different values of

ϵ

⁠, which are

ϵ

= 1.5h_w,

ϵ

= 2h_w,

ϵ

= 3h_w, and

ϵ

= 4h_w. In the four cases, the fractures all propagate horizontally to the right side and finally are connected with the cave. Our focus is to compare the crack opening displacement (COD) of the four cases. The COD along the fracture length when t = 2.7 s is computed by the following formula:

u_{y} (x, 2) = \frac{1}{2} \int_{- \infty}^{+ \infty} u (x, 2) \cdot \nabla v d y .

(62)

FIG. 8.

View large Download slide

The sensitivity analysis of $ϵ$ ⁠. (a) The geometric structure. (b) The crack opening displacement with different $ϵ$ ⁠.

The simulation result is shown in Fig. 8(b). We can find when $ϵ$ = 1.5, 2, and 3h_w, the COD under the three cases is approximately the same. However, when $ϵ$ = 4h_w, the computed COD is obviously different from others. We can also find there is a sudden increase in COD near the fracture–cave connection position. It is because the interaction is affected by the excessive diffusive interfacial thickness. From the simulation result, it is concluded that the reasonable choice of $ϵ$ is h_w < $ϵ$ < 4h_w. In other simulation examples of this paper, we make $ϵ$ = 2h_w.

V. NUMERICAL CASES AND DISCUSSION

The existence of natural cave can markedly impact the HF propagation. The interaction mechanism between HF and natural cave is very complex. In this section, we set the following numerical cases to investigate the interaction mechanism in depth.

A. Effect of the distance between fracture and cave

To investigate the effect of the distance between fracture and cave, four cases are conducted here. The geometric structure of these cases is shown in Fig. 9. The computation domain is 20 $\times$ 20 m in size, and an initial fracture is set at its center. The coordinates of point 0 are (0, 0). The injection point is in the fracture center, and its coordinates are (10, 10). The position of cave's center is (14, H_c), while its radius is 1.3 m. To simulate the different distances between fracture and cave, H_c is set as 10, 12, 13.2, and 15.5. It is obvious that when H_c increases, the distance between fracture and cave increases. The initial reservoir pressure is 0.1 MPa. The maximum and minimum principal stresses are applied to the four boundaries of the computation domain. Other parameters are shown in Table II.

FIG. 9.

View large Download slide

The geometric structure for numerical cases.

TABLE II.

Main parameters of the cases.

Parameters	Value	Unit
Young's modulus E	30	GPa
Poisson's ratio v	0.2	/
Initial fracture length L_f	0.6	m
Cave radius r_c	1.3	m
Biot coefficient $α$	1	/
Fluid injection rate Q	1 $\times$ 10^–4	m²/s
Maximum principal stress $σ_{H}$	0.5	MPa
Minimum principal stress $σ_{h}$	0.4	MPa
Porosity ϕ	0.05	/
Fluid viscosity μ	3 $\times$ 10^–4	Pa·s
Initial reservoir pressure p_r	0.10	MPa
Rock bulk modulus K_s	10	GPa
Fluid bulk modulus K_f	0.6	GPa
Reservoir permeability k_r	$1 \times 10^{- 16}$	m²
Critical energy release rate G_c	250	N/m

Parameters	Value	Unit
Young's modulus E	30	GPa
Poisson's ratio v	0.2	/
Initial fracture length L_f	0.6	m
Cave radius r_c	1.3	m
Biot coefficient $α$	1	/
Fluid injection rate Q	1 $\times$ 10^–4	m²/s
Maximum principal stress $σ_{H}$	0.5	MPa
Minimum principal stress $σ_{h}$	0.4	MPa
Porosity ϕ	0.05	/
Fluid viscosity μ	3 $\times$ 10^–4	Pa·s
Initial reservoir pressure p_r	0.10	MPa
Rock bulk modulus K_s	10	GPa
Fluid bulk modulus K_f	0.6	GPa
Reservoir permeability k_r	$1 \times 10^{- 16}$	m²
Critical energy release rate G_c	250	N/m

Figure 10 shows the simulation result of the HF propagation path. We can find the propagation path changes with the different distance. When H_c is 10, the HF propagates horizontally along the X direction (the maximum principal stress direction) and finally is connected with the cave. When H_c is 12, the HF first propagates horizontally and then deflects up near the location of cave. Finally, the HF is connected with the cave. When H_c is 13.2, the HF also deflects up and is connected with the cave. Yet, the location where HF begins to deflect is farther away from the cave compared with H_c is 12. When H_c is 15.5, the HF propagates along the X direction throughout. The cave does not affect the HF propagation path.

FIG. 10.

$The phase-field evolution with different distance between fracture and cave.$

View large Download slide

The phase-field evolution with different distance between fracture and cave.

Through the simulation result, it is found the stress concentration caused by cave can play a notable role in the HF propagation. It can induce the HF to propagate away from the direction of maximum principal stress. However, when the distance between fracture and cave is very large, this stress concentration is negligible to the HF propagation.

The pressure distribution with different distance is shown in Fig. 11. As seen from Fig. 11(d), the HF is not connected with the cave. The fluid pressure in the fracture domain is maximum. Due to the fluid leak off, the pressure increases in the reservoir area, which is close to the fracture domain. However, the fluid pressure in the cave is still equal to the initial reservoir pressure. It is because the cave is far away from the HF, so the injection fluid does not enter the cave. As seen from Figs. 11(a)–11(c), the pressure in the cave significantly increases. It is because the HF is connected with the cave, resulting in high-pressure injection fluid into the cave.

FIG. 11.

$The pressure distribution with different distances between fracture and cave.$

View large Download slide

The pressure distribution with different distances between fracture and cave.

The fluid pressure at injection point is a key parameter of hydraulic fracturing. It can be measured in engineering practice and has an important guiding significance for the optimization design of hydraulic fracturing. In this work, we analyze the pressure variation at injection point in detail. When H_c is 13.2, the fluid pressure at injection point vs time is shown in Fig. 12. The pressure rapidly increases before the fracture initiation. When the HF propagates, the pressure gradually declines. However, when time is 12.0 s, the fluid pressure steeply declines. According to the phase-field evolution contour, we find the HF is connected with cave at this time. After the connection, the injection fluid enters the cave. The fluid pressure at injection point begins to gradually increase.

FIG. 12.

View large Download slide

The fluid pressure at injection point vs time when H_c = 13.2.

Figure 13 presents the fluid pressure vs time at varied H_c. When H_c is 15.5, the HF is not connected with the cave. Hence, the pressure curve does not experience a sharp drop. Meanwhile, the gradual growth stage is not existent in the pressure curve. Comparing the curve of H_c is 13.2 and H_c is 12, we can find the pressure experiences a sharp drop earlier when H_c is 12. It is because the distance between fracture and cave is shorter when H_c is 12. Hence, the HF is connected with the cave earlier. Finally, the fluid pressure in the three cases all tends to be stable.

FIG. 13.

View large Download slide

The fluid pressure at injection point vs time with different H_c.

B. Effect of injection rate

Injection rate is a critical factor of hydraulic fracturing operation. Engineers can improve the efficacy of hydraulic fracturing by adjusting the injection rate. Here, our study pays attention to how the injection rate affects the interaction between HF and cave. The position of cave center is (14, 12.7). The injection rate Q is set as 1.0 $\times$ 10⁻⁴, 1.4 $\times$ 10⁻⁴, 1.8 $\times$ 10⁻⁴, and 2.2 $\times$ 10⁻⁴ m²/s. Other parameters conform to Table II.

Figure 14 presents the phase-field evolution contour. We can find the propagation path varies with different injection rate. The interaction pattern can be divided into two types. When Q is 1.0 $\times$ 10⁻⁴ and 1.4 $\times$ 10⁻⁴ m²/s, the HF is finally connected with the cave. When Q is 1.8 $\times$ 10⁻⁴ and 2.2 $\times$ 10⁻⁴ m²/s, the HF deflects along the cave, but not connected with it. For a clearer comparison of the four cases, we use v = 0.4 to display the HF profile on the right side of injection point, which is shown in Fig. 15. Here, we define the deflection point, which represents the point where the HF begins to deflect up. It can be seen the deflection point is farther away from injection point when Q increases. In addition, deflection angle of HF deviation from X direction decreases as Q increases. The result reveals that the cave effect on HF declines when the injection rate increases. Hence, the HF is more inclined to propagate along the maximum principal stress direction.

FIG. 14.

View large Download slide

The phase-field evolution with different injection rates.

FIG. 15.

$The hydraulic fracture propagation profile with different injection rates.$

View large Download slide

The hydraulic fracture propagation profile with different injection rates.

Figure 16 shows the fluid pressure at injection point with different injection rates. We can find the fracture initiation pressure raises as the injection rate increases. When Q is 1.0 $\times$ 10⁻⁴ and 1.4 $\times$ 10⁻⁴ m²/s, the fluid pressure experiences a sharp drop. It is because the HF is connected with cave. However, when Q is 1.8 $\times$ 10⁻⁴ m²/s, the fluid pressure does not experience a sharp drop since the HF is unable to connect with cave.

FIG. 16.

View large Download slide

The fluid pressure at injection point vs time with different injection rates.

C. Effect of cave radius

The cave radius r_c can also markedly impact the HF–cave interaction. In above section, we have found the stress concentration caused by cave can make the HF deflect. In this section, we devote more attention to investigate the impact on the propagation velocity of HF.

We conduct three cases here. In case 1, the cave radius is set as 0, which indicates there is no cave in the reservoir. In case 2 and case 3, the cave radius is set as 0.5 and 0.9 m, respectively. The coordinates of the cave center are (15, 10). Other parameters conform to Table II.

Figure 17 presents the phase-field evolution contour. When time is 6.0 s, we can find the fracture length is different in the three cases. In case 1, the fracture length is equal on both sides of injection point. However, in case 3, it is obvious that the fracture length on the right side of injection point is longer than that on the left side. In case 2, the length difference between the two sides is very small, with the right-side fracture length of 1.93 m and the left-side fracture length of 1.86 m. When time is 9.5 s, the HF is connected with the cave in case 3. The right-side fracture length is much longer than the left-side fracture length. When time is 19.0 s, the right-side fracture length does not increase in case 3, since the HF has been connected with cave, while the left-side fracture length has increased a little. In case 2, the right-side fracture length is 4.5 m, and the left-side fracture length is 3.9 m. There is a noticeable length difference between the two sides.

FIG. 17.

View large Download slide

The phase-field evolution with different cave radii.

Figure 18 presents the right-side fracture length vs time. The HF begins to propagate when time is 2.5 s in all three cases. Then, it can be found the HF propagation velocity is different among the three cases. The HF propagates faster as the cave radius increases. Contrasting case 1 and case 2, we can find that the existence of cave can make the HF propagate faster. Moreover, if there is no cave, it is apparent that the propagation velocity will gradually decrease with time. However, if the cave exists, the propagation velocity can suddenly increase near the cave.

FIG. 18.

$The fracture length on the right side of injection point vs time.$

View large Download slide

The fracture length on the right side of injection point vs time.

Some significant findings can be drawn from the simulation result. The stress concentration caused by cave is similar to an attraction force for the HF. If the cave location is far away from the maximum principal stress direction, it can make the HF deflect. If the cave location is in the direction of maximum principal stress, it can make the HF propagate faster. The propagation velocity raises as the cave radius increases. In addition, it can cause the asymmetric propagation of HF. The propagation velocity is greater on the side with cave.

D. Effect of in situ stress in three-dimensional example

The in situ stress can alter the stress concentration caused by cave, which further influences the HF–cave interaction. Here, we set a three-dimensional example to investigate the interaction under different in situ stress.

As discussed in Sec. III A, we use a structured eight-node element to discretize the spatial domain. Hence, our model is applicable to the three-dimensional simulation. The schematic of the three-dimensional example is shown in Fig. 19. There is a cubic computation domain, measuring 5 $\times$ 5 $\times$ 5 m in magnitude. Its center is placed a penny-shaped fracture, with 0.1 m radius. The coordinates of point 0 are (0, 0, 0). Hence, the coordinates of fracture center point are (2.5, 2.5, and 2.5). In addition, there is natural cave at the upper right of the penny-shaped fracture. The cave center is located at (3.5, 2.5, and 3.3), and the cave radius is 0.4 m. The penny-shaped fracture is pumped with fracturing fluid at a steady rate of $5 \times 10^{- 4} m^{3} / s$ ⁠. To investigate the effect of in situ stress, we set three cases here: case 1: σ_x = 1, σ_y = 1, and σ_z = 0.8 MPa; case 2: σ_x = 1, σ_y = 1, and σ_z = 0.6 MPa; case 3: σ_x = 1, σ_y = 1, and σ_z = 0.2 MPa. Other parameters are identical to those in Table II. In this simulation example, the element size is $0.05 \times 0.05 \times 0.05 m,$ and there are 1 030 301 nodes in total. We use the Intel Xeon Platinum 9242 CPU with 96 cores for parallel computing. The case 1 is computed within 12 h for 4000 time steps.

FIG. 19.

View large Download slide

Schematic of the cubic computation domain.

The propagation path of the penny-shaped fracture is presented in Fig. 20. We use the phase-field v = 0.25 to display the three-dimensional fracture shape and cave shape. In case 1, the in situ stress difference $σ_{d} = σ_{x} - σ_{z} = 0.2 MPa .$ In cases 2 and 3, the in situ stress difference is 0.4 and 0.8 MPa, respectively. It is generally believed that the penny-shaped fracture propagates in the direction orthogonal to the minimum in situ stress. However, in case 1, the penny-shaped fracture deflects up. It is because the stress concentration caused by cave is similar to an attraction force, which can change the fracture propagation path. In case 2, the penny-shaped fracture also deflects up. However, the deflection angle of penny-shaped fracture deviation from x–y plane decreases. In case 3, the penny-shaped fracture propagates horizontally, and it does not deflect. From the three cases, we can find the impact of cave stress concentration declines as the in situ stress difference increases. Hence, the penny-shaped fracture is more inclined to propagate in the direction orthogonal to the minimum in situ stress.

FIG. 20.

$The penny-shaped fracture propagation path in three-dimensional example.$

View large Download slide

The penny-shaped fracture propagation path in three-dimensional example.

Next, we test the capacity of our model in handling three-dimensional issues when multiple caves exist. We set two cases for evaluation. In case 3, there are two caves in the computation domain, with 0.4 m radius. Their locations are (3.5, 2.5, and 3.2) and (1.5, 2.5, and 1.8), respectively. In case 4, we add a new cave based on case 3. The new cave is located at (3.5, 2.5, and 1.8), and its radius is 0.25 m. The in situ stress of both case 3 and case 4 is the same as case 1. The fracture propagation path is presented in Fig. 21. When t = 100 s, the fracture predominantly propagates horizontally. Comparing Figs. 21(a) and 21(c), we can find the fracture propagates farther to the right side in case 4. It is because the two caves on the right side create larger attraction force in case 4. When t = 400 s, the fracture is connected with the left and right caves in case 3. However, in case 4, the fracture is not connected with the right cave. It only deflects up a little.

FIG. 21.

$The penny-shaped fracture propagation path in multiple cave cases.$

View large Download slide

The penny-shaped fracture propagation path in multiple cave cases.

The above example substantiates the capacity in our model to simulate the HF propagation with natural cave in three-dimensional issues. Moreover, we use the fixed mesh to characterize the natural cave and simulate the HF propagation. The proposed model does not need any special node treatment for characterizing the cave. In addition, it does not need any remeshing (the element size is smaller than $ϵ$ in the whole computation domain). These advantages fully prove the conciseness and efficiency of the proposed model.

VI. MULTIPLE CAVE SIMULATIONS

Sometimes, there are multiple natural caves in the geomaterials. These caves are distributed randomly and vary in size. The following example is presented to check the performance of our model in multiple cave issues. The geometric structure of this example is shown in Fig. 22. The computation domain is 20 × 20 m² in size, and its center is placed an initial fracture. The initial fracture length is 0.6 m, and the angle between the fracture and horizontal direction is 25°. The four outer boundaries of computation domain are fixed. The injection rate of fracturing fluid is 1 $\times$ 10⁻⁴ m²/s. Within the computation domain, there are six natural caves and their radii range from 0.85 to 1.1 m. In addition, the distances between the initial fracture and caves are 3.9–5.5 m. Other parameters are identical to those in Table II.

FIG. 22.

View large Download slide

The geometric structure for multiple cave issues.

Figure 23 shows the phase-field evolution contour. We can found the HF propagation path is significantly influenced by the multiple caves. When t is 2.5 s, the HF begins to extend. After 1.5 s, the HF has evidently grown. Because of the stress concentration caused by caves, the propagation path is not straight. When t is 8.1 s, the HF is connected with the lower-left cave. Since the lower two caves are larger than the up two caves, the attraction force of the lower two caves is larger. Therefore, the HF propagates faster along the lower-left direction. Then, the HF continues to extend and is connected with the up-right cave when t is 20.3 s.

FIG. 23.

View large Download slide

The phase-field evolution in multiple cave issues.

Figure 24 shows the fluid pressure variation at injection point over time. The pressure reaches its peak when fracture initiation. After that, the pressure gradually decreases. Subsequently, we can found the pressure experiences two sharp drops. The reason is that the HF is connected with the caves at these two time points. The finding is consistent with the study in Sec. V.

FIG. 24.

View large Download slide

The fluid pressure at injection point vs time in multiple cave issues.

This example demonstrates that our model can simulate the HF propagation in the geomaterials containing multiple caves. It also proves the wide applicability of our model.

VII. CONCLUSIONS

In this study, we establish a novel hydraulic fracturing model for the fluid-driven fracture propagation in poroelastic media containing the natural cave by the phase-field method. First, we derive the fracture-cave-reservoir flow governing equations. In reservoir domain, the fluid flow adheres to low-permeability Darcy flow; in cave domain, it follows free flow, while in fracture domain, it conforms to generalized Reynolds flow. The hydro-mechanical coupling is achieved by employing Biot poroelasticity theory and updating fracture width. The phase field is introduced to diffuse the sharp fracture and cave edge. Compared to existing study, our model does not need additional enrichment function to characterize the cave. The proposed model can naturally detect and simulate the hydraulic fracture (HF) propagation behavior in the existence of natural cave. The complicate calculation of stress intensity factor when HF approaches cave is avoided. The effectiveness of the proposed model in two-dimensional issue is fully displayed. Then, the model is extended to the three-dimensional issue and multiple cave issues.

The interaction mechanism between HF and natural cave is also investigated in accordance with the proposed model. Several findings are obtained: (1) The stress concentration caused by cave can affect the HF propagation. If the cave is far away from the maximum principal stress direction, it can make the HF deflect toward the cave. (2) The cave effect can change the HF propagation velocity. The fracture propagates faster on the side with cave. (3) Increasing the fracture–cave distance, injection rate, and in situ stress difference can all decline the cave effect. (4) The cave effect raises as the cave radius increases. (5) The fluid pressure at injection point will experience a sharp drop if the HF is connected with cave. After the connection, the fluid pressure will slowly increase.

ACKNOWLEDGMENTS

The authors gratefully acknowledge the funding support from Fourteenth Five-Year Major Science and Technology Project of CNPC (Grant No. 2021DJ1505).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Jie Jin: Data curation (lead); Investigation (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (lead). Xiaoqiang Wang: Formal analysis (lead); Methodology (lead); Software (lead). Xiaohua Liu: Funding acquisition (lead); Resources (equal). Yanmei Xu: Resources (equal); Supervision (lead). Detang Lu: Conceptualization (lead); Project administration (lead).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX: THE DISCRETIZING PROCESS OF THE FLOW EQUATIONS

Here, we define a indicator function $v_{c}$ to distinguish between cave domain and non-cave domain. Specifically, it has the same form as $v$ in Eq. (31). In non-cave domain, $v_{c} = (0, 1]$ ⁠. However, in cave domain, $v_{c} = 0$ ⁠. Let $Ω_{d} = Ω_{r} \cup Γ$ ⁠, it represents the non-cave domain (reservoir domain and fracture domain). First, the integral in $Γ$ and $Ω_{r}$ needs to be transformed to the integral in $Ω_{d}$ ⁠.

Chukwudozie et al.³⁴ proposed a regularized method to transform the integral in $Γ$ ⁠. It is realized through the integral term in fracture domain $Γ$ multiplied by $| \nabla v | .$ Following their method, the integral in $Ω_{r}$ is transformed through the integral term in $Ω_{r}$ multiplied by $v^{2}$ ⁠. To avoid the ill-conditioned system when $v = 0$ ⁠, they proposed that Biot modulus term can be used as a stabilizing term.

To sum up, the flow governing equation in non-cave domain becomes

\begin{matrix} \int_{Ω_{d}} ((\frac{ψ}{M} + v^{2} ψ \frac{α^{2}}{κ}) \frac{\partial p^{l + 1, j + 1}}{\partial t} + \frac{k_{r}}{μ} \nabla p^{l + 1, j + 1} \cdot \nabla ψ) d V + \int_{Ω_{d}} \frac{{(w^{j})}^{3}}{12 μ} (\nabla_{Γ} p^{l + 1, j + 1} \cdot \nabla_{Γ} ψ) | \nabla v | d V \\ = - \int_{Ω_{d}} ψ \frac{\partial w^{j}}{\partial t} | \nabla v | d V - \int_{Ω_{d}} ψ α \frac{\partial ε_{vol}^{j}}{\partial t} v^{2} d V + \int_{Ω_{d}} ψ \frac{α^{2}}{κ} \frac{\partial p^{j}}{\partial t} v^{2} d V - \int_{\partial_{N} Ω_{d}} ψ q_{n} d S \\ + \int_{Ω_{d}} ψ q_{r} v^{2} d V + \int_{Ω_{d}} ψ q_{f} | \nabla v | d V + \int_{Ω_{d}} (1 - v^{2}) \frac{ψ}{M} \frac{\partial p^{j}}{\partial t} d V + \int_{\partial Ω_{c}} ψ q_{c r}^{l} d S + \int_{\partial Γ_{c}} ψ w^{j} V_{f}^{l} \cdot n_{c} d s \end{matrix} .

(A1)

Since the fracture is diffused by PFM, the fluid exchange term between fracture and cave needs to be regularized. According to Chukwudozie's method,^34,80 the regularized form is as follows:

\int_{\partial Γ_{c}} ψ w V_{f} \cdot n_{c} d s ≃ \int_{\partial Ω_{c f}} ψ w V_{f} | \nabla v | \cdot n_{c} d S .

(A2)

The function inside the integral symbol on the right side of above formula can be seen as the “regularized” fluid exchange flux.

q_{c f} = w V_{f} | \nabla v | \cdot n_{c},

(A3)

where

\partial Ω_{c f}

is the diffused boundary between fracture and cave.

The pressure gradient over fracture face

\nabla_{Γ} p

is complicated. Thus, we need to approximate this term.

\nabla v

is a vector, which points toward the orientation where

v

changes fastest. The region with the fastest change of

v

is the transition region. Therefore, the fracture normal

n_{Γ}

can be approximated as

n_{Γ} \approx \frac{\nabla v}{| \nabla v |} = n_{Γ x} e_{x} + n_{Γ y} e_{y} + n_{Γ z} e_{z} .

(A4)

Then,

\nabla_{Γ} p

can be computed as

\nabla_{Γ} p \approx \nabla p - (\nabla p \cdot n_{Γ}) n_{Γ} \approx \nabla p - (\nabla p \cdot \frac{\nabla v}{| \nabla v |}) \frac{\nabla v}{| \nabla v |} .

(A5)

By substituting Eq. into Eq. , the specific form of

\nabla_{Γ} p

is obtained

\begin{matrix} \nabla_{Γ} p = [\frac{\partial p}{\partial x} (1 - n_{Γ x}^{2}) - \frac{\partial p}{\partial y} n_{Γ y} n_{Γ x} - \frac{\partial p}{\partial z} n_{Γ z} n_{Γ x}] e_{x} + [\frac{\partial p}{\partial y} (1 - n_{Γ y}^{2}) - \frac{\partial p}{\partial x} n_{Γ x} n_{Γ y} - \frac{\partial p}{\partial z} n_{Γ z} n_{Γ y}] e_{y} \\ + [\frac{\partial p}{\partial z} (1 - n_{Γ z}^{2}) - \frac{\partial p}{\partial x} n_{Γ x} n_{Γ z} - \frac{\partial p}{\partial y} n_{Γ y} n_{Γ z}] e_{z} \end{matrix} .

(A6)

Next, Eq. is discretized by the structured eight-node elements. Similar to Sec. III A, the test function

ψ

is rewritten as

ψ = N_{p} \hat{ψ}, \nabla ψ = B_{p} \hat{ψ},

(A7)

where

\hat{ψ}

is the vector composed of node values

ψ_{i}

⁠.

\hat{ψ} = {[\begin{matrix} ψ_{1} & ψ_{2} & \dots & ψ_{n} \end{matrix}]}^{T} .

The pressure gradient and test function gradient over fracture face are

\nabla_{Γ} p = B_{Γ p} \hat{p}, \nabla_{Γ} ψ = B_{Γ p} \hat{ψ},

(A8)

where

B_{Γ p} = [\begin{matrix} N_{1, x} (1 - n_{Γ x}^{2}) - N_{1, y} n_{Γ y} n_{Γ x} - N_{1, z} n_{Γ z} n_{Γ x} & \dots & N_{n, x} (1 - n_{Γ x}^{2}) - N_{n, y} n_{Γ y} n_{Γ x} - N_{n, z} n_{Γ z} n_{Γ x} \\ N_{1, y} (1 - n_{Γ y}^{2}) - N_{1, x} n_{Γ x} n_{Γ y} - N_{1, z} n_{Γ z} n_{Γ y} & \dots & N_{n, y} (1 - n_{Γ y}^{2}) - N_{n, x} n_{Γ x} n_{Γ y} - N_{n, z} n_{Γ z} n_{Γ y} \\ N_{1, z} (1 - n_{Γ z}^{2}) - N_{1, x} n_{Γ x} n_{Γ z} - N_{1, y} n_{Γ y} n_{Γ z} & \dots & N_{n, z} (1 - n_{Γ z}^{2}) - N_{n, x} n_{Γ x} n_{Γ z} - N_{n, y} n_{Γ y} n_{Γ z} \end{matrix}] .

The discrete form of Eq. is as follows:

R^{p} = [M_{p}] \frac{d {\hat{p}}^{l + 1, j + 1}}{d t} + [S_{p}] {\hat{p}}^{l + 1, j + 1} - [f_{p}] = 0,

(A9)

where

R^{p}

is the residuals,

M_{p}

is the compressibility matrix,

S_{p}

is the transportation matrix, and

f_{p}

is the coupling vector.

M_{p} = \int_{Ω_{d}} (\frac{N_{p}^{T} N_{p}}{M} + v^{2} N_{p}^{T} N_{p} \frac{α^{2}}{κ}) d V,

S_{p} = \int_{Ω_{d}} [\frac{k_{r}}{μ} B_{p}^{T} B_{p} + \frac{{(w^{j})}^{3}}{12 μ} | \nabla v | B_{Γ p}^{T} B_{Γ p}] d V,

\begin{matrix} f_{p} = \int_{Ω_{d}} N_{p}^{T} (- \frac{\partial w^{j}}{\partial t} | \nabla v | - α \frac{\partial ε_{vol}^{j}}{\partial t} v^{2} + \frac{α^{2}}{κ} \frac{\partial p^{j}}{\partial t} v^{2} + \frac{(1 - v^{2})}{M} \frac{\partial p^{j}}{\partial t} + q_{r} v^{2} + q_{f} | \nabla v |) d V \\ - \int_{\partial_{N} Ω_{d}} N_{p}^{T} q_{n} d S + \int_{\partial Ω_{c}} N_{p}^{T} q_{c r}^{l} d S + \int_{\partial Ω_{c f}} N_{p}^{T} q_{c f}^{l, j} d S \end{matrix} .

We use backward Euler scheme to discretize time term. Hence, the final discrete form is

R^{p} = [M_{p} + Δ t S_{p}] {[\hat{p}]}_{n + 1}^{l + 1, j + 1} - [M_{p}] {[\hat{p}]}_{n} - Δ t {[f_{p}]}_{n + 1} = 0 .

(A10)

In the cave domain, the flow governing equation is also discretized by backward Euler scheme.

p_{n + 1}^{l + 1} = p_{n} - \frac{Δ t}{V c_{f}} (\int_{\partial Ω_{c}} q_{c r}^{l + 1} d S + \int_{\partial Ω_{c f}} q_{c f}^{l + 1, j} d S) .

(A11)

REFERENCES

1.

A. C.

Correas

,

M.

Corrado

,

A.

Sapora

, and

P.

Cornetti

, “

Size-effect on the apparent tensile strength of brittle materials with spherical cavities

,”

Theor. Appl. Fract. Mech.

116

,

103120

(

2021

).

https://doi.org/10.1016/j.tafmec.2021.103120

Google Scholar

Crossref

2.

H. S.

Suh

and

W.

Sun

, “

An immersed phase field fracture model for microporomechanics with Darcy–Stokes flow

,”

Phys. Fluids

33

(

1

),

016603

(

2021

).

https://doi.org/10.1063/5.0035602

Google Scholar

Crossref

3.

M.

Rasoulzadeh

, “

Anomalous transport in a porous medium with randomly packed ellipse cavities

,”

Phys. Fluids

34

(

12

),

126607

(

2022

).

https://doi.org/10.1063/5.0123491

Google Scholar

Crossref

4.

X.

Du

,

Q.

Li

,

P.

Li

,

Y.

Xian

,

Y.

Zheng

, and

D.

Lu

, “

A novel pressure and rate transient analysis model for fracture-caved carbonate reservoirs

,”

J. Pet. Sci. Eng.

208

,

109609

(

2022

).

https://doi.org/10.1016/j.petrol.2021.109609

Google Scholar

Crossref

5.

Z.

Huang

,

H.

Wang

,

X.

Su

, and

W.

Liao

, “

A two-scale fractal permeability model for vuggy porous media

,”

Phys. Fluids

35

(

2

),

027116

(

2023

).

https://doi.org/10.1063/5.0135861

Google Scholar

Crossref

6.

F.

Zhang

,

M.

An

,

B.

Yan

,

Y.

Wang

, and

Y.

Han

, “

A novel hydro-mechanical coupled analysis for the fractured vuggy carbonate reservoirs

,”

Comput. Geotech.

106

,

68

–

82

(

2019

).

https://doi.org/10.1016/j.compgeo.2018.10.013

Google Scholar

Crossref

7.

Z.

Liu

,

X.

Tang

,

S.

Tao

,

G.

Zhang

, and

M.

Chen

, “

Mechanism of connecting natural caves and wells through hydraulic fracturing in fracture-cavity reservoirs

,”

Rock Mech. Rock Eng.

53

(

12

),

5511

–

5530

(

2020

).

https://doi.org/10.1007/s00603-020-02225-w

Google Scholar

Crossref

8.

P.

Liu

,

Y.

Ju

,

F.

Gao

,

P. G.

Ranjith

, and

Q.

Zhang

, “

CT identification and fractal characterization of 3-D propagation and distribution of hydrofracturing cracks in low-permeability heterogeneous rocks

,”

J. Geophys. Res.

123

(

3

),

2156

–

2173

, https://doi.org/10.1002/2017JB015048 (

2018

).

Google Scholar

Crossref

9.

S.

Zhou

,

X.

Zhuang

, and

T.

Rabczuk

, “

Phase-field modeling of fluid-driven dynamic cracking in porous media

,”

Comput. Methods Appl. Mech. Eng.

350

,

169

–

198

(

2019

).

https://doi.org/10.1016/j.cma.2019.03.001

Google Scholar

Crossref

10.

B.

Liu

,

Y.

Jin

, and

M.

Chen

, “

Influence of vugs in fractured-vuggy carbonate reservoirs on hydraulic fracture propagation based on laboratory experiments

,”

J. Struct. Geol.

124

,

143

–

150

(

2019

).

https://doi.org/10.1016/j.jsg.2019.04.007

Google Scholar

Crossref

11.

Y.

Wang

,

X.

Li

,

B.

Zhao

, and

Z.

Zhang

, “

3D numerical simulation of pulsed fracture in complex fracture-cavitied reservoir

,”

Comput. Geotech.

125

,

103665

(

2020

).

https://doi.org/10.1016/j.compgeo.2020.103665

Google Scholar

Crossref

12.

Z.

Zhou

,

W.

Wang

, and

T.

Xu

, “

Methodology of experiment to investigate connection between hydraulic fractures and vugs in fracture–vug carbonate formations

,”

J. Energy Resour. Technol.

143

(

11

),

113202

(

2021

).

https://doi.org/10.1115/1.4049779

Google Scholar

Crossref

13.

S.

Entrekin

,

M.

Evans-White

,

B.

Johnson

, and

E.

Hagenbuch

, “

Rapid expansion of natural gas development poses a threat to surface waters

,”

Front. Ecol. Environ.

9

(

9

),

503

–

511

(

2011

).

https://doi.org/10.1890/110053

Google Scholar

Crossref

14.

S. G.

Osborn

,

A.

Vengosh

,

N. R.

Warner

, and

R. B.

Jackson

, “

Methane contamination of drinking water accompanying gas-well drilling and hydraulic fracturing

,”

Proc. Natl. Acad. Sci. U. S. A.

108

(

20

),

8172

–

8176

(

2011

).

https://doi.org/10.1073/pnas.1100682108

Google Scholar

Crossref

PubMed

15.

T.

Belytschko

and

T.

Black

, “

Elastic crack growth in finite elements with minimal remeshing

,”

Int. J. Numer. Methods Eng.

45

(

5

),

601

–

620

(

1999

).

https://doi.org/10.1002/(SICI)1097-0207(19990620)45:5<601::AID-NME598>3.0.CO;2-S

Google Scholar

Crossref

16.

M.

Stolarska

,

D. L.

Chopp

,

N.

Moës

, and

T.

Belytschko

, “

Modelling crack growth by level sets in the extended finite element method

,”

Int. J. Numer. Methods Eng.

51

(

8

),

943

–

960

(

2001

).

https://doi.org/10.1002/nme.201

Google Scholar

Crossref

17.

E.

Gordeliy

and

A.

Peirce

, “

Coupling schemes for modeling hydraulic fracture propagation using the XFEM

,”

Comput. Methods Appl. Mech. Eng.

253

,

305

–

322

(

2013

).

https://doi.org/10.1016/j.cma.2012.08.017

Google Scholar

Crossref

18.

T.

Mohammadnejad

and

A. R.

Khoei

, “

An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model

,”

Finite Elem. Anal. Des.

73

,

77

–

95

(

2013

).

https://doi.org/10.1016/j.finel.2013.05.005

Google Scholar

Crossref

19.

O.

Nguyen

,

E. A.

Repetto

,

M.

Ortiz

, and

R. A.

Radovitzky

, “

A cohesive model of fatigue crack growth

,”

Int. J. Fract.

110

(

4

),

351

–

369

(

2001

).

https://doi.org/10.1023/A:1010839522926

Google Scholar

Crossref

20.

H.

Shimizu

,

S.

Murata

, and

T.

Ishida

, “

The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution

,”

Int. J. Rock Mech. Min. Sci.

48

(

5

),

712

–

727

(

2011

).

https://doi.org/10.1016/j.ijrmms.2011.04.013

Google Scholar

Crossref

21.

S.

Deng

,

H.

Li

,

G.

Ma

,

H.

Huang

, and

X.

Li

, “

Simulation of shale-proppant interaction in hydraulic fracturing by the discrete element method

,”

Int. J. Rock Mech. Min. Sci.

70

,

219

–

228

(

2014

).

https://doi.org/10.1016/j.ijrmms.2014.04.011

Google Scholar

Crossref

22.

K.

Wu

and

J. E.

Olson

, “

Simultaneous multifracture treatments: Fully coupled fluid flow and fracture mechanics for horizontal wells

,”

SPE J.

20

(

2

),

337

–

346

(

2015

).

https://doi.org/10.2118/167626-PA

Google Scholar

Crossref

23.

Z.-Z.

Yang

,

L.-P.

Yi

,

X.-G.

Li

, and

W.

He

, “

Pseudo-three-dimensional numerical model and investigation of multi-cluster fracturing within a stage in a horizontal well

,”

J. Pet. Sci. Eng.

162

,

190

–

213

(

2018

).

https://doi.org/10.1016/j.petrol.2017.12.034

Google Scholar

Crossref

24.

E.

Tanné

,

T.

Li

,

B.

Bourdin

,

J. J.

Marigo

, and

C.

Maurini

, “

Crack nucleation in variational phase-field models of brittle fracture

,”

J. Mech. Phys. Solids

110

,

80

–

99

(

2018

).

https://doi.org/10.1016/j.jmps.2017.09.006

Google Scholar

Crossref

25.

S.

Zhou

and

X.

Zhuang

, “

Phase field modeling of hydraulic fracture propagation in transversely isotropic poroelastic media

,”

Acta Geotech.

15

(

9

),

2599

–

2618

(

2020

).

https://doi.org/10.1007/s11440-020-00913-z

Google Scholar

Crossref

26.

G. A.

Francfort

and

J. J.

Marigo

, “

Revisiting brittle fracture as an energy minimization problem

,”

J. Mech. Phys. Solids

46

(

8

),

1319

–

1342

(

1998

).

https://doi.org/10.1016/S0022-5096(98)00034-9

Google Scholar

Crossref

27.

B.

Bourdin

,

G. A.

Francfort

, and

J. J.

Marigo

, “

Numerical experiments in revisited brittle fracture

,”

J. Mech. Phys. Solids

48

(

4

),

797

–

826

(

2000

).

https://doi.org/10.1016/S0022-5096(99)00028-9

Google Scholar

Crossref

28.

B.

Bourdin

,

G. A.

Francfort

, and

J.-J.

Marigo

, “

The variational approach to fracture

,”

J. Elasticity

91

(

1–3

),

5

–

148

(

2008

).

https://doi.org/10.1007/s10659-007-9107-3

Google Scholar

Crossref

29.

C.

Miehe

,

S.

Mauthe

, and

S.

Teichtmeister

, “

Minimization principles for the coupled problem of Darcy–Biot–type fluid transport in porous media linked to phase field modeling of fracture

,”

J. Mech. Phys. Solids

82

,

186

–

217

(

2015

).

https://doi.org/10.1016/j.jmps.2015.04.006

Google Scholar

Crossref

30.

S.

Lee

,

M. F.

Wheeler

, and

T.

Wick

, “

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

(

2016

).

https://doi.org/10.1016/j.cma.2016.02.037

Google Scholar

Crossref

31.

Z. A.

Wilson

and

C. M.

Landis

, “

Phase-field modeling of hydraulic fracture

,”

J. Mech. Phys. Solids

96

,

264

–

290

(

2016

).

https://doi.org/10.1016/j.jmps.2016.07.019

Google Scholar

Crossref

32.

Y.

Heider

and

B.

Markert

, “

A phase-field modeling approach of hydraulic fracture in saturated porous media

,”

Mech. Res. Commun.

80

,

38

–

46

(

2017

).

https://doi.org/10.1016/j.mechrescom.2016.07.002

Google Scholar

Crossref

33.

S.

Zhou

,

X.

Zhuang

, and

T.

Rabczuk

, “

A phase-field modeling approach of fracture propagation in poroelastic media

,”

Eng. Geol.

240

,

189

–

203

(

2018

).

https://doi.org/10.1016/j.enggeo.2018.04.008

Google Scholar

Crossref

34.

C.

Chukwudozie

,

B.

Bourdin

, and

K.

Yoshioka

, “

A variational phase-field model for hydraulic fracturing in porous media

,”

Comput. Methods Appl. Mech. Eng.

347

,

957

–

982

(

2019

).

https://doi.org/10.1016/j.cma.2018.12.037

Google Scholar

Crossref

35.

D. T.

Van

,

D.

Duc Hong

,

D.

Nguyen Dinh

, and

B.

Tinh Quoc

, “

Phase-field thermal buckling analysis for cracked functionally graded composite plates considering neutral surface

,”

Compos. Struct.

182

,

542

–

548

(

2017

).

https://doi.org/10.1016/j.compstruct.2017.09.059

Google Scholar

Crossref

36.

Y.

Xiao

,

Z.

Zeng

,

L.

Zhang

,

J.

Wang

,

Y.

Wang

,

H.

Liu

, and

C.

Huang

, “

A spectral element-based phase field method for incompressible two-phase flows

,”

Phys. Fluids

34

(

2

),

022114

(

2022

).

https://doi.org/10.1063/5.0077372

Google Scholar

Crossref

37.

X.

Zhuang

,

S.

Zhou

,

G. D.

Huynh

,

P.

Areias

, and

T.

Rabczuk

, “

Phase field modeling and computer implementation: A review

,”

Eng. Fract. Mech.

262

,

108234

(

2022

).

https://doi.org/10.1016/j.engfracmech.2022.108234

Google Scholar

Crossref

38.

D. J.

Youn

, “

Hydro-mechanical coupled simulation of hydraulic fracturing using the eXtended Finite Element Method (XFEM)

,” Ph. D. thesis (

Colorado School of Mines

,

2016

).

Google Scholar

39.

L.

Cheng

,

Z.

Luo

,

Y.

Yu

,

L.

Zhao

, and

C.

Zhou

, “

Study on the interaction mechanism between hydraulic fracture and natural karst cave with the extended finite element method

,”

Eng. Fract. Mech.

222

,

106680

(

2019

).

https://doi.org/10.1016/j.engfracmech.2019.106680

Google Scholar

Crossref

40.

Z.

Luo

,

N.

Zhang

,

L.

Zhao

,

J.

Zeng

,

P.

Liu

, and

N.

Li

, “

Interaction of a hydraulic fracture with a hole in poroelasticity medium based on extended finite element method

,”

Eng. Anal. Boundary Elem.

115

,

108

–

119

(

2020

).

https://doi.org/10.1016/j.enganabound.2020.03.011

Google Scholar

Crossref

41.

Z.

Liu

,

Q.

Lu

,

Y.

Sun

,

X.

Tang

,

Z.

Shao

, and

Z.

Weng

, “

Investigation of the influence of natural cavities on hydraulic fracturing using phase field method

,”

Arab. J. Sci. Eng.

44

(

12

),

10481

–

10501

(

2019

).

https://doi.org/10.1007/s13369-019-04122-z

Google Scholar

Crossref

42.

A.

Azarov

,

A.

Patutin

, and

S.

Serdyukov

, “

Hydraulic fracture propagation near the cavity in a poroelastic media

,”

Appl. Sci.

11

(

22

),

11004

(

2021

).

https://doi.org/10.3390/app112211004

Google Scholar

Crossref

43.

S.

Liu

,

Z.

Liu

, and

Z.

Zhang

, “

Numerical study on hydraulic fracture-cavity interaction in fractured-vuggy carbonate reservoir

,”

J. Pet. Sci. Eng.

213

,

110426

(

2022

).

https://doi.org/10.1016/j.petrol.2022.110426

Google Scholar

Crossref

44.

R.

Camacho-Velazquez

,

M.

Vasquez-Cruz

,

R.

Castrejon-Aivar

, and

V.

Arana-Ortiz

, “

Pressure-transient and decline-curve behavior in naturally fractured vuggy carbonate reservoirs

,”

SPE Res. Eval. Eng.

8

(

4

),

95

–

266

(

2005

).

https://doi.org/10.2118/77689-PA

Google Scholar

Crossref

45.

B.

Gao

,

Z.-Q.

Huang

,

J.

Yao

,

X.-R.

Lv

, and

Y.-S.

Wu

, “

Pressure transient analysis of a well penetrating a filled cavity in naturally fractured carbonate reservoirs

,”

J. Pet. Sci. Eng.

145

,

392

–

403

(

2016

).

https://doi.org/10.1016/j.petrol.2016.05.037

Google Scholar

Crossref

46.

X.

Liu

,

H.

Wang

, and

L.

Chang

, “

Equivalent permeability model of dual-porosity and bi-dispersed porous media based on the intermingled fractal units

,”

Phys. Fluids

35

(

3

),

033606

(

2023

).

https://doi.org/10.1063/5.0140041

Google Scholar

Crossref

47.

M. A.

Biot

, “

General theory of three-dimensional consolidation

,”

J. Appl. Phys.

12

(

2

),

155

–

164

(

1941

).

https://doi.org/10.1063/1.1712886

Google Scholar

Crossref

48.

W.

Xu

,

H.

Yu

,

J.

Zhang

,

C.

Lyu

,

Q.

Wang

,

M.

Micheal

, and

H.

Wu

, “

Phase-field method of crack branching during SC-CO₂ fracturing: A new energy release rate criterion coupling pore pressure gradient

,”

Comput. Methods Appl. Mech. Eng.

399

,

115366

(

2022

).

https://doi.org/10.1016/j.cma.2022.115366

Google Scholar

Crossref

49.

P. A.

Witherspoon

,

J. S. Y.

Wang

,

K.

Iwai

, and

J. E.

Gale

, “

Validity of cubic law for fluid-flow in a deformable rock fracture

,”

Water Resour. Res.

16

(

6

),

1016

–

1024

, https://doi.org/10.1029/WR016i006p01016 (

1980

).

Google Scholar

Crossref

50.

V.

Martin

,

J.

Jaffre

, and

J. E.

Roberts

, “

Modeling fractures and barriers as interfaces for flow in porous media

,”

SIAM J. Sci. Comput.

26

(

5

),

1667

–

1691

(

2005

).

https://doi.org/10.1137/S1064827503429363

Google Scholar

Crossref

51.

C.

Chukwudozie

, “

Application of the variational fracture model to hydraulic fracturing in poroelastic media

,” Doctoral dissertation (LSU, 2016), see https://digitalcommons.lsu.edu/gradschool_dissertations/788.

52.

L. H.

Liu

,

Y. S.

Ma

,

B.

Liu

, and

C. L.

Wang

, “

Hydrothermal dissolution of Ordovician carbonates rocks and its dissolution mechanism in Tarim Basin, China

,”

Carbonates Evaporites

32

(

4

),

525

–

537

(

2017

).

https://doi.org/10.1007/s13146-016-0309-2

Google Scholar

Crossref

53.

V.

Baques

,

E.

Ukar

,

S. E.

Laubach

,

S. R.

Forstner

, and

A.

Fall

, “

Fracture, dissolution, and cementation events in Ordovician carbonate reservoirs, Tarim Basin, NW China

,”

Geofluids

2020

,

9037429

.

https://doi.org/10.1155/2020/9037429

Crossref

54.

P.

Popov

,

Y.

Efendiev

, and

G.

Qin

, “

Multiscale modeling and simulations of flows in naturally fractured karst reservoirs

,”

Commun. Comput. Phys.

6

(

1

),

162

–

184

(

2009

), see http://www.global-sci.com/intro/article_detail.html?journal=undefined&article_id=7676.

Google Scholar

55.

P.

Popov

,

G.

Qin

,

L.

Bi

,

Y.

Efendiev

,

Z.

Kang

, and

J.

Li

, “

Multiphysics and multiscale methods for modeling fluid flow through naturally fractured carbonate karst reservoirs

,”

SPE Reservoir Eval. Eng.

12

(

2

),

218

–

231

(

2009

).

https://doi.org/10.2118/105378-PA

Google Scholar

Crossref

56.

A. F.

Gulbransen

,

V. L.

Hauge

, and

K.-A.

Lie

, “

A multiscale mixed finite-element method for vuggy and naturally fractured reservoirs

,”

SPE J.

15

(

2

),

395

–

403

(

2010

).

https://doi.org/10.2118/119104-PA

Google Scholar

Crossref

57.

J.

Yao

,

Z.

Huang

,

Y.

Li

,

C.

Wang

, and

X.

Lv

, “

Discrete fracture-vug network model for modeling fluid flow in fractured vuggy porous media

,” in

International Oil and Gas Conference and Exhibition in China

,

2010

.

Google Scholar

Crossref

58.

X.

Zhang

,

Z.

Huang

,

Q.

Lei

,

J.

Yao

,

L.

Gong

,

S.

Sun

, and

Y.

Li

, “

Connectivity, permeability and flow channelization in fractured karst reservoirs: A numerical investigation based on a two-dimensional discrete fracture-cave network model

,”

Adv. Water Resour.

161

,

104142

(

2022

).

https://doi.org/10.1016/j.advwatres.2022.104142

Google Scholar

Crossref

59.

P.

Chen

,

X.

Wang

,

H.

Liu

,

Y.

Huang

, and

H.

Zhang

, “

A pressure-transient model for a fractured-vuggy carbonate reservoir with large-scale cave

,”

Geosyst. Eng.

19

(

2

),

69

–

76

(

2016

).

https://doi.org/10.1080/12269328.2015.1093965

Google Scholar

Crossref

60.

X.

Du

,

Q.

Li

,

Z.

Lu

,

P.

Li

,

Y.

Xian

,

Y.

Xu

,

D.

Li

, and

D.

Lu

, “

Pressure transient analysis for multi-vug composite fractured vuggy carbonate reservoirs

,”

J. Pet. Sci. Eng.

193

,

107389

(

2020

).

https://doi.org/10.1016/j.petrol.2020.107389

Google Scholar

Crossref

61.

Y.

Li

,

Q.

Yu

,

C.

Jia

,

P.

Liu

,

Q.

Wang

, and

D.

Wang

, “

Rate transient analysis for coupling Darcy flow and free flow in bead-string fracture-caved carbonate reservoirs

,”

J. Pet. Sci. Eng.

195

,

107809

(

2020

).

https://doi.org/10.1016/j.petrol.2020.107809

Google Scholar

Crossref

62.

C.

Wei

,

S.

Cheng

,

J.

Song

,

D.

Shi

,

R.

Shang

,

L.

Zhu

, and

H.

Yu

, “

Pressure transient analysis for wells drilled into vertical beads-on-string caves in fracture-caved carbonate reservoirs: Field cases in Shunbei Oilfield

,”

J. Pet. Sci. Eng.

208

,

109280

(

2022

).

https://doi.org/10.1016/j.petrol.2021.109280

Google Scholar

Crossref

63.

C.

Xing

,

H.

Yin

,

H.

Yuan

,

J.

Fu

, and

G.

Xu

, “

Pressure transient analysis for fracture-cavity carbonate reservoirs with large-scale fractures-caves in series connection

,”

J. Energy Resour. Technol.

144

(

5

),

052901

(

2022

).

https://doi.org/10.1115/1.4051739

Google Scholar

Crossref

64.

X.

Li

,

X.

Huang

,

S.

Yu

,

D.

Lu

, and

X.

Du

, “

Numerical study on the transient pressure response of the vug in carbonate reservoirs

,”

Pet. Sci. Technol.

42

,

1

–

20

(

2023

).

Google Scholar

65.

W.

Shi

,

J.

Cheng

,

Y.

Liu

,

M.

Gao

,

L.

Tao

,

J.

Bai

, and

Q.

Zhu

, “

Pressure transient analysis of horizontal wells in multibranched fault-karst carbonate reservoirs: Model and application in SHB oilfield

,”

J. Pet. Sci. Eng.

220

,

111167

(

2023

).

https://doi.org/10.1016/j.petrol.2022.111167

Google Scholar

Crossref

66.

A. A.

Griffith

, “

The phenomena of rupture and flow in solids

,”

Philos. Trans. R. Soc., A

221

,

163

–

198

(

1921

).

https://doi.org/10.1098/rsta.1921.0006

Google Scholar

Crossref

67.

B.

Bourdin

,

C.

Chukwudozie

, and

K.

Yoshioka

, “

A variational approach to the numerical simulation of hydraulic fracturing

,” in

SPE Annual Technical Conference and Exhibition

(

OnePetro

,

2012

).

Google Scholar

Crossref

68.

L. C.

Evans

and

R. F.

Garzepy

,

Measure Theory and Fine Properties of Functions

(

Routledge

,

2018

).

Google Scholar

Crossref

69.

K.

Yoshioka

and

B.

Bourdin

, “

A variational hydraulic fracturing model coupled to a reservoir simulator

,”

Int. J. Rock Mech. Min. Sci.

88

,

137

–

150

(

2016

).

https://doi.org/10.1016/j.ijrmms.2016.07.020

Google Scholar

Crossref

70.

H.

Amor

,

J.-J.

Marigo

, and

C.

Maurini

, “

Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments

,”

J. Mech. Phys. Solids

57

(

8

),

1209

–

1229

(

2009

).

https://doi.org/10.1016/j.jmps.2009.04.011

Google Scholar

Crossref

71.

X.

Li

,

H.

Hofmann

,

K.

Yoshioka

,

Y.

Luo

, and

Y.

Liang

, “

Phase-field modelling of interactions between hydraulic fractures and natural fractures

,”

Rock Mech. Rock Eng.

55

(

10

),

6227

–

6247

(

2022

).

https://doi.org/10.1007/s00603-022-02970-0

Google Scholar

Crossref

72.

A.

Mikelic

and

M. F.

Wheeler

, “

Convergence of iterative coupling for coupled flow and geomechanics

,”

Comput. Geosci.

17

(

3

),

455

–

461

(

2013

).

https://doi.org/10.1007/s10596-012-9318-y

Google Scholar

Crossref

73.

T.

Wick

, “

Modified Newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation

,”

Comput. Methods Appl. Mech. Eng.

325

,

577

–

611

(

2017

).

https://doi.org/10.1016/j.cma.2017.07.026

Google Scholar

Crossref

74.

T.

Wick

, “

An error-oriented Newton/inexact augmented Lagrangian approach for fully monolithic phase-field fracture propagation

,”

SIAM J. Sci. Comput.

39

(

4

),

B589

–

B617

(

2017

).

https://doi.org/10.1137/16M1063873

Google Scholar

Crossref

75.

C.

Miehe

,

M.

Hofacker

, and

F.

Welschinger

, “

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

(

2010

).

https://doi.org/10.1016/j.cma.2010.04.011

Google Scholar

Crossref

76.

T.

Gerasimov

and

L.

De Lorenzis

, “

A line search assisted monolithic approach for phase-field computing of brittle fracture

,”

Comput. Methods Appl. Mech. Eng.

312

,

276

–

303

(

2016

).

https://doi.org/10.1016/j.cma.2015.12.017

Google Scholar

Crossref

77.

A. K.

Zheltov

, “

Formation of vertical fractures by means of highly viscous liquid

,” in

Proceedings of the 4th World Petroleum Congress

,

1955

.

Google Scholar

78.

J.

Geertsma

and

F. D.

Klerk

, “

A rapid method of predicting width and extent of hydraulically induced fractures

,”

J. Pet. Technol.

21

,

1571

–

1969

(

1969

).

https://doi.org/10.2118/2458-PA

Google Scholar

Crossref

79.

D. A.

Spence

and

P.

Sharp

, “

Self-similar solutions for elastohydrodynamic cavity flow

,”

Proc. R. Soc. A

400

(

1819

),

289

–

313

(

1985

).

https://doi.org/10.1098/rspa.1985.0081

Google Scholar

Crossref

80.

C.

Chukwudozie

,

B.

Bourdin

, and

K.

Yoshioka

, “

A variational approach to the modeling and numerical simulation of hydraulic fracturing under in situ stresses

,” in

Proceedings of the 38th Workshop on Geothermal Reservoir Engineering

(

Stanford Geothermal Program Stanford

,

California

,

2013

).

Google Scholar

2023

Author(s)

A novel hydraulic fracturing model for the fluid-driven fracture propagation in poroelastic media containing the natural cave

I. INTRODUCTION

II. MATHEMATICAL MODELS

A. The reservoir flow model

B. The fracture flow model

C. The cave flow model

D. Energy functional

E. Regularization by phase-field method

F. The decomposition of elastic strain energy

G. Minimization of energy functional

III. NUMERICAL IMPLEMENTATION

A. Governing equation discretization

B. Numerical implementation scheme

IV. MODEL VALIDATION

A. Comparison with KGD model

B. Comparison with Cheng's model

C. Sensitivity analysis of the regularization parameter

V. NUMERICAL CASES AND DISCUSSION

A. Effect of the distance between fracture and cave

B. Effect of injection rate

C. Effect of cave radius

D. Effect of in situ stress in three-dimensional example

VI. MULTIPLE CAVE SIMULATIONS

VII. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: THE DISCRETIZING PROCESS OF THE FLOW EQUATIONS

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

A novel hydraulic fracturing model for the fluid-driven fracture propagation in poroelastic media containing the natural cave

I. INTRODUCTION

II. MATHEMATICAL MODELS

A. The reservoir flow model

B. The fracture flow model

C. The cave flow model

D. Energy functional

E. Regularization by phase-field method

F. The decomposition of elastic strain energy

G. Minimization of energy functional

III. NUMERICAL IMPLEMENTATION

A. Governing equation discretization

B. Numerical implementation scheme

IV. MODEL VALIDATION

A. Comparison with KGD model

B. Comparison with Cheng's model

C. Sensitivity analysis of the regularization parameter

V. NUMERICAL CASES AND DISCUSSION

A. Effect of the distance between fracture and cave

B. Effect of injection rate

C. Effect of cave radius

D. Effect of in situ stress in three-dimensional example

VI. MULTIPLE CAVE SIMULATIONS

VII. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: THE DISCRETIZING PROCESS OF THE FLOW EQUATIONS

REFERENCES

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only