Phase-field modeling of rock fractures with roughness

Consider the solid body $\Omega\in\mathbb{R}^{\mathrm{dim}}$ ($\mathrm{dim}=2,3$) with boundary $\partial\Omega$. The boundary is suitably decomposed into the displacement (Dirichlet) boundary $\partial_{u}\Omega$ and the traction (Neumann) boundary $\partial_{t}\Omega$ such that $\partial_{u}\Omega\cap\partial_{t}\Omega=\emptyset$ and $\overline{\partial_{u}\Omega\cup\partial_{t}\Omega}=\partial\Omega$. The body may contain a set of discontinuities denoted by $\Gamma$, which may evolve in the time domain $\mathcal{T}$.

Figure 1 illustrates how phase-field modeling diffusely approximates the original discrete problem. The approximation begins by introducing the phase-field variable, $d\in[0,\,1]$, where $d=0$ indicates the undamaged (bulk) region and $d=1$ indicates the fully cracked (interface) region. We then introduce a crack density functional, $\Gamma_{d}(d,\operatorname{\nabla}d)$, that satisfies

A general form of $\Gamma_{d}(d,\operatorname{\nabla}d)$ is given by

Here, $L$ is the phase-field length parameter controlling the width of diffuse approximation. The specific form of $\Gamma_{d}(d,\operatorname{\nabla}d)$ is determined by the choice of $\alpha(d)$. Among a few forms of $\alpha(d)$ proposed in the literature [35, 36], here we choose $\alpha(d)=d$ for its relative simplicity and its relation to cohesive fracture [36]. We note that other choices of $\alpha(d)$ would equally be valid for the purpose of diffuse approximation of stationary cracks (see Fei and Choo [21] where $\alpha(d)=d^{2}$ is used), but for simulating evolving fractures, they would lead to different results. Substituting $\alpha(d)=d$ into Eq. (2) gives

The phase-field formulation gives rise to two governing equations—the linear momentum balance equation and the phase-field evolution equation—of which the primary variables are the displacement vector, $\bm{u}$, and the phase field, $d$. For brevity, we omit their derivations and refer to the relevant literature (e.g. [36, 22]). The governing equations may be written as

In the momentum balance equation (4), $\bm{\sigma}$ is the Cauchy stress tensor, $\rho$ is the mass density, and $\bm{g}$ is the gravitational acceleration vector. In the phase-field evolution equation (5), $\mathcal{G}_{c}$ is the critical fracture energy, $\mathcal{H}^{+}$ is the crack driving force, and $g(d)$ is the degradation function. It is noted that the specific forms of $\mathcal{H}^{+}$ and $g(d)$ should be chosen in accordance with the selected form of the crack density functional, $\Gamma_{d}(d,\operatorname{\nabla}d)$. For the particular form of $\Gamma_{d}(d,\operatorname{\nabla}d)$ in Eq. (3), $g(d)$ is given by

Here, $\mathcal{H}_{t}$ is defined as the threshold value of the crack driving force at the peak material strength. The specific expressions for $\mathcal{H}^{+}$ and $\mathcal{H}_{t}$ depend also on how the stress tensor is decomposed into its crack driving and non-driving parts. So they will be discussed after presenting the stress decomposition scheme.

To complete the problem statement, we write the boundary conditions as

where $\hat{\bm{u}}$ and $\hat{\bm{t}}$ are the prescribed displacement and traction vectors on the boundary, respectively, and $\bm{v}$ is the unit vector outward normal to the boundary. The initial conditions are given by

where $\bm{u}_{0}$ and $d_{0}$ are the initial displacement and phase-field solutions, respectively.

Given that our interest is in fracture under compressive loads, we adopt the stress decomposition scheme proposed by Fei and Choo [21], which is a verified way to incorporate frictional contact into phase-field-approximate discontinuities. It calculates the stress tensor $\bm{\sigma}$ as the weighted average of the stress in the rock matrix (bulk region), $\bm{\sigma}_{m}$, and the stress in the fracture (interface region), $\bm{\sigma}_{f}$,¹¹1$\bm{\sigma}_{m}$ and $\bm{\sigma}_{f}$ correspond to $\bm{\sigma}_{\mathrm{bulk}}$ and $\bm{\sigma}_{\mathrm{interface}}$ in Fei and Choo [21], respectively. with the weighting done according to the degradation function, $g(d)$. Specifically,

The matrix stress tensor can be calculated using a standard continuum constitutive relationship. Here we assume that the rock matrix is linear elastic, such that the constitutive relationship is given by

where $\bm{\varepsilon}$ is the infinitesimal strain tensor, and $\mathbb{C}^{\mathrm{e}}$ is the elastic stiffness tensor. The linear elastic stiffness tensor can be expressed using the bulk modulus, $K_{m}$, and the shear modulus, $G_{m}$ of the rock matrix as

with $\mathbb{I}$ and $\bm{1}$ denoting the fourth-order symmetric identity tensor and the second-order identity tensor, respectively.

The fracture stress tensor $\bm{\sigma}_{f}$, which is nonzero whenever $d>0$, is computed depending on the contact condition at the material point. To determine the contact condition, we introduce an interface-oriented coordinate system. In 2D, for example, the coordinate system is defined with the unit normal vector $\bm{n}$ and the unit tangent vector $\bm{m}$ of the discontinuity, see Fig. 2. In this interface-oriented coordinate system, the normal strain $\varepsilon_{\mathrm{N}}$ is calculated as

We can distinguish between open and closed contact conditions based on the sign of $\varepsilon_{\mathrm{N}}$: open if $\varepsilon_{\mathrm{N}}>0$, and closed otherwise. When the crack is open, $\bm{\sigma}_{f}=\bm{0}$. On the other hand, when the crack is closed, $\bm{\sigma}_{f}$ should be calculated such that it incorporates the contact behavior of the discontinuity.

Having decomposed the stress tensor as above, let us now determine the specific expression for the crack driving force, $\mathcal{H}^{+}$. Considering that rocks are quasi-brittle (rather than perfectly brittle) and cracks from rough discontinuities are mostly tensile [43], here we adopt the crack driving force for cohesive tensile fracture derived from a directionally decomposed stress tensor [17]. Employing the popular approach to crack irreversibility that uses the maximum $\mathcal{H}^{+}$ in history [44], the crack driving force is given by

Here, $M_{m}:=K_{m}+(4/3)G_{m}$ is the constrained modulus of the matrix, $\sigma_{1,m}$ is the maximum (tensile) principal stress in the matrix, and $\sigma_{p}$ is the peak tensile strength. The detailed derivation of Eq. (16) can be found from Fei and Choo [17]. Note that one may also arrive at the same expression for the crack driving force by extending the crack driving force in Steinke and Kaliske [45], where the same type of directional stress decomposition is used for brittle tensile fracture, to cohesive tensile fracture as done in Geelen et al. [36].

In previous works where the above stress decomposition scheme is employed [21, 22, 17, 24], the contact behavior is modeled using the standard Coulomb friction law. However, as explained in the Introduction, the Coulomb friction law does not incorporate a complete array of important features that emanate from the rough nature of rock fractures. In the following section, we develop a general constitutive framework that allows one to cast a standard model for rough rock fractures—originally formulated for discrete modeling of discontinuities—into $\bm{\sigma}_{f}$ in diffusive phase-field modeling.

In discrete approaches to modeling discontinuities inside continuous materials, the displacement field is often decomposed into its continous part, $\bm{u}_{c}$, and its discontinuous part—the displacement jump—$[\![\bm{u}]\!]:=\bm{u}^{+}-\bm{u}^{-}$. Mathematically, the decomposition can be written as

where $H(\bm{x})$ is the Heaviside function, defined as

with $\Omega_{+}$ and $\Omega_{-}$ denoting the subdomains that are separated by the discontinuity $\Gamma$.

To model the mechanical behavior of discontinuities, a constitutive law must be introduced relating the displacement jump $[\![\bm{u}]\!]$ to the mobilized traction $\bm{t}$ on the surface. This constitutive law is often expressed in a generic rate form as

where $\mathbb{D}$ is the instantaneous tangent stiffness of the fracture. It is convenient to develop such models using the surface aligned coordinate system in Fig. 2. We therefore decompose the displacement jump vector into its normal and tangential components as

where scalars $u_{f}$ and $v_{f}$ denote the magnitude of the normal closure and tangential slip, respectively.

Our first task is to devise a way to model the kinematics of fracture modeled by the phase-field method, which lacks an explicit representation of the displacement jump $[\![\bm{u}]\!]$. To begin, we calculate the strain tensor as the symmetric gradient of the decomposed displacement field (17)

where

are the continuous and discontinuous parts of the strain tensor, respectively, and $\delta_{\Gamma}(\bm{x})$ is the Dirac delta function emanating from the gradient of the Heaviside function, i.e. $\operatorname{\nabla}H_{\Gamma}(\bm{x})=\bm{n}\delta_{\Gamma}(\bm{x})$. Unfortunately, Eq. (23) is incompatible with the phase-field setting, because the Heaviside and Dirac delta functions cannot be defined when discontinuous geometry is approximated diffusely. To overcome this issue, we approximate Eq. (23) in the following two ways. First, we postulate that the displacement jump is constant along the discontinuity, such that $\operatorname{\nabla^{s}}[\![\bm{u}]\!]=\bm{0}$. The same postulate was also introduced in the assumed enhanced strain (AES) formulation of Regueiro and Borja [46]. Second, we replace the Dirac delta function by the crack density functional, $\Gamma_{d}(d,\operatorname{\nabla}d)$, which converges to the Dirac delta function as $L\rightarrow 0$ ($\Gamma$-convergence). The same approximation can be found from Verhoosel and de Borst [47]. As a result, $\bm{\varepsilon}_{f}$ simplifies to

Hereafter, we shall use Eq. (24) to calculate $\bm{\varepsilon}_{f}$ in the phase-field setting. Further, combining Eqs. (20) and (24), we can calculate the normal closure and tangential slip as

respectively, where

is the slip direction tensor.

It is also postulated that $\bm{\varepsilon}_{c}$ and $\bm{\varepsilon}_{f}$ can be additively decomposed into elastic and plastic (inelastic) parts as

where the superscripts $(\cdot)^{\mathrm{e}}$ and $(\cdot)^{\mathrm{p}}$ denote the elastic and plastic parts, respectively. Equivalently, the total elastic and plastic strains are

Although many frictional contact models (e.g. the standard Coulomb friction model) treat frictional slip as entirely inelastic, rough rock fracture can manifest significant elastic deformation due to asperity compressibility and contact effects. For this reason, elasticity is an integral part of modeling rock fracture with roughness.

Combining Eqs. (25), (26), and (29), we get the following approximations of the elastic and plastic components of normal and tangential displacement jumps:

A constitutive model for rock discontinuities can then be introduced to relate the kinematic quantities in Eqs. (32)–(35) to their corresponding stress components of $\bm{\sigma}_{f}$. For this purpose, we also decompose $\bm{\sigma}_{f}$ into its normal and tangential component as

where $\sigma_{\mathrm{N}}$ is the contact normal stress and $\tau$ is the shear stress in the fracture. Note that here $\sigma_{\mathrm{N}}$ (and thus the fracture displacement) is considered negative in compression, being consistent with the sign convention in standard phase-field modeling. In what follows, we formulate a constitutive framework for modeling the compression and shear behavior of rock fractures, from the elastic regime to the plastic regime.

We make use of a standard approach to rock joint elasticity that describes the fracture normal and shear response separately. For the normal traction, we adopt the hyperbolic function proposed by Bandis et al. [51]—perhaps the most popular approach in the rock mechanics community—given by

where $\kappa$ denotes the initial compressive stiffness, and $u^{\mathrm{e}}_{\max}$ the maximum closure of the fracture. The tangential traction may be described as

where $\mu$ denotes the shear modulus related to the frictional properties of the fracture. It is noted that $\mu$ may or may not be dependent on $\sigma_{\mathrm{N}}$ according to the modeling assumption.

We now reformulate Eqs. (38) and (39) for the phase-field formulation at hand to construct a constitutive relationship between $\bm{\sigma}_{f}$ and $\bm{\varepsilon}_{f}^{\mathrm{e}}$. In the present work, we consider the normal and tangential deformations of the fracture and neglect rotations and relative stretching, as in the classic AES method for strong discontinuities [52, 46]. Let us consider the normal deformation first. Substituting Eq. (32) into Eq. (38) and rearranging the resulting equation gives

Inserting Eq. (36) into Eq. (40), we can rewrite the above equation as

where $E_{f}$ is the Young’s modulus of fracture, defined as

Likewise, for the tangential deformation, we can use Eq. (33) to rewrite Eq. (39) as

where $G_{f}$ is the shear modulus of the fracture, defined as

To arrive at a constitutive relation between $\bm{\sigma}_{f}$ and $\bm{\varepsilon}_{f}^{\mathrm{e}}$, we multiply both sides of Eq. (41) by $\bm{n}\operatorname{\otimes}\bm{n}$ and get

Similarly, we multiply both sides of Eq. (43) by $(1/2)\bm{\alpha}$ and get

Adding Eqs. (45) and (46), we obtain the constitutive relation between $\bm{\sigma}_{f}$ and $\bm{\varepsilon}_{f}^{\mathrm{e}}$ as follows:

The continuous part of the elastic strain, $\bm{\varepsilon}_{c}^{\mathrm{e}}$, is related to $\bm{\sigma}_{f}$ as

Since $\bm{\varepsilon}^{\mathrm{e}}=\bm{\varepsilon}^{\mathrm{e}}_{c}+\bm{\varepsilon}^{\mathrm{e}}_{f}$, Eqs. (47) and (48) can be combined as

where

Note that $\mathbb{K}$ is not constant and depends on the stress state. This equation then allows us to derive the elastic tangent of the fracture stress tensor as

Noting that $E_{f}$ and $G_{f}$ depend on $\bm{\sigma}_{f}$ via Eqs. (42) and (44), we take derivatives of both sides of Eq. (49) with respect to $\bm{\varepsilon}^{\mathrm{e}}$ and obtain

Rearranging Eq. (53) gives $\mathbb{C}^{\mathrm{e}}_{f}$ as

where

with $\mu^{\prime}(\sigma_{\mathrm{N}})$ denoting the derivative of $\mu$ with respect to $\sigma_{\mathrm{N}}$. Note that $\mu^{\prime}(\sigma_{\mathrm{N}})=0$ if $\mu$ is considered stress-independent.

For modeling the inelastic deformation of rough rock fractures, we adopt a plasticity-like framework, which is standard for physically-motivated models for rock joints (e.g. [53, 34]). The general form of the yield function may be written as

where $\phi_{\mathrm{b}}$ is the basic friction angle, $\psi$ the dilation angle, and $\omega$ is an abrasion coefficient which accounts for the impact of asperity damage on mobilized shear strength. Here we consider $\phi_{\mathrm{b}}$ a constant, while viewing $\psi$ and $\omega$ as state variables. Next, we consider the potential function of the following general form [53]

which gives the non-associative flow rule as

with $\lambda$ denoting the plastic multiplier. It is noted that the form of potential function (57) accommodates a stress-dependent dilation angle, while keeping the definition of the dilation angle as $\tan\psi=\differential u^{\mathrm{p}}_{f}/\differential v^{\mathrm{p}}_{f}$.

We now specialize the general framework to the constitutive model proposed by White [34] to match experimentally observed behavior of rough fractures. Under monotonic (non-cyclic) loading, the dilation is given by

where $\psi_{\mathrm{r}}$ denotes a residual dilation angle, and $b$ is a characteristic slip length. In the seated position of the fracture ($v^{p}_{f}=0$) the initial dilation angle is zero. With accumulating slip ($v^{p}_{f}>0$) asperities ride up over one another and the dilation grows towards a residual dilation angle. The abrasion coefficient controls the mobilized shear strength via Eq. (56), and is given by

where $\omega_{0}$ is an abrasion parameter controlling the peak shear strength, $c$ is a softening constant, and $\Lambda$ denotes a history variable that quantifies the degradation of roughness. When the abrasion is assumed isotropic in all directions, it can be equated to the cumulative plastic slip, i.e.

It is noted that in the phase-field formulation, one can calculate $\dot{v}^{\mathrm{p}}_{f}$ by taking time derivatives on both sides of Eq. (35). This model allows for an initial peak in strength that then degrades as slip accumulates and asperities are worn away, as often seen in rough fracture experiments. See White [34] for further details.

The standard continuous-Galerkin finite element method can be used to solve the two-field governing equations, Eqs. (4) and (5), in which the unknown fields are the displacement vector field $\bm{u}$ and the phase field $d$. We introduce spaces for the trial solutions as

where $H^{1}$ denotes a Sobolev space of order one. Accordingly, spaces for the weighting functions are defined as

Through the standard weighted residual procedure, we arrive at the variational form of the governing equations,

The finite element discretization of Eqs. (66) and (67) is standard and its details are skipped for brevity. When modeling non-propagating rock fractures, Eq. (67) needs to be solved only once for the initialization of the phase field. To simulate propagating fractures, however, both Eqs. (66) and (67) should be solved in every load step. In either case, Eq. (67) can be solved separately in the same manner as existing phase-field models, because it is common to solve the two governing equations in a staggered manner [44]. Therefore, in what follows, we focus on the solution of Eq. (66).

Given that the constitutive relation is nonlinear, we use Newton’s method to solve the problem at hand. At each Newton iteration, we solve the Jacobian system given by

where $\bm{\mathcal{J}}_{u}$ denotes the Jacobian matrix, $\Delta\bm{U}$ the nodal displacement increment, and $\bm{\mathcal{R}}_{u}^{h}$ the discretized version of Eq. (66). The specific expression of $\bm{\mathcal{J}}_{u}$ is

where $\bm{\eta}^{h}$ is the discretized version of $\bm{\eta}$, and $\mathbb{C}$ denotes the stress-strain tangent operator. To evaluate Eq. (68) at each Newton iteration, one needs to update the stress tensor and the tangent operator at every material (quadrature) point. An algorithm for updating these two quantities is designed below.

Algorithm 1 presents a detailed procedure to update the stress tensor and tangent operator at every material point. Quantities at the previous time step $t_{n}$ are denoted with subscript $(\cdot)_{n}$, whereas quantities at the current time step $t_{n+1}$ are written without any subscript to avoid proliferation of subscripts. The algorithm is essentially the same as the predictor–corrector algorithm in the phase-field method for frictional interfaces [21], except the following three modifications. First, we have introduced a boolean flag open flag to keep the material point in the open mode when its contact condition is identified to be open in a Newton iteration. In every load step, the flag is initialized as false at the beginning and switched to true if the contact condition becomes open during a Newton iteration. We have experienced that the use of this flag makes the Newton iteration more robust. Second, due to the nonlinearity of elastic fracture deformation, we have designed a local Newton iteration to evaluate the interface stress tensor $\bm{\sigma}_{f}$ and the interface elastic tangent operator $\mathbb{C}^{\mathrm{e}}_{f}$—see Algorithm 2. Third, we evaluate the yield function in a semi-implicit manner, using the normal stress in the previous step. This semi-implicit approach greatly simplifies the calculation of the derivative of the yield function and thus the inelastic fracture tangent operator $\mathbb{C}_{f}$. Lastly, we have added a small tolerance $k_{\text{tol}}$ (e.g. $10^{-3}$) to the weight of the bulk stress tensor in the calculation of the overall stress tensor (cf. Eq. (12)), to avoid ill-posedness of the matrix system [10]. Although such a tolerance is often found to be unnecessary in standard phase-field models, we have found that its use is critical to numerical robustness of the phase-field formulation at hand, particularly for intersecting cracks.

In Algorithm 1, when the trial stress violates the requirement $F\leq 0$, we use a return mapping algorithm to correct the trial stress and strains such that they satisfy $F=0$. The return mapping is performed in the full strain space and described in the following.

The unknowns in the return mapping are the six independent components of the strain tensor and the discrete plastic multiplier. Using Kelvin notation for an elegant representation of tensor algebra [54], we write the unknown vector as

The equations to be satisfied can be written as residuals

It is again noted that we evaluate the value of $F$ in the residual vector (71) in a semi-implicit way, using $\sigma_{\mathrm{N},n}:=\bm{\sigma}_{f,n}:(\bm{n}\operatorname{\otimes}\bm{n})$.

Then we use Newton’s method to find a numerical solution. At each Newton iteration, we solve

for the search direction $\Delta\bm{x}$. The Jacobian matrix is given by

The derivatives required to complete the Jacobian matrix are provided in A. It is noted that the tensorial operations above have been written in Kelvin notation.

Once the return mapping is completed, we calculate the inelastic fracture tangent operator defined as

To evaluate this inelastic tangent operator, here we adopt the method used in Bryant and Sun [24]. To begin, we consider the trial elastic strain as a variable and re-linearize Eq. (71) following Eq. (7.127) in de Souza Neto et al. [55]. This gives

We then insert $\differential\bm{\varepsilon}^{\mathrm{e}}=(\mathbb{C}^{\mathrm{e}}_{f})^{-1}:\differential\bm{\sigma}_{f}$ into the above equation and obtain

Rearranging Eq. (76) gives [56, 57]

where

Differentiating Eq. (78) with respect to the trial elastic strain gives the inelastic interface tangent operator as

The last term in the bracket, i.e. $\differential\Delta\lambda/\differential\bm{\varepsilon}^{\mathrm{e},\tr}$, can be calculated by substituting Eq. (78) into Eq. (77). This operation gives

Inserting the above equation into Eq. (80) gives the final form of the inelastic fracture tangent operator as follows:

We begin by simulating shearing of a single discontinuity, which is the most straightforward setting for studying the behavior of a rough rock fracture. We adopt the configuration of a shear test performed by Lee et al. [66], in which an elastic block under a constant normal stress slides along a wider rigid block fixed to the bottom boundary. The specific geometry and boundary conditions are illustrated in Fig. 3. To investigate the sensitivity of the numerical results to the phase-field length parameter, we prepare three different phase-field distributions initialized with $L=1.6$ mm, $0.8$ mm and $0.4$ mm, as shown in Fig. 4. We then apply a prescribed horizontal displacement on the boundary of the upper block, except for the region where $d>0$ (i.e. the diffusely approximated discontinuity) for which a displacement boundary condition cannot be imposed (see Bryant and Sun [24] for a similar treatment). The rigid block is meshed but its displacement degrees of freedom are fixed. The XFEM solution to this problem is obtained by embedding the discontinuity into a rectangular domain, similar to how other XFEM studies have modeled similar problems (e.g. [67, 68]). The simulation proceeds with a uniform displacement rate of $0.1$ mm until the total horizontal displacement reaches $10$ mm ($100$ steps in total).

We first verify the phase-field formulation by comparing its results with those obtained by the XFEM. To ensure that the phase-field result is accurate enough, we assign a sufficiently small length parameter with quite fine discretization, namely $L=0.4$ mm and $L/h=20$. Figure 5 compares the results of the phase-field method and XFEM in terms of the shear stress and dilation in the discontinuity. One can see that the two methods provide nearly identical results. The results show typical shear stress and dilation responses of a rough fracture undergoing shearing.

Next, we examine the mesh sensitivity of the phase-field method by repeating the same problem with three different levels of refinement: $L/h=5$, $L/h=10$, and $L/h=20$. The results are presented in Fig. 6. We observe very little sensitivity to the value of $L/h$, which indicates that a rather coarse refinement level of $L/h=5$ would be good enough for practical purposes. The results also confirm that the proposed approximation of the discontinuous strain (cf. Eqs. (23) and (24)), which relies on the $\Gamma$-convergence of the crack density functional, provides mesh-insensitive solutions.

Lastly, in Fig. 7 we investigate the sensitivity of the method to the phase-field length parameter. Here, the problem is repeated with three different length parameters, $L=1.6$ mm, $L=0.8$ mm and $L=0.4$ mm (illustrated in Fig. 4), with a fixed refinement level of $L/h=20$. It can be seen that the shear stress results converge as $L$ decreases and that the dilation results are almost the same even when $L$ is fairly large at 1.6 mm.