<figure><img src="image_1.png" /></figure>

HAL
open science

# An adaptive phase field approach to 3D internal crack growth in rocks

Bin Xu, Tao Xu, Michael Heap, Alexandra R.L. Kushnir, Bo-Yi Su, Xiao-Cong Lan

► To cite this version:

Bin Xu, Tao Xu, Michael Heap, Alexandra R.L. Kushnir, Bo-Yi Su, et al. An adaptive phase field approach to 3D internal crack growth in rocks. Computers and Geotechnics, 2024, 173, pp.106551. 10.1016/j.compgeo.2024.106551. hal-04655075

HAL Id: hal-04655075

https://hal.science/hal-04655075v1

Submitted on 12 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# An adaptive phase field approach to 3D internal crack growth in rocks

Bin Xu<sup>a</sup>, Tao Xu<sup>a,*</sup>, Michael J. Heap<sup>b,c</sup>, Alexandra R.L. Kushnir<sup>d</sup>, Bo-yi Su<sup>a</sup>, Xiao-cong Lan<sup>a</sup>

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

**Keywords:**

Phase field method

3D crack

Crack growth

Wing crack

Adaptive mesh

## ABSTRACT

Modeling the 3D internal crack under compression entails complex fracture mechanics (mode I-II-III fracture), resulting in substantial computational costs and challenges in characterizing fracture morphology characterization for Phase Field Method (PFM) simulations. In the present paper, we developed a 3D adaptive crack propagation model integrated into the PFM to study the internal cracks in rocks. The proposed locally refined solution was implemented within the isogeometric analysis (IGA) framework, utilizing PHT-splines and employing the phase field variable as an error indicator to accurately capture non-smooth crack surfaces. The simulation begins with a coarser mesh, and the spatial discretization is adaptively refined in real-time as the computation progresses and the phase variable threshold is reached. By comparing our results with previous numerical tests, we validated the effectiveness of our model in simulating crack growth. Notably, a smaller length scale parameter led to a higher peak force in notched square plates under tension, and our adaptive PFM successfully reduced errors caused by length scale parameters. We further validated our model by comparing it with uniaxial compression experiments on 3D-printed resin samples with internal penny-shaped cracks, achieving good agreement with the experimental data. Using the 3D adaptive PFM, we predicted changes in crack morphology and fracturing mechanisms by varying the inclination angles of internal penny-shaped cracks. Our numerical results showed that as the inclination angle increased, wing cracks initiated closer to the end tips of the original crack, and their length decreased due to the wings wrapping effect. The corresponding fracturing mechanisms transition from mode I-III fracture (0°) to mode I-II-III fracture (30°, 45°, 60°), and eventually to mode I fracture (90°) as inclination angles increase. Comparing 2D and 3D models, we found that the growth of wing cracks in 3D was significantly influenced by the intermediate principal stress. We conclude that the extended 3D adaptive PFM-based model accurately predicted complex mechanical behavior and crack growth patterns in rock.

## 1. Introduction

Internal cracks are prevalent within the rock mass of various underground structures, such as those associated with hydrocarbon production (Chang et al., 2024), nuclear waste disposal (Li et al., 2017; Shao et al., 2024), tunnel stability (Li, 2024; Wang et al., 2022a), geothermal energy (Wasantha et al., 2021), mining (Lu et al., 2022; Xu et al., 2015), slope stability (Li et al., 2024; Yu et al., 2023; Yuan et al., 2021), and other rock engineering operations (Fu et al., 2020b; Xu et al., 2006). These internal cracks can trigger the formation and propagation of new cracks, leading to damage accumulation that can trigger spalling, failure, and collapse. In rock engineering projects, such collapses can result

in significant rock engineering disasters (Li, Xu, et al., 2024; Xu et al., 2014). Therefore, a comprehensive understanding of the behavior of internal cracks in rock masses is crucial.

Numerous experimental, theoretical, and numerical studies have been conducted to understand the failure of fractured rock under uniaxial compression. For example, previous research focused on the initiation and growth of 2D cracks has demonstrated that, in pre-cracked rock samples, wing cracks initiate and propagate towards the applied load, exhibiting virtually unrestrained and unregulated growth (Wong and Einstein, 2009; Xu et al., 2024; Xue et al., 2021; Yang, 2011). Under substantial compressive stresses, even the propagation of a single crack alone is sufficient to cause rock failure (Fig. 1(a) and Fig. 1(b)).

<footnote><sup>a</sup> Center for Rock Instability and Seismicity Research, Northeastern University, Shenyang 110819, China</footnote>

<footnote><sup>b</sup> Université de Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, 5 rue René Descartes, Strasbourg F-67084, France</footnote>

<footnote><sup>c</sup> Institut Universitaire de France (IUF), 1 rue Descartes, Paris 75231, France</footnote>

<footnote><sup>d</sup> Rock Physics and Geofluids Laboratory (RPGL), IIC, ENAC, École Polytechnique Fédérale de Lausanne, Station 18 CH-1015, Lausanne, Switzerland</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: xutao@mail.neu.edu.cn (T. Xu).</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 1. 2D wing crack pattern analysis in maximum principal strain field using (a) digital image correlation (Xu et al., 2022a) and (b) phase field method where a = 1 indicates the fully broken state in red, while a = 0 indicates the intact state in blue (Xu, Xu, et al., 2022a), compared with (c) the morphology of 3D experimental wing crack propagation (Dyskin, Sahouryeh, et al., 2003).</figcaption></figure>

However, the majority of existing literature has predominantly focused on the initiation and growth of internal 2D cracks, leaving a notable research gap in our understanding of internal 3D cracks.

Studying the growth of 3D internal cracks within rocks provides critical insights into the mechanical mechanism of crack expansion under true triaxial stress conditions (Wang et al., 2023a; Zhang et al., 2023a). For example, Gu et al. (2023) discovered that the 3D stress state determined the direction of rock crack growth, concluding that inhibiting tensile crack growth is crucial for controlling disasters in deep hard rock engineering. Jiang et al. (2024) investigated the failure mechanics of naturally jointed granite, finding that tensile cracks dominate the early stages of failure. Li and Fan (2024) analyzed the shear band formation and found that the inhomogeneity of the stress-strain relationship is a function of the confining pressure. Because of the difficulty observing 3D crack growth in opaque rock samples, researchers have used various transparent materials to better understand 3D crack growth behavior. A noteworthy investigation, conducted by Adams and Sines (1978), used PMMA (polymethyl methacrylate) to study crack

propagation patterns, identifying three distinct types of cracks: wing cracks, petal cracks, and fish fin cracks. Dyskin et al. (2003) conducted constant strain rate tests on transparent samples made of casting resin, cement, and mortar, each containing a single internal crack (Fig. 1 (c)). Their findings revealed significant differences in the growth of 3D cracks compared with that of 2D cracks under compression, specifically noting the constraints imposed by the length of the pre-existing cracks. Expanding on this research, Yin et al. (2014) and Zhou et al. (2018) conducted further experiments involving two 3D pre-existing cross-embedded flaws in PMMA. In their study, they observed the occurrence of anti-wing cracks and oblique secondary cracks, and emphasized the pivotal role of petal cracks in crack interactions. The investigation of 3D internal cracks using experimental methods continues to receive significant attention (Cui and Wong, 2023; Wang et al., 2022c). Concurrently, numerical methods have been employed as a complementary approach to enhance our understanding of internal crack mechanisms.

Numerical simulations have gained significant popularity for their ability to deliver accurate results quickly and predictions for various

<figure><img src="image_6.png" /><figcaption>Fig. 2. The approximate expression of the discrete crack in the phase field method (PFM). (a) Discrete crack. (b) Phase field crack. $a = 1$ indicates the fully broken state (i.e. red), while $a = 0$ indicates the intact state (i.e. blue).</figcaption></figure>
scenarios, reducing the dependence on time-consuming experiments, particularly in 3D scenarios (Carter et al., 2000). Different numerical methods address three-dimensional fracture problems, which can be categorized into mesh-based and meshfree methods. In mesh-based numerical methods, such as the Finite Element Method, and Finite Difference Method, re-meshing is normally required to simulate crack initiation, propagation, and coalescence. For example, Fu et al. (2017) employed an elastic-plastic model in FLAC3D to simulate 3D internal crack propagation, aiming to replicate experimental findings. However, these authors noted the requirement for superfine element divisions to achieve adequate resolution. Zhang et al. (2012) employed a 3D element partition method to simulate 3D cracks, addressing challenges associated with configuring joint elements and modifying the mesh scheme. Zhou et al., (2020a) proposed a new 3D numerical mesoscale model to describe the deformation and fracturing of brittle rocks that accounts for material heterogeneity and local material degradation. Zhang et al. (2020) used the FRANC3D software to simulate 3D cracks propagation and applied the M-integral method to calculate the mixed mode fracture parameters ($K_I, K_{II}, K_{III}$) at the crack tip. Conversely, meshfree methods (Zhou et al., 2014) offer advantages in modeling fracture problems by eliminating the need for re-meshing. For example, peridynamics (PD) (Wang et al., 2018a; Wang et al., 2017; Wang et al., 2018b) hold promise in addressing fundamental complexities in fracturing while reducing the reliance on phenomenological modeling assumptions and numerical tuning parameters. Shou et al. (2019) introduced an extended non-ordinary state-based PD for 3D crack simulation, accurately replicating failure patterns but yielding non-smooth crack surfaces. Tian and Zhou (2021) extended the PD to study 3D transversely isotropic and further proposed a general particle dynamics method (Yao et al., 2022). Compared to mesh-based numerical methods like the phase field method (PFM), determining suitable discretization parameters for PD, such as the optimal ratio between horizon and mesh size, can be challenging due to the sensitivity of the simulations. On the other hand, mesh-free simulations, including PD, require a computational framework that predicts both the initial position and direction of crack growth, as well as determines the crack surfaces in 3D conditions. Consequently, the numerical simulation of 3D internal crack propagation presents challenges in both mesh processing and computational efficiency. These difficulties have led to limited research on the numerical simulation of rock samples with 3D internal cracks subjected to compression.

In recent years, the phase field method (PFM) has gained increasing attention in the context of crack propagation (Liu et al., 2023; Wang et al., 2023b; Zhuang et al., 2022). PFM employs a scalar field, such as the phase field, to represent the area of material damage (Fig. 2), eliminating the need for an external failure criterion to predict crack growth. Fei et al., (2022a) incorporated the roughness of the crack surface into the PFM and further demonstrated its capability to transition from discontinuous to continuous modeling, with applications in fault rupture processes (Fei et al., 2022b). These findings highlight the potential of PFM for diverse geotechnical environments (Choo et al., 2021). To accurately model rock failure processes, different energy decomposition methods are employed within the PFM framework. For example, the volumetric and deviatoric strain split is considered in the work of Liu et al. (2022) and further used several classic examples to show the ability of the present framework in the uniaxial compression state. Additionally, factors such as cohesion and internal friction angles (Xu et al., 2022b) can be incorporated in the double PFM framework. Choo et al. (2023) introduced an innovative double-phase-field approach that uses different phase variables to characterize tensile or shear failure processes in rocks. However, most conventional PFMs are based on the assumptions of 2D plane stress or strain conditions, as highlighted in existing literatures (Cui et al., 2023; Zhuang, Zhou, et al., 2022). While these 2D models serve their purpose in specific scenarios (Xu, Xu, et al., 2024), they have limitations in accurately predicting non-planar fractures and the intricate interplay between cracks within

complex networks. As a result, they often yield underestimated or incorrect assessments of mechanical performance. There is a compelling need for 3D modeling to provide a more realistic representation of materials and structures responses to external loads. Such an approach offers a comprehensive understanding of behavior beyond the constraints of conventional 2D models. One significant advantage of PFM is its seamless extension from 2D to 3D within the standard finite element framework, making it a particularly attractive method for studying crack propagation in 3D. Francfort and Marigo (1998) introduced a variational model based on energy minimization principles, deviating from the conventional constraints associated with pre-existing cracks and well-defined crack paths. Subsequently, Bourdin et al. (2008) improved this model by incorporating a phase-field variable to regularize the sharp interface, further validated in subsequent studies (Borden et al., 2014; Goswami et al., 2020). These studies demonstrated convergence from the smooth crack surface to the sharp crack surface as the regularization parameter $l_0$ approaches zero. Recently, Nguyen-Thanh et al. (2020) formulated a groundbreaking PFM, using the isogeometric meshfree collocation method (IMCM) to simulate 3D polycrystalline material crack extension. Zhou et al., (2020b) used the PFM to model the hydraulic fracturing while considering the effect of the initial stress field. Wu et al. (2021) introduced a phase-field regularized cohesive zone model (PF-CZM), and the use of this innovative PF-CZM in the context of large-scale engineering structures is remarkably straightforward. Furthermore, Khalil et al. (2022) integrated AT1, AT2, and PF-CZM components into the phase field theory framework to model and predict the deformation and mechanical behavior of metals. Most of the fracturing mechanics discussed in current literature are driven by mode-I fracture, rather than more complex fracture mechanics, such as mixed-mode I-II-III fracture. In conclusion, the 3D PFM remains a challenging yet promising approach for analyzing rock deformation and cracking.

Although the PFM possesses several merits, a notable drawback stems from its high computational workload, requiring an exceptionally fine mesh discretization to accurately represent the phase-field profile. This requirement is not exclusive to phase field models but extends to various other damage or fracture models as well. Due to the length scale parameter employed in the PFM, the size of the damaged area becomes particularly significant relative to the grid size. Therefore, it is crucial to improve the computational efficiency for the practical application of PFM. The adoption of adaptive analysis has become a widely embraced and accepted method to address this challenge. In recent years, various adaptive analysis approaches for PFM in different discretization spaces have been proposed. The process of selecting a reasonable error indicator has become a crucial aspect in adaptive analysis approaches. It plays a vital role in ensuring the effectiveness and accuracy of the adaptive analysis approaches employed. For example, Borden et al. (2012) initially integrated PFM with local adaptive refinement, using the critical value of the phase-field variable as an indicator. Li et al. (2020) enhanced the refinement scheme by incorporating both the phase field and its gradient, specifically addressing crack resolution in a 3D polycrystalline material. Yue et al. (2022) employed a strategy of refining the meshes near the crack tips by adjusting the mesh size, thereby enhancing the precision and accuracy of the simulation in the vicinity of crack propagation. Freddi and Mingazzi (2022) proposed a novel approach wherein the strain energy is used to precisely approximate the topology of cracks, offering an advanced and accurate method for characterizing crack structures in the analyzed materials. Paul et al. (2020) presented an adaptive space-time PFM which is a discretized generalized-$\alpha$ method based on locally refinable non-uniform rational B-splines (LRNURBS). Hirshikesh et al. (2019) used the recovery-based error indicator in combination with quad-tree decomposition, revealing that this method achieved a lower computational cost compared with classical numerical examples. Hirshikesh et al. (2021) also used stability analysis in conjunction with quad-tree decomposition to refine the mesh, eliminating the need for post-processing techniques
in the refinement of fracture propagation. Xing et al. (2023) employed a quad-tree mesh refinement scheme, incorporating the combination of phase-field value and elemental strain energy density, to model fracture in hyperelastic composites. Walloth and Wollner (2022) highlighted the importance of the residual-type a posteriori error in classic test cases, such as the single edge notched tension test, shear test, and the L-shaped panel test. On one hand, most of these adaptive algorithms require prior knowledge of the crack path, which may not always be available due to the complexities involved in pre-processing and analysis, as well as the potential introduction of some degree of grid distortion. On the other hand, an appropriate and efficient adaptive program is crucial for accurately characterizing the complex fracture morphology (mode-I-II-III fracture), of a 3D internal crack under compressive stress. In contrast to the selection of other basis equations, we propose PHT-splines basis functions to approximate displacement and phase fields. PHT splines, as an extension of NURBS, provide the capability of local refinement and the presence of hanging nodes. PHT-splines also perform better in handling surfaces with irregular geometries, especially in 3D processing of complex crack surfaces generated by internal cracks (mode-I-II-III fracture) in the present paper, because they can dynamically adjust the level of refinement based on the surface characteristics. We also introduce an adaptive h-refinement scheme aimed at keeping refinement localized around singularities and regions of high gradients. Meanwhile, considering the various error indicators used in the PFM, it is challenging to definitively determine the superior choice. After careful deliberation, we have chosen the phase-field variable as the primary judgment criterion in our study. This selection is based on its inherent ability to provide a comprehensive representation of PFM systematic behavior, capturing intricate phase transitions and ensuring a nuanced assessment of model performance. The versatility and capacity of the phase-field variable to encapsulate complex phenomena make it a robust and reliable choice for guiding our evaluation process.

To accurately represent the non-smooth crack surfaces and twisting crack patterns resulting from mixed torsion and bending loads observed in internal crack tests (Fu et al., 2016; Wang, Li, et al., 2022c; Zhang et al., 2023b; Zhang, Yu, et al., 2020), comprehensive 3D simulations are necessary. In this paper, we aim to investigate the complex fracture behaviors of a 3D internal crack under uniaxial compression, along with the impact of inclination angles on 3D fracture morphology and fracturing mechanisms. Additionally, to address the efficiency challenges associated with 3D PFM computations, we propose an adaptive PFM for analyzing fracture in rock materials. This adaptive approach leverages an innovative finite element technique called Isogeometric Analysis (IGA) (Borden, Hughes, et al., 2014; Greco et al., 2024) to achieve both efficiency and accuracy. We leverage PHT-splines within the IGA framework to locally refine the crack path. PHT-splines offer an efficient and easy-to-implement local refinement algorithm, which provides a higher derivative of the field variable. Unlike traditional methods that require prior knowledge of the failure path, this method proposes an adaptive h-refinement scheme that dynamically refines the mesh at singularities and high-gradient regions based on error detection during computation. The governing equations for the displacement and phase fields are approximated using PHT-splines basis functions. To achieve a smooth representation of crack surfaces in the 3D internal local crack refinement strategy, a threshold is applied to the phase field variable. Using this method, we model the growth of internal cracks with different inclination angles. This innovation allows for precise capture of complex crack surfaces while significantly reducing computational costs. Our discussion encompasses exploring the correlations and differences between 2D and 3D crack propagation, including aspects such as wing crack initiation, propagation direction, and length. Subsequently, we analyze the factors contributing to the limitation of wing crack length and macroscopic failure in 3D scenarios.

The structure of the work is as follows: Section 2 provides detailed information about the theory and numerical implementation of the 3D adaptive PFM-based internal crack growth model, with a particular

emphasis on local adaptive refinement. In Section 3, we compare results from our adaptive PFM, in 2D and 3D, with previously-published results from notched square plate tension numerical tests using PFM. Section 4 presents a further validation of our adaptive PFM with laboratory experiments performed on 3D-printed resin samples, followed by models designed to understand the influence of internal crack angle on wing crack length and initiation locations. Finally, we provide in-depth discussions, including a comparison between the results of 2D and 3D numerical tests, and important conclusions at the end of the paper.

## 2. Phase field method (PFM) based 3D adaptive internal crack growth model

In this section, we provide a detailed exploration of the adaptive phase field constitutive theory. Within the framework of IGA, we use PHT-splines as the employed solution. To establish constraints on mesh refinement, we employ the phase field variable threshold as an error indicator. The chosen re-meshing procedure is based on an adaptive h-refinement scheme.

### 2.1. Phase field constitutive theory

A 3D computational solid domain, denoted as $\Omega$, is established, encompassing an external boundary, labeled as $\partial\Omega$, within which a specific region known as the crack zone $h$ is defined, as shown in Fig. 2. In the brittle field of rock deformation, the total energy $\mathfrak{E}$ of the 3D computational domain could be derived as (Ambati et al., 2014; Borden, Verhoosel, et al., 2012; Miehe and Mauthe, 2016):

$$ \mathfrak{E} = \psi_{\varepsilon} + \psi_{f} - \psi_{\text{ext}} \quad (1) $$

where

$$ \psi_{\varepsilon} = \int_{\Omega} \psi(\boldsymbol{\varepsilon}(\mathbf{u}), a) d\Omega \quad (2) $$

$$ \psi_{f} = \int_{h} G_{c} d\mathbf{h} \quad (3) $$

$$ \psi_{\text{ext}} = \int_{\Omega} \mathbf{b}(\mathbf{x}, t) \cdot \mathbf{u} d\Omega + \int_{\partial\Omega_t} \mathbf{t}(\mathbf{x}, t) \cdot \mathbf{u} d\partial\Omega_t \quad (4) $$

where $\mathfrak{E}$ is the total potential energy, $\psi_{\varepsilon}$ is the elastic strain energy, $\mathbf{u}$ is the displacement field, $\boldsymbol{\varepsilon}(\mathbf{u})$ is the elastic strain, $a$ is the phase field variable, $\psi_{f}$ is the fracture energy, $G_{c}$ is the critical energy release rate, $\psi_{\text{ext}}$ is energy made by external boundary forces, and $\mathbf{x}$ and $t$ are related to the space and time of the 3D computational domain, respectively.

To distinguish the tension-compression asymmetry in 3D PFM, the stored elastic strain energy is separated into two distinct components: tensile and compressive. This division involves the decomposition of the energy portions of the strain tensor, which are further refined through the strain spectral decomposition of Borden, Verhoosel, et al. (2012):

$$ \boldsymbol{\varepsilon}^{\pm} = \sum_{i=1}^{3} \langle \boldsymbol{\varepsilon}_{i} \rangle_{\pm} \mathbf{n}_{i} \otimes \mathbf{n}_{i} \quad (5) $$

where, for the purpose of differentiating the asymmetric behavior of crack propagation under stress, $\boldsymbol{\varepsilon}^{+}$ is the positive strain tensor reflecting tensile stress and $\boldsymbol{\varepsilon}^{-}$ is the negative strain tensor reflecting compressive stress. The operators are defined by $\langle \mathbf{x} \rangle_{\pm} = (\mathbf{x} \pm |\mathbf{x}|)/2$.

The relationship between strain and elastic strain energy density in the anisotropic phase field model can be derived from Miehe et al. (2010):

$$ \psi_{\pm}(\boldsymbol{\varepsilon}) = \frac{\lambda}{2}(\mathrm{tr}(\boldsymbol{\varepsilon}))_{\pm}^{2} + \mu(\mathrm{tr}(\boldsymbol{\varepsilon})_{\pm}^{2}) \quad (6) $$

where $\psi_{\pm}(\boldsymbol{\varepsilon})$ indicates the standard free energy density of undamaged
<figure><img src="image_7.png" /><figcaption>Fig. 3. In a three-dimensional scenario, the mesh represents the discretized structure formed by interconnected elements, with nodes indicating the discrete points within this mesh where relevant properties or values are defined or evaluated. (a) The schematic of a 3D sample. (b), (c), and (d) represent a cross-section of a three-dimensional sample at the first, second, and third refinement step, respectively.</figcaption></figure>

rock matrix, and $\lambda$ and $\mu$ are Lamé constants. We also have:

$$ \psi_{\varepsilon} = \int_{\Omega} \psi(\boldsymbol{\varepsilon}(\mathbf{u}), a) d\Omega = \int_{\Omega} [(1-a)^2 + k] \psi_{+}(\boldsymbol{\varepsilon}) + \psi_{-}(\boldsymbol{\varepsilon}) d\Omega \quad (7) $$

where $k$ is a very small, and positive, integer ($k = 1 \times 10^{-9}$ by default) that keeps the numeric calculation stability.

The fracture energy required for the formation of macrocrack is expressed as:

$$ \psi_f = \int_{\Omega} G_c \cdot \left( \frac{a^2}{2l_0} + \frac{l_0}{2} (\nabla a)^2 \right) d\Omega \quad (8) $$

where the length scale parameter, denoted as $l_0$ (Fig. 2), is a crucial material property employed to characterize the extent of crack propagation within a rock material, specifically describing the width of the discrete crack in a physical context.

The formula for crack growth can be derived through the minimization of the total potential energy. In order to find the equilibrium solution, the variation of the total energy functional must be equal to zero for all permissible values of variational operators. Consequently, the strong form of the governing equation and boundary conditions for the phase field can be determined by computing the first variation of Eq (1):

$$ 2\psi_{0}^{+}(1-a) - \frac{G_{c}}{l_{0}}(a - l_{0}^{2}\Delta a) = 0 \quad (9) $$

$$ \nabla a \cdot \mathbf{n} = 0 \quad (10) $$

where $\mathbf{n}$ is the normal vector of the boundary.

Similarly, the results of the displacement field are given by:

$$ \nabla \boldsymbol{\sigma} + \mathbf{b}(\mathbf{x}, t) = 0 \quad (11) $$

$$ \boldsymbol{\sigma} \cdot \mathbf{n} = t(\mathbf{x}, t) \quad (12) $$

The positive elastic energy density (Miehe, Hofacker, et al., 2010) serves as the driving force for phase field evolution and is introduced to enforce irreversibility conditions in crack growth:

$$ H(\mathbf{x}, t) = \max_{i \in [0, T]} \psi_{0}^{+}((\boldsymbol{\varepsilon}(\mathbf{u}), a)) \quad (13) $$

where $H(\mathbf{x}, t)$ represents a monotonic function that ensures the irreversibility of damage.

For 3D internal pre-existing cracks in rock, the pre-existing crack must satisfy the following initial boundary conditions (Borden, Verhoosel, et al., 2012):

$$ H_0(\mathbf{x}) = \begin{cases} \frac{BG_c}{2l_0} - \frac{BG_c d(\mathbf{x}, l)}{l_0^2}, & d(\mathbf{x}, l) \le \frac{l_0}{2} \\ 0, & d(\mathbf{x}, l) > \frac{l_0}{2} \end{cases} \quad (14) $$

where $d(\mathbf{x}, l)$ is defined as a function representing the shortest distance from any point $\mathbf{x}$ in the domain to the line $l$ that represents a pre-existing crack. $B = \frac{1}{1-a}$ and, when $a$ approaches unity, $B$ tends to infinity. In the present paper, we followed the approach of Borden, Verhoosel, et al. (2012) and selected $B = 1 \times 10^3$ to generate the pre-existing crack.

## 2.2. PFM implement through adaptive IGA with PHT-splines

In the pursuit of discrete crack modeling and accurate crack growth using the PFM, the quality of computational mesh plays a crucial role. An excessively coarse mesh can compromise calculation accuracy, while an overly refined mesh consumes excessive computational resources. In the solution adopted here, within the framework of isogeometric analysis, PHT-splines are used as basis functions. Basis functions are chosen over other methods such as B-splines and NURBS (Non-Uniform Rational B-Splines) (Li et al., 2022) due to their advantageous properties. PHT splines are an extension of NURBS and offer the capability of local refinement and the existence of hanging nodes. They are based on polynomial splines, which provide both smoothness and accuracy. PHT-splines excel in capturing the local features of surfaces, particularly in representing complex crack patterns involved in mode-I-II-III fracture simulations presented in this study. This allows for a more flexible and efficient representation of complex geometries and enables adaptive mesh refinement. The properties of PHT splines not only meet the requirements of design and analysis, but also provide new perspectives for isogeometric methods.

The basic theory of the adaptive IGA with PHT-splines is explored in this section. The basis functions of PHT-splines are constructed on a hierarchical T-mesh, starting at level 0 with a tensor-product mesh. The mesh is then advanced to higher levels, denoted by $i$, using the “cross-insertion” technique. The cross-insertion technique involves subdividing specific elements by connecting midpoints of opposing edges, resulting in the creation of 4 sub-elements for surfaces or 8 sub-elements for volumes. This meticulous subdivision process enhances the level of structural detail and refinement, as illustrated in Fig. 3. During the transition from refinement level $i$ to $i + 1$, the mesh undergoes subdivision into sub-elements. However, at the $i + 1$ level, the initial meshes from level $i$ remain undivided, and only the specific sub-element undergoes further division at this subsequent level. Based on a given surface T-mesh $\mathcal{T}$, $F$ represents all the set of all the cells in $\mathcal{T}$ and $\Omega$ denotes the region occupied by the cells $\phi$ (Kang et al., 2015). The polynomial
spline space is defined as:

$$P(m, n, \alpha, \beta, \mathbb{T}) := \{p(x, y) \in C^{\alpha, \beta}(\Omega) | p(x, y)|_{\phi} \in Q_{mn} \text{ for any } \phi \in F\} \quad (15)$$

where $Q_{mn}$ denotes the space of all polynomials with bi-degree $(m, n)$, and $C^{\alpha, \beta}(\Omega)$ represents the space comprising all bivariate functions that are continuous in $\Omega$, with order $\alpha$ in the x-direction and $\beta$ in the y-direction. When $m \ge 2\alpha + 1$ and $n \ge 2\beta + 1$, the dimension of the spline space $P(m, n, \alpha, \beta, T)$ is provided in Deng et al. (2006). For cubic splines with $C^1$ continuity, the dimension formula can be elegantly expressed as follows:

$$P(3, 3, 1, 1, \mathbb{T}) = 4(V^b + V^+) \quad (16)$$

where $V^b$ and $V^+$ denote the boundary and crossing vertices, respectively, in the T-mesh $\mathbb{T}$. Similarly, for a given 3D T-mesh $\mathbb{T}$, the dimension formula can be defined as:

$$P(3, 3, 3, 1, 1, \mathbb{T}) = 8(V^b + V^+) \quad (17)$$

At level 0, the surfaces and volumes of the PHT-splines are characterized by a precise formulation, which entails expressing them as a linear combination of the basis functions of PHT-splines along with their associated control points $P_i$. We have:

$$\hat{S}_i(s, t) = \sum_{i=1}^{n} \hat{R}_i(s, t) P_i, \quad \hat{R}_i(s, t) = \frac{Q_i(s, t) w_i}{\sum_{\zeta=1}^{n} Q_\zeta(s, t) w_\zeta} \quad (18)$$

$$\hat{V}(s, t, r) = \sum_{i=1}^{n} \hat{R}_i(s, t, r) P_i, \quad \hat{R}_i(s, t, r) = \frac{Q_i(s, t, r) w_i}{\sum_{\zeta=1}^{n} Q_\zeta(s, t, r) w_\zeta} \quad (19)$$

where $\hat{R}_i(s, t)$ and $R_i(s, t, r)$ represent the bivariate and trivariate NURBS basis functions, respectively, $w_i$ and $w_\zeta$ are weights, and $\hat{S}_i(s, t)$ and $\hat{V}(s, t, r)$ are the surface and volume of the PHT-splines, respectively.

At any refinement level $i$, the formulation of the surface and volume of the PHT-splines is derived as:

$$\hat{S}^i(s, t) = \sum_{i=1}^{n^i} \hat{R}_i^i(s, t) P_i^i \quad (20)$$

$$\hat{V}^i(s, t, r) = \sum_{i=1}^{n^i} \hat{R}_i^i(s, t, r) P_i^i \quad (21)$$

Continuing the iterative refinement process with the $(i + 1)$-th step for the surface and volume of the PHT-splines involves further formulating the $(i + 1)$-th step to enhance the precision and complexity of the model:

$$\hat{S}^{i+1}(s, t) = \sum_{i=1}^{n_{i+1}} \hat{R}_i^{i+1}(s, t) P_i^{i+1}, \quad \hat{R}_i^{i+1}(s, t) = \frac{Q_i^{i+1}(s, t) w_i^{i+1}}{\sum_{\zeta=1}^{n_{i+1}} Q_\zeta^{i+1}(s, t) w_\zeta^{i+1}} \quad (22)$$

$$\hat{V}^{i+1}(s, t, r) = \sum_{i=1}^{n_{i+1}} \hat{R}_i^{i+1}(s, t, r) P_i^{i+1}, \quad \hat{R}_i^{i+1}(s, t, r) = \frac{Q_i^{i+1}(s, t, r) w_i^{i+1}}{\sum_{\zeta=1}^{n_{i+1}} Q_\zeta^{i+1}(s, t, r) w_\zeta^{i+1}} \quad (23)$$

PHT-splines, regarded as a generalization of B-spline basis functions, find extensive applications in computational modeling owing to their hierarchical mesh structure. The basis functions of PHT-splines inherit the desirable properties of B-splines and NURBS basis functions, facilitating efficient local refinement processes.

## 2.3. Adaptive h-refinement scheme using PHT-splines

The iterative re-meshing procedure in the context of adaptive mesh refinement uses the h-refinement scheme, as discussed by Kiran et al. (2023). This scheme leverages the capabilities of PHT-splines to

<figure><img src="image_8.png" /><figcaption>Fig. 4. Geometry of the 2D square plates (side length of 1 mm) containing a single-edge crack (0.5 mm in length) positioned centrally along the left side of the plate (red line).</figcaption></figure>

facilitate mesh refinement while minimizing computational costs associated with the refinement process. The h-refinement scheme does not require the pre-specification of crack refinement paths, thus avoiding the drawbacks associated with pre-defined crack propagation paths. These drawbacks include dependency on prior knowledge, risks of grid distortion, increased computational costs, and limitations in accommodating various types of crack morphologies.

The main objective of adaptive mesh refinement is to progressively optimize the mesh throughout the simulation, thereby incrementally adding detail to the mesh where necessary, while reducing to a maximum the computational time needed. However, when dealing with geomaterials and accounting for factors such as mineral or grain size, selecting a suitable refinement criterion can be challenging. Therefore, external criteria are employed to define the extent of the mesh refinement. In this context, we explore the phase field damage threshold in T-junctions as a guiding parameter for mesh refinement limitations. The phase field damage threshold plays a crucial role in determining when and where mesh refinement occurs. Although the influence of the phase field damage threshold on computational efficiency within the model has not been explicitly addressed, established methods for determining this threshold have been employed. In this study, a phase field damage threshold value of 0.5, as suggested by Borden, Verhoosel, et al. (2012) and Goswami, Anitescu, et al. (2020), has been selected. Notably, the chosen damage threshold of 0.5 is strategically set below the threshold for crack initiation, aiding in accurately identifying the region surrounding the crack. When the phase field damage variable is lower than the damage threshold, no mesh refinement occurs. Conversely, if the phase field damage variable exceeds the damage threshold, a maximum refinement level is established in conjunction with the length scale parameter. This process ensures that the mesh refinement is concentrated in the vicinity of the crack, where it is most needed for capturing accurate crack propagation behavior.

## 3. Validation of the adaptive PFM-based internal crack growth model

An adaptive isogeometric approach, grounded in the PFM, is
<table><caption>**Table 1**

Input parameters for phase field model to validate and then explore the effect of the notch length scale in a 2D notched square plate.</caption><thead><tr><td>Elastic modulus E (GPa)</td><td>Critical energy release rate G<sub>c</sub> (N/mm)</td><td>Poisson's ratio</td><td>Length scale parameter l<sub>0</sub> (mm)</td></tr></thead><tbody><tr><td>210</td><td>2.7</td><td>0.3</td><td>0.0075, 0.01, 0.0125</td></tr></tbody></table>

employed to simulate the three-dimensional progression of internal crack growth. The refinement criterion for local adaptive mesh refinement is determined by a threshold value of the phase-field, serving as the error indicator in the computational scheme. To validate the accuracy and robustness of the proposed framework, we assess its effectiveness against several classical examples. Specifically, we apply tension to square plates containing a single-edged crack to look at the propagation of tensile (mode-I) cracks in 2D (Section 3.1) and 3D (Section 3.2). If the observed crack growth patterns and mechanical properties align well with the results of previously-published examples using PFM, it will serve as an effective indicator to validate and affirm the effectiveness of our adaptive model. Additionally, we conducted a comparative analysis using locally-refined meshes and explored the sensitiveness of length scale parameters on the peak force of notched square plates. This validation ensures that the model is reliable and capable of simulating realistic crack propagation phenomena, which is crucial for the subsequent exploration of 3D internal crack growth in rock in Section 4.

### 3.1. Validation of the 2D notched square plate tension tests

2D square plates with a single notch positioned centrally along the left side were analyzed (the geometry of the plate is shown in Fig. 4). The bottom surface of the model was fixed, while the top surface was subjected to tensile stress. The normal direction of the remaining boundaries was constrained. An upward incremental displacement step of

$Δu = 5 × 10⁻⁵mm$ was applied to the top surface, resulting in tensile stress conditions in the model. The input mechanical parameters for the notched square plate are shown in Table 1. Various length scale parameters for the notch, $l_0 = 0.0075$ mm, 0.01 mm, and 0.0125 mm, were analyzed (Table 1). Simulation here are compared to the solution for the notched plate provided using the approach outlined in Miehe, Hofacker, et al. (2010).

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Load-displacement curves determined using the adaptive phase field method for 2D square plates containing notches with different length scale parameters of $l_0 = 0.0075$ mm, 0.01 mm, 0.0125 mm (solid lines). Also plotted is a load-displacement curve from Li et al. (2020) who used $l_0 = 0.0075$ mm (dashed line).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Mesh refinement (top, blue squares) and progress of crack propagation (bottom, red lines) in a 2D square plate are investigated using the proposed PFM framework ($l_0 = 0.0075$ mm). The images on the left show the start of crack propagation in the plate, the images in the middle show the plate after the crack has propagated a certain distance (from left to right), and the images of the right show the plate after the crack has reached the other side of the plate.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Geometry of the 3D square plates (side length of 1 mm and thickness of 0.2 mm) containing a single-notch crack (0.5 mm in length) positioned centrally along the left side of the plate (dark brown plane).</figcaption></figure>

Fig. 5 illustrates the progressive development of the 2D crack and the localized mesh refinement at different stages. The refinement process focuses on enhancing the mesh resolution in the vicinity of the crack, gradually transitioning from a finely refined region around the crack to a coarser mesh in the surrounding areas. Additionally, as the crack propagates, this refinement strategy is systematically implemented to adapt to the evolving crack geometry, ensuring accurate representation of the behavior of the crack. Fig. 6 illustrates the load-displacement curves for each of the three crack length scale parameters ($l_0$ = 0.0075 mm, 0.01 mm, and 0.0125 mm) compared with previously-published model results obtained by Li, Nguyen-Thanh, et al. (2020) for a length scale parameter $l_0$ = 0.0075 mm. It is observed that a smaller length scale parameter results in a higher peak force in the adaptive mesh in Fig. 6, which is in agreement with the results obtained by Li, Nguyen-Thanh, et al. (2020). Despite the variation in length scale parameter (i.e., ranging from 0.0075 to 0.0125 mm), the observed changes in peak force in our results are minimal compared with those seen in other results using PFM, such as those by Li, Nguyen-Thanh, et al. (2020). This can be attributed to the error reduction achieved by the adaptive method used in our approach, as it selectively refines the mesh in the damage area rather than across the entire computational domain. In theory, PFM is expected to reduce errors caused by length scale parameters, and the proposed adaptive PFM has been shown to effectively reduce these errors.

## 3.2. Validation of the 3D notched square plate tension tests

We have extended our phase field model to investigate tension tests

<caption>**Table 2**</caption>

Input parameters for phase field model to validate and then explore the effect of the notch length scale in a 3D notched square plate.

<table><thead><tr><td>Elastic modulus E (GPa)</td><td>Critical energy release rate G<sub>c</sub> (N/mm)</td><td>Poisson's ratio</td><td>Length scale parameter l<sub>0</sub>(mm)</td></tr></thead><tbody><tr><td>210</td><td>2.7</td><td>0.3</td><td>0.075, 0.0625, 0.05, 0.0313</td></tr></tbody></table>

on notched square plates in 3D, as illustrated in Fig. 7. Notably, in 3D, the plate is assigned a thickness of 0.2 mm in the 3D configuration. The normal direction of all boundaries in the 3D notched square plate is constrained, maintaining consistency with the 2D case described above. The material properties (elastic modulus, critical energy release rate, and Poisson's ratio) of the 3D notched square plate remained unchanged from those used in the 2D case (compare Tables 1 and 2). For the 3D simulations, however, length scale parameters of 0.075 mm, 0.0625 mm, 0.05 mm and 0.0313 mm were employed, as detailed in Table 2.

Fig. 8 depicts the 3D refinement specifically tailored for local cracks. Remarkably, a smoothly integrated transition zone is observed between the crack surface and the undamaged region. The adaptive refinement is strategically applied to the plate area susceptible to failure, thus affirming the efficacy of the adaptive method in seamlessly transitioning from 2D to 3D scenarios. Furthermore, notable mesh refinement is observed in the crack width region while, in other areas, the mesh size remains relatively larger. Fig. 9 shows the load-displacement curves obtained using our adaptive PFM approach compared with the adaptive PFM results obtained by Goswami, Anitescu, et al. (2020), who used a length scale parameter $l_0$ = 0.0313 mm. Fig. 9 shows that a decrease in the length scale parameter from 0.0313 to 0.075 mm results in an decrease in the peak force, as also observed for the 2D models (Fig. 6). Importantly, our numerical results exhibit commendable consistency with the findings of Goswami, Anitescu, et al. (2020) (Fig. 9).

## 4. Modeling of 3D internal crack growth in rock

In this section, we aim to further validate the effectiveness of the 3D adaptive PFM-based crack growth model. To do so, we compare results from our model with those obtained from previously-published uniaxial compression tests on resin-based artificial rocks fabricated using 3D printing (3DP) (from Zhou et al., 2019). These artificial rocks contained a single penny-shaped crack inclined at an angle of $\beta$ = 30° (Zhou et al., 2019). Following this validation, we used the 3D adaptive PFM-based crack growth model to predict the rock cracking process for various crack inclination angles of $\beta$ = 0°, 45°, 60°, 90° to understand the relationship between wing crack morphology and crack inclination angle.

The numerical samples used in the PFM were rectangular cuboids with a height of 60 mm and a square base of 30 mm (Fig. 10), the same geometry as used in the experiments of Zhou et al. (2019). In the center of the numerical samples, an inclined penny-shaped crack with a radius of 5 mm and a negligible thickness was positioned. The bottom of the numerical rock sample was fixed, and we applied downward incremental displacement steps $\Delta u$ = $1 \times 10^{-5}$ mm to the top surface, while leaving the remaining surfaces of the model unrestricted. The physical input parameters for the PFM are given in Table 3, based on the experimental results of Zhou et al. (2019). To ensure that the diffuse crack surface modeled by the PFM to converge to a sharp crack surface, or to model materials with high initiation stresses at crack tips, a small value of $l_0$ is required. The authors (Miehe, Hofacker, et al., 2010; Miehe and Mauthe, 2016) have proposed that the phase field model parameter $l_0$ be related to the mesh size, suggesting that $l_0$ should be approximately twice the minimum mesh size: $l_0$ = $2h_{size} \approx 0.5mm$.

Fig. 11 shows the numerical results of crack initiation and growth morphology obtained using the 3D adaptive PFM-based internal crack growth model for sample containing a penny-shaped crack inclined at
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Crack propagation and corresponding mesh (red squares) refinement in a 3D square plate are investigated using the proposed PFM framework ($l_0$ = 0.05 mm). The images on the left show the start of crack propagation in the plate, the images in the middle show the plate after the crack has propagated a certain distance (from left to right), and the images of the right show the plate after the crack has reached the other side of the plate.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Load-displacement curves determined using the adaptive phase field method for 3D square plates containing notches with different length scale parameters of $l_0$ = 0.075 mm, 0.0625 mm, 0.05 mm and 0.0313 mm (solid lines). Also plotted is a load-displacement curve from Goswami et al. (2020) who used $l_0$ = 0.0313 mm (dashed and dotted line).</figcaption></figure>

30°. During the deformation of the rock sample, wing cracks gradually originated near the tips of the pre-existing crack (Fig. 11 (b)). As axial loading increased, these wing cracks extended in length, accompanied by substantial crack wrapping (Fig. 11 (c)). Eventually, when the wing cracks reached a certain length, equivalent to the radius of the pre-existing penny-shaped crack, they ceased to propagate further (Fig. 11 (d)). Comparing these numerical findings with the corresponding experimental results (from Zhou et al. (2019) and Fu, Chen, et al. (2016)), as illustrated in the Fig. 12, it is evident that the numerical findings closely align with the experimental results. However, in the experiments of Zhou et al. (2019), the growth of the wing cracks in the lower-left and upper-right regions of the sample did not occur simultaneously (Fig. 12 (a) and Fig. 12 (b)). Additionally, an anti-wing crack emerged in the lower-right region of the pre-existing crack, as depicted in Fig. 12 (c). This phenomenon can be attributed to the heterogeneity of rock, introducing variations in its physical properties: consequently, these variations in physical properties result in different propagation rates for the wing crack and the occurrence of new cracks. Compared with the side and front views of the specimen's fracture pattern (Fu et al., 2016), as shown in Fig. 12 (d), the 3D adaptive PFM simulations in this study agree well with the experimental observations, specifically in capturing the frequently observed wing crack pattern. Additionally, the 3D PFM-based internal crack growth model was not capable of depicting the eventual failure pattern. According to published studies, despite the absence of any subsequent cracks beyond the wing crack, the rock sample continued to experience extended axial deformation (Dyskin, Sahouryeh, et al., 2003; Wang et al., 2022b). In these published studies, in contrast to the results of our modeling, the ultimate failure of the
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** 3D Geometry conditions of a numerical sample (in light green) containing a single inclined penny-shaped crack (with a radius of 5 mm) in the central region (in dark green).</figcaption></figure>

<table><caption>**Table 3**

Input parameters for phase field model to explore the effect of the crack inclination angle on the crack modes of rock.</caption><thead><tr><td>Elastic modulus E (GPa)</td><td>Critical energy release rate G<sub>c</sub> (N/m)</td><td>Poisson's ratio</td><td>Length scale parameter l<sub>0</sub> (mm)</td></tr></thead><tbody><tr><td>4.3</td><td>6</td><td>0.27</td><td>0.5</td></tr></tbody></table>

samples occurred abruptly, culminating in a sudden and decisive breakage. As observed in experiments, this unstable crack propagation often results in dynamic phenomena, such as disintegration (Zhou and Zhang, 2020). This phenomenon necessitates a transition from continuous to discontinuous simulation methods. Therefore, the phase field method based on continuous method fails to capture the collapse of the specimen, and exhibits a sudden stress drop in the stress–strain curve (Fig. 12 (e)). In Fig. 12 (e), stress–strain curves generated by the 3D adaptive PFM demonstrate consistency with experimental findings reported by Zhou et al. (2019). It becomes apparent that the 3D adaptive internal crack growth model adeptly captures the stress–strain relationship, overcoming the limitations of the continuous phase field method in capturing the collapse and exhibiting a smooth stress–strain response.

Following this validation, we conducted numerical uniaxial compression tests on rock samples containing pre-existing penny-shaped cracks. The penny-shaped crack was positioned at the center of each sample and inclined at various angles $\beta$ with respect to the horizontal

direction. To investigate the behavior of these pre-cracked rock samples, four test groups were performed, each corresponding to a different inclination angle for the penny-shaped crack. The specified inclination angles were as follows: $\beta = 0^{\circ}$, $45^{\circ}$, $60^{\circ}$, and $90^{\circ}$. These numerical rock samples were cuboid prismatic in shape with dimensions of 30 mm × 30 mm × 60 mm (length × width × height). We used the validated values in Table 3 as the physical input parameters for the model. In each time step of the testing process, we carefully selected the magnitude of compressional displacement increment, $\Delta u = 1 \times 10^{-5}$ mm, to ensure accurate data collection and analysis. The different inclination angles of the numerical results were also further compared with the experimental observations conducted by Xu et al. (2023) and Fu, Chen, et al. (2016).

Fig. 13 shows snapshots of the front, top, and side views of the numerical rock samples to highlight the progression of crack initiation and growth for each of the crack inclination angles ($0^{\circ}$, $45^{\circ}$, $60^{\circ}$, and $90^{\circ}$). For rock samples with $\beta = 0^{\circ}$ (Fig. 13 (a)), the initiation of wing cracks occurred at the midpoint of the penny-shaped crack (Fig. 13 (a-II)). As axial loading increased, the width of the wing cracks progressively increased. Simultaneously, the length of these wing cracks increased, propagating parallel to the direction of the applied load. The propagation of wing cracks ceased when their length reached 1.2 times the radius of the original penny-shaped crack (Fig. 13 (a-III)). Subsequently, wing wrapping emerged perpendicular to the existing wing cracks. Similarly, as axial loading increased, the surface area of the wing wrapping gradually increased (Fig. 13 (a-IV)). For rock samples with $\beta = 45^{\circ}$ and $60^{\circ}$ (Fig. 13 (b) and (c)), the initiation of wing cracks was notable near the tips of penny-shaped cracks. These wing cracks exhibited a gradual reduction in length. Wing wrapping, originating around the original penny-shaped cracks, actively resisted the propagation of these wing cracks. For rock samples with a $\beta$ of $90^{\circ}$ (Fig. 13 (d)), wing cracks initiated specifically at the end tips of the original penny-shaped crack. These cracks then propagated parallel to the applied loading. The observed length of these wing cracks is significantly smaller, approximately 0.1 times the size of the original penny-shaped crack (Fig. 13 (d-IV)). Importantly, no additional cracks were observed to initiate at the later stages of deformation. We further compared the numerical results with the experimental 3D inner crack progressive growth pattern, as illustrated in Fig. 14. It is evident that numerical findings closely align with the experimental results. However, the laboratory results also exhibit a distinct failure pattern known as splitting failure. These authors (Fu, Chen, et al., 2016; Xu, Li, et al., 2023) suggest that this characteristic arises as an artificial phenomenon due to the processing technology employed for creating internal cracks or the significant difference in stiffness between the samples and the bearing platform of the testing machine. The presence of friction during the compression of internal defects exacerbates this effect. For samples such as glass or 3D-printed resin-based artificial rocks, which are ideal brittle materials with high stiffness, stress concentration and localized failure occur at specific points due to imperfect alignment of the sample's end with the bearing platform of the testing machine.

## 5. Discussion

### 5.1. The advantage of the 3D adaptive PFM-based crack growth model in rock

Investigating the growth of rock cracks is essential due to its profound implications for understanding the fatigue and fracture characteristics of rocks. By understanding the mechanisms of crack propagation, we can better predict and prevent potential damage to rock structures, ensuring the safety of various engineering projects. Moreover, this deeper insight into crack growth contributes significantly to our understanding of geological processes, such as deformation associated with tectonic or volcanic activity. This, in turn, can help improve early warning systems and mitigate the risks associated with disasters
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** 3D snapshots of crack initiation and growth (a-d), with (d) illustrating the front, top, and side views captured through PFM on numerical samples containing a penny-shaped crack inclined at an angle of 30°. The adaptive mesh is indicated by the blue boxes and the crack is shown in beige.</figcaption></figure>

such as the collapse of engineered underground spaces. Therefore, delving into the growth of rock cracks not only enhances rock engineering practices, but also positively influences sustainable development.

Until now, numerical methods have been widely employed to analyze and simulate fracture mechanism in rocks. Accurately simulating crack propagation in rock, especially in 3D, has been a paramount challenge within the field of geomechanics. Various numerical methods require suitable failure criteria for modeling crack propagation. Some common approaches involve using microelements to represent macroscopic properties. For instance, the Rock Failure Process Analysis (RFPA) method, developed by Tang et al. (2000), incorporates a damage variable to reduce the elastic modulus based on the maximum tensile stress and Mohr-Coulomb failure criteria within a finite element framework. In the Discrete Element Method (DEM) (Fu et al., 2020a), rocks are treated as assemblies of rigid particles, and bond breakage signifies the formation of microcracks. The Discrete Fracture Network (DFN) approach (Ye et al., 2020) is also used within volumetric finite element elements to model and simulate crack behavior. However, these

methods often require a large number of elements to adequately represent cracks at the microscale and achieve spontaneous crack growth. Understanding the behavior of individual cracks remains a challenge using these methods. In addition to the methods mentioned above, the majority of numerical methods rely on stress intensity factors (SIFs) to determine the crack initiation and propagation direction. The failure criteria such as maximum tensile stress (Hobbs, 1964), minimum strain energy density (Sih, 1974), and maximum energy release rate (Nuismer, 1975) are commonly used. However, extending these methods to 3D cases becomes challenging due to the consideration of mode-III cracks in fracture mechanics.

In PFM, the criteria for crack growth are inherently incorporated as a fundamental component of the solution, eliminating the requirement for external criteria. This is particularly advantageous when dealing with the complex process of 3D crack propagation, as it minimizes the substantial computational resources required for the explicit tracking of crack growth. Fig. 15 depicts the 3D internal crack growth in rock based on existing studies. While these methods can serve as references, it is evident that replicating the effect of wing wrapping has proven
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Crack growth snapshots from experimental samples (3D-printed resin samples) that contained a pre-existing crack with an inclination angle of 30° (a–c) (Zhou, Zhu, et al., 2019). Side and front views of the fracture pattern of the specimen, which contains a pre-existing crack with an inclination angle of 30°, are depicted in (d) (Fu, Chen, et al., 2016). Comparisons are made between the axial stress–strain curves of rock samples containing internal cracks at a 30° angle and the experimental results (Zhou, Zhu, et al., 2019), utilizing the proposed 3D adaptive internal crack growth model (e).</figcaption></figure>

challenging. The 3D-based PFM internal crack growth model is better suited to capture accurate test information and reflect the details of 3D wing crack growth. As a result, we combined the findings obtained from deforming 3D-printed resin-based artificial rock samples (Zhou, Zhu, et al., 2019) and subsequently developed a 3D adaptive PFM-based model for simulating internal crack growth, as illustrated in Fig. 11 and Fig. 13. The outcomes of the developed 3D adaptive PFM faithfully replicate the mechanisms underlying crack initiation and growth. Fig. 16 summarizes the mechanics of the shift and propagation of wing cracks. As the inclination angle increases, the location of wing crack initiation undergoes a gradual shift from the central position of the penny-shaped crack towards the tips (Fig. 15). Fig. 16 (b) also illustrates the variation in wing crack initiation stress concerning crack inclination angles, showcasing an initial decrease followed by an increase in stress levels. Despite this shift, the propagation direction remains consistent, moving towards the direction of the applied loading. Subsequently, the wing wrappings relocate from the center of the wing crack to both sides, ultimately encircling and wrapping around the wing cracks (Fig. 16).

Fig. 17 displays the damage volume of the rock sample along with the corresponding refined mesh size, which is determined by the phase field variable. It is evident that the mesh size is smaller in the vicinity of the penny-shaped crack compared to volumes farther away. The damage area, represented by the phase field variable, indicates that the mesh near the wing crack is also refined, while other volumes remain unfined. In Fig. 18, the external refinement mesh size is presented, revealing that it is noticeably larger than the internal mesh size shown in

Fig. 17. This observation suggests that the new crack initiated from the internal penny-shaped crack does not reach the boundary of the numerical rock sample. However, for various inclination angles, the refinement of the mesh size varies, as illustrated, for example, in the left lower regions of (b) and (c) in Fig. 18. This discrepancy indicates that different inclination angles result in different crack propagation lengths. Importantly, using local mesh refinement, rather than uniformly-refined meshes, has significantly reduced computation time. This adaptive approach allows for a more efficient allocation of computational resources, focusing on the regions of interest where crack propagation occurs, while maintaining coarser meshes in less critical areas.

## 5.2. 3D internal crack initiation and growth mechanics in rock under compression

As mentioned in the introduction, the behavior of wing crack differs between 2D and 3D cases. In this study, we conducted simulations to analyze the growth of 3D internal penny-shaped cracks with various inclination angles under uniaxial compression. The dynamics of crack evolution, including initiation, growth, and ultimate fracture can vary at different stages. Crack initiation locations, readily observable and quantifiable, serve as a key indicator for characterizing the crack initiation process. We will now discuss differences between crack initiation locations in 2D and 3D using PFM. To do so, we defined the distance between the wing crack initiation location from the crack tips as $a_1$. Additionally, the half of the internal penny-shaped crack length is
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** The numerical rock samples containing penny-shaped cracks at various inclination angles: (a) 0°, (b) 45°, (c) 60°, and (d) 90° subjected to uniaxial compression using the 3D adaptive PFM-based crack growth model, revealing stages of crack initiation and growth (I-IV). In (IV), the front, top and side views of the sample are shown.</figcaption></figure>

denoted as $a_2$. We then calculated the average ratio of $a_1/a_2$, as illustrated in Fig. 19 (a). The results of 2D PFM are taken from Xu, Xu, et al., (2022a). Fig. 19 (a) shows a comparison between crack initiation locations in 2D and 3D. The data of Fig. 19 (a) show that there is no discernible difference in the initiation location of cracks between the 2D and 3D scenarios.

The propagation direction of the 3D internal wing crack aligns with that of the 2D cases, both extending toward the applied loading. However, a significant difference arises in the length of the wing crack, as shown in Fig. 11 and Fig. 13. We noticed a correlation between the inclination angle of the original penny-shaped cracks and the length of the wing crack. To standardize the underlying relationships across different cases, we calculated the average ratio of the length of the wing cracks to the radius of the penny-shaped cracks, as shown in Fig. 19 (b).

In this context, we define the perpendicular distance between the final tip of the wing crack propagation morphology and the internal penny-shaped crack plane as $a_3$. Furthermore, the radius of the penny-shaped crack is denoted as $a_2$. The average ratio of $a_3/a_2$ was then calculated, as illustrated in Fig. 19 (b). The results of 2D PFM are again taken from the Xu et al. (2022). In the 2D case, the propagation of wing cracks is unconstrained, resulting in consistent lengths across various inclination angle scenarios, defaulting to an average ratio of 2 (Fig. 18 (b)). In other words, wing crack lengths remain consistent regardless of the crack inclination angle in the 2D case. However, in the 3D case, where the inclination of the pre-existing penny-shaped crack plays a crucial role, the length of the wing cracks decreases with increasing crack inclination, ranging from 1.2 to 0.1 times the radius of the initial penny-shaped cracks (Fig. 18 (b)). Similar to the initiation position of wing
<figure><img src="image_18.png" /><figcaption>**Fig. 13.** (continued).</figcaption></figure>

cracks, the length of 3D wing cracks is also dependent on the crack inclination angle. Dyskin, Sahouryeh, et al. (2003) observed that during uniaxial compressive strength tests, the growth of a wing crack ceased when its length reached approximately 1.0 to 1.5 times the radius of the pre-existing crack, persisting in this state until the final rupture of the specimen. Tang et al. (2015) similarly found that when the length of wing cracks reached 1.5 times the radius of the original penny-shaped crack, the propagation of the wing crack halted. Subsequently, numerous microcracks emerged within the rock samples, leading to the coalescence of these microcracks and ultimately causing damage to the rock specimens.

Fig. 20 illustrates the fracturing mechanisms of 3D inner crack at various angles (0°, 30°-60°, to 90°). Notably, 3D numerical rock models exhibit constrained crack propagation within both the vertical and horizontal planes, as shown in Fig. 20. In contrast to 2D rock cases where cracks propagate freely within a vertical plane, the complex stress state in 3D scenarios leads to irregular crack propagation paths. This dual

influence from both vertical and horizontal planes introduces intricacies in the growth dynamics of 3D wing cracks. According to fracture mechanics classification, it is evident that the crack transitions from mode I-III fracture to mode I-II-III fracture, and finally to mode I fracture, as the inclination angles increase. In the case of a specimen containing a penny-shaped crack, loading typically occurs radially, resulting in a stress field with both radial and circumferential symmetry within a homogeneous specimen. Consequently, the direction of crack propagation is often influenced by this symmetry, and the newly formed cracks at both ends of the crack exhibit a distinct symmetry. In this scenario, the crack may tend to propagate along the radial direction of the specimen due to the highest stress experienced at the center of the specimen. Propagating in this direction reduces the degree of stress concentration, thereby satisfying the material's deformation and fracture requirements. For rock samples with $\beta = 0°$, under uniaxial compression loading, the stress field of the specimen typically generates radially and circumferentially. This means that the center of the specimen experiences
<figure><img src="image_19.png" /><figcaption>Fig. 13. (continued).</figcaption></figure>

radial compressive stress, while the crack surface experiences circumferential tensile stress. The formation of the wing crack leads to stress redistribution, causing an increase in circumferential tensile stress around the center. Consequently, the formation of circumferential cracks becomes more likely, as cracks tend to propagate along the direction of maximum tensile stress, which is perpendicular to the primary crack. Similar to the 2D scenarios, cracks initially exhibit mode-I fracture. However, unlike in 2D situations, they later exhibit mode-III fracture due to the influence of out-of-plane shear stress, which is neglected in 2D analyses. For crack angles ranging from 30°, 45°, and 60°, the ultimate fracture surface of the sample exhibits a predominantly three-dimensional curved profile. Some areas of the crack surface appeared twisted, indicating a typical Mode I-II-III fracture pattern. For rock samples with $\beta$ = 90°, the ultimate fracture surface of the sample exhibits a flat and smooth profile, aligning with the initial internal crack. For vertical pre-existing cracks parallel to the loading direction, since the loading is applied radially in the specimen, these cracks are not directly exposed to the circumferential stress field. Instead, they

primarily experience compressive stresses in the radial direction. Thus, under uniaxial compression loading, vertical circular pre-existing cracks typically only undergo the influence of the radial stress field and are not directly affected by the circumferential stress field. The expansion path of these cracks is constrained by the radial stress field, resulting in propagation predominantly along the radial direction with minimal expansion in the circumferential direction. Their failure mode resembles that of the noneccentric internal crack under three-point bending, displaying a characteristic Mode I fracture pattern. However, a notable discrepancy is observed: the propagation velocity of vertical internal crack consistently remains below the critical speed during uniaxial compression. Moreover, no discernible traits of dynamic crack propagation manifest on the fracture surface. Therefore, the simplification of 3D internal crack propagation to a 2D plane is deemed unreasonable as it overlooks the crucial role played by mode-III cracks, which are of significant importance in fracture mechanics (Zhang and Tong, 1996). It is precisely the influence of these mode-III cracks that initiates the wing wrapping, thereby impeding the propagation of the wing crack.
<figure><img src="image_20.png" /><figcaption>Fig. 13. (continued).</figcaption></figure>

While our 3D PFM-based crack growth model successfully simulates the initiation and propagation behavior of wing cracks, it does not capture the ultimate failure pattern, as illustrated in Fig. 13. In 3D uniaxial compressive strength tests involving a centrally positioned penny-shaped crack within rock samples, lateral propagation of wing wrapping encompasses both the wing cracks and the central penny-shaped crack. This lateral expansion effectively restrains further growth of the wing cracks, preventing the macroscopic axial fracture of the rock sample. Despite the substantial length of the wing cracks for rock samples with $\beta = 0°$ (Fig. 13 (a)), their ultimate size proves insufficient when compared to the overall dimensions of the rock sample, thus preventing its ultimate splitting. Even though some split failure occurred, it was primarily attributed to the processing technology applied to internal cracks or the significant difference in stiffness between the samples and the bearing platform of the testing machine. Based on the findings of published studies, although no additional cracks develop beyond the initial wing crack, the rock sample undergoes prolonged axial deformation (Dyskin, Sahouryeh, et al., 2003; Wang, Dyskin, et al., 2022b). The ultimate failure of the samples in these studies is characterized by a sudden, decisive breakage. Even in uniaxial compressive strength tests on rock samples with multiple pre-existing

cracks, the coalescence of wing cracks induced by these pre-existing cracks does not significantly contribute to the macroscopic failure of the rock sample (Dyskin, Sahouryeh, et al., 2003; Yang et al., 2010; Yin, Wong, et al., 2014; Zhou, Zhang, et al., 2018). The presence of petal cracks plays a pivotal role in the interaction and coalescence of multiple cracks, ultimately leading to the macroscopic failure of the rock. Consequently, the growth of 3D wing cracks within a shared plane is less pronounced compared to the coalescence observed in 2D cracks.

However, studies have shown that the expansion of wing cracks parallel to the applied loading can lead to macroscopic failure in biaxial compression ($\sigma_1 > \sigma_2 > \sigma_3 = 0$) tests (Sahouryeh et al., 2002; Wang, Dyskin, et al., 2022b). These investigations revealed that the growth of wing cracks is unincumbered, primarily due to the inconspicuous presence of wing wrappings. This phenomenon is attributed to the intermediate principal stress, which hinders the growth of wing wrappings while simultaneously enveloping and constraining the expansion of wing cracks. This implies that the intermediate principal stress, often overlooked in 3D rock cracking processes, plays a crucial role. Taylor et al. (1972) found that the intermediate principal stress exhibits a threshold effect: when below this threshold, wing wrappings occurs around the initial cracks; however, upon exceeding this threshold, wing
<figure><img src="image_21.png" /><figcaption>Fig. 14. Experimental observation of wing cracks progressive growth patterns from cuboid samples that contain a pre-existing shaped crack with an inclination angle of 60° (a), 90° (b) (Xu, Li, et al., 2023) and 45° (c) (Fu, Chen, et al., 2016).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 15. Comparative study of various numerical methods for 3D internal crack propagation (a) the non-ordinary state-based peridynamics by Shou, Zhou, et al. (2019), (b) the finite difference method by Fu, Chen, et al. (2016), and (c) the meshless method by Zhang et al. (2016).</figcaption></figure>

cracks experience unrestricted growth, resembling the behavior observed in 2D scenarios. In this study, we simulated the 3D growth of wing cracks under uniaxial compression conditions ($\sigma_1 > \sigma_2 = \sigma_3 = 0$),

signifying the absence of lateral constraints. This implies that as long as the intermediate principal stress remains below a certain threshold, the expansion of wing cracks is constrained, preventing macroscopic failure induced by wing cracks.

## 6. Conclusions

In this study, we successfully implemented and integrated a 3D adaptive crack growth model into PFM to provide valuable insights into the mechanics of internal cracks in rock samples. The solution was executed within the isogeometric analysis framework, using PHT-splines and employing the phase field variable as an error indicator. Validation of the phase field-based 3D crack growth model was conducted through a comparison with previously-published numerical tests performed on notched plates. The incorporation of adaptive PFM significantly improves computational efficiency in 3D scenarios. Subsequent to this validation, we conducted a comprehensive analysis on 3D numerical rock samples containing a pre-existing penny-shaped crack at their center subjected to uniaxial compression. Through our investigations, we observed that increasing the inclination angle of the internal penny-shaped crack – from perpendicular to the loading direction to parallel to the loading direction – resulted in a gradual shift of the wing crack initiation location from the middle towards the end tips of the penny-shaped crack.
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 16.** Sketch depicting the initiation positions (a) of wing cracks with changing pre-existing, penny-shaped crack inclination angle, and the corresponding crack initiation stress (b).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 17. The numerical grid size of the rock specimen and the corresponding damage conditions. (a) 0°, (b) 30°, (c) 45°, (d) 60°, and (e) 90°. *a* = 1 indicates the fully broken state (i.e. red), while *a* = 0 indicates the intact state (i.e. blue).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 18. The external profile chart illustrates the adaptive mesh refinement for numerical 3D rock samples with a cubic shape containing penny-shaped cracks with different inclination angles: (a) 0°, (b) 30°, (c) 45°, (d) 60°, and (e) 90°.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 19. The numerical results from Phase Field Modeling (PFM) comparing 2D and 3D simulations. 2D data are from Xu et al. (2022). (a) illustrates the initiation locations of wing cracks at different penny-shaped crack inclinations, while (b) presents the length of wing cracks at different penny-shaped crack inclinations.</figcaption></figure>

shaped crack, accompanied by a corresponding reduction in the length of the wing cracks. The corresponding fracturing mechanisms transition from mode I-III fracture (0°) to mode I-II-III fracture (30°, 45°, 60°), and

eventually to mode I fracture (90°) as inclination angles increase. Additionally, comparing the results obtained from 3D PFM with those from 2D simulations (from Xu et al., 2022), we discovered that
<figure><img src="image_28.png" /><figcaption>Fig. 20. Fracturing mechanisms of 3D internal crack at various angles from 0°, 30°-60°, to 90°. ($\sigma_1$ is axial compressive stress.).</figcaption></figure>

macroscopic failure of rock in 3D is significantly influenced by the in-
termediate principal stress. This crucial factor, often overlooked in 3D
rock cracking processes, plays a substantial role in restraining the
growth of wing cracks and preventing macroscopic failure induced by
these cracks. Overall, our study provides valuable insights into the
behavior of cracks in 3D rock samples under loading. The integration of
adaptive PFM and the consideration of the intermediate principal stress
contribute to a comprehensive understanding of crack growth mecha-
nisms and the factors influencing macroscopic failure in rock samples.
These findings offer practical implications for fracture mechanics and
can aid in the development of effective strategies for mitigating rock
failure in engineering and geotechnical applications.

CRediT authorship contribution statement

**Bin Xu:** Writing – original draft, Validation, Formal analysis, Data curation. **Tao Xu:** Writing – review & editing, Supervision, Resources, Investigation, Funding acquisition. **Michael J. Heap:** Writing – review & editing, Data curation. **Alexandra R.L. Kushnir:** Writing – review & editing, Data curation. **Bo-yi Su:** Software, Investigation. **Xiao-cong Lan:** Software.

Declaration of competing interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

This work was supported by the National Natural Science Foundation
of China (grant numbers 42172312, 52211540395). M.J. Heap ac-
knowledges support from the Institut Universitaire de France (IUF), and
A.R.L. Kushnir acknowledges the support of the Swiss National Science
Foundation (SNSF) Prima grant ReSiDue (PR00P2-201611).

References

Adams, M., Sines, G. 1978. Crack extension from flaws in a brittle material subjected to compression. Tectonophysics 49(1-2). Doi: 10.1016/0040-1951(78)90099-9.

Ambati, M., Gerasimov, T., De Lorenzis, L., 2014. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput. Mech. 55, 383-405. https://doi.org/10.1007/s00466-014-1109-y.

Borden, M.J., Verhoosel, C.V., Scott, M.A., et al., 2012. A phase-field description of dynamic brittle fracture. Comput. Methods Appl. Mech. Eng. 217-220, 77-95. https://doi.org/10.1016/j.cma.2012.01.008.

Borden, M.J., Hughes, T.J., Landis, C.M., et al., 2014. A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Comput. Methods Appl. Mech. Eng. 273, 100-118.

Bourdin, B., Francfort, G.A., Marigo, J.J., 2008. The variational approach to fracture. J. Elast. 91, 5-148. https://doi.org/10.1007/s10659-007-9107-3.

Carter, B., Wawrzynek, P., Ingraffea, A., 2000. Automated 3-D crack growth simulation. Int. J. Numer. Meth. Eng. 47, 229-253. https://doi.org/10.1002/(SICI)1097-0207 (20000110)30:47:1<3:229>:AID-NME769>3.0.CO;2-2.

Chang, X., Wang, X., Yang, C., et al., 2024. Simulation and optimization of fracture pattern in temporary plugging fracturing of horizontal shale gas wells. Fuel 359, 130378. https://doi.org/10.1016/j.fuel.2023.130378.

Choo, J., Sohail, A., Fei, F., et al., 2021. Shear fracture energies of stiff clays and shales. Acta Geotech. 16, 2291-2299. https://doi.org/10.1007/s11440-021-01145-5.

Choo, J., Sun, Y., Fei, F., 2023. Size effects on the strength and cracking behavior of flawed rocks under uniaxial compression: from laboratory scale to field scale. Acta Geotech. 18, 3451-3468. https://doi.org/10.1007/s11440-023-01806-7.

Cui, H., Du, C., Zhang, H., 2023. Applications of phase field methods in modeling fatigue fracture and performance improvement strategies: A review. Metals 13. https://doi.org/10.3390/met13040714.

Cui, X., Wong, L.N.Y., 2023. Extensive propagation of 3D wing cracks under compression. Rock Mech. Rock Eng. 56, 6135-6142. https://doi.org/10.1007/s00603-023-03358-4.

Deng, J., Chen, F., Feng, Y., 2006. Dimensions of spline spaces over T-meshes. J. Comput. Appl. Math. 194, 267-283. https://doi.org/10.1016/j.cam.2005.07.009.

Dyskin, A.V., Sahouryeh, E., Jewell, R.J., et al., 2003. Influence of shape and locations of initial 3-D cracks on their growth in uniaxial compression. Eng. Fract. Mech. 70, 2115-2136. https://doi.org/10.1016/s0013-7944(02)00240-0.

Fei, F., Choo, J., Liu, C., et al., 2022a. Phase-field modeling of rock fractures with roughness. Int. J. Numer. Anal. Meth. Geomech. 46, 841-868. https://doi.org/10.1002/nag.3317.

Fei, F., Mia, M.S., Elbanna, A.E., et al., 2022b. A phase-field model for quasi-dynamic nucleation, growth, and propagation of rate-and-state faults. Int. J. Numer. Anal. Meth. Geomech. 47, 187-211. https://doi.org/10.1002/nag.3465.

Francfort, G.A., Marigo, J.J., 1998. Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids 46, 1319-1342.

Freddi, F., Mingazzi, L., 2022. Mesh refinement procedures for the phase field approach to brittle fracture. Comput. Methods Appl. Mech. Eng. 388 https://doi.org/10.1016/j.cma.2021.114214.

Fu, J.W., Chen, K., Zhu, W.S., et al., 2016. Progressive failure of new modelling material with a single internal crack under biaxial compression and the 3-D numerical simulation. Eng. Fract. Mech. 165, 140-152. https://doi.org/10.1016/j.engfracmech.2016.08.002.

Fu, T.-F., Xu, T., Heap, M.J., et al., 2020a. Mesoscopic time-dependent behavior of rocks based on three-dimensional discrete element grain-based model. Comput. Geotech. 121, 103472 https://doi.org/10.1016/j.compgeo.2020.103472.

Fu, T.-F., Xu, T., Wasantha, P.L.P., et al., 2020b. Time-dependent deformation and fracture evolution around underground excavations. Geomat. Nat. Haz. Risk 11, 2615-2633. https://doi.org/10.1080/19475705.2020.1856202.

Fu, J.-W., Zhang, X.-Z., Zhu, W.-S., et al., 2017. Simulating progressive failure in brittle jointed rock masses using a modified elastic-brittle model and the application. Eng. Fract. Mech. 178, 212-230. https://doi.org/10.1016/j.engfracmech.2017.04.037.
Goswami, S., Anitescu, C., Rabczuk, T., 2020. Adaptive fourth-order phase field analysis for brittle fracture. Comput. Methods Appl. Mech. Eng. 361 https://doi.org/10.1016/j.cma.2019.112808.

Greco, L., Patton, A., Negri, M., et al., 2024. Higher order phase-field modeling of brittle fracture via isogeometric analysis. Eng. Comput. 1–20 https://doi.org/10.1007/s00366-024-01949-5.

Gu, L., Feng, X.-T., Kong, R., et al., 2023. Influence of excavation stress paths on failure feature of deep hard rocks. Acta Geotech. 19, 1107–1128. https://doi.org/10.1007/s11440-023-02003-2.

Hirshikesh, Jansari, C., Kannan, K., et al., 2019. Adaptive phase field method for quasi-static brittle fracture using a recovery based error indicator and quadtree decomposition. Engineering Fracture Mechanics 220: 106599. Doi: 10.1016/j.engfracmech.2019.106599.

Hirshikesh, H., Pramod, A., Waisman, H., et al., 2021. Adaptive phase field method using novel physics based refinement criteria. Comput. Methods Appl. Mech. Eng. 383 https://doi.org/10.1016/j.cma.2021.113874.

Hobbs, D.W. 1964. The tensile strength of rocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 1: 385-396. Doi: 10.1016/0148-9062(64)90005-1.

Jiang, M., Feng, X.-T., Zhao, J., et al., 2024. Time-dependent failure characteristics of natural jointed granite of deep tunnel under different dip angles conditions. Acta Geotech. 19, 3075–3093. https://doi.org/10.1007/s11440-023-02201-y.

Kang, H., Xu, J., Chen, F., et al., 2015. A new basis for PHT-splines. Graph. Model. 82, 149–159. https://doi.org/10.1016/j.gmod.2015.06.011.

Khaliil, Z., Elghazouli, A.Y., Martínez-Paneda, E., 2022. A generalised phase field model for fatigue crack growth in elastic-plastic solids with an efficient monolithic solver. Comput. Methods Appl. Mech. Eng. 388 https://doi.org/10.1016/j.cma.2021.114286.

Kiran, R., Nguyen-Thanh, N., Yu, H., et al., 2024. Adaptive isogeometric analysis-based phase-field modeling of interfacial fracture in piezoelectric composites. Eng. Fract. Mech. 288 https://doi.org/10.1016/j.engfracmech.2023.109181.

Li, C.C., 2024. Field observations and interpretation of extensional fracture in hard rock surrounding deep underground openings. Deep Resour. Eng. 1, 100006 https://doi.org/10.1016/j.deepre.2024.100006.

Li, X., Fan, G., 2024. On strain localization of aeolian sand in true triaxial apparatus. Acta Geotech. 19, 3115–3128. https://doi.org/10.1007/s11440-024-02273-4.

Li, Y., Hu, S., Sun, X., et al., 2017. A review: applications of the phase field method in predicting microstructure and property evolution of irradiated nuclear materials. npj Comput. Mater. 3 https://doi.org/10.1038/s41524-017-0018-y.

Li, W., Nguyen-Thanh, N., Zhou, K., 2020. Phase-field modeling of brittle fracture in a 3D polycrystalline material via an adaptive isogeometric-meshfreeapproach. Int. J. Numer. Meth. Eng. 121, 5042–5065. https://doi.org/10.1002/nme.6509.

Li, Z., Xu, T., Zhao, L., et al., 2024. Enhancing stability analysis of open-pit slopes via integrated 3D numerical modeling and data monitoring. Eng. Fail. Anal. 163, 108495 https://doi.org/10.1016/j.engfailanal.2024.108495.

Li, Y., Yu, T., Natarajan, S., 2022. An adaptive isogeometric phase-field method for brittle fracture in rock-like materials. Eng. Fract. Mech. 263 https://doi.org/10.1016/j.engfracmech.2022.108298.

Li, S., Wang, Y., Peng, C., et al., 2022. A thermodynamically consistent phase field model for mixed-mode fracture in rock-like materials. Comput. Methods Appl. Mech. Eng. 392 https://doi.org/10.1016/j.cma.2022.114642.

Li, S., Wang, Y., Wu, W., 2023. A modified phase-field model for cohesive interface failure in quasi-brittle solids. Int. J. Mech. Sci. 252 https://doi.org/10.1016/j.ijmecsci.2023.108368.

Lu, J., Xu, T., Tang, X., et al., 2022. Nanoindentation-based characterization of micromechanical properties of greenish mudstone from deep Fushun West open-pit mine (Fushun city, China). Geomech. Geophys. Geo-Energy Geo-Resour. 8 https://doi.org/10.1007/s40948-022-00371-9.

Miehe, C., Hofacker, M., Welschinger, F., 2010. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput. Methods Appl. Mech. Eng. 199, 2765–2778. https://doi.org/10.1016/j.cma.2010.04.011.

Miehe, C., Mauthe, S., 2016. Phase field modeling of fracture in multi-physics problems. Part III. Crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media. Comput. Methods Appl. Mech. Eng. 304, 619–655. https://doi.org/10.1016/j.cma.2015.09.021.

Nguyen-Thanh, N., Li, W., Huang, J., et al., 2020. Adaptive higher-order phase-field modeling of anisotropic brittle fracture in 3D polycrystalline materials. Comput. Methods Appl. Mech. Eng. 372 https://doi.org/10.1016/j.cma.2020.113434.

Nuismer, R., 1975. An energy release rate criterion for mixed mode fracture. Int. J. Fract. 11, 245–250. https://doi.org/10.1007/BF00038891.

Paul, K., Zimmermann, C., Mandadapu, K.K., et al., 2020. An adaptive space-time phase field formulation for dynamic fracture of brittle shells based on LR NURBS. Comput. Mech. 65, 1039–1062. https://doi.org/10.1007/s00466-019-01807-y.

Sahouryeh, E., Dyskin, A., Germanovich, L., 2002. Crack growth under biaxial compression. Eng. Fract. Mech. 69, 2187–2198. https://doi.org/10.1016/S0013-7944(02)00015-2.

Shao, J.-F., Yu, Z., Liu, Z., et al., 2024. Numerical analysis of thermo-hyromechanical process related to deep geological radioactive repository. Deep Resour. Eng. 1, 100001 https://doi.org/10.1016/j.deepre.2024.100001.

Shou, Y., Zhou, X., Berto, F., 2019. 3D numerical simulation of initiation, propagation and coalescence of cracks using the extended non-ordinary state-based peridynamics. Theor. Appl. Fract. Mech. 101, 254–268. https://doi.org/10.1016/j.tafmec.2019.03.006.

Sih, G.C., 1974. Strain-energy-density factor applied to mixed mode crack problems. Int. J. Fract. 10, 305–321. https://doi.org/10.1007/BF00035493.

Tang, C., Liu, H., Lee, P., et al., 2000. Numerical studies of the influence of microstructure on rock failure in uniaxial compression — Part I: effect of heterogeneity. Int. J. Rock Mech. Min. Sci. 37, 555–569. https://doi.org/10.1016/S1365-1609(99)00122-7.

Tang, H.-D., Zhu, Z.-D., Zhu, M.-L., et al., 2015. Mechanical behavior of 3D crack growth in transparent rock-like material containing preexisting flaws under compression. Adv. Mater. Sci. Eng. 2015, 1–10. https://doi.org/10.1155/2015/193721.

Taylor, M.A., Jain, A.K., Ramey, M.R., 1972. Path dependent biaxial compressive testing of an all-lightweight aggregate concrete. J. Am. Concr. Inst. 69, 758–764. https://doi.org/10.14359/11282.

Tian, D.-L., Zhou, X.-P., 2021. A continuum-kinematics-inspired peridynamic model of anisotropic continua: Elasticity, damage, and fracture. Int. J. Mech. Sci. 199 https://doi.org/10.1016/j.ijmecsci.2021.106413.

Walloth, M., Wollner, W., 2022. A posteriori estimator for the adaptive solution of a quasi-static fracture phase-field model with irreversibility constraints. SIAM J. Comput. 44, B479–B505. https://doi.org/10.48550/arXiv:2106.09469.

Wang, Y., Zhou, X., Shou, Y., 2017. The modeling of crack propagation and coalescence in rocks under uniaxial compression using the novel conjugated bond-based peridynamics. Int. J. Mech. Sci. 128–129, 614–643. https://doi.org/10.1016/j.ijmecsci.2017.05.019.

Wang, Y., Zhou, X., Wang, Y., et al., 2018b. A 3-D conjugated bond-pair-based peridynamic formulation for initiation and propagation of cracks in brittle solids. Int. J. Solids Struct. 134, 89–115. https://doi.org/10.1016/j.ijsolstr.2017.10.022.

Wang, Y., Borja, R.I., Wu, W., 2023b. Dynamic strain localization into a compaction band via a phase-field approach. J. Mech. Phys. Solids 173. https://doi.org/10.1016/j.jmps.2023.105228.

Wang, H., Dyskin, A., Pasternak, E., et al., 2022b. Possible mechanism of spallation in rock samples under uniaxial compression. Eng. Fract. Mech. 269 https://doi.org/10.1016/j.engfracmech.2022.108577.

Wang, F., Huang, H., Zhang, D., et al., 2022a. Cracking feature and mechanical behavior of shield tunnel lining simulated by a phase-field modeling method based on spectral decomposition. Tunn. Undergr. Space Technol. 119 https://doi.org/10.1016/j.tust.2021.104246.

Wang, H., Li, H., Tang, L., et al., 2022c. Fracturing behavior of brittle solids containing 3D internal crack of different depths under ultrasonic fracturing. Int. J. Min. Sci. Technol. 32, 1245–1257. https://doi.org/10.1016/j.ijmstc.2022.09.008.

Wang, J., Wu, S., Chang, X., et al., 2023a. A novel three-dimensional nonlinear unified failure criterion for rock materials. Acta Geotech. https://doi.org/10.1007/s11440-023-11424-w.

Wang, Y.-T., Zhou, X.-P., Kou, M.-M., 2018a. Three-dimensional numerical study on the failure characteristics of intermittent fissures under compressive-shear loads. Acta Geotech. 14, 1161–1193. https://doi.org/10.1007/s11440-018-0709-7.

Wasantha, P.L.P., Xu, T., Xu, B., 2021. Effect of fluid injection rate on the hydraulic fracture propagation characteristics. IOP Conf. Ser.: Earth Environ. Sci. 861 https://doi.org/10.1088/1755-1315/861/3/032095.

Wong, L.N.Y., Einstein, H.H., 2009. Crack coalescence in molded gypsum and carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech. Rock Eng. 42, 475–511. https://doi.org/10.1007/s00603-008-0002-4.

Wu, J.Y., Huang, Y., Zhou, H., et al., 2021. Three-dimensional phase-field modeling of mode I + II/III failure in solids. Comput. Methods Appl. Mech. Eng. 373 https://doi.org/10.1016/j.cma.2020.113537.

Xing, C., Yu, T., Sun, Y., et al., 2023. An adaptive phase-field model with variable-node elements for fracture of hyperelastic materials at large deformations. Eng. Fract. Mech. 281 https://doi.org/10.1016/j.engfracmech.2023.109115.

Xu, J., Li, H., Wang, H., et al., 2023. Experimental study on 3D internal penny-shaped crack propagation in brittle materials under uniaxial compression. Deep Undergr. Sci. Eng. 2, 37–51. https://doi.org/10.1002/dug2.12037.

Xu, T., Tang, C.A., Yang, T.H., et al., 2006. Numerical investigation of coal and gas outbursts in underground collieries. Int. J. Rock Mech. Min. Sci. 43, 905–919.

Xu, T., Xu, Q., Deng, M., et al., 2014. A numerical analysis of rock creep-induced slide: a case study from Jiweishan Mountain, China. Environ. Earth Sci. 72, 2111–2128. https://doi.org/10.1007/s12665-014-3119-7.

Xu, B., Xu, T., Xue, Y., et al., 2022a. Phase-field modeling of crack growth and interaction in rock. Geomech. Geophys. Geo-Energy Geo-Resour. 8 https://doi.org/10.1007/s40948-022-00497-w.

Xu, B., Xu, T., Xue, Y., et al., 2024. Phase field modeling of mixed-mode crack in rocks incorporating heterogeneity and frictional damage. Eng. Fract. Mech., 109936 https://doi.org/10.1016/j.engfracmech.2024.109936.

Xu, T., Yang, T.-H., Chen, C.-F., et al., 2015. Mining induced strata movement and rock behavior in underground coal mine. Geomech. Geophys. Geo-Energy Geo-Resour. 1, 79–89. https://doi.org/10.1007/s40948-015-01002-0.

Xu, Y., Zhou, S., Xia, C., et al., 2022b. A new phase field model for mixed-mode brittle fractures in rocks modified from triple shear energy criterion. Acta Geotech. https://doi.org/10.1007/s11440-022-01589-3.

Xue, Y.C., Xu, T., Zhu, W.C., et al., 2021. Full-field quantification of time-dependent and independent deformation and fracturing of double-notch flawed rock using digital image correlation. Geomech. Geophys. Geo-Energy Geo-Resour. 7 https://doi.org/10.1007/s40948-021-00302-0.

Yang, S.Q., 2011. Crack coalescence behavior of brittle sandstone samples containing two coplanar fissures in the process of deformation failure. Eng. Fract. Mech. 78, 3059–3081. https://doi.org/10.1016/j.ijmecsci.2011.09.002.

Yang, L., Li, S.C., Guo, Y.S., et al., 2010. Influence of space between 3D-cracks on mechanical feature of mortar material. Yantu Lixue/Rock Soil Mech. 31, 32–38. https://doi.org/10.16285/j.rsm.2010.01.034.

Yao, W.-W., Zhou, X.-P., Qian, Q.-H., 2022. From statistical mechanics to nonlocal theory. Acta Mech. 233, 869–887. https://doi.org/10.1007/s00707-021-03123-0.
Ye, Z., Wang, J.G., Zhao, Y., 2020. A thermal-hydraulic-mechanical coupling study of heat extraction from the geothermal reservoir with a discrete fracture network. *Geofluids* 2020, 1–18. https://doi.org/10.1155/2020/8875918.

Yin, P., Wong, R.H.C., Chau, K.T., 2014. Coalescence of two parallel pre-existing surface cracks in granite. Int. J. Rock Mech. Min. Sci. 68, 66–84. https://doi.org/10.1016/j.ijrmms.2014.02.011.

Yu, X.-Y., Xu, T., Heap, M.J., et al., 2023. Time-dependent virtual crack model of rock with application to slope stability. Eng. Anal. Bound. Elem. 154, 172–185. https://doi.org/10.1016/j.enganabound.2023.05.043.

Yuan, Y., Xu, T., Heap, M.J., et al., 2021. A three-dimensional mesoscale model for progressive time-dependent deformation and fracturing of brittle rock with application to slope stability. Comput. Geotech. 135, 104160. https://doi.org/10.1016/j.compgeo.2021.104160.

Yue, Q., Zhou, W., Wang, Q., et al., 2022. An adaptive phase-field model based on bilinear elements for tensile-compressive-shear fracture. Comput. Math. Appl. 105, 112–135. https://doi.org/10.1016/j.camwa.2021.11.010.

Zhang, Y., Feng, X.-T., Yao, Z., et al., 2023a. Failure characteristics and development mechanism of fault rockburst in a deep TBM tunnel: A case study. Acta Geotech. 18, 5575–5596. https://doi.org/10.1007/s11440-023-01883-8.

Zhang, T.Y., Tong, P., 1996. Fracture mechanics for a mode III crack in a piezoelectric material. Int. J. Solids Struct. 33, 343–359. https://doi.org/10.1016/0020-7683(95)00046-D.

Zhang, Z.N., Wang, D.Y., Zheng, H., et al., 2012. Interactions of 3D embedded parallel vertically inclined cracks subjected to uniaxial compression. Theor. Appl. Fract. Mech. 61, 1–11. https://doi.org/10.1016/j.tafmec.2012.08.001.

Zhang, Z., Yu, S., Yang, E., 2020. Theoretical and numerical study on three-dimensional internal crack specimen under uniaxial compression. Geotech. Geol. Eng. 38, 4059–4072. https://doi.org/10.1007/s10706-020-01278-6.

Zhang, Z., Wang, H., Zhong, L., et al., 2023b. Experimental and numerical research of the behavior of 3D internal cracks under uniaxial compression using 3D-ILC technology. Theor. Appl. Fract. Mech. 124 https://doi.org/10.1016/j.tafmec.2023.103819.

Zhang, D., Zhang, B., Wang, W., et al., 2016. Study on the built-in ellipse 3D crack propagation under uniaxial compression by meshless method. Chin. J. Appl. Mech. 33, 483–489. https://doi.org/10.11776/cjam.33.03.B026.

Zhou, X.P., Bi, J., Qian, Q.H., 2014. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech. Rock Eng. 48, 1097–1114. https://doi.org/10.1007/s00603-014-0627-4.

Zhou, G., Xu, T., Heap, M.J., et al., 2020a. A three-dimensional numerical meso-approach to modeling time-independent deformation and fracturing of brittle rocks. Comput. Geotech. 117, 103274 https://doi.org/10.1016/j.compgeo.2019.103274.

Zhou, X.-P., Zhang, J.-Z., 2020. AE event rate characteristics of flawed granite: from damage stress to ultimate failure. Geophys. J. Int. 222, 795–814. https://doi.org/10.1093/gji/ggaa207.

Zhou, X.-P., Zhang, J.-Z., Wong, L.N.Y., 2018. Experimental study on the growth, coalescence and wrapping behaviors of 3D cross-embedded flaws under uniaxial compression. Rock Mech. Rock Eng. 51, 1379–1400. https://doi.org/10.1007/s00603-018-1406-4.

Zhou, T., Zhu, J.B., Ju, Y., et al., 2019. Volumetric fracturing behavior of 3D printed artificial rocks containing single and double 3D internal flaws under static uniaxial compression. Eng. Fract. Mech. 205, 190–204. https://doi.org/10.1016/j.engfracmech.2018.11.030.

Zhou, S., Zhuang, X., Rabczuk, T., 2020b. Phase field method for quasi-static hydro-fracture in porous media under stress boundary condition considering the effect of initial stress field. Theor. Appl. Fract. Mech. 107 https://doi.org/10.1016/j.tafmec.2020.102523.

Zhuang, X., Zhou, S., Huynh, G.D., et al., 2022. Phase field modeling and computer implementation: A review. Eng. Fract. Mech. 262 https://doi.org/10.1016/j.engfracmech.2022.108234.