Phase field modelling of crack propagation, branching and coalescence in rocks

TL;DR

This paper introduces a phase field model implemented in COMSOL to simulate complex crack behaviors in rocks, including propagation, branching, and coalescence, validated against experimental data.

Contribution

The paper presents a novel phase field modeling approach for rock crack analysis, incorporating strain decomposition and complex 3D crack pattern simulations.

Findings

Model accurately predicts crack propagation and coalescence.

Simulations match experimental results.

Applicable to complex 3D crack scenarios.

Abstract

We present a phase field model (PFM) for simulating complex crack patterns including crack propagation, branching and coalescence in rock. The phase field model is implemented in COMSOL and is based on the strain decomposition for the elastic energy, which drives the evolution of the phase field. Then, numerical simulations of notched semi-circular bend (NSCB) tests and Brazil splitting tests are performed. Subsequently, crack propagation and coalescence in rock plates with multiple echelon flaws and twenty parallel flaws are studied. Finally, complex crack patterns are presented for a plate subjected to increasing internal pressure, the (3D) Pertersson beam and a 3D NSCB test. All results are in good agreement with previous experimental and numerical results.

Equations

$$L={\Psi }_{kin}(\dot{\mathbf{u}})-\underset{{\Psi }_{\epsilon }}{\underbrace{{\int }_{\Omega }{\psi }_{\epsilon }(\mathbf{\epsilon })\mathrm{d}\Omega }}-\underset{{\Psi }_f}{\underbrace{{\int }_{\Gamma }{G}_c\mathrm{d}S}}+\underset{W_{ext}}{\underbrace{{\int }_{\Omega }\mathbf{b}\cdot \mathbf{u}\mathrm{d}\Omega +{\int }_{\partial {\Omega }_{h_i}}\mathbf{f}\cdot \mathbf{u}\mathrm{d}S}}$$

$$\mathbf{\epsilon }=\frac{1}{2}\left[\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\mathrm{T}}\right]$$

$${\Psi }_{kin}(\dot{\mathbf{u}})=\frac{1}{2}{\int }_{\Omega }\rho {\dot{\mathbf{u}}}^2\mathrm{d}\Omega$$

$$\varphi =\{\begin{array}{l}0,\text{if material is intact}\\ 1,\text{if material is cracked}\end{array}$$

$$\varphi (x)=e^{-|x|/l_0}$$

$$\gamma (\varphi ,▽\varphi )=\frac{{\varphi }^2}{2l_0}+\frac{l_0}{2}|\nabla \varphi {|}^2$$

$${\int }_{\Gamma }{G}_c\mathrm{d}S={\int }_{\Omega }{G}_c\left(\frac{{\varphi }^2}{2l_0}+\frac{l_0}{2}|\nabla \varphi {|}^2\right)\mathrm{d}\Omega$$

$${\mathbf{\epsilon }}_{\pm }=\sum _{a=1}^d⟨{\epsilon }_a{⟩}_{\pm }{\mathbf{n}}_a\otimes {\mathbf{n}}_a$$

$${\psi }_{\epsilon }^{\pm }(\mathbf{\epsilon })=\frac{\lambda }{2}⟨\mathrm{tr}(\mathbf{\epsilon }){⟩}_{\pm }^2+\mu \mathrm{tr}\left({\mathbf{\epsilon }}_{\pm }^2\right)$$

$$\{\begin{array}{rl}\lambda & =\frac{E\nu }{(1+\nu )(1-2\nu )}\\ \mu & =\frac{E}{2(1+\nu )}\end{array}$$

$${\psi }_{\epsilon }(\mathbf{\epsilon })=\left[(1-k)(1-\varphi {)}^2+k\right]{\psi }_{\epsilon }^{+}(\mathbf{\epsilon })+{\psi }_{\epsilon }^{-}(\mathbf{\epsilon })$$

$$\begin{gathered} L=\frac{1}{2}{\int }_{\Omega }\rho {\dot{\mathbf{u}}}^2\mathrm{d}\Omega -{\int }_{\Omega }\left\{\left[(1-k)(1-\varphi {)}^2+k\right]{\psi }_{\epsilon }^{+}(\mathbf{\epsilon })+{\psi }_{\epsilon }^{-}(\mathbf{\epsilon })\right\}\mathrm{d}\Omega - \\ {\int }_{\Omega }{G}_c\left(\frac{{\varphi }^2}{2l_0}+\frac{l_0}{2}|\nabla \varphi {|}^2\right)\mathrm{d}\Omega +{\int }_{\Omega }\mathbf{b}\cdot \mathbf{u}\mathrm{d}\Omega +{\int }_{\partial {\Omega }_{h_i}}\mathbf{f}\cdot \mathbf{u}\mathrm{d}S \end{gathered}$$

$$\{\begin{array}{r}\text{Div}(\mathbf{\sigma })+\mathbf{b}=\rho \ddot{\mathbf{u}}\\ \left[\frac{2l_0(1-k){\psi }_{\epsilon }^{+}}{G_c}+1\right]\varphi -l_0^2{\nabla }^2\varphi =\frac{2l_0(1-k){\psi }_{\epsilon }^{+}}{G_c}\end{array}$$

$$\mathbf{\sigma }$$

$$=\left[(1-k)(1-\varphi {)}^2+k\right]\left[\lambda ⟨tr(\mathbf{\epsilon }){⟩}_{+}\mathbf{I}+2\mu {\mathbf{\epsilon }}_{+}\right]+\lambda ⟨tr(\mathbf{\epsilon }){⟩}_{-}\mathbf{I}+2\mu {\mathbf{\epsilon }}_{-}$$

$$H(\mathbf{x},t)=\underset{x\in [0,t]}{\max }{\psi }_{\epsilon }^{+}\left(\mathbf{\epsilon }(\mathbf{x},s)\right)$$

$${\psi }_{\epsilon }^{+}-H\le 0,\dot{H}\ge 0,\dot{H}({\psi }_{\epsilon }^{+}-H)=0$$

$$\{\begin{array}{r}\text{Div}(\mathbf{\sigma })+\mathbf{b}=\rho \ddot{\mathbf{u}}\\ \left[\frac{2l_0(1-k)H}{G_c}+1\right]\varphi -l_0^2{\nabla }^2\varphi =\frac{2l_0(1-k)H}{G_c}\end{array}$$

$$\{\begin{array}{rlr}& \mathbf{u}=\mathbf{g}& \mathrm{on}\partial {\Omega }_{g_i}\\ & \mathbf{\sigma }\cdot \mathbf{m}=\mathbf{f}& \mathrm{on}\partial {\Omega }_{h_i}\\ & \nabla \varphi \cdot \mathbf{m}=0& \mathrm{on}\partial \Omega \end{array}$$

$$\{\begin{array}{rlr}& \mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x})& \mathbf{x}\in \Omega \\ & \mathbf{v}(\mathbf{x},0)={\mathbf{v}}_0(\mathbf{x})& \mathbf{x}\in \Omega \\ & \varphi (\mathbf{x},0)={\varphi }_0(\mathbf{x})& \mathbf{x}\in \Omega \end{array}$$

$${\sigma }_{cr}=\frac{9}{16}\sqrt{\frac{E{G}_c}{3l_0}}$$

$$l_0=\frac{27E{G}_c}{256{\sigma }_{cr}^2}$$

$${\int }_{\Omega }\left(-\rho \ddot{\mathbf{u}}\cdot \delta \mathbf{u}-\mathbf{\sigma }:\delta \mathbf{\epsilon }\right)\mathrm{d}\Omega +{\int }_{\Omega }\mathbf{b}\cdot \delta \mathbf{u}\mathrm{d}\Omega +{\int }_{{\Omega }_{h_i}}\mathbf{f}\cdot \delta \mathbf{u}\mathrm{d}S=0$$

$${\int }_{\Omega }-2(1-k)H(1-\varphi )\delta \varphi \mathrm{d}\Omega +{\int }_{\Omega }{G}_c\left(l_0\nabla \varphi \cdot \nabla \delta \varphi +\frac{1}{l_0}\varphi \delta \varphi \right)\mathrm{d}\Omega =0$$

$$\mathbf{u}={\mathbf{N}}_u\mathbf{d},\varphi ={\mathbf{N}}_{\varphi }\hat{\mathbf{\varphi }}$$

$$\bm{N}_{u}=\left[\begin{array}[]{ccccccc}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{array}\right],\hskip 14.22636pt\bm{N}_{\phi}=\left[\begin{array}[]{cccc}N_{1}&N_{2}&\dots&N_{n}\end{array}\right]$$

$$\delta \mathbf{u}={\mathbf{N}}_u\delta \mathbf{d},\delta \varphi ={\mathbf{N}}_{\varphi }\delta \hat{\mathbf{\varphi }}$$

$$\mathbf{\epsilon }={\mathbf{B}}_u\mathbf{d},\nabla \varphi ={\mathbf{B}}_{\varphi }\hat{\mathbf{\varphi }},\mathbf{\delta }\epsilon ={\mathbf{B}}_u\delta \mathbf{d},\nabla \varphi ={\mathbf{B}}_{\varphi }\delta \hat{\mathbf{\varphi }}$$

$$\bm{B}_{u}=\left[\begin{array}[]{ccccccc}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{array}\right],\hskip 5.69046pt\bm{B}_{\phi}=\left[\begin{array}[]{cccc}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{array}\right]$$

$$-(\delta \mathbf{d}{)}^{\mathrm{T}}\left[{\int }_{\Omega }\rho {\mathbf{N}}_u^{\mathrm{T}}{\mathbf{N}}_u\mathrm{d}\Omega \ddot{\mathbf{d}}+{\int }_{\Omega }{\mathbf{B}}_u^{\mathrm{T}}{\mathbf{D}}_e{\mathbf{B}}_u\mathrm{d}\Omega \mathbf{d}\right]+(\delta \mathbf{d}{)}^{\mathrm{T}}\left[{\int }_{\Omega }{\mathbf{N}}_u^{\mathrm{T}}\mathbf{b}\mathrm{d}\Omega +{\int }_{{\Omega }_{h_i}}{\mathbf{N}}_u^{\mathrm{T}}\mathbf{f}\mathrm{d}S\right]=0$$

$$-(\delta \hat{\mathbf{\varphi }}{)}^{\mathrm{T}}{\int }_{\Omega }\left\{{\mathbf{B}}_{\varphi }^{\mathrm{T}}{G}_c{l}_0{\mathbf{B}}_{\varphi }+{\mathbf{N}}_{\varphi }^{\mathrm{T}}\left[\frac{G_c}{l_0}+2(1-k)H\right]{\mathbf{N}}_{\varphi }\right\}\mathrm{d}\Omega \hat{\mathbf{\varphi }}+(\delta \hat{\mathbf{\varphi }}{)}^{\mathrm{T}}{\int }_{\Omega }2(1-k)H{\mathbf{N}}_{\varphi }^{\mathrm{T}}\mathrm{d}\Omega =0$$

Full text

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

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

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

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

5 Institute of Continuum Mechanics, Leibniz University Hannover, Hannover 30167, Germany

• Corresponding author: timon.rabczuk@tdt.edu.vn

Keywords: Phase field, Rock, COMSOL, Crack propagation, Crack branching

1 Introduction

Fracture-induced failure has gained extensive concern in engineering because of the huge threat to engineering safety (Anderson 2005). The prediction of fracture in rock is challenging. Rock masses have many pre-existing flaws, such as micro cracks, voids and soft minerals. Many efforts have been made to study crack propagation in rock, see for instance the contributions in Bobet and Einstein (1998), Wong et al. (2001), Sagong and Bobet (2002), Wong and Einstein (2009), Park and Bobet (2009), Park and Bobet (2010), Lee and Jeon (2011), and Zhou et al. (2014). However, many studies focus on uniaxial compressive loads since tensile loads or more complicated load cases, which are more difficult to perform in practical tests.

Numerical methods are a good alternative to study fracture problems. They are less expensive than experimental tests and can provide physical insight difficult to gain through ’pure’ experimental testing. Computational methods for fracture can be classified in discrete and continuous approaches. Efficient remeshing techniques (Areias and Rabczuk 2017, Areias et al. 2013, Areias and Rabczuk 2013), multiscale method (Budarapu et al. 2014a, b, Yang et al. 2015), strain-softening element (Areias et al. 2014), the extended finite element method (Nanthakumar et al. 2014, Moës and Belytschko 2002), the phantom node method (Rabczuk et al. 2008a, Chau-Dinh et al. 2012, Vu-Bac et al. 2013) and specific meshfree methods (Rabczuk et al. 2007a, Rabczuk and Zi 2007, Rabczuk et al. 2007b, 2008a, Rabczuk and Samaniego 2008, Rabczuk et al. 2008b, Amiri et al. 2014a) are classical representatitves of the first class. The cracking particles method (CPM) (Rabczuk and Belytschko 2004, 2007, Rabczuk et al. 2010), Peridynamics (Rabczuk and Ren 2017) and dual-horizon peridynamics (Ren et al. 2016, 2017) are also discrete crack approaches but they share the simplicity of continuous approaches to fracture as they also do not require any explicit representation of the crack surface and any crack tracking algorithms. Element-erosion (Belytschko and Lin 1987, Johnson and Stryk 1987) directly sets the stresses of the elements to zero when the elements fulfill the fracture criterion. However, the element-erosion method cannot simulate crack branching correctly (Song et al. 2008). Gradient models (Thai et al. 2016), non-local models (Pijaudier-Cabot et al. 2004), models based on the screend-poisson equation (Areias et al. 2016a) and also phase field models are typical continuous approaches to fracture.

In this paper, we pursue the phase field model (PFM) (Bourdin et al. 2008, Miehe et al. 2010a, b, Hesch and Weinberg 2014, Borden et al. 2012) to model crack propagation, branching and coalescence in rock. The origins of the PFM can be traced back to Bourdin et al. (2008), but a thermodynamic consistent framework was first presented by Miehe et al. (2010a). Considerable attention has been paid to PFMs due to their ease in implementation and applicability to multi-physics problems. The PFM does not treat the crack as a physical discontinuity but uses a scalar field (the phase field) to smoothly transit the intact material to the broken one. Thus, the sharp crack is represented by a ’damage-like’ zone. The shape of the crack is controlled by a length scale parameter and propagation of the crack is obtained through the solution of a differential equation. Thus, the PFM does not require any external criterion for fracture and additional work to track the fracture surface algorithmically (Borden et al. 2012). It is believed that for this reason, the phase field is therefore has some advantage over other approaches in modeling branching and merging of multiple cracks.

Phase field models have been discretized in the context of the finite element method (Areias et al. 2016b), meshfree methods (Amiri et al. 2014b) and isogeometric analysis (Borden et al. 2012); the latter two approaches use a fourth-order differential equation for the phase field exploiting the higher continuity of the meshfree and isogeometric approximation. The PFM for brittle cracks has also been implemented in commercial software such as ABAQUS (Msekh et al. 2015, Liu et al. 2016). However, the extension of the implementation in ABAQUS to problems with more fields – as hydraulic fracturing – is difficult. Hence, we present an implementation of the phase field model in COMSOL Multiphysics, a software particularly dedicated to multi-field modeling.

This paper is organized as follows. The phase field model for brittle fractures is presented in Section 2. Subsequently, the numerical implementation of the phase field model in COMSOL is described in Section 3. Then, simulations of initiation, propagation, branching, and coalescence of cracks in rock are shown in Section 4 before Section 5 concludes our manuscript.

2 Theory of phase field modeling

2.1 Theory of brittle fracture

Consider an elastic body $\Omega\subset\mathbb{R}^{d}$ ($d\in\{1,2,3\}$) as shown in Figure 1, whose external boundary and internal discontinuity boundary are denoted as $\partial\Omega$ and $\Gamma$, respectively; $\bm{x}$ is the position vector and $\bm{u}(\bm{x},t)\subset\mathbb{R}^{d}$ the displacement vector at time $t$. In Fig. 1, the body $\Omega$ 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 also the time-dependent Neumann conditions on $\partial\Omega_{h_{i}}\in\Omega$; $\bm{b}(\bm{x},t)\subset\mathbb{R}^{d}$ is the body force and $\bm{f}(\bm{x},t)$ the traction on boundary $\partial\Omega_{h_{i}}$.

Given that the stored elastic energy can be transformed into dissipative forms of energy, the classical Griffith’s theory (Anderson 2005) for brittle fracture states that the crack starts to propagate when the stored energy is sufficient to overcome the fracture resistance of the material. Therefore, the crack propagation is regarded as a process to minimize a free energy $L$ that consists of the kinetic energy $\Psi_{kin}(\bm{\dot{u}})$, elastic energy $\Psi_{\varepsilon}$, fracture energy $\Psi_{f}$ and external work $W_{ext}$:

[TABLE]

where $\bm{\dot{u}}=\frac{\partial{\bm{u}}}{\partial t}$, $\psi_{\varepsilon}$ is the elastic energy density, and $G_{c}$ is the critical energy release rate. The linear strain tensor $\bm{\varepsilon}=\bm{\varepsilon}(\bm{u})$ is given by

[TABLE]

The kinetic energy is given by

[TABLE]

where $\rho$ indicates the density.

2.2 Phase filed approximation for the fracture energy

The phase field method (Miehe et al. 2010a, b, Borden et al. 2012) uses a scalar field, i.e. the phase field, to smear out the crack surface (see Fig. 1) over the domain $\Omega$. The phase field $\phi(\bm{x},t)\in[0,1]$ has to satisfy the following conditions:

[TABLE]

A typical one dimensional phase field approximated by the exponential function is given by (Miehe et al. 2010a)

[TABLE]

$l_{0}$ denoting the length scale parameter, which controls the transition region of the phase field and thereby reflects the width of the crack. The distribution of the one dimensional phase field is shown in Fig. 2. The crack region will have a larger width as $l_{0}$ increases and the phase field will represent a sharp crack when $l_{0}$ tends to zero.

It can be shown that the crack surface density per unit volume of the solid is given by (Miehe et al. 2010a)

[TABLE]

Thus, the fracture energy is approximated by

[TABLE]

The variational approach (Bourdin et al. 2000) states that the crack surface energy is transformed from the elastic energy, which drives the evolution of the phase field. To capture cracks only under tension, the elastic energy is decomposed into tensile and compressive parts (Miehe et al. 2010b):

[TABLE]

where $\bm{\varepsilon}_{+}$ and $\bm{\varepsilon}_{-}$ are the tensile and compressive strain tensors, respectively. In addition, $\varepsilon_{a}$ is the principal strain and $\bm{n}_{a}$ is the direction of the principal strain. The operators $\langle\centerdot\rangle_{\pm}$ are defined as : $\langle\centerdot\rangle_{\pm}=(\centerdot\pm|\centerdot|)/2$. Consequently, the positive and negative elastic energy densities are expressed as

[TABLE]

where $\lambda>0$ and $\mu>0$ are the Lamé constants. The Lamé constants are related to the Young’s modulus $E$ and Poisson’s ratio $\nu$ of the solid through the well known relation:

[TABLE]

The phase field is assumed to affect only the positive elastic energy density, which introduces a stiffness reduction as (Borden et al. 2012)

[TABLE]

where $0<k\ll 1$ is a stability parameter for avoiding numerical singularities because the positive elastic energy density disappears as the phase field $\phi$ tends to 1.

2.3 Governing equations

By substituting Eqs. (3), (7), and (11) into Eq. (1), the energy functional $L$ is rewritten as

[TABLE]

Employ the first variation of the functional $\delta L=0$, it can be shown that the strong form of the governing equations are given by (Borden et al. 2012)

[TABLE]

where $\ddot{\bm{u}}=\frac{\partial^{2}{\bm{u}}}{\partial t^{2}}$ and $\bm{\sigma}$ is Cauchy stress tensor given by

[TABLE]

with unit tensor $\bm{I}$ $\in\mathbb{R}^{d\times d}$.

The phase field requires the irreversibility condition $\Gamma(\bm{x},s)\in\Gamma(\bm{x},t)(s<t)$ during compression or unloading, i.e. the crack cannot be healed. Therefore, we introduce a strain-history field $H(\bm{x},t)$ (Miehe et al. 2010a, b) to ensure a monotonically increasing phase field:

[TABLE]

The history field $H$ satisfies the Kuhn-Tucker conditions for loading and unloading (Miehe et al. 2010a):

[TABLE]

By replacing $\psi_{\varepsilon}^{+}$ by $H(\bm{x},t)$ in Eq. (13), the strong form is rewritten as

[TABLE]

We denote $\bm{m}$ as the outward-pointing normal vector to the boundaries, and the governing equations are subjected to the Dirichlet and Neumann boundary conditions

[TABLE]

along with the initial conditions

[TABLE]

2.4 The choice of $l_{0}$

Borden et al. (2012) and Zhang et al. (2017) proposed an analytical solution for the critical tensile stress $\sigma_{cr}$ that a one-dimensional bar can sustain:

[TABLE]

There is an apparent singularity when $l_{0}$ tends to zero, i.e. in case of a sharp crack, which is phyiscally meaningless. However, assuming all other parameters except of $l_{0}$ are known, eq. (20) can be solved for $l_{0}$:

[TABLE]

In Eq. (21), the critical energy release rate $G_{c}$ and Young’s modulus $E$ can be obtained by conducting regular experimental tests, while the critical stress $\sigma_{cr}$ can be approximated by the tensile strength $\sigma_{t}$ through the standard tensile test. Hence, Eq. (21) can estimate the length scale. Note that the accuracy is unknown for more complex cases.

3 Numerical implementation

3.1 Finite element method

We use the finite element method to solve the governing equations (17) given in weak form by

[TABLE]

and

[TABLE]

We use the standard vector-matrix notation and denote the nodal values of the displacement and phase field as $\bm{u}_{i}$ and $\phi_{i}$. The discretization is thereby given by

[TABLE]

where $\bm{d}$ and $\hat{\bm{\phi}}$ are the vectors consisting of node values $\bm{u}_{i}$ and $\phi_{i}$. $\bm{N}_{u}$ and $\bm{N}_{\phi}$ are shape function matrices given by

[TABLE]

where $n$ is the node number in one element and $N_{i}$ is the shape function of node $i$. The same discretization is applied to the test functions and we obtain

[TABLE]

where $\delta\bm{d}$ and $\delta\hat{\bm{\phi}}$ are the vectors consisting of node values of the test functions.

The gradients are thereby calculated by

[TABLE]

where $\bm{B}_{u}$ and $\bm{B}_{\phi}$ are the derivatives of the shape functions defined by

[TABLE]

By applying the finite element approximation, the equations of weak form (22) and (23) are then written as

[TABLE]

[TABLE]

where $\bm{D}_{e}$ is the degraded stiffness matrix. $\bm{D}_{e}$ can be calculated from the fourth order elasticity tensor $\bm{D}$:

[TABLE]

where $J_{ijkl}=\delta_{ij}\delta_{kl}$, $\delta_{ij}$ being the Kronecker and $P_{ijkl}^{\pm}=\sum_{a=1}^{3}\sum_{b=1}^{3}H_{\varepsilon}(\varepsilon_{a})\delta_{ab}n_{ai}n_{aj}n_{bk}n_{bl}+\sum_{a=1}^{3}\sum_{b\neq a}^{3}\frac{1}{2}\frac{\langle\varepsilon_{a}\rangle_{\pm}-\langle\varepsilon_{b}\rangle_{\pm}}{\varepsilon_{a}-\varepsilon_{b}}n_{ai}n_{bj}(n_{ak}n_{bl}+n_{bk}n_{al})$ with $n_{ai}$ the $i$-th component of vector $\bm{n}_{a}$; $H_{\varepsilon}\langle x\rangle$ is the Heaviside function:

[TABLE]

Since $P_{ijkl}^{\pm}$ cannot be computed when $\varepsilon_{a}=\varepsilon_{b}$, we apply a “perturbation” technology for the principal strains (Miehe 1993) with an unchanged $\varepsilon_{2}$:

[TABLE]

with the perturbation $\delta=1\times 10^{-9}$.

For admissible arbitrary test functions, Eqs. (29) and (30) always hold, thereby producing the discretized weak form as

[TABLE]

[TABLE]

where $\bm{F}_{u}^{ine}$, $\bm{F}_{u}^{int}$, and $\bm{F}_{u}^{ext}$ are the inertial, internal, and external forces for the displacement field and $\bm{F}_{\phi}^{int}$ and $\bm{F}_{\phi}^{ext}$ are the internal and external force terms of the phase field. In addition, the mass and stiffness matrices follow

[TABLE]

In this paper, we use a staggered scheme to solve the displacement and phase fields. Thus, the Newton-Raphson approach is adopted to obtain the residual of the discrete equations $\bm{R}_{u}=\bm{F}_{u}^{ext}-\bm{F}_{u}^{ine}-\bm{F}_{u}^{int}=0$ and $\bm{R}_{\phi}=\bm{F}_{\phi}^{ext}-\bm{F}_{\phi}^{int}=0$, respectively.

3.2 COMSOL implementation

We implemented the phase field approach into COMSOL Multiphysics, which can simulate mathematical and physical problems easily by adding application-specific modules. Therefore, it is suitable for multi-field modeling. In this paper, we construct three main modules, namely, the Solid Mechanics, History-strain and Phase Field Modules, which employ the standard finite element discretization in space as described in Subsection 3.1. In addition, a pre-set Storage Module is employed to evaluate and store the intermediate field variables in a time step, such as the positive elastic energy and principal strains.

Based on a linear elastic material library, the Solid Mechanics Module is used for the displacement $\bm{u}$. The boundary and initial conditions shown in Section 2 are added to the Solid Mechanics Module, while a non-linear stress-strain relationship is considered. The stiffness matrix $D_{e}$ is modified in a time step as follows

[TABLE]

The governing equation (17) is presented for dynamic crack problems. However, for a quasi-static problem, the inertia term must be neglected in the Solid Mechanics Module. The Phase Field Module is established for the phase field $\phi$ by revising a pre-defined module, which is governed by the Helmholtz equation. The boundary condition Eq. (18) and initial condition (19) are also implemented in this module. For the history strain field $H$, the Distributed ODEs and DAEs Interfaces are used to construct the History-strain Module, where the history strain field is not solved directly. We use a “previous solution” solver to record the results in the previous time steps and obtain the field $H$ by applying the following format in COMSOL:

[TABLE]

Additionally, the initial condition $H_{0}(\bm{x})=0$ is used for the History-strain Module.

Figure 3 shows the relationship among all the established modules. The mechanical responses, such as the principal strains, the directions of principal strain, and the elastic energy, are naturally exported from the Solid Mechanics Module and stored in the Storage Module. The History-strain Module then call the positive elastic energy and update the local history strain field $H$. The Phase Field Module employs the updated $H$ to solve the phase field. In a time step, the updated phase field and the stored principal strains as well as the directions are used to modify the stiffness matrix of the Solid Mechanics Module and subsequently to obtain the mechanical responses.

The detailed procedure of the staggered (segregated) scheme is depicted in Fig. 4. The equations of displacement, history strain and phase-field are solved independently. An implicit Generalized-$\alpha$ method (Borden et al. 2012) is used for time integration. The Generalized-$\alpha$ method is unconditionally stable and requires a prediction step and a correction step. When a new time step starts, an initial guess from the linear extrapolation of the previous solution is used in the prediction stage. Newton-Raphson iterations are used to solve the residuals for each module in the correction stage. That is, in the iteration step $j+1$ of a given time step $i$, the displacement $\bm{u}_{i}^{j+1}$ is first solved by using the results ($\bm{u}_{i}^{j}$, $H_{i}^{j}$ and $\phi_{i}^{j}$) of the previous iteration step $j$. Following the updated displacements $\bm{u}_{i}^{j+1}$, the history strain is updated. Subsequently, another Newton-Raphson iteration is applied to solve the phase-field $\phi_{i}^{j+1}$. The total relative error is estimated and the iterations continue until the tolerance requirement is met. The maximum number of iteration in one time step is set as 50 in our simulations. We accelerate the convergence by using the Anderson acceleration (Toth and Kelley 2015) where the dimension of the iteration space field is chosen as more than 50. A flow chart of our implementation is shown in Fig. 5. The source code can be found in “https://sourceforge.net/projects/phasefieldmodelingcomsol/”.

4 Numerical examples of crack propagation, branching and coalescence in rocks

4.1 Simulation of notched semi-circular bend (NSCB) tests

Let us consider the notched semi-circular bend (NSCB) test first. Geometry and boundary condition of the rock specimen are shown in Fig.6. The mechanical properties of the rock are taken from the Laurentian granite (LG) from Grenville province of Canada (Gao et al. 2015). The rock density $\rho$ is 2630 kg/m3, while the Young’s modulus $E$ and Poisson’s ratio $\nu$ are 92 GPa and 0.21, respectively. We follow the 1D solution of Borden et al. (2012) for the critical stress of material softening, and choose $G_{c}$ = 7.6 J/m2 and the length scale $l_{0}=4.5\times 10^{-4}$ m. This produces a critical stress close to the quasi-static tensile strength of the specimen (12.8 MPa) (Gao et al. 2015).

The phase field modeling is performed by using 93587 6-node quadratic triangular elements and the maximum element size is $h=2.25\times 10^{-4}$ m. We apply a vertical displacement on the top of the specimen to drive crack propagation from the tip of the notch. During the simulation, we apply the displacement increment $\Delta u=5\times 10^{-7}$ mm in each time step.

Figure 7 shows the crack initiation and propagation in the rock specimen by using the phase field model. When the applied displacement $u$ reaches to $6.72\times 10^{-3}$ mm, the crack initiates from the tip of the notch. Subsequently, the crack propagates straightly in the vertical direction when $u$ accumulates to $6.74\times 10^{-3}$ mm and $6.77\times 10^{-3}$ mm. When the displacement reaches to $6.86\times 10^{-3}$ mm, the tip of the propagating crack is close to the upper boundary of the specimen and failure of the semi-circular rock specimen occurs. The crack patterns obtained by the phase field simulation are in good agreement with the experimental results in Gao et al. (2015).

We test the crack patterns for $G_{c}=5.6$, 6.6, 8.6, and 9.6 J/m2. The results show that $G_{c}$ has little influence on the crack patterns. Different $G_{c}$ have the same crack paths. Figure 8 gives the curves of the reaction force on the upper boundary of the specimen versus the displacement for different $G_{c}$. A sudden drop of the load is observed after the crack initiation. In addition, all the curves have the same slope before the crack initiation under different $G_{c}$, while the maximum load increases as $G_{c}$ increases .

We also test the influence of the maximum mesh size $h$ under a fixed $G_{c}$ = 7.6 J/m2 and the length scale $l_{0}=4.5\times 10^{-4}$ m. We choose $h=2.25\times 10^{-4}$ m, $1.13\times 10^{-4}$ m, and $5.63\times 10^{-5}$ m in the tests. The resulting load-displacement curve is shown in Fig. 9. As expected, the load-displacement curve converges with the mesh refinement. To test the influence of the length scale $l_{0}$, we fix the mesh size $h=2.25\times 10^{-4}$ m and $G_{c}$ = 7.6 J/m2, and then present the load-displacement curve for different $l_{0}$ in Fig. 10. Figure 10 shows a decreasing peak load and displacement range when the length scale increases.

4.2 Simulation of Brazil splitting tests

Brazil splitting tests are commonly used to obtain the tensile strength of rocks and many researchers simulated crack propagation in a Brazilian disc under compression, such as Cai (2013) and Zhou and Wang (2016). Figure 11 presents the geometry of the Brazilian disc along with the boundary conditions.

These material parameters are adopted: $\rho=2630$ kg/m3, $E$ = 31.5 GPa, and $\nu$ = 0.25. The length scale parameter $l_{0}$ is fixed to 1 mm. 26700 linear triangular elements (base mesh) are used to discretize the disc with the maximum element size $h=0.5$ mm, $k$ is set to $1\times 10^{-9}$. Finally, we conduct the simulation by using five different $G_{c}$: 50, 75, 100, 125, and 150 J/m2.

Figure 12 shows the crack initiation and propagation in the Brazil splitting tests for $G_{c}$ = 100 J/m2. When the displacement $u$ approaches a value of 0.476 mm, the crack occurs in the center of the disc where the maximum tensile stress occurs which is in good agreement with the experimental and analytical results (Atkinson et al. 1982, Entacher et al. 2015). When $u$ reaches a value of 0.477 mm, the crack propagates with an increasing width. Then, the crack continues to propagate and the crack tips move close to the upper and bottom ends of the disc when $u$ reaches to 0.478 mm. The crack branching is observed when $u$ is close to 0.480 mm. In addition, the crack cannot penetrate deeply into the ends of the disc because of the locally compressed area around both ends.

We compare the curve of the reaction force on the upper end of the Brazilian disc versus the displacement $u$ with the experimental result in Fig. 13. The experimental curve is originated from the work of Erarslan and Williams (2012). Figure 13 shows that the phase field model can reproduce the experimental results well. The main reason for the difference in Fig. 13 is that there are contact issues in the actual Brazilian tests and a compaction stage is commonly observed before the elastic stage.

Figure 14 presents the curves of the reaction force on the upper end of the Brazilian disc versus the displacement $u$ for different $G_{c}$. Similar to the NSCB tests, a sudden drop of the load is observed after the phase field increases to 1. In addition, the peak load increases with the increase in $G_{c}$. We then test the influence of the length scale $l_{0}$ under a fixed $G_{c}=100$ J/m2 and mesh size $h=0.5$ mm. The tested length scale parameters are $l_{0}=$ 0.5, 1, 2, and 3 mm, respectively. The load-displacement curves under different length scale $l_{0}$ are shown in Fig. 15. As shown in Fig. 15, the peak load of the specimen decreases with the increase in the length scale $l_{0}$.

We also test the influence of mesh size $h$ under a fixed $G_{c}=100$ J/m2 and $l_{0}$ = 1 mm. The maximum mesh size $h$ is set as 1, 0.5, 0.25, and 0.125 m, respectively. The resulting load-displacement curve is shown in Fig. 16. The load-displacement curve converges with mesh refinement as expected. In addition, in the Brazil splitting test, the tensile strength $\sigma_{t}$ of the rock specimen is given by

[TABLE]

where $P_{peak}$ is the peak load, and $D$ and $L$ are the diameter and length of the rock specimen. Thus, we compare the tensile strength by the phase field simulation and the critical stress for 1D by Borden et al. (2012) in Fig. 17. The tensile strength increases at a decreasing rate with the increase in $G_{c}$. The critical stress has the same trend as the tensile strength. However, the critical stress is far larger than the tensile strength. This observation indicates that the critical stress for 1D analysis cannot be applied directly to 2D or 3D.

4.3 Propagation of multiple echelon flaws

We consider a 50 mm $\times$ 50 mm square rock sample subjected to tension. The rock sample has three pre-existing flaws, whose position and geometry are shown in Fig. 18. All the flaws have the same length, spacing, and inclination angle of $45^{\circ}$. These parameters are adopted: the rock density $\rho=2500$ kg/m3, the Young’s modulus $E=30$ GPa, the Poisson’s ratio $\nu=0.333$, $G_{c}=3$ J/m2, $k=1\times 10^{-5}$, and the length scale $l_{0}=0.25$ mm. Vertical displacements are applied on the top and bottom boundaries of the rock sample as shown in Fig. 18. The simulation is performed by using 106852 6-node quadratic triangular elements where the maximum element size $h$ is 0.25 mm. In each time step, the displacement increment is $\Delta u=5\times 10^{-7}$ mm.

Figure 19(a)-(f) shows the propagation and coalescence of the three pre-existing flaws. We also calculate the reaction force on the upper boundary of the rock sample, and present the load-displacement curves in Fig. 20. As the displacement $u$ increases, both the load and the phase field around the tips of the flaws increase. When $u$ reaches to $2.6\times 10^{-3}$ mm, the phase field increases close to 1 and the load approaches to the maximum value. When $u$ reaches to $2.63\times 10^{-3}$ mm, the first tensile cracks occurs around the left and right tips of the flaw $②$. These two cracks are perpendicular to the direction of the applied displacement. In addition, the load reaches to the maximum as shown in Fig. 20. When the displacement $u$ reaches to $2.64\times 10^{-3}$ mm, the cracks from the tips of the flaw $②$ propagates perpendicular to $u$ while the load starts to drop after the peak load. Additionally, new cracks occur from the right tip of the flaw $①$ and the left tip of the flaw $③$.

Figure 19(d) shows the coalescences of the cracks from the tips of the flaws $①$, $②$, and $③$ when $u=2.65\times 10^{-3}$ mm. Thus, the flaw $②$ links the flaws $①$ and $③$, which is accompanied by a sudden drop of the load after the maximum value in Fig. 20. At the same time, new cracks from the left tip of the flaw $①$ and the right tip of the flaw $③$ are observed. When $u=2.66\times 10^{-3}$ mm, the cracks from the flaws $①$ and $③$ continue to propagate and the load decreases. When the displacement increases to $2.69\times 10^{-3}$ mm, the cracks initiating from the left tip of the flaw $①$ and the right tip of the flaw $③$ reach the boundaries of the rock sample. This indicates the rock loses its load-bearing capacity, which can also be verified by Fig. 20.

We now consider the same square rock sample subjected to tension with nine pre-existing flaws in Fig. 21. These flaws have the same inclination angle of $45^{\circ}$, while the lengths and spacing are not fixed. The position and geometry of the flaws are shown in Fig. 21. The parameters are the same as those in the example of three flaws. A total of 107982 6-node quadratic elements are used to discretize the rock sample and the maximum element size $h$ is 0.25 mm. The displacement increment $\Delta u=1\times 10^{-6}$ mm is applied in each time step.

Figure 22(a)-(f) shows the propagation and coalescence process of the nine pre-existing flaws, and the reaction force on the upper boundary of the rock sample is depicted in Fig. 23. As the displacement $u$ increases, the phase field around the tips of the flaws and the load both increase. The first tensile cracks initiate from the left tips of the flaws $⑧$ and $⑨$ when the displacement $u$ reaches to $3.02\times 10^{-3}$ mm. At this time, the load achieves the maximum value. As the displacement increases to $3.04\times 10^{-3}$ mm, the crack from the left tip of the flaw $⑧$ propagates and links up the flaw $⑦$. The crack from the left tip of the flaw $⑨$ continues to propagate while new cracks initiate from the left tip of the flaw $⑦$ and the right tips of the flaws $⑧$ and $⑨$. The load then has a drop after the peak pint (a). When $u$ reaches to $3.06\times 10^{-3}$ mm and $3.08\times 10^{-3}$ mm, the cracks initiating from the left tip of the flaw $⑨$ and the right tip of the flaw $⑧$ continue to propagate at a decreasing rate and at a small angle with the horizontal. However, the cracks from the left tip of the flaw $⑦$ and the right tip of the flaw $⑨$ propagate at relatively large velocity nearly along the horizontal direction. The load decreases sharply as the applied displacement increases.

Figure 22(e) shows that the crack from the left tip of the flaw $⑦$ propagates close to the left boundary of the rock sample when $u=3.12\times 10^{-3}$ mm. The crack initiating from the right tip of the flaw $⑨$ propagates at a relatively small rate because of a larger distance from the boundary where the applied displacement is applied. When $u=3.14\times 10^{-3}$ mm, the cracks from the left tip of the flaw $⑦$ and the right tip of the flaw $⑨$ both reaches the left and right boundaries of the rock sample, indicating the failure of rock sample. The rock loses its load-bearing capacity and the load drops to near 0 in Fig. 23. In addition, Fig. 22 shows no cracks initiation from the tips of the flaws $①$-$⑥$. The reason is the stress shielding and amplification effects due to the interaction of the flaws Zhou et al. (2015).

4.4 Propagation and coalescence of twenty parallel flaws

A 2D square rock sample with twenty parallel pre-existing flaws subjected to tension is tested. All flaws have the same length, spacing, and inclination angle of $0^{\circ}$. We consider the flaws in doubly periodic rectangular and diamond-shaped arrays, respectively. The arrangement and geometry of the flaws are depicted in Fig. 24. The rock sample are 50 mm $\times$ 50 mm. These parameters are adopted: the rock density $\rho=2450$ kg/m3, the Young’s modulus $E=30$ GPa, the Poisson’s ratio $\nu=0.3$, $G_{c}=100$ J/m2, $k=1\times 10^{-9}$, and the length scale $l_{0}=0.4$ mm. The rock sample is discretized by using uniform 8-node quadratic elements with the element size $h=0.2$ mm. We adopt the displacement increment $\Delta u=5\times 10^{-6}$ mm for each time step.

Figure 25 presents the load-displacement curves for the square rock sample with twenty parallel flaws. A sudden drop of the load is also observed after the maximum value is obtained. Figure 26 presents the propagation and coalescence of the doubly periodic rectangular array of twenty pre-existing flaws. As the displacement $u$ increases, the phase field $\phi$ concentrates at the tip of each flaw and the load achieves the maximum when $u=2.185\times 10^{-2}$ mm. When the displacement $u$ reaches to $2.19\times 10^{-2}$ mm, the first tensile cracks initiate from the left and right tips of the flaws $①$, $⑤$, $\footnotesize{17}⃝$, and $\footnotesize{20}⃝$. In addition, the load has a small drop from the peak. When $u=2.195\times 10^{-2}$ mm, the first cracks continue to propagate perpendicular to the direction of the applied displacement and the load decreases. However, when $u=2.2\times 10^{-2}$ mm, the cracks fully connect the flaws $①$, $⑤$, $\footnotesize{17}⃝$, and $\footnotesize{20}⃝$ with the flaws $②$, $③$, $\footnotesize{18}⃝$, and $\footnotesize{19}⃝$. Moreover, some new cracks initiate from the right tips of the flaws $②$ and $\footnotesize{18}⃝$ as well as the left tips of the flaws $③$ and $\footnotesize{19}⃝$. The load then drops to approximately $60\%$ of the maximum load. When $u=2.205\times 10^{-2}$ mm, the cracks initiating from the right tips of the flaws $②$ and $\footnotesize{18}⃝$ and the cracks from the left tips of the flaws $③$ and $\footnotesize{19}⃝$ coalesce. Finally, the cracks from the left tips of the flaws $①$ and $\footnotesize{17}⃝$ as well as the cracks from the right tips of the flaws $④$ and $\footnotesize{20}⃝$ reach the left and right sides of the rock sample when $u=2.21\times 10^{-2}$ mm, indicating that the rock loses the load-bearing capacity. In addition, no cracks initiate from the tips of the flaws $⑤$ - $\footnotesize{16}⃝$.

Figure 27 presents the propagation and coalescence of the twenty pre-existing flaws in the diamond-shaped array. The crack patterns for the diamond-shaped array are different from those for the doubly periodic rectangular array. As the displacement $u$ increases, the phase field $\phi$ rapidly increases around the tips closest to the left and right boundaries of the rock sample. When the displacement $u$ reaches to $2.11\times 10^{-2}$ mm, the load increases to the maximum. When $u=2.115\times 10^{-2}$ mm, the tensile cracks first emanate from the left tips of the flaws $⑤$ and $\footnotesize{13}⃝$ as well as the right tips of the flaws $④$ and $\footnotesize{20}⃝$. When $u=2.12\times 10^{-2}$ mm, another two cracks initiate from the left tips of the flaws $④$ and $\footnotesize{20}⃝$. The first four cracks continue to propagate along the horizontal direction, while the load at $u=2.12\times 10^{-2}$ mm is close to that at $u=2.11\times 10^{-2}$ mm and $u=2.115\times 10^{-2}$ mm.

The flaws close to the top and bottom edges of the sample are placed perpendicular to the direction of loading. The arrangement of the internal flaws has little effect on the propagation of the edge flaws. Therefore, cracks propagating from those edge flaws will show typical features of Mode-I fracture. As shown in Fig. 27, when $u=2.130\times 10^{-2}$ mm, the originally initiating cracks propagate as expected, while new cracks emanate from the right tips of the flaws $③$, $⑤$, $\footnotesize{13}⃝$, and $\footnotesize{19}⃝$. The load then starts to drop from the peak region. However, when $u=2.135\times 10^{-2}$ mm, the cracks fully connect the flaws $③$ and $\footnotesize{19}⃝$ with the flaws $④$ and $\footnotesize{20}⃝$. In addition, two new cracks initiate from the left tips of the flaws $③$ and $\footnotesize{19}⃝$. At this time, the load decreases to approximately half of the peak load. As the displacement $u$ increases to $2.145\times 10^{-2}$ mm, the propagating cracks and the flaws $②$, $③$, $⑤$, $⑥$, $\footnotesize{13}⃝$, $\footnotesize{14}⃝$, $\footnotesize{18}⃝$, and $\footnotesize{19}⃝$ coalesce. At the same time, the cracks from the left tips of the flaws $⑤$ and $\footnotesize{13}⃝$ as well as the cracks from the right tips of the flaws $④$ and $\footnotesize{20}⃝$ reach the left and right sides of the rock sample. The load at $u=2.145$ mm is less than zero followed by new cracks emanating from the left tip of the flaw $②$ and the right tip of the flaw $⑥$.

When the displacement $u$ reaches to $2.195\times 10^{-2}$ mm and $2.3\times 10^{-2}$ mm, the crack initiating from the right tip of the flaw $\footnotesize{14}⃝$ slowly propagates; however, the cracks from the left tip of the flaw $②$ and the right tip of the flaw $⑥$ continue to propagate obliquely at a relatively large velocity. The difference in the crack patterns from Figs. 26 and 27 implies that initiation, propagation and coalescence of the cracks are significantly affected by the arrangement of the flaws.

4.5 Crack branching in a plate subjected to internal pressure

This example is a square plate subjected to internal pressure with geometry and boundary conditions shown in Fig. 28. The internal pressure is applied on the upper and lower boundaries of the notch with $\bar{p}=1$ MPa/s, while the outer boundaries of the plate are traction-free. In this example, dynamic cracks are considered, and these parameters are adopted: the rock density $\rho=2450$ kg/m3, the Poisson’s ratio $\nu=0.3$, $G_{c}=1$ J/m2, $k=1\times 10^{-9}$, and the length scale $l_{0}=0.4$ mm. The plate is discretized by using uniform Q4 elements with the element size $h=l_{0}$, and we adopt the time step size $\Delta t=0.01$ $\mu$s.

We consider the plate as a heterogeneous material and apply a Weibull distribution to the Young’s modulus:

[TABLE]

where $\varphi$ is the probability density function and coefficient $m$ determines the shape of $\varphi$. $m$ also reflects the homogeneity of the material. As $m$ increases, the material becomes more homogeneous and vice versa. In this paper, we consider $m=1$, 3, 5, 7, and 9 along with $m=\infty$ representing a homogeneous plate. We use $E_{0}=30$ GPa and Fig. 29 shows the distribution of Young’s modulus for different $m$.

Figure 30 presents the crack patterns of the plate subjected to internal pressure for different $m$. Complex crack patterns, such as many crack branching, are observed. These observations are different from the previous examples. When the time $t$ reaches 20 $\mu$s, cracks initiate from the left and right tips of the notch and then start to branch. When the time $t$ increases to 30 $\mu$s, the branching cracks propagate. When $t=$ 50$\mu$s, large “damage” region are observed around the notch where $\phi$ is large. At $t=60$ $\mu$s, some cracks initiate from the upper and lower boundaries of the plate, while the cracks from the notch keep propagating. Finally, when $t=70$ $\mu$s, many new cracks occur inside the plate and the cracks from the upper and lower boundaries of the plate start to branch. For a smaller $m$, more crack branching is observed during the crack propagation. Therefore, a smaller $m$ produces more complex crack patterns.

4.6 3D Petersson beam

In this example, the phase field model is applied to a single edge notched beam subjected to three-point bending (the so-called Petersson beam). A 3D simulation is conducted in this example. The geometry and boundary conditions are depicted in Fig. 31 according to Petersson (1981). The thickness of the beam is 50 mm. The following mechanical properties of the beam are chosen Petersson (1981): Young’s modulus $E$ = 27 GPa, Poisson’s ratio $\nu$ = 0.21, critical energy release rate $G_{c}$ = 56 J/m2 and length scale $l_{0}$ = 1 mm. We choose the maximum element size $h$ = 5 mm in most of the beam but $h$ = 1 mm in the region where the crack is expected to propagate. A displacement increment of $\Delta u=2.5\times 10^{-4}$ mm is used.

Figure 32 compares the load-displacement curves obtained of the phase field model with the experimental test (Petersson 1981). The results obtained by the phase field model are in good agreement with the experimental test. Figure 33 presents the crack propagation in the beam at the displacements $u=0.34$ mm, $0.355$ mm, and $0.46$ mm.

4.7 3D NSCB tests

In the last example, we simulate the crack propagation in a 3D NSCB specimen. The geometry and boundary conditions are similar to the 2D tests except a thickness of 16 mm in 3D. The same parameters are used as those in the 2D tests and $G_{c}=7.6$ J/m2. We refine the elements with the maximum size $h=4.5\times 10^{-4}$ m in the region where the crack is expected to propagate, while in the rest region $h=1.8\times 10^{-3}$ m. In addition, the displacement increment $\Delta u=5\times 10^{-7}$ mm is applied.

Figure 34 presents the crack propagation at the displacements $u=6.67\times 10^{-3}$ mm, $6.71\times 10^{-3}$ mm, and $6.745\times 10^{-3}$ mm. Only the domain with $\phi>0.95$ are displayed for the crack shape. The crack patterns are the same as those in the 2D simulation, and also in good agreement with the results of the experimental tests (Gao et al. 2015). The example of 3D NSCB test shows the ability and practicability of the phase field method in modeling crack propagation of rocks in 3D.

5 Conclusions

The phase field theory for fracture is applied to study the crack propagation, branching and coalescence in rocks. Implementation details of the phase field modeling in COMSOL are presented with the consideration of cracks only due to tension. The numerical simulations of the 2D notched semi-circular bend (NSCB) tests and Brazil splitting tests are then performed. The presented results are in good agreement with those of the previous experimental tests. Subsequently, the crack propagation and coalescence in plates with multiple echelon flaws and twenty parallel flaws are studied. We also present the complex crack patterns in a plate subjected to internal pressure, the increase of which produces crack propagation and branching. Finally, the simulation of a 3D Petersson beam and a 3D NSCB test are performed to show the practicability of phase field modeling in 3D rocks.

All the numerical examples presented by this work show that the initiation, propagation, coalescence, and branching of cracks are autonomous, while the phase field modeling does not require external criterion for fracture and setting propagation path in advance. These observations highlight the advantages of the phase field method over other numerical methods in modeling complex crack propagation in rocks. Therefore, the phase field modeling approach will be useful and practicable for other crack problems in rock engineering in future research. In addition, the presented phase field model cannot predict the shear cracks when a rock reaches its shear strength. The reason is that the shear strength is not involved in the formulation of the phase field method and the crack propagation is only driven by the elastic energy. In this sense, a modified phase field model coupled with the shear model of rocks will be also attractive in the future.

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.