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, Γ ( d , ∇ 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.

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. ∇ ⋅ ∂ ψ ϵ , d I , ∇ d I , d I I , ∇ d I I ∂ ∇ d I − ∂ ψ ϵ , d I , ∇ d I , d I I , ∇ d I I ∂ d I = 0 , (1) ∇ ⋅ ∂ ψ ϵ , d I , ∇ d I , d I I , ∇ d I I ∂ ∇ d I I − ∂ ψ ϵ , d I , ∇ d I , d I I , ∇ d I I ∂ d I I = 0 , (2) where d I and d I I are the phase fields representing tensile and shear fractures, respectively, ϵ is the (infinitesimal) strain tensor, and ψ is the potential energy density function. The potential energy density function is comprised of three parts as ψ = ψ e + ψ f + ψ d , (3) where ψ e is the elastic strain energy, ψ f is the energy dissipation due to frictional slip, and ψ 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 ′ d I H I − G I π L 2 L 2 ∇ ⋅ ∇ d I − 2 + 2 d I = 0 , (4) − g I I ′ d I I H I I − G I I π L 2 L 2 ∇ ⋅ ∇ d I I − 2 + 2 d I I = 0 . (5) Here, g I ′ ( d I ) and g I I ′ ( d I I ) are the degradation functions, G I and G I I are the fracture energies, and H I and H I I are the crack driving forces for tensile and shear fractures, respectively. Considering that G I and G I I of geomaterials can be significantly different, the crack driving forces are calculated based on the F -criterion (Shen and Stephansson, 1994) as θ = arg max θ F θ | ϵ , where F θ ≔ H I ϵ , θ G I + H I I ϵ , θ G I I , (6) where θ is the angle between the crack normal direction and the major principal stress direction on the slip plane. It is noted that the F criterion lets the mixed-mode fracture develop such that it minimizes energy dissipation. The specific expressions for H I ( ϵ , θ ) and H I I ( ϵ , θ ) 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 ′ ( d I ) and g I I ′ ( d I I ) , 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.

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.

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 G I = 98 J / m 2 and G I I = 1500 J / 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 ≥ 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).

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.

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, σ 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 σ 2 with respect to the flaw and the magnitude of σ 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 σ 2 are shown in Figure 4. When σ 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 σ 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 σ 1 is parallel to the flaw plane, the two flaws did not coalesce. These differences indicate that the orientation of σ 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 σ 2 on the 3D mixed-mode cracking behavior of rocks—see Sun et al. (2024) for details.

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.