Phase-field modeling of geologic fractures

Geologic fractures such as joints, faults, and slip surfaces govern the stability and performance of many subsurface systems in the built environment. As such, a variety of approaches have been developed for computational modeling of geologic fractures. Yet none of them lends itself to a straightforward utilization with the classical finite element method widely used in practice. Over the past decade, phase-field modeling has become a popular approach for simulating fracture, because it can be implemented simply with the standard finite element method without any surface-tracking algorithms. However, the standard phase-field formulations do not incorporate several critical features of geologic fractures, including frictional contact, pressure-dependence, quasi-brittleness, mode-mixity, and their combined impacts on cracking. This article provides a brief report of a novel phase-field model that incorporates these features of geologic fractures in a well-verified and validated manner. Remarkably, the phase-field model allows one to simulate the combination of cohesive tensile fracture and frictional shear fracture without any algorithms for surface tracking and contact constraints. It is also demonstrated how phase-field modeling enables us to gain insights into geologic fractures that are challenging to investigate experimentally.

1 Introduction

Subsurface systems are replete with fractures of various kinds. Examples include joints and faults, which are tensile and shear fractures, respectively, from the viewpoint of fracture mechanics (Pollard and Fletcher, 2005; Schultz, 2019). Slip surfaces can also be viewed as shear discontinuities, of which propagation can be described by fracture mechanics (Palmer and Rice, 1973; Puzrin and Germanovich, 2005). The deformation and growth of these fractures govern the stability and performance of a wide range of geotechnical systems from slopes and earth retaining systems to tunnels and caverns. Also, modern subsurface energy technologies such as unconventional resource recovery and enhanced geothermal systems purposely generate fractures in deep underground deposits. Unlike other types of fractures, the geologic fractures in these problems often involve frictional contact and show significant pressure dependence. Also, fractures in geomaterials are quasi-brittle instead of purely brittle, and their critical fracture energies in tensile and shear usually differ by an order of magnitude. Computational modeling is essential to addressing such complex behavior of geologic fractures.

A variety of approaches have been developed for computational modeling of geologic fractures. Following the classical approaches in computational fracture mechanics, earlier works treated geologic fractures as sharp discontinuities across which displacement and/or strain fields are discontinuous. These approaches honor the sharpness of crack geometry, being consistent with fracture mechanics principles derived from sharp cracks. However, it is challenging to represent sharp crack geometry in numerical methods. For example, the standard finite element method can only model cracks aligned with element boundaries. To overcome this limitation, embedded finite elements (e.g., Regueiro and Borja, 1999; Foster et al., 2007) and extended/generalized finite elements (e.g., Liu and Borja, 2009; Sanborn and Prévost, 2011) have been proposed for modeling geologic fractures that pass through inside elements. However, these non-standard finite element methods require not only additional shape functions for enrichment but also sophisticated algorithms (e.g., the level-set method) to track crack surfaces.

Over the past decade, the phase-field method has become a popular approach for computational modeling of fractures (e.g., Bourdin et al., 2008; Miehe et al., 2010b; Borden et al., 2012). The phase-field method approximates a sharp crack surface as a diffuse interface and describes its evolution by a partial differential equation formulated from fracture mechanics concepts. In this way, the phase-field method can simulate complex cracking processes such as kinking, branching, and coalescence without any surface-tracking algorithms, being solvable by the standard finite elements widely used in practice. While a phase-field model demands a significant computational cost, it has been increasingly affordable and widely applied to many kinds of fracture problems in various materials. Nevertheless, the standard phase-field model of fracture and its variants do not incorporate several distinct features of geologic fractures: frictional contact, pressure-dependence, quasi-brittleness, mode-mixity, and roughness, among others. Therefore, although these phase-field models have been employed in many geomechanical applications (e.g., Lee et al., 2016; Zhang et al., 2017; Choo and Sun, 2018a; b; Ha et al., 2018; Zhou et al., 2018), there is a need for new phase-field models tailored to geologic fractures.

Recently, the author and his coworkers have developed a class of novel phase-field models for geologic fractures, which incorporate frictional contact (Fei and Choo, 2020a), quasi-brittle shear fracture with friction and pressure effects (Fei and Choo, 2020b), mixed-mode fracture (Fei and Choo, 2021), roughness effects (Fei et al., 2022), and inertia effects with rate- and state-dependent friction Fei et al. (2023). This report reviews these phase-field models with particular attention to the double-phase-field model (Fei and Choo, 2021) which allows one to simulate the combination of cohesive tensile fracture and frictional shear fracture without any algorithms for surface tracking and contact constraints. The double-phase-field model has not only been well-validated against various experimental data (Fei et al., 2021) but has also been applied to gain insights into geologic fractures that are challenging to investigate experimentally (Choo et al., 2023; Sun et al., 2024). The remainder of this paper introduces the formulation, validation, and applications of the phase-field model.

2 Phase-field model

2.1 Geometry approximation

The departure point of the phase-field modeling of fracture is to diffusely approximate the sharp geometry of a crack (e.g., a line crack in a 2D domain, and a surface crack in a 3D domain) by a phase field, $d$, whose value ranges from 0 to 1. Usually, the phase field is defined such that $d=0$ in an intact (undamaged) region, $d=1$ in a fully cracked region, and $0<d<1$ in a transition region which can also be viewed as a partially damaged region. The distribution of the phase-field value is determined by a crack density functional, $\mathrm{\Gamma }(d,\nabla d)$. An important parameter of the crack density functional is the so-called length regularization parameter, $L$, which determines the thickness of the diffuse approximation. By definition, the phase-field approximation converges to sharp geometry as $L$ approaches zero, and hence $L$ should be chosen to be sufficiently small to mimic a sharp crack. At the same time, as $L$ becomes smaller, a finer discretization is necessary to obtain a sufficiently accurate numerical solution. As such, one must determine the value of $L$ considering the balance between solution accuracy and computational efficiency.

Mixed-mode fractures—the combination of tensile and shear fractures—are common in geomaterials. Standard phase-field models approximate both tensile and shear fractures with the same phase-field variable. This approach works well when the characteristics of tensile and shear fractures are similar, but it is undesirable when the two types of fractures exhibit contrasting characteristics. The latter is the case for geologic fractures: Shear fractures in geomaterials have much higher fracture energies than tensile fractures and involve frictional effects. Therefore, Fei and Choo (2021) have proposed a double-phase-field model that approximates tensile and shear fractures separately by two distinct phase fields. Figure 1 illustrates the double-phase-field approximation of mixed-mode fractures.

Figure 1

Figure 1. Phase-field approximation of mixed-mode (tensile and shear) fractures. Source: Fei and Choo (2021).

2.2 Formulation

Phase-field modeling describes the cracking process as the evolution of the phase-field variable over time. The governing equation for the phase-field evolution can be derived either from a variational principle or microforce theory. The latter approach allows one to easily incorporate complex features such as quasi-brittleness, contact dependence, and friction effects. From microforce theory, the governing equations for tensile and shear fractures, respectively, are derived as.

$\nabla \cdot \left(\frac{\partial \psi \left(\mathit{ϵ},d_I,\nabla d_I,d_{II},\nabla d_{II}\right)}{\partial \nabla d_I}\right)-\frac{\partial \psi \left(\mathit{ϵ},d_I,\nabla d_I,d_{II},\nabla d_{II}\right)}{\partial d_I}=0,(1)$

$\nabla \cdot \left(\frac{\partial \psi \left(\mathit{ϵ},d_I,\nabla d_I,d_{II},\nabla d_{II}\right)}{\partial \nabla d_{II}}\right)-\frac{\partial \psi \left(\mathit{ϵ},d_I,\nabla d_I,d_{II},\nabla d_{II}\right)}{\partial d_{II}}=0,(2)$

where $d_I$ and $d_{II}$ are the phase fields representing tensile and shear fractures, respectively, $\mathit{ϵ}$ is the (infinitesimal) strain tensor, and $\psi$ is the potential energy density function. The potential energy density function is comprised of three parts as

$\psi ={\psi }^{\text{e}}+{\psi }^{\text{f}}+{\psi }^{\text{d}},(3)$

where ${\psi }^{\text{e}}$ is the elastic strain energy, ${\psi }^{\text{f}}$ is the energy dissipation due to frictional slip, and ${\psi }^{\text{d}}$ is the energy dissipation due to crack surface generation. Substituting the specific expressions for these terms (see Fei and Choo (2021) for details) into Eqs 1, 2, the governing equations become.

$-g_I^{\prime }\left(d_I\right){\mathcal{H}}_I-\frac{{\mathcal{G}}_I}{\pi L}\left(2L^2\nabla \cdot \nabla d_I-2+2d_I\right)=0,(4)$

$-g_{II}^{\prime }\left(d_{II}\right){\mathcal{H}}_{II}-\frac{{\mathcal{G}}_{II}}{\pi L}\left(2L^2\nabla \cdot \nabla d_{II}-2+2d_{II}\right)=0.(5)$

Here, $g_I^{\prime }(d_I)$ and $g_{II}^{\prime }(d_{II})$ are the degradation functions, ${\mathcal{G}}_I$ and ${\mathcal{G}}_{II}$ are the fracture energies, and ${\mathcal{H}}_I$ and ${\mathcal{H}}_{II}$ are the crack driving forces for tensile and shear fractures, respectively. Considering that ${\mathcal{G}}_I$ and ${\mathcal{G}}_{II}$ of geomaterials can be significantly different, the crack driving forces are calculated based on the $\mathcal{F}$-criterion (Shen and Stephansson, 1994) as

$\theta =\arg {\max }_{\theta }\left[\mathcal{F}\left(\theta \right)\right]{|}_{\mathit{ϵ}},\text{where}\mathcal{F}\left(\theta \right)≔\frac{{\mathcal{H}}_I\left(\mathit{ϵ},\theta \right)}{{\mathcal{G}}_I}+\frac{{\mathcal{H}}_{II}\left(\mathit{ϵ},\theta \right)}{{\mathcal{G}}_{II}},(6)$

where $\theta$ is the angle between the crack normal direction and the major principal stress direction on the slip plane. It is noted that the $\mathcal{F}$ criterion lets the mixed-mode fracture develop such that it minimizes energy dissipation. The specific expressions for ${\mathcal{H}}_I(\mathit{ϵ},\theta )$ and ${\mathcal{H}}_{II}(\mathit{ϵ},\theta )$ depend on the contact condition of the material point, and they are adapted from phase-field models for cohesive tensile fracture (Wu, 2017) and frictional shear fracture (Fei and Choo, 2020b), respectively. See Fei and Choo (2021) for details.

The degradation functions, $g_I^{\prime }(d_I)$ and $g_{II}^{\prime }(d_{II})$, play a critical role in the stress-strain response of the phase-field model. The degradation functions of the standard phase-field models give rise to stress-strain responses that are dependent on the length parameter, $L$. This is because the standard phase-field models regularize brittle fracture having stress singularity at the crack tip and the length parameter controls the degree of the regularization of stress singularity. In other words, as the length parameter becomes larger, the “strength” of the crack tip—regularized stress singularity due to the phase-field approximation—becomes smaller. However, for cohesive (quasi-brittle) fracture, the stress-strain response should be independent of the length parameter, $L$; otherwise, the fracture energy is dependent on the length parameter, violating the premise of cohesive fracture.

To make the stress-strain response of a phase-field model independent of the length parameter, one must employ a particular degradation function derived to provide a length-insensitive stress-strain response. The double-phase-field model employs two such degradation functions, one derived for cohesive tensile fracture (Wu, 2017) and the other derived for frictional shear fracture (Fei and Choo, 2020b). In this way, it is possible to directly use the tensile and shear strength properties measured from experiments, without concerning the chosen value of $L$. This feature greatly helps the calibration of the phase-field model with laboratory data.

Last but not least, it should be emphasized that the double-phase-field formulation can be solved well with the standard finite element method. (See Fei and Choo (2021) for the full finite element formulation.) Specifically, one can use a three-field finite element discretization in which the displacement field and the two phase fields are the unknown variables, with a material (integration) point update algorithm to calculate contact-dependent stresses. The displacement fields and phase fields can be solved in a staggered manner for robustness (Miehe et al., 2010a). In doing so, no algorithm is necessary for tracking the crack geometry or imposing the contact constraints. This ease of implementation with the standard finite element method is a distinct advantage of the phase-field method.

2.3 Material parameters

The material parameters of the double-phase-field model are summarized as follows.

$•$ Elasticity parameters: Assuming isotropic linear elasticity, the model takes two elasticity moduli. The elasticity moduli can be measured through standard laboratory experiments.

$•$ Tensile fracture parameters: For modeling cohesive tensile fracture, the two parameters of traction–separation law, namely, the tensile strength and tensile (Mode I) fracture energy, are required. The fracture energy can be calculated as the area of the softening region (Wu, 2017).

$•$ Shear fracture parameters: Adopting the Mohr-Coulomb criterion, the model takes the cohesion strength and peak friction angle to represent pressure-dependent shear strength. For simplicity, the residual friction angle is assumed to be equal to the peak friction angle, while this assumption can be relaxed for rough fracture involving shear-induced dilation, as in Fei et al. (2022). Lastly, the shear (Mode II) fracture energy is necessary to furnish slip-weakening law. The shear fracture energy can be estimated from the post-peak response in the ring shear test or inferred from those in the triaxial compression test (Choo et al., 2021).

It is noted that when fracture test data are available, the tensile and shear fracture energies can also be calibrated to match the fracture data.

2.4 Validation

The double-phase-field model has been validated with several experimental data on mixed-mode fracture in rock specimens with preexisting flaws under compression. Fei and Choo (2021) and Fei et al. (2021) present qualitative and quantitative validations with data from uniaxial compression tests on rock specimens with various flaw configurations. These papers have shown that the double-phase-field model can simulate complex mixed-mode cracking patterns and stress-strain responses similar to the experimental data. Yet these validations were restricted to 2D (planar) fracture under uniaxial compression. More recently, in Sun et al. (2024), the double-phase-field model has been validated with fractures in two types of more general conditions, namely, 1) 2D fracture under true triaxial compression, and 2) 3D fracture under uniaxial compression. It is noted that although 3D fractures under true triaxial compression are the most general condition, they cannot be characterized experimentally through the existing techniques. In the following, the validation of 2D fractures under true triaxial compression is reported briefly.

To validate the double-phase-field method for cracking behavior under true triaxial compression, it was applied to simulate a true triaxial compression test on a fine sandstone specimen with two coplanar, fully-penetrating 2D flaws (Zhou et al., 2021). The width, length, and ligament length of each flaw were 1 mm wide, 16 mm, and 8 mm, respectively, and the flaw inclination angle was $45°$ from the horizontal. More details about the experimental setup can be found in Zhou et al. (2021). The material parameters of the double-phase-field model were adopted directly from Zhou et al. (2021) if they were measured. Other material parameters that had not been measured by Zhou et al. (2021) (Young’s modulus and the tensile and shear fracture energies) were calibrated based on the stress-strain data in the reference. It is noted that the tensile and shear fracture energies were calibrated to be ${\mathcal{G}}_I=98$ ${\mathrm{J}/\mathrm{m}}^2$ and ${\mathcal{G}}_{II}=1500$ ${\mathrm{J}/\mathrm{m}}^2$, respectively, indicating that the shear fracture energy was about 15 times greater than the tensile fracture energy. This ratio of the fracture energies is typical for rocks. For numerical simulation, the phase-field length parameter was set as $L=0.3$ mm. To ensure the accuracy of the solution, the elements near the flaw were locally refined such that the element size $h$ satisfies $L/h\ge 2$. Then, the loading protocol of the true triaxial experiment was replicated through a combination of stress and displacement controls.

Figure 2 compares the simulation and experimental results of mixed-mode fracture from two coplanar flaws under true triaxial compression. One can see that the crack type, crack geometry, and crack coalescence pattern of the simulation and experimental results are quite similar. From the outer tips of the flaws, primary and secondary tensile cracks emerged, and then oblique shear cracks developed from both the inner and outer tips of the flaws. Also, from the shear crack fronts, anti-wing cracks grew toward the direction of the major principal stress. Eventually, the two flaws coalesced into a shear crack, and coplanar shear cracks nucleated from the outer tips of the flaws and grew along the flaw planes. Apart from this qualitative agreement in terms of the cracking pattern, the simulation and experimental results show good quantitative agreement in terms of the stress-strain curves. It is noted that the same conclusion has been reached from other validation studies (Fei and Choo, 2021; Fei et al., 2022; Sun et al., 2024).

Figure 2

Figure 2. Validation of the double-phase-field model with experimental data on cracking from two coplanar flaws under true triaxial compression (Zhou et al., 2021). (A) Cracking patterns from the simulation and experiments. (B) 3D mixed-mode processes in the simulation (Points A–E correspond to those in the stress-strain curve below). (C) Difference between the major principal stress $({\sigma }_1)$ and $({\sigma }_2)$ versus strain from the simulation and experiments. (Note that the experimental stress-strain curve is redrawn without the initial nonlinear portion related to the closure of preexisting defects.) Source: Sun et al. (2024).

3 Applications of phase-field model

This section introduces two examples of how the phase-field model has been applied to investigate mixed-mode fracture in rocks which are extremely challenging to characterize by the existing experimental methods. The two examples are: 1) size effect on mixed-mode fracture in rocks with preexisting flaws (Choo et al., 2023), and 2) intermediate principal stress effect on the 3D cracking behavior of rocks under true triaxial compression (Sun et al., 2024).

3.1 Size effect on mixed-mode fracture in rocks

Laboratory specimens with preexisting flaws, such as the specimen simulated in the validation example, have been conventionally used as small-scale analogs of rock masses, and their failure behavior under compression has been extensively studied through experimental and numerical methods. However, no study has been concerned with the energetic size effect—determined by the relative size between the fracture process zone (FPZ) and the structure size—on the failure behavior of rock specimens with preexisting flaws. In quasi-brittle materials like rocks, the size of FPZ is significant in relatively small structures like laboratory specimens, and it gives rise to rather ductile failure behavior that deviates from the description of linear elastic fracture mechanics (LEFM). However, in large structures like field-scale rock masses, the size of FPZ is negligible compared with the structural size; in this case, the structure fails like a purely brittle material described by LEFM. This energetic size effect was first revealed in Bažant (1984) in the context of tensile failure in quasi-brittle materials and has been investigated for various types of quasi-brittle failures. Nevertheless, the energetic size effect on mixed-mode fracture in rocks under compression has not been investigated systematically. One primary reason may be that it is highly challenging to examine the size effect on mixed-mode fracture in flawed rocks under compression with the existing experimental methods. For example, it is extremely difficult to prepare field-scale (meter-scale) rock specimens with preexisting flaws and characterize the mixed-mode fracture under compression.

The double-phase-field model can be an ideal tool to investigate the energetic size effect on mixed-mode fracture in flawed rocks. The reasons are: 1) it can simulate complex crack patterns without geometry tracking algorithms, 2) it can distinguish between tensile and shear cracks naturally, and 3) it combines phase-field models of cohesive tensile and shear fractures in quasi-brittle materials. As such, in Choo et al. (2023), the double-phase-field model is leveraged to perform the first systematic investigation of mixed-mode fracture in rocks under compression. It should be noted that unlike tensile fracture where size affects the strength but not crack morphology, mixed-mode fracture in rocks under compression may be subjected to two types of size effects in terms of crack morphology as well as strength. Indeed, experimental studies have found different cracking patterns in rocks and rock-like materials with the same specimen and flaw geometry (e.g., Lee and Jeon, 2011), and the difference may be attributed to the different brittleness/ductility of the materials. For this reason, it is intriguing to explore whether structure size also affects the mixed-mode cracking pattern in flawed rocks under compression.

To investigate the size effect on the mixed-mode cracking behavior of rock masses under compression, the double-phase-field model is used to simulate a series of uniaxial compression tests on single- and double-flawed rock specimens of seven different sizes with geometrical similarity. The sizes of the seven specimens were determined such that they range from millimeter-scale specimens to meter-scale specimens. Specifically, setting a 76.2 mm wide and 152.4 mm tall laboratory specimen as a reference, seven values of the scaling factor $(D)$ were considered: $D$ = 0.25, 0.5, 1, 2, 4, 8, and 16. ($D$ = 1 corresponds to the reference size.). The material parameters of the double-phase-field model were assigned from those calibrated to molded gypsum in Fei and Choo (2021), which have been shown to reproduce laboratory test data on the mixed-mode cracking behavior of molded gypsum (Bobet and Einstein, 1998) in both qualitative and quantitative manners. Therefore, when $D$ is 1, the simulation result must be analogous to the laboratory test data. It is of interest how the result would become different when the value of $D$ is much greater or smaller than 1.

Figure 3 presents the simulated mixed-mode cracking patterns in double-flawed specimens of seven sizes ($D$ = 0.25, 0.5, 1, 2, 4, 8, and 16). It can be seen that the cracking patterns in the specimens of $D$ = 1, 2, and 4 are the same as the cracking pattern observed in the laboratory-scale specimen made of molded gypsum (Bobet and Einstein, 1998). This is expected because the material parameters are calibrated to the molded gypsum. However, when the specimen size is smaller than the laboratory scale (i.e., when $D$ = 0.25 and 0.5), the mixed-mode fracture exhibits a different pattern, and this pattern resembles the cracking pattern observed in Carrara marble at the laboratory scale (Wong, 2008). Notably, the FPZ size of Carrara marble is known to be greater than that of molded gypsum (Wong and Einstein, 2009), so this transition in the cracking pattern can be interpreted as the consequence of the fact that the FPZ size becomes relatively more significant as the specimen size becomes smaller. Conversely, when the specimen size is significantly larger than the laboratory scale (i.e., when $D$ = 8 and 16), another kind of cracking pattern emerges, in which shear fracture is nearly absent. Remarkably, more or less the same cracking pattern has been observed in laboratory specimens made of PMMA (Ingraffea and Heuze, 1980), which is a highly brittle material. This is also consistent with the energetic size effect: When the structure size is large, the FPZ size becomes negligible, making the structure highly brittle. These results thus suggest that the energetic size effect exists in the failure behavior of rock masses under compression, and they provide insight into how to bridge laboratory-scale observations and field-scale processes. For example, for a laboratory-scale investigation of the failure behavior of a field-scale rock mass, it would be desirable to use a more brittle material (e.g., PMMA) with the same structural geometry.

Figure 3

Figure 3. Mixed-mode cracking patterns in double-flawed specimens with different scaling factors: $D$ = 0.25, 0.5, 1, 2, 4, 8, and 16. Also presented are cracking patterns observed in the laboratory-scale specimens of Carrara marble (Wong, 2008), molded gypsum (Bobet and Einstein, 1998), and PMMA (Ingraffea and Heuze, 1980). Source: Choo et al. (2023).

3.2 Intermediate principal stress effect on the 3D cracking behavior of rocks

So far, mixed-mode fracture in rocks has mostly been studied under uniaxial compression. However, uniaxial compression is far from in-situ stress conditions of rocks. Underground rocks are usually under true triaxial stress conditions in which all three (major, intermediate, and minor) principal stresses are compressive and distinct. Even at surfaces such as excavation boundaries, the stress condition is biaxial in which only the minor principal stress is zero. Therefore, the intermediate principal stress, ${\sigma }_2$, is usually significant in in-situ stress states, and hence its effect should be considered properly.

Besides the stress condition, the 3D geometrical feature of the preexisting flaw(s) is a critical factor of cracking behavior under compression. Most existing studies have studied cracks emanating from 2D (planar, penetrating) preexisting flaws. However, a few studies have shown that when the preexisting flaw is 3D (internally embedded), new types of cracking patterns such as petal cracks and crack wrapping emerge (e.g., Dyskin et al., 2003; Yin et al., 2014; Lu et al., 2015). Yet such 3D cracking behavior has only been studied under uniaxial and biaxial compression regimes. Therefore, the 3D cracking behavior of flawed rocks under a wide range of true triaxial stress conditions remains elusive.

Indeed, it is virtually impossible to experimentally characterize 3D cracking processes in rocks under true triaxial stress conditions. This is because while a high-speed imaging system is necessary to identify a mixed-mode fracture process faithfully, it cannot be applied to a rock specimen under a true triaxial cell. Alternatively, high-fidelity numerical simulations based on sound physical principles can be employed to investigate 3D cracking processes under true triaxial compression. To this end, the double-phase-field method can again serve as an ideal tool because of its ability to simulate complex 3D tensile and shear fractures individually.

In Sun et al. (2024), a series of numerical true triaxial compression tests were performed on single-flawed and double-flawed cubic specimens. The flaws were internally embedded such that fully 3D cracking patterns could develop. The test protocol was designed to investigate how the orientation of ${\sigma }_2$ with respect to the flaw and the magnitude of ${\sigma }_2$ affect the cracking pattern and the peak stress. Among many simulation data in Sun et al. (2024), cracking patterns in the double-flawed specimen (flaw inclination angle = 45$°$) under different orientations of ${\sigma }_2$ are shown in Figure 4. When ${\sigma }_2$ was parallel to the flaw plane, tensile wing cracks developed from the flaw tips, and planar shear cracks grew in the flaw direction. When ${\sigma }_3$ was parallel to the flaw plane, tensile wing cracks grew from the flaw edge. Also the so-called fish-fin cracks and several conjugate pairs of shear cracks developed. In these two cases, the two preexisting flaws coalesced. However, when ${\sigma }_1$ is parallel to the flaw plane, the two flaws did not coalesce. These differences indicate that the orientation of ${\sigma }_2$ has a profound effect on the mixed-mode cracking behavior of rocks. Apart from this, there are several new observations and findings about the control of ${\sigma }_2$ on the 3D mixed-mode cracking behavior of rocks—see Sun et al. (2024) for details.

Figure 4

Figure 4. 3D cracking patterns in a double-flawed specimen (inclination angle: 45$°$. (A) The intermediate principal stress $({\sigma }_2)$ is parallel to the flaw plane. (B) The minor principal stress $({\sigma }_3)$ is parallel to the flaw plane. (C) The major principal stress $({\sigma }_1)$ is parallel to the flaw plane. Source: Sun et al. (2024).

Remarkably, based on the observations made from the numerical true triaxial compression tests, three mechanisms of the cracking behavior of 3D flawed rocks under true triaxial conditions were proposed. First, the normal stress on the flaw controls the tensile fracture. Second, the Coulomb stress on the flaw controls the shear fracture. Third, the Coulomb-to-normal stress ratio affects the mixed-mode cracking pattern which controls the peak stress. These mechanisms would not have been uncovered without the use of the double-phase-field model.

4 Closure

This paper has introduced a novel phase-field approach for computational modeling of geologic fractures and its application to investigating the mechanics of complex rock fractures. Particular emphasis has been placed on the double-phase-field model developed for mixed-mode fracture in rocks and similar quasi-brittle materials. The double-phase-field model has two standout features: 1) it can simulate complex fractures without sophisticated algorithms for crack geometry tracking and contact constraints, and 2) it can naturally distinguish between tensile and shear fractures. The combination of these two features makes the model an ideal tool for studying the mechanics of complex mixed-mode fractures that cannot be well characterized by the existing experimental methods. The model is however not without drawbacks. The most critical drawback is that phase-field modeling entails a significant computational cost because an extremely fine discretization is necessary around the crack. Another drawback is that it is less straightforward to extract discrete quantities (e.g., aperture) from phase-field approximated fractures. Work is underway to overcome these drawbacks through the combined use of novel formulations, efficient algorithms, machine learning methods, and modern computing platforms.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

JC: Writing–original draft, Conceptualization, Funding acquisition, Methodology.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. RS-2023-00209799).

Acknowledgments

The author wishes to express his deep thanks to his former students Fan Fei and Yuan Sun for their collaboration on the works reported in this paper.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Bažant, Z. P. (1984). Size effect in blunt fracture: concrete, rock, metal. J. Eng. Mech. 110, 518–535. doi:10.1061/(asce)0733-9399(1984)110:4(518)

CrossRef Full Text | Google Scholar

Bobet, A., and Einstein, H. (1998). Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int. J. Rock Mech. Min. Sci. 35, 863–888. doi:10.1016/s0148-9062(98)00005-9

CrossRef Full Text | Google Scholar

Borden, M. J., Verhoosel, C. V., Scott, M. A., Hughes, T. J., and Landis, C. M. (2012). A phase-field description of dynamic brittle fracture. Comput. Methods Appl. Mech. Eng. 217, 77–95. doi:10.1016/j.cma.2012.01.008

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Choo, J., Sohail, A., Fei, F., and Wong, T.-f. (2021). Shear fracture energies of stiff clays and shales. Acta Geotech. 16, 2291–2299. doi:10.1007/s11440-021-01145-5

CrossRef Full Text | Google Scholar

Choo, J., and Sun, W. (2018a). Coupled phase-field and plasticity modeling of geological materials: from brittle fracture to ductile flow. Comput. Methods Appl. Mech. Eng. 330, 1–32. doi:10.1016/j.cma.2017.10.009

CrossRef Full Text | Google Scholar

Choo, J., and Sun, W. (2018b). Cracking and damage from crystallization in pores: coupled chemo-hydro-mechanics and phase-field modeling. Comput. Methods Appl. Mech. Eng. 335, 347–379. doi:10.1016/j.cma.2018.01.044

CrossRef Full Text | Google Scholar

Choo, J., Sun, Y., and 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. doi:10.1007/s11440-023-01806-7

CrossRef Full Text | Google Scholar

Dyskin, A., Sahouryeh, E., Jewell, R., Joer, H., and Ustinov, K. (2003). Influence of shape and locations of initial 3-D cracks on their growth in uniaxial compression. Eng. Fract. Mech. 70, 2115–2136. doi:10.1016/s0013-7944(02)00240-0

CrossRef Full Text | Google Scholar

Fei, F., and Choo, J. (2020a). A phase-field method for modeling cracks with frictional contact. Int. J. Numer. Methods Eng. 121, 740–762. doi:10.1002/nme.6242

CrossRef Full Text | Google Scholar

Fei, F., and Choo, J. (2020b). A phase-field model of frictional shear fracture in geologic materials. Comput. Methods Appl. Mech. Eng. 369, 113265. doi:10.1016/j.cma.2020.113265

CrossRef Full Text | Google Scholar

Fei, F., and Choo, J. (2021). Double-phase-field formulation for mixed-mode fracture in rocks. Comput. Methods Appl. Mech. Eng. 376, 113655. doi:10.1016/j.cma.2020.113655

CrossRef Full Text | Google Scholar

Fei, F., Choo, J., Liu, C., and White, J. A. (2022). Phase-field modeling of rock fractures with roughness. Int. J. Numer. Anal. Methods Geomechanics 46, 841–868. doi:10.1002/nag.3317

CrossRef Full Text | Google Scholar

Fei, F., Mia, M. S., Elbanna, A. E., and Choo, J. (2023). A phase-field model for quasi-dynamic nucleation, growth, and propagation of rate-and-state faults. Int. J. Numer. Anal. Methods Geomechanics 47, 187–211. doi:10.1002/nag.3465

CrossRef Full Text | Google Scholar

Fei, F., Sun, Y., Wong, L. N. Y., and Choo, J. (2021). “Phase-field modeling of mixed-mode fracture in rocks with discontinuities: from laboratory scale to field scale,” in 55th US rock mechanics/geomechanics symposium. ARMA–21–1223.

Google Scholar

Foster, C. D., Borja, R. I., and Regueiro, R. A. (2007). Embedded strong discontinuity finite elements for fractured geomaterials with variable friction. Int. J. Numer. Methods Eng. 72, 549–581. doi:10.1002/nme.2020

CrossRef Full Text | Google Scholar

Ha, S. J., Choo, J., and Yun, T. S. (2018). Liquid CO₂ fracturing: effect of fluid permeation on the breakdown pressure and cracking behavior. Rock Mech. Rock Eng. 51, 3407–3420. doi:10.1007/s00603-018-1542-x

CrossRef Full Text | Google Scholar

Ingraffea, A. R., and Heuze, F. E. (1980). Finite element models for rock fracture mechanics. Int. J. Numer. Anal. Methods Geomechanics 4, 25–43. doi:10.1002/nag.1610040103

CrossRef Full Text | Google Scholar

Lee, H., and Jeon, S. (2011). An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int. J. Solids Struct. 48, 979–999. doi:10.1016/j.ijsolstr.2010.12.001

CrossRef Full Text | Google Scholar

Lee, S., Wheeler, M. F., and Wick, T. (2016). Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model. Comput. Methods Appl. Mech. Eng. 305, 111–132. doi:10.1016/j.cma.2016.02.037

CrossRef Full Text | Google Scholar

Liu, F., and Borja, R. I. (2009). An extended finite element framework for slow-rate frictional faulting with bulk plasticity and variable friction. Int. J. Numer. Anal. Methods Geomechanics 33, 1535–1560. doi:10.1002/nag.777

CrossRef Full Text | Google Scholar

Lu, Y., Wang, L., and Elsworth, D. (2015). Uniaxial strength and failure in sandstone containing a pre-existing 3-D surface flaw. Int. J. Fract. 194, 59–79. doi:10.1007/s10704-015-0032-3

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Miehe, C., Welschinger, F., and Hofacker, M. (2010b). Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int. J. Numer. methods Eng. 83, 1273–1311. doi:10.1002/nme.2861

CrossRef Full Text | Google Scholar

Palmer, A. C., and Rice, J. R. (1973). The growth of slip surfaces in the progressive failure of over-consolidated clay. Proc. R. Soc. Lond. A. Math. Phys. Sci. 332, 527–548. doi:10.1098/rspa.1973.0040

CrossRef Full Text | Google Scholar

Pollard, D. D., and Fletcher, R. C. (2005). Fundamentals of structural geology. Cambridge University Press.

Google Scholar

Puzrin, A. M., and Germanovich, L. (2005). The growth of shear bands in the catastrophic failure of soils. Proc. R. Soc. A Math. Phys. Eng. Sci. 461, 1199–1228. doi:10.1098/rspa.2004.1378

CrossRef Full Text | Google Scholar

Regueiro, R. A., and Borja, R. I. (1999). A finite element model of localized deformation in frictional materials taking a strong discontinuity approach. Finite Elem. Analysis Des. 33, 283–315. doi:10.1016/s0168-874x(99)00050-5

CrossRef Full Text | Google Scholar

Sanborn, S. E., and Prévost, J. H. (2011). Frictional slip plane growth by localization detection and the extended finite element method (xfem). Int. J. Numer. Anal. Methods Geomechanics 35, 1278–1298. doi:10.1002/nag.958

CrossRef Full Text | Google Scholar

Schultz, R. A. (2019). Geologic fracture mechanics. Cambridge University Press.

Google Scholar

Shen, B., and Stephansson, O. (1994). Modification of the G-criterion for crack propagation subjected to compression. Eng. Fract. Mech. 47, 177–189. doi:10.1016/0013-7944(94)90219-4

CrossRef Full Text | Google Scholar

Sun, Y., Fei, F., Wong, L. N. Y., and Choo, J. (2024). Intermediate principal stress effects on the 3D cracking behavior of flawed rocks under true triaxial compression. Rock Mech. Rock Eng. 57, 4607–4634. doi:10.1007/s00603-024-03777-x

CrossRef Full Text | Google Scholar

Wong, L. N. Y., and Einstein, H. H. (2009). Crack coalescence in molded gypsum and carrara marble: part 2—microscopic observations and interpretation. Rock Mech. Rock Eng. 42, 513–545. doi:10.1007/s00603-008-0003-3

CrossRef Full Text | Google Scholar

Wong, N. Y. (2008). Crack coalescence in molded gypsum and Carrara marble. Ph.D. thesis, Massachusetts Institute of Technology.

Google Scholar

Wu, J.-Y. (2017). A unified phase-field theory for the mechanics of damage and quasi-brittle failure. J. Mech. Phys. Solids 103, 72–99. doi:10.1016/j.jmps.2017.03.015

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Zhang, X., Sloan, S. W., Vignes, C., and Sheng, D. (2017). A modification of the phase-field model for mixed mode crack propagation in rock-like materials. Comput. Methods Appl. Mech. Eng. 322, 123–136. doi:10.1016/j.cma.2017.04.028

CrossRef Full Text | Google Scholar

Zhou, S., Zhuang, X., Zhu, H., and Rabczuk, T. (2018). Phase field modelling of crack propagation, branching and coalescence in rocks. Theor. Appl. Fract. Mech. 96, 174–192. doi:10.1016/j.tafmec.2018.04.011

CrossRef Full Text | Google Scholar

Zhou, X.-P., Peng, S.-L., Zhang, J.-Z., and Berto, F. (2021). Failure characteristics of coarse and fine sandstone containing two parallel fissures subjected to true triaxial stresses. Theor. Appl. Fract. Mech. 112, 102932. doi:10.1016/j.tafmec.2021.102932

CrossRef Full Text | Google Scholar