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

Shuwei Zhou; Xiaoying Zhuang; and Timon Rabczuk

arXiv:1902.09954·physics.geo-ph·February 27, 2019

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

Shuwei Zhou, Xiaoying Zhuang, and Timon Rabczuk

PDF

TL;DR

This paper introduces a phase-field model for simulating fracture propagation in poroelastic media, integrating Biot's theory with a fluid-interpolating phase field, verified through numerical examples and analytical comparisons.

Contribution

It presents a novel phase-field approach for fracture in poroelastic materials, coupling fluid flow and fracture mechanics within a unified framework.

Findings

01

Model accurately predicts fracture in poroelastic media

02

Successful implementation in Comsol Multiphysics

03

Demonstrates capability with 2D and 3D fluid injection simulations

Abstract

This paper proposes a phase field model for fracture in poroelastic media. The porous medium is modeled based on the classical Biot poroelasticity theory and the fracture behavior is controlled by the phase field model. Moreover, the fracture propagation is driven by the elastic energy where the phase field is used as an interpolation function to transit fluid property from the intact medium to the fully broken one. We use a segregated (staggered) scheme and implement our approach in Comsol Multiphysics. The proposed model is verified by a single-phase solid subjected to tension and a 2D specimen subjected to an increasing internal pressure. We also compare our results with analytical solutions. Finally, we show 2D and 3D examples of internal fluid injection to illustrate the capability of the proposed approach.

Tables2

Table 1. Table 1: Parameters for the single-phasic solid subjected to tension

$μ$	80.77 GPa	$λ$	121.15 GPa	$G_{c}$	2700 N/m
$k$	$1 \times 10^{- 9}$	$l_{0}$	$7.5 \times 10^{- 3}$ mm, $1.5 \times 10^{- 2}$ mm

Table 2. Table 2: Parameters for the notched specimen subjected to internal pressure

$μ$	80.77 GPa	$λ$	121.15 GPa	$G_{c}$	2700 N/m	$k$	$1 \times 10^{- 9}$	$l_{0}$	$4 \times 10^{- 3}$ m
$c_{1}$	0.4	$c_{2}$	1.0	$ε_{p r}$	0.002	$ρ_{r}$ , $ρ_{f}$	$1.0 \times 10^{3}$ kg/m³	$α_{r}$	0.002
$q_{r}$	0	$q_{f}$	0	$K_{r}$	$1 \times 10^{- 15}$ m²	$K_{f}$	$1.333 \times 10^{- 6}$ m²
$c_{r}$	$1 \times 10^{- 8}$ 1/Pa	$c_{f}$	$1 \times 10^{- 8}$ 1/Pa	$μ_{r}$	$1 \times 10^{- 3}$ Pa $\cdot$ s	$μ_{f}$	$1 \times 10^{- 3}$ Pa $\cdot$ s

Equations119

Ψ (u, Γ) = \int_{Ω} ψ_{ε} (ε) d Ω - \int_{Ω} α p \cdot (\nabla \cdot u) d Ω + \int_{Γ} G_{c} d S - \int_{Ω} b \cdot u d Ω - \int_{\partial Ω_{h_{i}}} f \cdot u d S

Ψ (u, Γ) = \int_{Ω} ψ_{ε} (ε) d Ω - \int_{Ω} α p \cdot (\nabla \cdot u) d Ω + \int_{Γ} G_{c} d S - \int_{Ω} b \cdot u d Ω - \int_{\partial Ω_{h_{i}}} f \cdot u d S

ε_{ij} = \frac{1}{2} (\frac{\partial u _{i}}{\partial x _{j}} + \frac{\partial u _{j}}{\partial x _{i}})

ε_{ij} = \frac{1}{2} (\frac{\partial u _{i}}{\partial x _{j}} + \frac{\partial u _{j}}{\partial x _{i}})

ψ_{ε} (ε) = \frac{1}{2} λ ε_{ii} ε_{j j} + μ ε_{ij} ε_{ij}

ψ_{ε} (ε) = \frac{1}{2} λ ε_{ii} ε_{j j} + μ ε_{ij} ε_{ij}

γ (ϕ, ▽ ϕ) = \frac{ϕ ^{2}}{2 l _{0}} + \frac{l _{0}}{2} \frac{\partial ϕ}{\partial x _{i}} \frac{\partial ϕ}{\partial x _{i}}

γ (ϕ, ▽ ϕ) = \frac{ϕ ^{2}}{2 l _{0}} + \frac{l _{0}}{2} \frac{\partial ϕ}{\partial x _{i}} \frac{\partial ϕ}{\partial x _{i}}

\int_{Γ} G_{c} d S = \int_{Ω} G_{c} [\frac{ϕ ^{2}}{2 l _{0}} + \frac{l _{0}}{2} \frac{\partial ϕ}{\partial x _{i}} \frac{\partial ϕ}{\partial x _{i}}] d Ω

\int_{Γ} G_{c} d S = \int_{Ω} G_{c} [\frac{ϕ ^{2}}{2 l _{0}} + \frac{l _{0}}{2} \frac{\partial ϕ}{\partial x _{i}} \frac{\partial ϕ}{\partial x _{i}}] d Ω

⎩ ⎨ ⎧ ε^{+} = a = 1 \sum d ⟨ ε_{a} ⟩^{+} n_{a} \otimes n_{a} ε^{-} = a = 1 \sum d ⟨ ε_{a} ⟩^{-} n_{a} \otimes n_{a}

⎩ ⎨ ⎧ ε^{+} = a = 1 \sum d ⟨ ε_{a} ⟩^{+} n_{a} \otimes n_{a} ε^{-} = a = 1 \sum d ⟨ ε_{a} ⟩^{-} n_{a} \otimes n_{a}

⎩ ⎨ ⎧ ψ_{ε}^{+} (ε) = \frac{λ}{2} ⟨ t r (ε) ⟩^{+ 2} + μ t r (ε^{+ 2}) ψ_{ε}^{-} (ε) = \frac{λ}{2} ⟨ t r (ε) ⟩^{- 2} + μ t r (ε^{- 2})

⎩ ⎨ ⎧ ψ_{ε}^{+} (ε) = \frac{λ}{2} ⟨ t r (ε) ⟩^{+ 2} + μ t r (ε^{+ 2}) ψ_{ε}^{-} (ε) = \frac{λ}{2} ⟨ t r (ε) ⟩^{- 2} + μ t r (ε^{- 2})

ψ_{ε} (ε) = [(1 - k) (1 - ϕ)^{2} + k] ψ_{ε}^{+} (ε) + ψ_{ε}^{-} (ε)

ψ_{ε} (ε) = [(1 - k) (1 - ϕ)^{2} + k] ψ_{ε}^{+} (ε) + ψ_{ε}^{-} (ε)

L = \int_{Ω} {[(1 - k) (1 - ϕ)^{2} + k] ψ_{ε}^{+} (ε) + ψ_{ε}^{-} (ε)} d Ω - \int_{Ω} α p \cdot (\nabla \cdot u) d Ω + \int_{Ω} G_{c} [\frac{ϕ ^{2}}{2 l _{0}} + \frac{l _{0}}{2} \frac{\partial ϕ}{\partial x _{i}} \frac{\partial ϕ}{\partial x _{i}}] d Ω - \int_{Ω} b_{i} u_{i} d Ω - \int_{\partial Ω_{h_{i}}} f_{i} u_{i} d S

L = \int_{Ω} {[(1 - k) (1 - ϕ)^{2} + k] ψ_{ε}^{+} (ε) + ψ_{ε}^{-} (ε)} d Ω - \int_{Ω} α p \cdot (\nabla \cdot u) d Ω + \int_{Ω} G_{c} [\frac{ϕ ^{2}}{2 l _{0}} + \frac{l _{0}}{2} \frac{\partial ϕ}{\partial x _{i}} \frac{\partial ϕ}{\partial x _{i}}] d Ω - \int_{Ω} b_{i} u_{i} d Ω - \int_{\partial Ω_{h_{i}}} f_{i} u_{i} d S

δ L = ① \int_{\partial Ω_{h_{i}}} [(σ_{ij} - α p δ_{ij}) m_{j} - f_{i}] δ u_{i} d S - ② \int_{Ω} [(σ_{ij} - α p δ_{ij})_{, j} + b_{i}] δ u_{i} d Ω - ③ \int_{Ω} [2 (ϕ - 1) (1 - k) ψ_{ε}^{+} + \frac{G _{c} ϕ}{l _{0}} - G_{c} l_{0} \frac{\partial ^{2} ϕ}{\partial x _{i}^{2}}] δ ϕ d Ω + ④ \int_{\partial Ω} (\frac{ϕ}{x _{i}} m_{i}) δ ϕ d S = 0

δ L = ① \int_{\partial Ω_{h_{i}}} [(σ_{ij} - α p δ_{ij}) m_{j} - f_{i}] δ u_{i} d S - ② \int_{Ω} [(σ_{ij} - α p δ_{ij})_{, j} + b_{i}] δ u_{i} d Ω - ③ \int_{Ω} [2 (ϕ - 1) (1 - k) ψ_{ε}^{+} + \frac{G _{c} ϕ}{l _{0}} - G_{c} l_{0} \frac{\partial ^{2} ϕ}{\partial x _{i}^{2}}] δ ϕ d Ω + ④ \int_{\partial Ω} (\frac{ϕ}{x _{i}} m_{i}) δ ϕ d S = 0

σ_{ij} = [(1 - k) (1 - ϕ)^{2} + k] \frac{\partial ψ _{ε}^{+}}{\partial ε _{ij}} + \frac{\partial ψ _{ε}^{-}}{\partial ε _{ij}}

σ_{ij} = [(1 - k) (1 - ϕ)^{2} + k] \frac{\partial ψ _{ε}^{+}}{\partial ε _{ij}} + \frac{\partial ψ _{ε}^{-}}{\partial ε _{ij}}

σ = [(1 - k) (1 - ϕ)^{2} + k] [λ ⟨ t r (ε) ⟩^{+} I + 2 μ ε^{+}] + λ ⟨ t r (ε) ⟩^{-} I + 2 μ ε^{-}

σ = [(1 - k) (1 - ϕ)^{2} + k] [λ ⟨ t r (ε) ⟩^{+} I + 2 μ ε^{+}] + λ ⟨ t r (ε) ⟩^{-} I + 2 μ ε^{-}

σ^{p or} (ε) = σ (ε) - α p I, in Ω \times (0, T]

σ^{p or} (ε) = σ (ε) - α p I, in Ω \times (0, T]

⎩ ⎨ ⎧ \frac{\partial σ _{ij}^{p or}}{\partial x _{j}} + b_{i} = 0, [\frac{2 l _{0} ( 1 - k ) ψ _{ε}^{+}}{G _{c}} + 1] ϕ - l_{0}^{2} \frac{\partial ^{2} ϕ}{\partial x _{i}^{2}} = \frac{2 l _{0} ( 1 - k ) ψ _{ε}^{+}}{G _{c}}, in Ω \times (0, T] in Ω \times (0, T]

⎩ ⎨ ⎧ \frac{\partial σ _{ij}^{p or}}{\partial x _{j}} + b_{i} = 0, [\frac{2 l _{0} ( 1 - k ) ψ _{ε}^{+}}{G _{c}} + 1] ϕ - l_{0}^{2} \frac{\partial ^{2} ϕ}{\partial x _{i}^{2}} = \frac{2 l _{0} ( 1 - k ) ψ _{ε}^{+}}{G _{c}}, in Ω \times (0, T] in Ω \times (0, T]

H (x, t) = s \in [0, t] max ψ_{ε}^{+} (ε (x, s)), in Ω \times (0, T]

H (x, t) = s \in [0, t] max ψ_{ε}^{+} (ε (x, s)), in Ω \times (0, T]

ψ_{ε}^{+} - H \leq 0, \dot{H} \geq 0, \dot{H} (ψ_{ε}^{+} - H) = 0

ψ_{ε}^{+} - H \leq 0, \dot{H} \geq 0, \dot{H} (ψ_{ε}^{+} - H) = 0

⎩ ⎨ ⎧ \frac{\partial σ _{ij}^{p or}}{\partial x _{i}} + b_{i} = 0, [\frac{2 l _{0} ( 1 - k ) H}{G _{c}} + 1] ϕ - l_{0}^{2} \frac{\partial ^{2} ϕ}{\partial x _{i}^{2}} = \frac{2 l _{0} ( 1 - k ) H}{G _{c}}, in Ω \times (0, T] in Ω \times (0, T]

⎩ ⎨ ⎧ \frac{\partial σ _{ij}^{p or}}{\partial x _{i}} + b_{i} = 0, [\frac{2 l _{0} ( 1 - k ) H}{G _{c}} + 1] ϕ - l_{0}^{2} \frac{\partial ^{2} ϕ}{\partial x _{i}^{2}} = \frac{2 l _{0} ( 1 - k ) H}{G _{c}}, in Ω \times (0, T] in Ω \times (0, T]

⎩ ⎨ ⎧ σ_{ij}^{p or} m_{j} = f_{i}, \frac{\partial ϕ}{\partial x _{i}} m_{i} = 0, o n \partial Ω_{h_{i}} \times (0, T] o n \partial Ω \times (0, T]

⎩ ⎨ ⎧ σ_{ij}^{p or} m_{j} = f_{i}, \frac{\partial ϕ}{\partial x _{i}} m_{i} = 0, o n \partial Ω_{h_{i}} \times (0, T] o n \partial Ω \times (0, T]

\frac{\partial}{\partial t} (ε_{p r} ρ_{r}) + \nabla \cdot (ρ_{r} v_{r}) = q_{r} - ρ_{r} α_{r} \frac{\partial ε _{v o l}}{\partial t}

\frac{\partial}{\partial t} (ε_{p r} ρ_{r}) + \nabla \cdot (ρ_{r} v_{r}) = q_{r} - ρ_{r} α_{r} \frac{\partial ε _{v o l}}{\partial t}

v_{r} = - \frac{K _{r}}{μ _{r}} (\nabla p + ρ_{r} g)

v_{r} = - \frac{K _{r}}{μ _{r}} (\nabla p + ρ_{r} g)

\frac{\partial}{\partial t} (ε_{p r} ρ_{r}) = ρ_{r} S_{r} \frac{\partial p}{\partial t}

\frac{\partial}{\partial t} (ε_{p r} ρ_{r}) = ρ_{r} S_{r} \frac{\partial p}{\partial t}

S_{r} = ε_{p r} c_{r} + \frac{( α _{r} - ε _{p r} ) ( 1 - α _{r} )}{K _{V r}}

S_{r} = ε_{p r} c_{r} + \frac{( α _{r} - ε _{p r} ) ( 1 - α _{r} )}{K _{V r}}

ρ_{r} S_{r} \frac{\partial p}{\partial t} + \nabla \cdot (ρ_{r} v_{r}) = q_{r} - ρ_{r} α_{r} \frac{\partial ε _{v o l}}{\partial t}

ρ_{r} S_{r} \frac{\partial p}{\partial t} + \nabla \cdot (ρ_{r} v_{r}) = q_{r} - ρ_{r} α_{r} \frac{\partial ε _{v o l}}{\partial t}

ρ_{f} S_{f} \frac{\partial p}{\partial t} + \nabla \cdot (ρ_{f} v_{f}) = q_{f}

ρ_{f} S_{f} \frac{\partial p}{\partial t} + \nabla \cdot (ρ_{f} v_{f}) = q_{f}

v_{f} = - \frac{K _{f}}{μ _{f}} (\nabla p + ρ_{f} g)

v_{f} = - \frac{K _{f}}{μ _{f}} (\nabla p + ρ_{f} g)

χ_{r} (\cdot, ϕ) := χ_{r} (x, t, ϕ) = 1 in Ω_{r} (t), an d χ_{r} (\cdot, ϕ) = 0 in Ω_{f} (t)

χ_{r} (\cdot, ϕ) := χ_{r} (x, t, ϕ) = 1 in Ω_{r} (t), an d χ_{r} (\cdot, ϕ) = 0 in Ω_{f} (t)

χ_{f} (\cdot, ϕ) := χ_{f} (x, t, ϕ) = 1 in Ω_{f} (t), an d χ_{f} (\cdot, ϕ) = 0 in Ω_{r} (t)

χ_{f} (\cdot, ϕ) := χ_{f} (x, t, ϕ) = 1 in Ω_{f} (t), an d χ_{f} (\cdot, ϕ) = 0 in Ω_{r} (t)

χ_{r} (\cdot, ϕ) = \frac{c _{2} - ϕ}{c _{2} - c _{1}}, an d χ_{f} (\cdot, ϕ) = \frac{ϕ - c _{1}}{c _{2} - c _{1}}

χ_{r} (\cdot, ϕ) = \frac{c _{2} - ϕ}{c _{2} - c _{1}}, an d χ_{f} (\cdot, ϕ) = \frac{ϕ - c _{1}}{c _{2} - c _{1}}

ρS \frac{\partial p}{\partial t} + \nabla \cdot (ρ v) = q_{m} - ρ α χ_{r} \frac{\partial ε _{v o l}}{\partial t}

ρS \frac{\partial p}{\partial t} + \nabla \cdot (ρ v) = q_{m} - ρ α χ_{r} \frac{\partial ε _{v o l}}{\partial t}

S = ε_{p} c + \frac{( α - ε _{p} ) ( 1 - α )}{K _{V r}}

S = ε_{p} c + \frac{( α - ε _{p} ) ( 1 - α )}{K _{V r}}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

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

Shuwei Zhou1,2, Xiaoying Zhuang2, Timon Rabczuk3,4∗

Abstract

This paper proposes a phase field model for fracture in poroelastic media. The porous medium is modeled based on the classical Biot poroelasticity theory and the fracture behavior is controlled by the phase field model. Moreover, the fracture propagation is driven by the elastic energy where the phase field is used as an interpolation function to transit fluid property from the intact medium to the fully broken one. We use a segregated (staggered) scheme and implement our approach in Comsol Multiphysics. The proposed model is verified by a single-phase solid subjected to tension and a 2D specimen subjected to an increasing internal pressure. We also compare our results with analytical solutions. Finally, we show 2D and 3D examples of internal fluid injection to illustrate the capability of the proposed approach.

1 Institute of Structural Mechanics, Bauhaus-University Weimar, Weimar 99423, Germany

2 Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, P.R. China

3 Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam

4 Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam

Corresponding author: Timon Rabczuk

Keywords: Phase field, Hydraulic fractures, Poroelasticity, Comsol

1 Introduction

Fracture in poroelastic media is of major importance in mechanical, energy and environmental engineering [1, 2]. In particular, predicting fracture propagation during water injection is critical for hydraulic fracturing (fracking), an interesting technique used to extract petroleum and natural gas (e.g. shale gas). Fractures are created by fracking to connect wellbores with expected petroleum or gas. Thus, vast amounts of resources, which are inaccessible in previous decades, will be extracted promising great economic benefits. However, fracking may also result in leakage of fracturing fluid or gas, which might contaminate the groundwater because of the unintentionally created channels [3, 4]. Other reasons for objecting fracking are risks to air quality loss and surface contamination from spills [1]. Therefore, the vast economic benefits and continued controversies require accurate mathematical models and proper numerical simulation tools to treat complex problems on fracture propagation in poroelastic media [5].

Fracture models for single-phasic solids can be classified into two categories, i.e. discrete and continuous approaches. Discrete approaches introduce discontinuities which capture the jump in the displacement field. Among the most popular computational methods for fracture are the extended finite element method (XFEM) [6], cohesive elements, element-erosion techniques [7] and remeshing techniques [8]. Discrete crack methods are usually computationally challenging and require complex crack tracking procedures. Furthermore, the lack of reliable cracking criteria reduces the robustness and reliability of discrete crack propagation, at least when they are based on a continuous crack topology. The cracking particle method [9], peridynamics [10] and dual-horizon peridynamics [11] are discrete crack approaches which do not exhibit such drawbacks and are therefore well suited for dynamic fracture including complex crack interactions. They do not model the crack as continuous surface which also significantly facilitates the implementation. However, fluid flow through the opening of discrete cracks cannot be captured in such method due to the absence of crack path continuity. There are a few contributions modeling fracture through discrete cracks in extended finite element and meshfree methods, see e.g. the contributions in [12]. However, they are commonly restricted to a few number of cracks. Different from the discrete approaches are continuous approaches to fracture which do not introduce discontinuities in the displacement field. Popular approaches include gradient damage models [13], screened-poisson models [14] and phase field models [15]. All these models introduce an intrinsic length scale into the discretization and smear the fracture over a localization band of finite width drastically facilitating their implementation.

To date, there are a few approaches of discrete and continuous approaches for hydraulic fracturing (HF). For example, Wu et al. [16] modeled 2D crack propagation using the BEM. Lecampion [17] applied the extended finite element method to hydraulic fracturing. The fluid pressure was applied along a line fracture in a 2D impermeable medium to drive the crack. However, the complexity of discrete approaches for HF remains high and studies are limited to fairly simple cases, mostly in 2D. PFMs seem to be a valuable alternative, especially since they can be easily extended to coupled problems. PFMs for HF have been proposed for instance in [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Wheeler et al. [19] rewrote the energy functional by including poroelastic terms and succeeded in extending the phase field model to porous media. However, the variation in the reservoir and fracture domains with time is treated as a moving boundary problem and thereby extra work is needed for implementation. Later, Mikelić et al. [21] fully coupled the three fields: elasticity, phase field, and pressure. They modified the energy functional from their previous work, the flow in their entire system is governed by Biot equations. The permeability tensor was also modified to consider a higher permeability along the fracture. The implementation approaches of Mikelić et al. [21] were then enhanced using adaptive element schemes by [23]. Yoshioka and Bourdin [26] proposed approaches to couple the phase field model to reservoir simulators. Miehe et al. [20, 27] proposed new minimization and saddle point principles for Darcy-Biot-type flow in fractured poroelastic media coupled with phase field modeling. The evolution of the phase field was driven by the effective stress in the solid skeleton and a stress threshold was set. Moreover, the flow in the fractures was set as Poiseuille-type by modeling Darcy flow with an anisotropic permeability tensor. Recently, Ehlers and Luo [28] embedded a phase-field approach in the theory of porous media to model dynamic hydraulic fracturing. Santillán et al. [29] proposed an immersed-fracture formulation for impermeable porous media. However, these recently developed approaches for extending PFM to HF are quite complicated and computationally expensive.

In this paper another attempt for phase field modeling of fluid-driven fractures is made. The porous medium is modeled based on the classical Biot poroelasticity theory and the fracture behavior is controlled by the phase field model. We add an additional pressure-related term in the energy functional and derive the governing equations of the strong form. Different from Miehe et al. [20, 27], the fracture propagation is driven by the elastic energy and we set no stress threshold in the formulation. In addition, the phase field is used as an interpolation function to transit fluid property from the intact medium to the fully broken one. We use the FEM software Comsol Multiphysics to implement the proposed approach. A segregated (staggered) scheme is used and thereby the displacement, pressure, and phase field are solved independently. We verify our approach for two cases: a single-phasic solid and a 2D specimen subjected to an increasing internal pressure. In the second case, the proposed approach is compared with analytical solutions. Subsequently, more complex examples are studied.

The paper is organized as follows. We give the mathematical models for fracture in Section 2 followed by Section 3 which shows the numerical implementation of the phase field model in Comsol. In Section 4, we then verify the numerical simulation by two cases: a single-phasic solid subjected to tension and a 2D specimen subjected to an increasing internal pressure. Afterwards, we give 2D and 3D examples of internal fluid injection in Section 5. Finally, we end with conclusions regarding our work in Section 6.

2 Mathematical models for fracture in porous media

2.1 Theory of brittle fracture

Let us consider a bounded domain $\Omega\subset\mathbb{R}^{d}$ ( $d\in\{2,3\}$ ) with external boundary $\partial\Omega$ and an internal discontinuity boundary $\Gamma$ . Furthermore, we assume:

•

The intrinsic length scale parameter of the phase field is larger than the pore size.

•

The porous medium is linear elastic, homogeneous, and isotropic.

•

The fluid in the medium is compressible and viscous.

Let $T>0$ be the computational time interval and the displacement of body $\Omega$ at time $t\in[0,T]$ is denoted by $\bm{u}(\bm{x},t)\subset\mathbb{R}^{d}$ where $\bm{x}$ is the position vector. The displacement field satisfies the time-dependent Dirichlet boundary conditions, $u_{i}(\bm{x},t)=g_{i}(\bm{x},t)$ , on $\partial\Omega_{g_{i}}\in\Omega$ , and the time-dependent Neumann conditions on $\partial\Omega_{h_{i}}\in\Omega$ . The body $\Omega$ is subjected to a body force $\bm{b}(\bm{x},t)\subset\mathbb{R}^{d}$ and a traction $\bm{f}(\bm{x},t)$ acts on the the boundary $\partial\Omega_{h_{i}}$ .

Griffith’s theory [30] assumes the energy to create a fracture surface per unit area is equal to the critical fracture energy density $G_{c}$ , which is also commonly referred to as the critical energy release rate. Thus, for a pure single-phasic solid, the energy functional for the solid $\Psi(\bm{u},\Gamma)$ can be expressed in terms of the elastic energy $\psi_{\varepsilon}(\bm{\varepsilon})$ , fracture energy and work from external forces. However, for a porous solid, the fluid pressure must be considered for the solid and an additional pressure-related term is used [21, 22, 24]. In this paper, we follow the idea of Lee et al. [24] and rewrite the energy functional as

[TABLE]

where $p:\Omega\times[0,T]\rightarrow\mathbb{R}$ is the fluid pressure, $\alpha\in[\epsilon_{p},1]$ is the Biot coefficient and $\epsilon_{p}$ and

[TABLE]

is the linear strain tensor. For an isotropic and linear elastic solid, the elastic energy density $\psi_{\varepsilon}(\bm{\varepsilon})$ is given by [31]

[TABLE]

where $\lambda,\mu>0$ are the Lamé constants.

2.2 Phase field approximation for fracture

The phase field $\phi(\bm{x},t)\in[0,1]$ diffuses the fracture surface $\Gamma$ over a volume which avoids complex crack tracking procedures. $\phi=1$ indicates the ’fully’ cracked region while $\phi=0$ denotes the unbroken material. A length scale parameter $l_{0}$ is used to control the transition region where the phase field varies from 0 to 1, thus, reflecting the shape of a diffusive crack (see Fig. 1 for details). The width of the crack region will increase with increasing $l_{0}$ . It can be shown that the crack surface density per unit volume of the solid is given by [31]

[TABLE]

By using Eq. (4), the fracture energy in Eq. (1) is calculated by

[TABLE]

The variational approach for fracture [30] states that the elastic energy drives the evolution of the phase field. To ensure cracks only under tension, the elastic energy is decomposed into tensile and compressive parts [32]. We follow the decomposition of Miehe et al. [31] where the strain tensor $\bm{\varepsilon}$ is decomposed as follows:

[TABLE]

where $\bm{\varepsilon}^{+}$ and $\bm{\varepsilon}^{-}$ are the tensile and compressive strain tensors, respectively; $\varepsilon_{a}$ and $\bm{n}_{a}$ are the principal strains and their directions. The operators $\langle\centerdot\rangle^{+}$ and $\langle\centerdot\rangle^{-}$ are defined as $\langle\centerdot\rangle^{+}=(\centerdot+|\centerdot|)/2$ , and $\langle\centerdot\rangle^{-}=(\centerdot-|\centerdot|)/2$ . Using the decomposed strain tensor, the tensile and compressive parts of the elastic energy density are given by

[TABLE]

We follow Borden et al. [32] and assume that the phase field affects only the tensile part of the elastic energy density. Next, the stiffness reduction is modeled by the following equation:

[TABLE]

where $0<k\ll 1$ is a model parameter that prevents the tensile part of the elastic energy density from disappearing and avoids numerical singularity when the phase field $\phi$ tends to 1.

2.3 Governing equations for evolution of the phase field

Based on the phase field approximation for the fracture energy (5) and decomposition of the elastic energy (8), the energy functional (1) can be rewritten as

[TABLE]

The variational approach [30] stresses that initiation, propagation and branching of the crack $\Gamma(\bm{x},t)$ at time $t\in[0,T]$ for ${\bm{x}}\in\Omega$ occur when the potential achieves a minimum value. Hence, we calculate the first variation of the functional $L$ and set it zero, yielding

[TABLE]

where $m_{j}$ is the component of the outward-pointing normal vector of the boundary and $\sigma_{ij}$ is the component of the effective trial stress tensor $\bm{\sigma}(\bm{\varepsilon})$ given by:

[TABLE]

where $\bm{I}$ is the identity tensor $\in\mathbb{R}^{d\times d}$ .

We recall the Cauchy stress tensor $\bm{\sigma}^{por}$ :

[TABLE]

Because of the arbitrariness of the variation $\delta u$ and $\delta\phi$ , Eq. (10) $②$ and $③$ give rise to the governing equations:

[TABLE]

The variational approach [30] also emphasizes that the irreversibility condition $\Gamma(\bm{x},s)\in\Gamma(\bm{x},t)(s<t)$ must be satisfied, which ensures that a crack cannot heal. To ensure a monotonically increasing phase field, the irreversibility condition during compression or unloading is accomplished by introducing a strain-history field $H(\bm{x},t)$ [31, 15, 32] defined by

[TABLE]

The history field $H$ satisfies the Kuhn-Tucker conditions for loading and unloading

[TABLE]

More details can be found in [31]. Replacing $\psi_{\varepsilon}^{+}$ by $H(\bm{x},t)$ in Eq. (14), the strong form is obtained by

[TABLE]

In addition, Eq. (10) $①$ and $④$ give rise to the Neumann conditions of the displacement and phase field,

[TABLE]

2.4 Governing equations for the flow field

To formulate the flow equations in the domain $\Omega$ , we divide the whole domain into three parts: $\Omega_{r}(t)$ , $\Omega_{f}(t)$ and $\Omega_{t}(t)$ . $\Omega_{r}(t)$ represents the unbroken domain (reservoir domain) and $\Omega_{f}(t)$ is the fracture domain. $\Omega_{t}(t)$ is the transition domain betwen $\Omega_{r}(t)$ and $\Omega_{f}(t)$ . Here, we follow Lee et al. [24] and use the phase field as an indicator function. Two thresholds $c_{1}$ and $c_{2}$ are set to separate the three flow domains. We consider a subdomain as the fracture domain $\Omega_{f}(t)$ if $\phi\geq c_{2}$ and as the reservoir domain $\Omega_{r}(t)$ if $\phi\leq c_{1}$ . In the transition domain, $c_{1}<\phi<c_{2}$ .

We use Darcy’s law to describe the flow field in the poroelastic medium. For the reservoir domain $\Omega_{r}(t)$ , the mass conservation leads to

[TABLE]

where $\rho_{r}$ , $q_{r}$ , $\varepsilon_{pr}$ , and $\alpha_{r}$ are the density of fluid, source term, porosity, and Biot coefficient in the reservoir domain, respectively; $\varepsilon_{vol}=\nabla\cdot\bm{u}$ is the volumetric strain of $\Omega_{r}(t)$ .

Darcy’s velocity in $\Omega_{r}(t)$ is given by

[TABLE]

where $K_{r}$ and $\mu_{r}$ are the permeability and fluid viscosity of $\Omega_{r}(t)$ , respectively; $\bm{g}$ is the gravity.

Based on the storage model in [33], we have

[TABLE]

where $S_{r}$ , the storage coefficient of $\Omega_{r}$ , is given by

[TABLE]

with $c_{r}$ the fluid compressibility and $K_{Vr}$ the bulk modulus of the reservoir domain. Thus, the equation of mass conservation (19) translates to

[TABLE]

For the fracture domain $\Omega_{f}(t)$ , the volumetric strain vanishes from the equation of mass conservation:

[TABLE]

where $\rho_{f}$ , $S_{f}$ , and $q_{f}$ are the fluid density, storage coefficient, and source term in the fracture domain $\Omega_{f}(t)$ .

The storage coefficient $S_{f}$ is equal to the fluid compressibility $c_{f}$ and the Darcy’s velocity in $\Omega_{f}(t)$ is given by

[TABLE]

with $K_{f}$ and $\mu_{f}$ the permeability and fluid viscosity of $\Omega_{f}(t)$ .

We follow Lee et al. [24] and define two linear indicator functions to connect the transition domain with the reservoir and fracture domains: $\chi_{r}$ and $\chi_{f}$ , which satisfies

[TABLE]

In the transition domain, the two linear functions are defined as

[TABLE]

Figure 2 shows more details of the linear indicator functions $\chi_{r}$ and $\chi_{f}$ as well as illustration of the fracture and reservoir domains. The fluid and solid properties of the transition domain $\Omega_{t}$ are obtained by interpolation of the reservoir and fracture domains with the indicator functions $\chi_{r}$ and $\chi_{f}$ . Thus, the mass conservation in the transition domain is rewritten as

[TABLE]

where $\rho=\rho_{r}\chi_{r}+\rho_{f}\chi_{f}$ , $\alpha=\alpha_{r}\chi_{r}+\alpha_{f}\chi_{f}$ , $q_{m}$ is the source term, and the storage coefficient $S$ is represented by

[TABLE]

with $c=c_{r}\chi_{r}+c_{f}\chi_{f}$ . Note that $\varepsilon_{p}=0$ and $\alpha=1$ for the fracture domain and thereby $\varepsilon_{p}=\varepsilon_{pr}\chi_{r}$ and $\alpha=\alpha_{r}\chi_{r}+\chi_{f}$ .

Eq. (29) is also applicable for the whole domain $\Omega$ because it can be transformed into Eq. (23) for $\Omega_{r}$ and Eq. (24) for $\Omega_{f}$ . The Darcy’s velocity $\bm{v}$ is then calculated by

[TABLE]

where $K=K_{r}\chi_{r}+K_{f}\chi_{f}$ is the effective permeability and $\mu=\mu_{r}\chi_{r}+\mu_{f}\chi_{f}$ is the effective fluid viscosity. Finally, the governing equation for the flow field in the porous domain $\Omega$ is rewritten in terms of the fluid pressure $p$ :

[TABLE]

Note that the determination of the coefficients in the phase-field model will be tackled in future research. Within Comsol, it is easily possible to include an indirect coupling such as the dependence of the permeability of the fracture domain on the crack opening. Relevant contributions to such an indirect coupling are for instance proposed in [34, 35, 36, 37]. In this context, we also mention the experimental work by Yeo et al. [34] who explored the effect of shear deformation on the aperture and permeability of rock fracture. However, the crack opening cannot be extracted in a straight forward manner in the PFM because the sharp crack is represented by a smooth phase field. Therefore, for simplicity, we use an unchanged fluid property in the fracture domain in our presented examples, which also shows favorable results.

2.5 Initial and boundary conditions

For time-transient problems, the following initial conditions must be fulfilled:

[TABLE]

Here, the initial phase field $\phi_{0}=1$ in a local domain can be used to model a pre-existing crack [32]. Boundary conditions for the displacement and phase field are introduced in Subsections 2.1 and 2.3. For the fluid pressure, the Dirichlet condition on $\partial\Omega_{D}$ and Neumann condition on $\partial\Omega_{N}$ with $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$ are prescribed:

[TABLE]

with $p_{D}$ the prescribed pressure on the Dirichlet boundary and $M_{N}$ the mass flux on the Newman boundary.

3 Numerical implementation

3.1 Finite element discretization

The proposed phase field approach is applied in the framework of finite element method. Thus, we first derive the weak form of all the governing equations as

[TABLE]

,

[TABLE]

and

[TABLE]

The nodal values for the three fields ( $\bm{u}$ , $\phi$ , and $p$ ) are defined as $\bm{u}_{i}$ , $\phi_{i}$ , and $p_{i}$ . Afterwards, the fields are discretized by the standard 3D notation as

[TABLE]

where $n$ is the total number of nodes in each element and $N_{i}$ is the shape function of node $i$ . The gradients of the three fields are thereby given by

[TABLE]

where $\bm{B}_{i}^{u}$ , $\bm{B}_{i}^{\phi}$ , and $\bm{B}_{i}^{p}$ are derivatives of the shape functions:

[TABLE]

For 2D, the above equations (41) can be simplified after a little adjustment. Due to the arbitrariness of the test functions, the external force $\bm{F}_{i}^{u,ext}$ and inner force $\bm{F}_{i}^{u,int}$ of the mechanical field are described by

[TABLE]

The inner force term of the phase field are also obtained by

[TABLE]

Finally, when the gravity is neglected, the inner force $F_{i}^{p,int}$ , viscous force $F_{i}^{p,vis}$ , and external force $F_{i}^{p,ext}$ of the pressure field is given by

[TABLE]

Thus, according to Eqs. (36), (37) and (38), contribution of node $i$ to the residual of the discrete equations of stress equilibrium, evolution of phase field, and mass conservation of pressure field is given as

[TABLE]

We use the segregated scheme to solve the displacement, phase field and fluid pressure. The segregated scheme ensures that the energy functional of each field remains convex. The Newton-Raphson approach is adopted to obtain $\bm{R}_{i}^{u}=0$ , $R_{i}^{\phi}=0$ , and $R_{i}^{p}=0$ , respectively. The tangents on the element level are thereby calculated by

[TABLE]

where $\bm{D}$ is the elasticity tensor of four order given by $\bm{D}={\partial\bm{\sigma}}/{\partial\bm{\varepsilon}}$ .

3.2 Comsol implementation

The phase field model is implemented into Comsol. We establish three main modules first: Poroelasticity Module, History-strain Module and Phase Field Module. The Poroelasticity Module is used to solve the displacement field and fluid pressure. The History-strain Module and Phase Field Module are established to solve the other two fields $H$ and $\phi$ . These modules are all written in strong forms and solved based on standard finite element discretization in space domain and finite difference discretization in time domain. The implementation in Comsol needs many intermediate variables, such as the elastic energy, principal strain and direction vectors, and thereby we also establish a pre-set Storage Module to evaluate and store these intermediate variables during each time step.

3.3 Module setup

The Poroelasticity Module is set up based on a linear elastic material library and a transient formulation of Darcy’s law. The boundary and initial conditions in the Poroelasticity Module are implemented as shown in Section 2. However, the elasticity matrix in a time step requires modification. The elasticity matrix comes from the elasticity tensor of fourth order $\bm{D}$ given by

[TABLE]

where $H_{\varepsilon}\langle x\rangle$ is a Heaviside function: $H_{\varepsilon}\langle x\rangle=1$ if $x>0$ and $H_{\varepsilon}\langle x\rangle=0$ if $x\leq 0$ , and $J_{ijkl}=\delta_{ij}\delta_{kl}$ with $\delta_{ij}$ and $\delta_{kl}$ the Kronecker symbols. Finally, the stiffness matrix $\bm{G}$ is rewritten in a standard material data ordering:

[TABLE]

with $G_{ijkl}=\bar{G}_{ijkl}+\tilde{G}_{ijkl}$ .

$\bar{G}_{ijkl}$ is related to trace of the strain $tr(\bm{\varepsilon})$ :

[TABLE]

Based on the algorithm for fourth-order isotropic tensor [38], $\tilde{G}_{ijkl}$ is expressed as

[TABLE]

where

[TABLE]

and

[TABLE]

with $n_{ai}$ the $i$ -th component of vector $\bm{n}_{a}$ .

Because of the fractional form given in Eqs. (51) and (52), Eq. (50) cannot be calculated if $\varepsilon_{a}=\varepsilon_{b}$ . To better apply Eq.(50) in Comsol, we therefore refer to [39] and use a “perturbation” for the principal strains:

[TABLE]

with the perturbation $\delta=1\times 10^{-9}$ for this paper. The second principal strain and volumetric strain remain unchanged.

We establish the Phase Field Module by revising the Helmholtz equation in Comsol. Coefficients of the Helmholtz equation must have the same form as the governing equation (17). The boundary condition in Eq. (18) and initial condition (33) are also implemented in this module.

Comsol has stable Distributed ODEs and DAEs Interface to solve distributed ordinary differential equations (ODE) and differential-algebraic equations (DAE). Thus, we use the Distributed ODEs and DAEs Interfaces to establish the History-strain Module, in which we do not really solve the history-strain field but call a “previous solution” function of Comsol solver to record and update $H(\bm{x},t)$ . The format of the equations written into the History-strain Module is as follows,

[TABLE]

The History-strain Module also needs initial conditions, commonly, $H_{0}(\bm{x})=0$ . However, pre-existing cracks can be generated by introducing the following initial conditions [32]:

[TABLE]

When $d=0$ , substituting $H_{0}$ into (17)2 will give rise to

[TABLE]

If $\phi$ becomes quite close to 1, $B$ will become quite large. Therefore, we follow [32] and adopt $B=1\times 10^{3}$ to artificially generate the initial cracks if necessary.

3.4 Segregated scheme

Figure 3 shows all the established modules and their relationship. The Storage Module stores the principal strains and the direction of principal strain from the Poroelasticity Module. Some intermediate variables such as the elastic energy and the component $G_{ijkl}$ are also calculated in the Storage Module. The History-strain module then uses the positive part of the elastic energy $\psi_{\varepsilon}^{+}$ calculated from the Storage Module to solve and updated the local history strain field. Afterwards, the Phase Field Module utilizes the updated history strain $H$ to solve the phase field. Because of the variation in the phase field, the mechanical and hydraulic properties of the domain need to be updated. Therefore, the stiffness matrix of the Poroelasticity Module is revised according to the previously stored intermediate variables and the updated phase field. The three types of domain $\Omega_{R}(t)$ , $\Omega_{F}(t)$ , and $\Omega_{T}(t)$ are also distinguished and assigned the updated hydraulic property. Finally, the displacement and pressure in the Poroelasticity Module are coupled to be solved. Because Fig. 3 shows strong coupling of all the established modules, we use a segregated scheme to solve the coupled system of equations as indicated in Fig. 4. The segregated scheme will make the simulation much easier to converge.

The segregated scheme means that all the fields are solved independently in one time step. We place the displacement and pressure in one segregated step and solve them together. Afterwards, the history strain and phase field are in another two segregated steps. Thus, the equations of poroelasticity (displacement and pressure), history strain and phase field are solved independently. The implicit Generalized- $\alpha$ method [32] is used to ensure unconditional stability for the calculation. When the time reaches $t_{i}$ , the linear extrapolation of the previous solution provides the initial guess for the three segregated steps ( $(\bm{u},p)_{i}^{j=0}$ , $H_{i}^{j=0}$ and $\phi_{i}^{j=0}$ ). Then, the Newton-Raphson method is used to solve each segregated step. For a new iteration step $j+1$ in the given time step $i$ , the displacement and pressure $(\bm{u},p)_{i}^{j+1}$ is first solved by using the results from previous iteration step $j$ ( $(\bm{u},p)_{i}^{j}$ , $H_{i}^{j}$ and $\phi_{i}^{j}$ ) after a Newton-Raphson iteration. Based on the updated displacement $\bm{u}_{i}^{j+1}$ , the history strain field is then updated. The Newton-Raphson iteration is subsequently applied to achieve a new phase field $\phi_{i}^{j+1}$ by $\bm{u}_{i}^{j+1}$ , $H_{i}^{j+1}$ . Finally, the total relative error of the solutions from the previous and present iteration steps $\varepsilon_{r}$ is calculated. If $\varepsilon_{r}$ is less than the tolerance $\varepsilon_{t}$ , the calculation for the time step $i$ is finished and the time will increase to $t_{i+1}$ . If $\varepsilon_{r}\geq\varepsilon_{t}$ , a new iteration step will be used until $\varepsilon_{r}<\varepsilon_{t}$ . The tolerance $\varepsilon_{t}=1\times 10^{-4}$ is chosen in this paper. Thus, we will obtain all the solutions at all time points by the segregated scheme.

Since the three fields are highly coupled, the iteration is difficult to converge when the porous medium starts to fracture. The iteration steps will increase to a huge number. In Comsol, the Anderson acceleration method is used for nonlinear convergence acceleration by the information from previous Newton iterations [40]. Therefore, we use the Anderson acceleration method to accelerate convergence. In this paper, the dimension of iteration space field is set as more than 50 to control the number of iteration increments. The maximum number of iteration in one time step is set as 50 in our simulations. Finally, our implementation procedure of the phase field modeling for the fluid-driven crack propagation in Comsol is shown in Fig. 5.

4 Verification of the proposed approach

4.1 A single-phasic solid subjected to tension

The first example is a single-phasic solid subjected to tension. This example has been tested by Miehe et al. [31, 15], Liu et al. [41], Hesch and Weinberg [42], Mikelic et al. [1], and Ehlers and Luo [28]. The geometry and boundary conditions in Fig. 6 are used and the calculation parameters are listed in Table 1. These parameters are taken from [41].

Note that the gravity in all the examples is neglected. We test the notched plate without a pressure ( $p=0$ ). Thus, calculation for the pressure is removed from the segregated steps. The vertical displacement $u(t)=vt$ is prescribed on the top edge of the plate with $v=1\times 10^{-6}$ mm/s. We use 64516 uniform Q4 elements to discretize the solid with the element size $h$ around $3.96\times 10^{-3}$ mm. We use the same element size as in Hesch and Weinberg [42] so that the results are directly comparable. We apply a time step $\Delta t=10$ s for the first 4500 s, afterwards, $\Delta t=1$ s for the remaining time. This example has been solved within 4 h 20 min for $l_{0}=1.5\times 10^{-2}$ mm and 4 h 11 min for $l_{0}=7.5\times 10^{-3}$ mm on two I5-6200U CPUs. As shown in Fig. 7, a horizontal crack propagates through the solid as described in [31, 15, 41, 42, 1, 28]. In addition, Fig. 8 gives the load-displacement curve on the top edge of the solid. The calculated curves are in good agreement with the results by Hesch and Weinberg [42]. Only marginal differences are observed because different algorithms are used.

4.2 A 2D specimen subjected to an increasing internal pressure

The second example simulates the fluid-driven crack propagation of a saturated notched specimen subjected to an increasing internal pressure along the notch. The geometry and boundary conditions are depicted in Fig. 9. We use this example to verify the feasibility of the phase filed model for fracture in poroelastic media. Due to the double symmetry in Fig. 9a, only the top-right quarter of the specimen is modeled, cf. Fig. 9b. AIl the parameters are listed in Table 2.

We discretize the domain with uniformly spaced Q4 elements. Two mesh levels are used: Mesh 1 with the element size $h=2\times 10^{-3}$ m and Mesh 2 with $h=1\times 10^{-3}$ m. A total of 2000 time steps are used with the time increment $\Delta t=1$ s. For Meshes 1 and 2, the time cost on two I5-6200 CPUs is 19 h 18 min and 53 h 23 min, respectively. Analytical solutions are used to verify the numerical results. Considering a single pre-existing crack with a length of $2l_{c}$ in $y=0$ plane, the displacement driven by the pressure in the $y$ direction under the plain strain assumption is given by [43] as

[TABLE]

where $E_{p}=E/(1-\nu^{2})$ is the plane strain Young’s modulus and $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio, respectively.

According to Yoshioka and Bourdin [26] and Mikelic et al. [1], the fracture volume is $V=2\pi pl_{c}/E_{p}$ . Thus, the dissipated energy $E_{b}$ are calculated by

[TABLE]

The energy release rate $G_{c}$ are calculated as $E_{b}/(2l_{c})$ according to Griffith’s theory. Therefore, the critical pressure for fracture propagation is

[TABLE]

Figure 10 compares the displacement along the notch between the present results and the analytical solution. Both the numerical results and analytical solution are in good agreement. Moreover, the displacements given by the analytical solution are larger than those by the phase field model. The critical pressure for fracture initiation between the present results and analytical solution is shown in Fig. 11. The critical pressure increases with increasing $G_{c}$ and the numerical results are larger than those obtained by the analytical solution. In addition, a finer mesh (Mesh 2) will lead to a closer critical pressure to the analytical solution. The slight difference between the simulated critical pressure and the analytical solution occurs because the analytical solution is derived based on a sharp crack shape assumption with fully impermeable crack surfaces. On the other hand, the pre-existing notch in the PFM is a smooth representation of the newly formed crack around the notch tip and the notch surface is slightly permeable (with very low permeability and porosity).

With the increase in the applied pressure, the width of the notch starts to increase and the fracture propagates along the horizontal direction. Figure 12 shows the evolution of the phase field $\phi$ in (a)-(c) and the fluid pressure in (d)-(f). No significant variation in the width of the crack is observed as the time $t$ goes. In addition, the maximum pressure occurs inside the crack although radial distribution of the fluid pressure along the crack is also observed. The pressure gradient is much smaller along the crack than perpendicular to the crack. Therefore, as observed, the crack almost has a uniformly distributed fluid pressure.

Figure 13 gives the streamlines of the fluid flow at $t=1550$ s for Mesh 1. Two directions of fluid flow are observed. The major one is the flow along the propagating crack. The other flow is initially perpendicular to the notch and the crack and finally perpendicular to the boundaries.

5 2D and 3D examples of internal fluid injection

In this section, 2D and 3D examples of specimens subjected to internal fluid injection are presented to demonstrate capability of the proposed modeling approach. Pre-existing cracks are modeled by introducing an initial history field and the source term $q_{F}=$ 100 kg/(m ${}^{3}\cdot\textrm{s}$ ) is used for all examples. The crack propagation is therefore driven by the fluid injection in the pre-existing cracks. The parameters for all calculations can be found in Table 2.

5.1 Two perpendicular propagating cracks in 2D

Let us consider a crack propagation example of two perpendicular pre-existing cracks. Fluid is injected into both cracks and the geometry and boundary conditions are shown in Fig. 14. The length of the cracks is 1 m. The element size is $h=$ $2\times 10^{-2}$ m. A total of 1160 time steps are used with the time increment $\Delta t=0.01$ s. The time cost on two I5-6200 CPUs is 12 h 49 min for this example.

Figure 15 shows the evolution of the phase field. As the time increases, only the horizontal crack propagates in the early period while there is no propagation for the vertical crack. When $t=$ 6.8 s, the horizontal and vertical cracks coalesce. Subsquently, the horizontal crack propagates to the left boundary of the specimen and the vertical crack propagates at a large angle to the horizontal direction. Figure 16 presents the pressure distribution for different times and Fig. 17 depicts the pressure-time curve for the two perpendicular pre-existing cracks. The pressure of the two cracks is almost the same and only little fluctuations occur when the cracks coalesce.

5.2 Propagtion of two parallel penny-shaped cracks in 3D

In the last example, we demonstrate the capability of the proposed approach for 3D problems. We set two parallel penny-shaped cracks in the domain of (-2 m, 2 m)3. The two initial cracks are centered at (0, 0, 0.3 m) and (0, 0, -0.3 m) and on the planes of $y=0.3$ m and $y=-0.3$ m, respectively. Both cracks have the same initial radius of 0.5 m. We set the length scale parameter to $l_{0}=0.06$ m. The maximum element size $h$ is set as 0.06 m and a total of 774,480 6-node prism elements are used to discretize the domain. The time increment is set as $\Delta t=0.01$ s and 1330 time steps are totally used with the time cost of 175 h 30 min on four CPUs.

Crack propagation patterns of the two parallel penny-shape cracks in 3D are shown in Fig. 18. The radii of the cracks become more than 0.5 m when $t=13.2$ s. However, when $t=13.3$ s, the cracks propagate to bowl-shaped volumes. The 3D crack patterns obtained by the proposed approach are in good agreement with those reported in Lee et al. [24]. However, only a coarse mesh used in the 3D simulation and no extra re-meshing and adaptive techniques are adopted in this paper. This fully indicates the strong capability of phase field models to simulate hydraulic fracturing.

6 Conclusions

A phase field model for fracture in poroelastic media is proposed in this paper. The porous medium is modeled based on the classical Biot poroelasticity theory and the fracture behavior is captured by the phase field model. The fracture propagation is driven by the elastic energy and the phase field is used as an interpolation function to transit fluid property from the intact medium to the fully broken one. The model is implemented in Comsol Multiphysics and solved by a segregated scheme. The proposed approach is verified for two benchmark problems: a single-phasic solid subjected to tension and a 2D specimen subjected to an increasing internal pressure. Results are compared to analytical solutions. Finally, 2D and 3D examples of internal fluid injection are presented to show capability of the proposed approach in hydraulic fracturing.

The proposed phase field approach overcomes the disadvantages of previous phase field approaches in HF. Taking the computational cost as an example, Yoshioka and Bourdin [26] used up to 100 CPUs to run a general 3D example within several days while we run the 3D example in this paper in a week only by 4 CPUs. Moreover, we give full implementation details of the approach and make the implementation fully available in Comsol Multiphysics. The results presented in this work can be also extended to inelastic, partially saturated, or heterogeneous porous media and multiple complex fractures.

Acknowledgement

The financial support provided by the Sino-German (CSC-DAAD) Postdoc Scholarship Program 2016, the Natural Science Foundation of China (51474157), and RISE-project BESTOFRAC (734370) is gratefully acknowledged.

Bibliography43

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Mikelic et al. [2013] Andro Mikelic, Mary F Wheeler, and Thomas Wick. A phase field approach to the fluid filled fracture surrounded by a poroelastic medium. ICES report , 1315, 2013.
2Liu et al. [2014] Yaqun Liu, Haibo Li, Chaowen Luo, and Xuechao Wang. In situ stress measurements by hydraulic fracturing in the western route of south to north water transfer project in china. Engineering Geology , 168:114–119, 2014.
3Vidic et al. [2013] Radisav D Vidic, Susan L Brantley, Julie M Vandenbossche, David Yoxtheimer, and Jorge D Abad. Impact of shale gas development on regional water quality. Science , 340(6134):1235009, 2013.
4Ren et al. [2017] Feng Ren, Guowei Ma, Lifeng Fan, Yang Wang, and Hehua Zhu. Equivalent discrete fracture networks for modelling fluid flow in highly fractured rock mass. Engineering Geology , 229:21–30, 2017.
5Figueiredo et al. [2017] Bruno Figueiredo, Chin-Fu Tsang, Jonny Rutqvist, and Auli Niemi. The effects of nearby fractures on hydraulically induced fracture propagation and permeability changes. Engineering Geology , 228:197–213, 2017.
6Moës and Belytschko [2002] Nicolas Moës and Ted Belytschko. Extended finite element method for cohesive crack growth. Engineering fracture mechanics , 69(7):813–833, 2002.
7Johnson and Stryk [1987] Gordon R Johnson and Robert A Stryk. Eroding interface and improved tetrahedral element algorithms for high-velocity impact computations in three dimensions. International Journal of Impact Engineering , 5(1-4):411–421, 1987.
8Areias and Rabczuk [2017] P Areias and T Rabczuk. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design , 132:27–41, 2017.