A modified damage and fracture phase field model considering heterogeneity for rock‐like materials

Abstract

Damage and fracture are the most extensive failure modes of rock materials, which may easily induce disaster and instability of engineering structures. This study developed a nonlocal damage fracture phase field model for rocks considering the heterogeneity of rocks. The modified phase field model introduced the heterogeneity of fracture parameters and modified the governing equations. Meanwhile, the free energy was constructed by the elastic strain energy sphere‐bias decomposition and the plastic strain energy. As for the numerical implementation, the three layers finite elements method structure was used in the frame of the finite element method. The ability of the modified phase field model has been illustrated by reproducing the experiment results of rock samples with pre‐existing cracks under compression.

INTRODUCTION

The damage fracture phase field model, which establishes a connection between damage and fracture mechanics, is gradually becoming a hot spot in solid fracture mechanics. Due to its advantages in analyzing discontinuous mechanical behavior, the traditional damage fracture phase field model is mainly used to study the fracture behavior of materials dominated by tensile cracks. However, it is less capable of analyzing the damage fracture behavior of nonhomogeneous rock materials under complex external loads. In geotechnical structures, rocks are mostly under complex compressive stresses and cracks show a composite extensional pattern of tension–compression–shear. Therefore, to describe the damage fracture behavior of rock materials, relevant preliminary studies have been conducted on the damage fracture phase field model of rock materials (Bryant & Sun 2018, 2021; Fei & Choo, 2020a, 2020b; You, Waisman, & Zhu, 2021; Zhang et al., 2017; Zhou & Zhuang, 2018; Zhou et al., 2018). Choo and Sun (2018) proposed a coupling model of damage phase field and plastic deformation based on the Drucker–Prager plastic model with compression cap, considering the influence of confining pressure and strain rate on the elastic-brittle and plastic flow of rock materials. Based on the Palmer–Rice fracture theory, Fei and Choo (2020a) proposed a phase field model of frictional shear fracture of geological materials by considering the frictional slip of geological shear fractures and deducing the threshold of fracture shear energy. Rock materials have a low resistance to tensile stress, and rock cracks present complex fractures mixed with tensile and shear fractures. Fan and Choo (2021) proposed a two-phase field fracture model, including cohesive tensile fracture and frictional shear fracture, to improve the rock fracturing mode of strain energy decomposition, potential energy calculation, and energy control. Mollaali et al. (2019) established a CO₂ fracture phase field model in the rock mass, which was improved by the exponential function relationship between rock permeability and the phase field value. Wang et al. (2020) improved the unified phase field model based on the unified tensile fracture criterion, enabling the improved phase field model to simultaneously simulate the tensile fracture and shear expansion of rock materials. According to their findings, the fracture mode conversion was realized by changing the material property variation. Schuler et al. (2020) analyzed the characteristics of rock crack growth under the chemical–mechanical coupling of rock porous media, improved the phase field fracture model by defining the coupling damage variables of chemical damage and mechanical damage, and quantified the influence of chemical environment on rock fracture behavior. Cao et al. (2020) established a three-dimensional numerical model of damage and fracture of heterogeneous rock materials. The homogenization method of fast Fourier transform technology was used to analyze the influence of heterogeneous rock materials on damage and fracture. Zhou et al. (2019) proposed a phase field model for complex cracking in rocks based on the elastic strain decomposition principle. The results of rock notch semicircle and Brazilian splitting tests were in good agreement with the results of phase field simulation analysis. To simulate quasi-static crack extension in rocks, an adaptive phase field method was proposed to consider crack extension due to compressive strain, and the results showed that the phase field model does not require a predetermined fracture criterion and the reduced global mesh refinement is highly adaptable to crack modeling. You et al. (2020) proposed an irreversible thermodynamic framework for a plastic-damage coupled phase-field fracture model and incorporated tensile damage variables and compressive damage variables into the modeling framework to achieve plastic-damage coupling by introducing the damage phase field into the yield function. Guével et al. (2020) developed a modified phase-field fracture model based on the full dissipative thermodynamic theory to characterize the effect of microstructure on the micromechanical response of rock materials and established a viscous phase field model for microstructure evolution under coupled chemical–mechanical effects. Wang et al. (2021) proposed a phase field simulation method for damage and fracture of rock-like materials considering materials' heterogeneity and spatial variability. The damage evolution characteristics of the cavity under the action of the temperature field, deformation field, and fluid pressure field were analyzed by establishing a phase field model of the rock body considering tensile cracking and compression-shear cracking.

In this study, a modified rock damage fracture phase field model was developed, and the plastic-free energy and heterogeneity were introduced to make the model more consistent with the characteristics of rock materials. The paper is structured as follows: Section 2 reviews the theoretical background of the damage fracture phase field model briefly; Section 3 presents the model modification and numerical implementation; Section 4 demonstrates the fracture characteristics of rocks with prefabricated cracks under compressive loading; and the last section concludes the whole study.

BASIC THEORIES

Phase field model

Francfort and Marigo (1998) topologized discontinuous sharp cracks within solid materials as crack bands with a certain width based on the elliptic regularization method of the Mumford–Shah generalization, using damage variables to define the damage phase field variables for regional gradient variations within the crack band. Expressed as a scalar d(x) with a value domain of [0,1], the undamaged state of the material is characterized using (d = 0) and the fully damaged state is characterized with (d = 1), as shown in Figure 1. $\mathit{\Omega }$ is solid configuration, $\partial \Omega$ is the boundary of $\mathit{\Omega }$, $\mathit{\Gamma }$ is the sharp crack, and l is the width of the dispersion crack.

[IMAGE OMITTED. SEE PDF]

According to the Dirichlet-type boundary condition, the damage phase field is the solution of the χ² differential equation, which takes the form shown in the following equation: 1$d(x)-l^{2}d″(x)=0.$

According to the Euler equation of the variational principle, the expression for the damage phase field can be obtained as 2$d=\mathrm{Arg}\left\{\inf \left[\frac{1}{2}{\int }_{\Omega }(d^2+l^2{d\prime }^2)\mathrm{d}V\right]\right\}.$

Miehe et al. (2010) defined the crack surface density function γ(d). The crack density function $\mathrm{\Gamma }(d)$ is related to the damage phase field d and the spatial gradient $\nabla d$, as shown in the following equation: 3$\Gamma (d)={\int }_{\mathit{\Omega }}\gamma (d,\nabla d)\mathrm{d}V=\frac{1}{2}{\int }_{\mathit{\Omega }}(d^2/l+l{\mid \nabla d\mid }^2)\mathrm{d}V,$where γ(d,∇d) denotes the crack surface density function; d and ∇d represent the crack damage phase field and its spatial gradient, respectively.

Specification of the free energy

The total free energy of the system can be expressed as the sum of elastic strain energy, crack surface energy, and plastic strain energy on the premise of ignoring the inertia effect of the system. 4$\psi (\mathit{u},d,\nabla d)={\int }_{\Omega }{\psi }_{\mathrm{ela}}(\mathit{u},d)\mathrm{d}V+{\int }_{\Omega }{\psi }_{\mathrm{sur}}(d,\nabla d)\mathrm{d}V+{\int }_{\Omega }{\psi }_{\mathrm{pla}}(\mathit{u},d)\mathrm{d}V,$where ${\psi }_{\mathrm{ela}}$ is elastic strain energy, ${\psi }_{\mathrm{sur}}$ is crack surface energy, and ${\psi }_{\mathrm{pla}}$ is plastic strain energy.

Elastic strain energy

The initial elastic strain energy of a rock material can be expressed as Gurtin (2008). 5${\psi }_{\mathrm{ela}}^0\left(\epsilon \right)=\frac{1}{2}\mathit{\epsilon }\left(u\right):\mathit{E}:\mathit{\epsilon }\left(u\right),$where $\mathit{\epsilon }\left(u\right)$ is the strain tensor and E is the elastic stiffness tensor.

Considering the effect of cracking on the free energy of the system, as the damage develops and the crack expands, the increase in the damage phase field leads to a decrease in the Helmholtz free energy function, and the energy degradation function g(d) is introduced to quantitatively characterize the effect of damage on the Helmholtz free energy (Miehe et al., 2015). 6$\psi (\epsilon ,d)=\psi (\epsilon ,g\left(d\right))=g\left(d\right){\psi }_{\mathrm{ela}}^0\left(\epsilon \right).$

To characterize the elastic strain energy under tensile and compressive loading, the elastic strain energy is decomposed into tensile and compressive strain energy. Defined as positive for compression and negative for tension, the elastic strain energy is shown in the following equation: 7${\psi }_{\mathrm{ela}}^0\left(\epsilon \right)={\psi }_{\mathrm{ela}}^0{\left(\epsilon \right)}^{+}+{\psi }_{\mathrm{ela}}^0{\left(\epsilon \right)}^{-}.$

The strain energy is decomposed using the sphere-bias decomposition method proposed by Amor et al. (2009). The tensile and compressive portions of the elastic strain energy density are shown in the following equation: ${\psi }_{\mathrm{ela}}^0{\left(\epsilon \right)}^{-}=\left(\frac{1}{2}{K}^{\mathrm{v}}{\epsilon }^{\mathrm{vol}}:{\epsilon }^{\mathrm{vol}}\right)({⟨{\epsilon }^{\mathrm{vol}}⟩}_{-}/{\epsilon }^{\mathrm{vol}}),$ 8$\begin{array}{ccl}{\psi }_{\mathrm{ela}}^0{\left(\epsilon \right)}^{+}& =& \left(\frac{1}{2}{K}^{\mathrm{v}}{\epsilon }^{\mathrm{vol}}:{\epsilon }^{\mathrm{vol}}\right)({⟨{\epsilon }^{\mathrm{vol}}⟩}_{-}/{\epsilon }^{\mathrm{vol}})\\ & & +\mu {\epsilon }^{\text{dev}}:{\epsilon }^{\text{dev}},\end{array}$where K^v is volume modulus, ε^vol is volume strain, ε^dev is deviator strain, and μ is Poisson's ratio.

The operators $⟨\cdot ⟩$ are defined as ${⟨x⟩}_{+}=\frac{\mid x\mid +x}{2},{⟨x⟩}_{-}=\frac{\mid x\mid -x}{2}.$

Crack surface energy

The energy required to produce diffuse cracks in a sharp crack topology is defined as the crack surface energy, expressed as a function of Griffith energy release rate and microcrack density as shown in the following equation: 9${\psi }_{\mathrm{sur}}(d,\nabla d)={\int }_{\Omega }{G}_{\mathrm{c}}\gamma (d,\nabla d)\mathrm{d}V,$where G_c is the critical energy release rate.

The key to the crack surface energy is the critical energy release rate and the local energy term, which can be obtained from the crack surface energy expression. The critical energy release rate is an inherent property of the material and varies significantly between crack types. For the local energy term, the classical phase field theory of the crack surface density function is shown in the following equation: 10$\gamma (d,\nabla d)=\frac{1}{2}\left(\frac{d^2}{l}+l{\mid \nabla d\mid }^2\right).$

Plastic strain energy

Based on the damage fracture phase field method, Molnár et al. (2020) studied the influence of material plasticity characteristics on instantaneous fracture toughness and proposed an expression for the plastic strain energy density associated with plastic deformation of the material; You, Waisman, Chen, et al. (2021) studied an elastic-plastic damage phase field model for quasibrittle geological materials considering the asymmetric properties of tension-compression. The plastic strain energy density includes the following hardening and each homogeneous hardening part, while the damage phase field parameter d is considered to affect only the plastic strain energy, independent of the variation of the elastic strain energy.

The specific equation for the plastic strain energy density is expressed as 11$\begin{array}{ccl}{\psi }_{\mathrm{pla}}^0\left({\epsilon }_{\mathrm{eq}}^{\mathrm{pl}},d\right)& =& g\left(d\right){\psi }_{\mathrm{pla}}^0\left({\epsilon }_{\mathrm{eq}}^{\mathrm{pl}}\right)\\ & =& g\left(d\right)\left(\frac{1}{2}H{\epsilon }_{\mathrm{eq}}^{\mathrm{pl}}:{\epsilon }_{\mathrm{eq}}^{\mathrm{pl}}+\frac{1}{2}h{\varsigma }^2\right),\end{array}$where H is the hardening modulus, ${\epsilon }_{\mathrm{eq}}^{\mathrm{pl}}$ is equivalent plastic strain, h denotes a nonnegative parameter related to each homogeneous hardening, and the internal variable $\varsigma$ is the cumulative plastic bias strain matrix parametrization.

Substituting Equations (8), (10), and (11) into (4), the following form can be obtained: $\begin{array}{l}\psi (\mathit{u},d,\nabla d)\\ =g\left(d\right)[\left(\frac{K^{\mathrm{v}}{\epsilon }^{\mathrm{vol}}:{\epsilon }^{\mathrm{vol}}}{2}\right)\left(\frac{{⟨{\epsilon }^{\mathrm{vol}}⟩}_{-}}{{\epsilon }^{\mathrm{vol}}}\right)(1+\mu {\epsilon }^{\text{dev}}:{\epsilon }^{\text{dev}})\\ +\left(\frac{H{\epsilon }_{\mathrm{eq}}^{\mathrm{pl}}:{\epsilon }_{\mathrm{eq}}^{\mathrm{pl}}+h{\varsigma }^2}{2}\right)]+{\int }_{\Omega }{G}_{\mathrm{c}}\gamma (d,\nabla d)\mathrm{d}V.\end{array}$

Governing equations

Based on the variational form of the total energy of the system, the strong form of the side value problem can be obtained by applying Gauss's theorem to obtain the macroscopic equilibrium equation as 12$\mathrm{div}\sigma +\mathit{b}=0,$ $\mathit{\sigma }\mathit{\cdot }\mathit{n}=\mathit{t},$where $\mathit{t}$ is external force, $\mathit{b}$ is body force, $\mathit{n}$ is external normal direction.

Miehe et al. (2010) introduced the concept of a history variable or history field, establishing a correspondence between the historical maximum energy and the phase field, specifically in the control equation, using the history variable H instead of the strain energy in the control equation, that is, 13$\mathit{H}=\max \left[{\psi }_{\mathrm{ela}}(\epsilon ,t)\right]t\in [0,t].$

The resulting governing equation for the damage fracture phase field is $\nabla \cdot \sigma +\mathit{b}=0,$ 14$\frac{2G_{\mathrm{c}}l}{\mathrm{\pi }}\Delta d-\frac{\partial g\left(d\right)}{\partial d}\mathit{H}-\frac{2G_{\mathrm{c}}}{\mathrm{\pi }l}(1-d)=0.$

Based on a critical energy release rate criterion, Shen and Stephansson (1994) proposed a linear damage criterion, namely the F-criterion, as shown in the following equation: 15$F=\frac{G_{\mathrm{I}}\left(\theta \right)}{G_{\mathrm{Ic}}}+\frac{G_{\mathrm{II}}\left(\theta \right)}{G_{\mathrm{IIc}}}.$

Based on the form of the F-criterion, the historical variable $\mathit{H}$ was divided into ${\mathit{H}}_1$ and ${\mathit{H}}_2$ (Fan & Choo, 2021). So the form of the governing equations for rock-like materials is obtained as follows: $\nabla \cdot \sigma +\mathit{b}=0,$ 16$d-l_{\mathrm{c}}^2\Delta d=l_{\mathrm{c}}\frac{\partial g\left(d\right)}{\partial d}\left(\frac{{\mathit{H}}_1}{G_{\mathrm{Ic}}}+\frac{{\mathit{H}}_2}{G_{\mathrm{IIc}}}\right).$

The critical energy release rate, G_c, is a key fracture mechanics parameter for crack extension and a determining factor for crack extension in the damage phase field of a solid material. Therefore, the assumption that the critical energy release rate of a rock satisfies the Weibull distribution is an effective means of characterizing the nonhomogeneity of rock. It is usually assumed that the parameters of rock mechanical properties conform to the Weibull distribution law in the form of a probability density function, as shown in the following equation (Gao et al., 2020): 17$P\left(C\right)=\frac{m}{\lambda }{\left(\frac{C}{\lambda }\right)}^{m-1}\exp \left[-{\left(\frac{C}{\lambda }\right)}^m\right],$where m denotes the shape parameter or shape factor, λ is the scale parameter, and C is a random variable.

The shape parameter m reflects the nonhomogeneous nature of the fine structure of the material. The smaller the value of m, the greater the dispersion of the fine structure of the material, and the greater the degree of nonuniformity of the material. Cracks in rocks under compression loading present a complex composite extension, including both tensile and compression-shear cracks. Therefore, the tensile crack inhomogeneity and compression-shear crack inhomogeneity are defined as the shape parameters (m_t and m_c) of the probability density function of the Weibull distribution, respectively. Therefore, the tensile crack nonhomogeneity m_t and compression-shear crack nonhomogeneity m_c are introduced to modify the governing equations, as shown in the following equation: $\nabla \cdot \sigma +\mathit{b}=0,$ 18$d-l_{\mathrm{c}}^2\Delta d=l_{\mathrm{c}}\frac{\partial g\left(d\right)}{\partial d}\left[{\left(\frac{{\mathit{H}}_1}{G_{\mathrm{Ic}}}\right)}^{1+\frac{1}{m_{\mathrm{t}}}}+{\left(\frac{{\mathit{H}}_2}{G_{\mathrm{IIc}}}\right)}^{1+\frac{1}{m_{\mathrm{c}}}}\right].$

According to the above equation, ${\mathit{H}}_1/G_{\mathrm{Ic}}$ and ${\mathit{H}}_2/G_{\mathrm{IIc}}$ are the driving terms of tensile crack propagation and compression-shear crack propagation, respectively. The nonhomogeneity parameters of m_t and m_c are introduced in governing equations, and their main effect is to increase the driving force of crack propagation to a certain extent.

NUMERICAL IMPLEMENTATION

Finite element discrete scheme

The corresponding finite element discretization format is derived for the weak form of the solid damage phase field control equation, which includes the displacement field and the damage phase field, and the corresponding variational form is approximated by the Gallione method. The displacement field $u(x)$ and damage phase field $d(x)$ approximation are estimated as $u\left(x\right)=N_u\hat{u},\epsilon \left(x\right)={\mathit{B}}_u\hat{u},$ 19$d\left(x\right)=N_d\hat{d},\nabla d\left(x\right)={\mathit{B}}_d\hat{d},$where $\hat{u}$ and $\hat{d}$ are the displacement field and damage phase field vectors at the nodes; N_u and ${\mathit{B}}_u$ are the shape functions and its derivative of displacement field; N_d and ${\mathit{B}}_d$ are the shape functions and its derivative of damage phase field.

The discrete format of the finite element corresponding to the weak form of the displacement field residual $r^u$ and damage phase field residual $r^d$ are 20$r^u={\int }_{\Omega }{N}_u^T\mathit{b}\mathrm{d}\Omega +{\int }_{\Omega }{N}_u^T\mathit{t}\mathrm{d}S-{\int }_{\Omega }{\mathit{B}}_u^T\sigma \mathrm{d}v,$ 21$r^d={\int }_{\Omega }\left[N_u^T\left(\frac{G_{\mathrm{c}}}{l_{\mathrm{c}}}d-2\frac{\partial g\left(d\right)}{\partial d}\mathit{H}\right)+4l_{\mathrm{c}}{\mathit{B}}_d^T\nabla d\right]\mathrm{d}\Omega \mathit{,}$

To solve for the displacement and damage phase fields at each operator step, the Newton–Raphson method is used to solve the nonlinear equations. 22$\left[\begin{array}{ll}{\mathit{K}}^{\mathit{uu}}& {\mathit{K}}^{\mathit{ud}}\\ {\mathit{K}}^{\mathit{du}}& {\mathit{K}}^{\mathit{dd}}\end{array}\right]\left[\begin{array}{l}d\\ u\end{array}\right]=-\left[\begin{array}{l}{r}^u\\ r^d\end{array}\right],$where stiffness matrix ${\mathit{K}}^{\mathit{uu}}, {\mathit{K}}^{\mathit{ud}}, {\mathit{K}}^{\mathit{du}}, {\mathit{K}}^{\mathit{dd}}$ are expressed, respectively, as ${\mathit{K}}^{\mathit{uu}}={\int }_{\Omega }{\mathit{B}}^T\left(\frac{\partial \sigma }{\partial \epsilon }\right)\mathit{B}\mathrm{d}\Omega \mathit{,}\mathit{ }{\mathit{K}}^{\mathit{ud}}={\int }_{\Omega }{\mathit{B}}^T\left(\frac{\partial \sigma }{\partial d}\right)N\mathrm{d}\Omega \mathit{,}$ $\begin{array}{l}{\mathit{K}}^{\mathit{du}}={\int }_{\Omega }{N}^T\left(\frac{\partial g\left(d\right)}{\partial d}\frac{\partial \psi }{\partial \epsilon }\right)\mathit{B}\mathrm{d}V,\\ {\mathit{K}}^{dd}={\int }_{\Omega }{\mathit{N}}^T\mathit{N}\left(1+4\frac{l_{\mathrm{c}}}{G_{\mathrm{c}}}H\right)+4l_{\mathrm{c}}^2\mathit{B}{}^T\mathit{B}\mathrm{d}V.\end{array}$

TL-FEMS numerical calculation procedure

The coupled equations of the displacement field and damage phase field in the damage phase field model for solid materials are solved numerically using the finite element method. The open source and application of the damage fracture phase field algorithm for solids based on ABAQUS software developments have been promoted by scholars such as Wu et al. (2020), Seleš et al. (2019), and Martínez-Pañeda et al. (2018). In terms of ABAQUS software development, numerical algorithms such as the two-layer finite element structure proposed by Fang et al. (2019) and the three-layer finite element structure proposed by Jeong et al. (2018) and Molnár et al. (2020) have shown that the layered finite element structure has a stable and powerful computational capability. Based on the advantages of the layered finite element structure, a finite element structure containing multiple layers of user-defined subroutines is used to provide a base program for subsequent rock damage fracture phase field models. The ABAQUS-based user-defined element subroutine (UEL) solves the displacement and damage phase fields of the solid phase-field damage model using the three layers finite elements method structure (TL-FEMS) numerical calculation program, which includes a displacement field (UEL subroutine), a damage phase field (UEL subroutine) and a layer for the display of numerical results only (UMAT subroutine). The FEMS numerical calculation program consists of a first layer for solving the displacement field of the solid structure, a second layer for analyzing the damage phase field of the solid structure, and a third layer for displaying the numerical calculation results. This is shown in Figure 2.

[IMAGE OMITTED. SEE PDF]

NUMERICAL EXAMPLES

Single crack condition

To study the fracture characteristics of rock damage under different prefabricated crack inclination conditions, numerical models of rock containing a single crack are established. The size of numerical samples is 152.4 mm × 76.2 mm, the crack length is 12.7 mm, and the crack inclination is set to 30°, 45°, and 60°, as shown in Figure 3. The model meshes with a 0.1 mm grid encapsulation at the tip of the crack and a feature-length l_c of two times the grid encapsulation size. The model material parameters are determined according to Bobet and Einstein (1998), with Young's modulus of 4.26 GPa, Poisson's ratio of 0.21, the tensile strength of 3.2 MPa, and compressive strength of 32 MPa, m_t = m_c = 3 and critical energy release rates of 16 and 205 J/m² for type I and II cracks, respectively.

[IMAGE OMITTED. SEE PDF]

Figure 4 shows the results of the numerical analysis of single prefabricated crack expansion under compression loading. With the application of displacement at the top of the model, the crack tip first appears as a tensioned wing-type inflection fold crack, with the direction of expansion parallel to the direction of the maximum main compression load. The length of the tensioned wing-type inflection fold crack gradually decreases as the prefabricated crack inclination angle increases, and the crack expansion has a significant self-similarity feature when the prefabricated crack inclination angle is 60°. Continuous displacement loading enlarges the initial inflection crack area, and a local damage area, namely, the fracture process area, appears in the crack tip area; the secondary cracks are tension cracks, which sprout in the fracture process area at the crack tip, and the expansion direction is identical with the maximum compression load direction.

[IMAGE OMITTED. SEE PDF]

According to the results of indoor mechanical tests on single prefabricated cracked rocks (Wong & Einstein, 2009a, 2009b), local rock spalling occurs at the crack tip under compression shear, and the fracture process zone is formed at the crack tip, as shown in Figure 5, which becomes the breeding area for secondary cracks. Compared with the results of rock mechanics tests, the rock damage fracture phase field model can numerically reproduce the crack extension characteristics of the rock under compressive loading, mainly including tension wing crack extension, fracture process zone formation, and secondary tension crack development.

[IMAGE OMITTED. SEE PDF]

The crack propagation behaviors of double cracks were studied. Figure 6 shows the position features of double cracks.

[IMAGE OMITTED. SEE PDF]

Colinear double cracks

The numerical model geometry of the colinear double crack is 152.4 mm × 76.2 mm, the crack dimension is 2a = 12.7 mm and the rock bridge lengths are w = 0 and c = a = 6.35 mm. The crack inclination is set to 45° and the material parameters are the same as in the case of single precast cracked rock. Figure 7 shows the results of the numerical analysis of the colinear double crack extension under compression loading with displacements of 0.38, 0.62, and 0.74 mm.

[IMAGE OMITTED. SEE PDF]

With the gradual application of axial displacement, the tensioned wing inflection crack initiates in the prefabricated crack tip. There is a significant difference in the expansion law of prefabricated outer end crack and inner end crack. The expansion law of the outer end crack is identical to that of a single crack, with the expansion direction parallel to the direction of the maximum main compression load. Crack expansion at the inner end of the crack is inhibited, and there is some kind of resistance to lock the crack tip wing rupture, known as the “interlocking effect.” As the loading continues, crack expansion would break through the interlock, resulting in crack expansion through the rock bridge, and eventually causing the overall destruction of the rock. The results of the damage fracture phase field method are in line with the test results, indicating that the damage fracture phase field model can describe the prefabricated crack expansion law. The common line double crack expansion law in the initial rupture process is the same as the single crack expansion law, but the crack interaction between the inner ends of the common line double crack is complex; the interlocking effect between the inner ends leads to short-time expansion stagnation.

Noncollinear double cracks

To investigate the fracture characteristics of noncoincident double-cracked rock damage, numerical models were established with the w = c = a = 6.35 mm and α = 45° and 60°. The material parameters were consistent with the previous model. Figures 8 and 9 show the numerical and experimental results of the nonconcurrent double crack extension under compression loading (Bobet & Einstein, 1998; Wong, 2009a, 2009b).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Winged-inflection cracks appear first at the nonconcurrent double-cracked ends. They are in the same direction as the maximum main compressive load, with the winged cracks gradually expanding as the applied displacement increases. Fracture process zones and compression-shear cracks develop in the 45° nonconcurrently cracked rock bridge area under compression loading, with both inner ends developing in both directions and eventually penetrating. Cracks in the area of the 60° noncongruent break rock bridge under compression loading, with the fracture process zone and compression-shear cracks developing mainly in the inner end of the upper crack and eventually merging with and penetrating the inner end of the lower crack.

Based on the results of the numerical analysis of 45° nonconcurrent double cracks under compression loading, the influence of different rock bridge distances on the expansion characteristics of nonconcurrent double cracks was analyzed. With the working conditions as follows: w = a = 6.35 mm, c = 2a = 12.7 mm, and w = 2a = 12.7 mm, c = 2a = 12.7 mm, the analytical results are shown in Figures 10 and 11.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Under the compressive displacement load, the upper and lower crack tips show tensioned wing inflection crack extension. Still, the development of the wing cracks is weaker than the outer tip wing cracks of the two crack tips in the rock bridge region, and the fracture process zone. Compression shear damages are developed in the rock bridge region. The cracks eventually converge and penetrate the rock bridge region. By comparing with the test results, the rock damage fracture phase field model can describe the crack extension characteristics well and the results are in good agreement. Compared with the classical phase field method, the modified model can better simulate the fracture process zone of the crack tip, especially the crack in the rock bridge region.

CONCLUSIONS

The phase field method was applied to study crack propagation and coalescence in rock-like materials. In this study, a modified phase field model was proposed for simulating crack propagation under compression. The proposed model introduced the heterogeneity parameters on the critical energy release rates for mode I and mode II cracks. Further, the contribution of plastic strain energy to the shear crack was also taken into account.

The modified phase field model was implemented in the framework of the finite element method, and developed by TL-FEMS. The ability of the modified phase field model has been illustrated by reproducing the experiment results of rock samples with pre-existing cracks under compression. The modified phase field model successfully reproduces the wing cracks, the fracture process zone at the tip of the crack, and the typical type of crack coalescence at the rock bridge region.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.

References

Amor H, Marigo JJ, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids. 2009;57(8):1209‐1229.

Bobet A, Einstein HH. Fracture coalescence in rock‐type materials under uniaxial and biaxial compression. Int J Rock Mech Min Sci. 1998;35(7):863‐888.

Bryant EC, Sun W. A mixed‐mode phase field fracture model in anisotropic rocks with consistent kinematics. Comput Methods Appl Mech Eng. 2018;342:561‐584.

Bryant EC, Sun W. Phase field modeling of frictional slip with slip weakening/strengthening under non‐isothermal conditions. Comput Methods Appl Mech Eng. 2021;375: [eLocator: 113557].

Cao YJ, Shen WQ, Shao JF, Wang W. A novel FFT‐based phase field model for damage and cracking behavior of heterogeneous materials. Int J Plast. 2020;133(10): [eLocator: 102786].

Choo J, Sun WC. Coupled phase‐field and plasticity modeling of geological materials: from brittle fracture to ductile flow. Comput Methods Appl Mech Eng. 2018;330(3):1‐32.

Fan F, Choo J. Double‐phase‐field formulation for mixed‐mode fracture in rocks. Comput Methods Appl Mech Eng. 2021;376(1): [eLocator: 113655].

Fang J, Wu C, Rabczuk T, et al. Phase field fracture in elasto‐plastic solids: Abaqus implementation and case studies. Theor Appl Fract Mech. 2019;103: [eLocator: 102252].

Fei F, Choo J. A phase‐field method for modeling cracks with frictional contact. Int J Numer Methods Eng. 2020a;121(4):740‐762.

Fei F, Choo J. A phase‐field model of frictional shear fracture in geologic materials. Comput Methods Appl Mech Eng. 2020b;369: [eLocator: 113265].

Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids. 1998;46(8):1319‐1342.

Gao C, Zhou Z, Li Z, Li L, Cheng S. Peridynamics simulation of surrounding rock damage characteristics during tunnel excavation. Tunn Undergr Space Technol. 2020;97: [eLocator: 103289].

Guével A, Rattez H, Veveakis E. Viscous phase‐field modeling for chemo‐mechanical microstructural evolution: application to geomaterials and pressure solution. Int J Solids Struct. 2020;207:230‐249.

Gurtin ME. A theory of grain boundaries that accounts automatically for grain misorientation and grain‐boundary orientation. J Mech Phys Solids. 2008;56(2):640‐662.

Jeong H, Signetti S, Han TS, Ryu S. Phase field modeling of crack propagation under combined shear and tensile loading with hybrid formulation. Comput Mater Sci. 2018;155:483‐492.

Martínez‐Pañeda E, Golahmar A, Niordson CF. A phase field formulation for hydrogen assisted cracking. Comput Methods Appl Mech Eng. 2018;342(12):742‐761.

Miehe C, Mauthe S, Teichtmeister S. Minimization principles for the coupled problem of darcy–biot‐type fluid transport in porous media linked to phase field modeling of fracture. J Mech Phys Solids. 2015;82:186‐217.

Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase‐field models of fracture: variational principles and multi‐field FE implementations. Int J Numer Methods Eng. 2010;83(10):1273‐1311.

Mollaali M, Ziaei‐Rad V, Shen Y. Numerical modeling of CO₂ fracturing by the phase field approach. J Nat Gas Sci Eng. 2019;70: [eLocator: 102905].

Molnár G, Gravouil A, Seghir R, Réthoré J. An open‐source Abaqus implementation of the phase‐field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation. Comput Methods Appl Mech Eng. 2020;365(6): [eLocator: 113004].

Schuler L, Ilgen AG, Newell P. Chemo‐mechanical phase‐field modeling of dissolution‐assisted fracture. Comput Methods Appl Mech Eng. 2020;362: [eLocator: 112838].

Seleš K, Lesičar T, Tonković Z, Sorić J. A residual control staggered solution scheme for the phase‐field modeling of brittle fracture. Eng Fract Mech. 2019;205:370‐386.

Shen B, Stephansson O. Modification of the G‐criterion for crack propagation subjected to compression. Eng Fract Mech. 1994;47(2):177‐189.

Wang M, Yu Z, Jin Y, Shao J. Modeling of damage and cracking in heterogeneous rock‐like materials by phase‐field method. Mech Res Commun. 2021;114(6): [eLocator: 103612].

Wang Q, Feng YT, Zhou W, Cheng Y, Ma G. A phase‐field model for mixed‐mode fracture based on a unified tensile fracture criterion. Comput Methods Appl Mech Eng. 2020;370(15): [eLocator: 113270].

Wong L, Einstein HH. Systematic evaluation of cracking behavior in specimens containing single flaws under uniaxial compression. Int J Rock Mech Min Sci. 2009a;46(2):239‐249.

Wong L, Einstein HH. Crack coalescence in molded gypsum and carrara marble: part 1. macroscopic observations and interpretation. Rock Mech Rock Eng. 2009b;42(3):475‐511.

Wu JY, Huang Y, Nguyen VP. On the BFGS monolithic algorithm for the unified phase field damage theory. Comput Methods Appl Mech Eng. 2020;360(3): [eLocator: 112704].

You T, Waisman H, Chen WZ, Shao JF, Zhu QZ. A novel micromechanics‐enhanced phase‐field model for frictional damage and fracture of quasi‐brittle geomaterials. Comput Methods Appl Mech Eng. 2021;385: [eLocator: 114060].

You T, Waisman H, Zhu QZ. Brittle‐ductile failure transition in geomaterials modeled by a modified phase‐field method with a varying damage‐driving energy coefficient. Int J Plast. 2021;136: [eLocator: 102836].

You T, Zhu QZ, Li PF, Shao JF. Incorporation of tension‐compression asymmetry into plastic damage phase‐field modeling of quasi brittle geomaterials. Int J Plast. 2020;124:71‐95.

Zhang X, Sloan SW, Vignes C, Sheng D. A modification of the phase‐field model for mixed mode crack propagation in rock‐like materials. Comput Methods Appl Mech Eng. 2017;322(Aug 1):123‐136.

Zhou S, Zhuang X. Adaptive phase field simulation of quasi‐static crack propagation in rocks. Undergr Space. 2018;3(3):190‐205.

Zhou S, Zhuang X, Rabczuk T. Phase field modeling of brittle compressive‐shear fractures in rock‐like materials: a new driving force and a hybrid formulation. Comput Methods Appl Mech Eng. 2019;355(10):729‐752.

Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theor Appl Fract Mech. 2018;96:174‐192.