^{1}Key Laboratory of In-situ Property Improving Mining of Ministry of Education, Taiyuan University of Technology, Taiyuan, China^{2}Shanxi Institute of Energy, Yuci, China

The paper presents a framework for accelerating the phase field modeling of compressive failure of rocks. In this study, the Drucker-Prager failure surface is taken into account in the phase field model to characterize the tension-compression asymmetry of fractures in rocks. The degradation function that decouples the phase-field and physical length scales is employed, in order to reduce the mesh density in large structures. To evaluate the proposed approach, four numerical examples are given. The results of the numerical experiments demonstrate the accuracy and efficiency of the proposed approach in tracking crack propagation paths in rock materials under Drucker-Prager criterion.

## 1 Introduction

The phase-field fracture model gains great popularity in computational fracture mechanics in recent years due to its capability of capturing complex fracture patterns including crack initiation, propagation, bifurcation and coalescence. Because the conventional form of phase field is based on the assumption of tension and compression symmetry, it is not applicable to rock materials [1], whose tensile and compressive strengths show significant differences [2]. To reproduce the fracture behaviors exhibiting asymmetric tension–compression characteristics, Zhou et al. [3]and Wang et al. [4] developed new driving force formulations, whereby Mohr–Coulomb criterion can be introduced to phase field fracture modeling. Navidtehrani et al [5–7] proposed a general framework for decomposing the strain energy density under multi-axial loading, which enables us to simulate compressive failure in rocks under Drucker-Prager criterion.

The phase-field fracture model has limitations for simulating large-scale rock models [1]. In the traditional formulation, the phase field length scale is linked with the physical process zone length scale for a given material strength [8–10]. In the analysis of the structures whose sizes are orders of magnitude larger than their physical length scales, the mesh density can be prohibitively high, leading to an unaffordable computational cost in practice. To address this issue, Wu et al. proposed a new type of degradation function which is insensitive to the length scale [11–20]. After that, Lo et al [18, 20] presented a degradation function that decouples the phase field length scale from the physical length scale, which reduces the mesh density and thus enables one to simulate crack propagation in large-scale rock masses with phase field methods [21–28].

In this paper, we combine the work of Navidtehrani et al. [5–7] and Lo et al [18, 20] to accelerate the fracture phase field modeling of Drucker–Prager failure by using the degradation functions decoupling the phase field and physical length scales. The remainder of the paper is organized as follows. Section 2 introduces the phase field fracture model for Drucker-Prager failure. Section 3 explains the degradation function that separates phase field length scales from physical length scales. The numerical experiments were conducted in Section 4, followed by the conclusions in Section 5. To demonstrate the accuracy of this method in capturing the crack patterns of rock materials.

## 2 Phase field fracture model

According to [29], the total potential energy of an elastic body is composed of the elastic energy of the elastic body and the crack surface energy

where:

Here we use the fracture variational criterion inherited and developed from the traditional Griffith theory, which is still based on the elastic strain energy and energy release rate. Griffith believes that there are many small cracks or defects in actual materials. Under the action of external forces, stress concentration will occur near these cracks and defects. When the stress reaches a certain level, the cracks will start to expand and cause fracture. However, Griffith theory [29] has the defect that it cannot solve the problems of crack generation, propagation angle and instability bifurcation, so the fracture of materials is further studied by using the fracture variational criterion.

### 2.1 Estimation of fracture surface energy using phase field variables

B. Bourdin et al. [30] realized the fracture variational criterion numerically for the first time by introducing phase field variables. In this paper, we define a scalar variable that changes in the interval of [0, 1] to be a phase field variable, and use

L0∈R+ is an important model parameter to control the range of crack diffusion fracture transition zone (0<

From formula (2), we can express the total crack surface energy in the elastic body with the following equation

### 2.2 Elastic strain energy decomposition

According to the above fracture variational criteria, the crack surface energy and elastic strain energy are tied closely. If the elastic strain energy is not decomposed, the pseudo bifurcation of the crack will occur. To solve this problem, this paper decomposes the elastic strain energy in tension and compression based on the method proposed by C. Miehe et al. [32], so that the tensile part of the elastic strain energy drives the evolution of the phase field. To this end, the strain tensor is first spectral decomposed [33]:

In the formula,

The strain after spectral decomposition can decompose the elastic strain energy density:

Here

Where:

### 2.3 Characterization of Drucker-Prager fracture surface

According to Navidtehrani [5–7] Drucker-Prager fracture criterion applicable to brittle or quasi-brittle materials, such as rock or concrete, is shown as follows

A is the equation about uniaxial tension (

For quasi-brittle materials, the mechanical properties within the Drucker – Prager failure range can still be regarded as linear elastic. Only when the stress reaches the failure surface, the linear elastic behavior will be transformed into non-linear action. At the same time, the failure of the material will lead to the weakening of the tensile strength and compressive strength. We can get the relationship between the material strength and the phase field damage by using the degradation function

The phase field parameter (𝜙 = 0) represents the integrity of the material. (𝜙 = 1) represents the complete destruction of the material. At this time, the Drucker – Prager parameters are as follows

Governing equations of phase field fracture model.

Under the condition of considering kinetic energy *T* at the same time, according to the crack surface energy expressed in Eq.3 and the elastic strain energy density expressed in Eq. 6, the following formula can be obtained:

In the formula: ν is the material density. Thus, the expression of Lagrange function can be written:

According to Hamilton’s principle and ignoring the physical force, the Lagrangian function L takes the variation of {

In the formula:

Where:

The state variable *H* represents the maximum tensile elastic strain energy. The maximum value of the *,*

Where:

### 2.4 Principle of virtual work. Balance of forces

Now, we will use the principle of virtual work to derive the equilibrium equations of the coupled deformation fracture system. Cauchy stress is introduced

Among δ Represents a virtual quantity. This equation must be applicable to any field Ω and any kinematically permissible change of virtual quantity. Therefore, by applying the Gauss divergence theorem, the local force balance is given by the following formula:

boundary conditions:

### 2.5 Constitutive theory

Based on the phase field crack modeling and theoretical framework, we first discuss the relationship between the change of phase field variables around the crack and the length scale by modifying the degradation function. It is feasible to modify the degradation function under our phase-field model to affect the peak stress of the phase-field model. The following is our equilibrium equation, free energy and related constitutive equation. The balance equation is as follows:

Phase field parameters usually introduced by phase field method

It can be seen from the second law of thermodynamics in isothermal form that for the system, the sum of the work done by mechanical traction and body force and the work done by external surface and body force is greater than or equal to Helmholtz free energy ψ Total change of. The integral form of the second law is

Assume that for any volume ψ depends on

According to the standard Coleman and Noll procedure [35], if,

The remaining items lead to unequal reduced dissipation,

The attenuation of this dissipative inequality always satisfies the following conditions

In the fracture phase field method of brittle materials similar to rock, the value of the

The following formula introduces a form of Helmholtz free energy to apply the frame to brittle fracture

Among

Rooted in the variational principle of linear elastic fracture energy proposed by Francfort and Marigo [30], and referring to the elliptic regularization method of Mumford-Shah functional in computer image segmentation, it is called AT2 model [36–42], and its manifestation in the phase-field fracture model is

According to our degradation function, it is easy to obtain the peak stress:

The final governing equation is as follows:

And

## 3 The degradation function

For the sake of completeness, we introduce the degradation function proposed by Lo et al [18, 20]. In this section. Eq. 30 shows that the peak stress is not only related to the properties of the material itself, such as material strength

As shown in Figure 2, the function curve under different s values. The characteristic of this degradation function is that a parameter s is introduced to control the specific peak stress through different values of s. We can easily calculate the peak stress of different s values by the following formula.

**FIGURE 2**. The curvature of the degradation function is controlled by the parameter s, corresponding to different change rates. For large s values, the change of degradation function is similar to that of traditional AT2 model.

we separate the phase field length scale

At this time, the ratio of peak stress is also different. We use the new peak stress symbol

For the value of s, we can find that when s ≈ 1.0148,

## 4 Numerical examples

### 4.1 Tension square plate with unilateral crack

Consider the tension square plate containing initial cracks as shown in Figure 3Awith a side length of 50 mm, which is divided into 46118 quadrilateral elements as shown in Figure 3B, and locally densified at the predicted crack growth path, where the grid size is about 0.1 mm. Let the problem be a plane strain problem, with an elastic modulus of 250 kN/mm^{2}, a Poisson’s ratio of 0.2, and a fracture toughness of Gc = 2.5 × 10−^{3} kN/mm^{2}. Using displacement loading method, t = 1s, ΔU = 0.15 mm until complete failure. During the calculation, l_{0} = 0.2 mm is taken. First, the geometric discontinuity method is used to preset the initial crack, that is, the upper and lower elements are geometrically separated at the crack. The calculated crack growth process is shown in Figure 3.

**FIGURE 3**. **(A)** Square plate under tension/m. **(B)** Meshing diagram (Tight zone dimensions = 0.2, other = 1) **(C)** Crack propagation path (S = 5). **(D)** Crack propagation path (S = 200). **(E)** Crack propagation path (traditional degradation).

Based on the numerical simulation results presented in Figures 3, 4, it can be observed that the model employing our proposed phase-field fracture method generates a crack propagation path that is more complex and closely resembles a rock fracture section during the simulation of the tensile square plate. Conversely, the crack growth path generated by the traditional degradation function appears more linear. These findings indicate that our proposed phase-field fracture model possesses distinctive characteristics in accurately capturing crack propagation under complex stress-displacement boundary conditions, which enables it to fit large-scale rock models while maintaining crack tracking precision. Consequently, further numerical simulations will be conducted to comprehensively investigate, analyze and compare the performance of our proposed model.

**FIGURE 4**. **(A)** Meshing diagram = 0.4. **(B)** Crack propagation path (S = 1.107). **(C)** Crack propagation path (S = 1.01).

### 4.2 Shear square plate with unilateral crack

The same material parameters as in Example 4.1 are used to study the crack growth process of square plates under shear. The calculation model is shown in Figure 5A, and the vertical displacement on all boundaries is constrained to zero. The displacement loading method is adopted, and the displacement increment of each loading step is Δ U = 2.5mm, the non-uniform grid is adopted, and the grid is densified at the position where the crack growth is expected. The minimum grid size is about 0.2 mm. The grid model is shown in Figure 5B.

**FIGURE 5**. **(A)** Sheared square plates. **(B)** Meshing diagram (Tight zone dimensions = 0.2, other = 1) **(C)** Crack propagation path (S = 200). **(D)** Crack propagation path (traditional degradation).

In the experiment involving a shear square plate with a unilateral crack, numerical simulations were conducted by varying the S value (S = 200), and the results obtained are shown in Figures 5D, E. Comparison of these results with those obtained by the traditional degradation function, such as Figure 5C, reveals that the tracked crack propagation path remains the same. This observation is consistent with our theoretical calculations of the relationship between S value and peak stress.

For the following research, we plan to investigate the effect of varying the S value on the physical and phase field scales. Specifically, we will choose values of S = 1.01 and 1.1 to expand the physical scale and the multiple of the phase field scale. The choice of S = 1.01 will result in a phase field scale that is 16 times the physical scale, while S = 1.107 will yield a phase field scale that is 4 times the physical scale. The model cell mesh size will be set to 0.4 mm. Through this verification process, we aim to demonstrate that our phase-field model has lower demands for mesh accuracy than the traditional phase-field model.

Based on the results presented in Figure 6, it can be observed that a gradual decrease of parameter S leads to an increase in the proportion of phase field length to physical length. This results in the phase field length being several times or even ten times greater than the physical length. As the proportional relationship between the mesh density and phase field length is a key factor affecting fracture tracking accuracy in the phase-field method, it is possible to reduce the mesh density without compromising the accuracy of crack tracking, as long as the proportion of mesh density and phase field length is maintained. This outcome is highly valuable for large-scale rock numerical simulations, as it results in significant computational resource savings and enhances the efficacy of the phase-field method for such applications.

**FIGURE 6**. **(A)** Meshing diagram (Tight zone dimensions = 0.4). **(B)** Crack propagation path (S = 1.107). **(C)** Crack propagation path (S = 1.01).

### 4.3 Direct shear test of cracked plates

Next, direct shear tests are simulated to evaluate the degradation function changes in the test configuration. The geometry and boundary conditions of the model are shown in Figure 7A. Apply a lateral displacement on the top edge, which is equal to 0.15 mm. The boundary conditions are vertical constraints at the top and fixed constraints at the bottom. Modulus of elasticity E = 25 GPa and Poisson’s ratio *ν* = 0.2, fracture parameters are given by Gc = 0.15 kJ/m^{2} and l = 0.2 mm. In order to save computing resources, the predicted crack propagation area is encrypted, and the cell size of the encrypted part is at least half of the phase field scale, so as to save computer resources without affecting the computing accuracy.

**FIGURE 7**. **(A)** Direct shear test (DST). **(B)** Meshing diagram (Tight zone dimensions = 0.1, other = 1) and the calculated crack growth process. **(C)** Crack propagation path (S = 200). **(D)** Crack propagation path (traditional degradation).

The crack initiation and propagation mode observed in the simulation is in line with the rock shear fracture experiment, where the crack initiates from the edge and progresses towards the center. The simulation results presented in the figure above compare the crack growth path obtained using the specific degradation function with S = 200 and the traditional degradation function. The result display that the crack propagation paths traced by the traditional degradation function are similar when S = 200.

As the value of S approaches 1, the phase field characteristic length and the physical characteristic length become decoupled, resulting in a change in the quantitative relationship between the phase field length and physical scale. This leads to a significant increase in the phase field scale, which can become several times or even dozens of times larger than the physical scale. For instance, when S = 1.107, the phase field scale is observed to be four times larger than the physical scale. To accurately track the crack propagation path using the phase field method, a smaller mesh size than the phase field scale is required. Therefore, with the increase in phase field scale, excessively dense grids can be discarded. By maintaining the proportion between the phase field scale and mesh size, a mesh size of 0.4 mm can be set to achieve equally accurate crack propagation path tracking.

As shown in Figure 8, when using a smaller S value, increase the cell size of the model and use a 0.04 mm mesh size. At this time, very accurate crack tracking curve of rock direct shear test can still be obtained. It is verified that this phase field model can more effectively simulate large-scale rock fracture.

**FIGURE 8**. **(A)** Meshing diagram (Tight zone dimensions = 0.5, other = 2). **(B)** Crack propagation path (S = 1.107).

### 4.4 Notched plate with eccentric holes

In this case, we model the rock plate with holes, as shown in the figure. This is an eccentric plate with three holes with a length of 120 mm and a width of 70 mm. The size and distribution of the holes are indicated in the Figure 9 below. The boundary condition is that the bottom is fixed vertically and horizontally, and 0.15 mm vertical displacement is applied to the top. The material property parameter in this simulation is Young’s modulus E = 25 GPa, Poisson’s ratio *ν* = 0.20, Phase field scale ^{2}. At the same time, we compare the different degradation functions affected by multiple S and the cracks tracked by the traditional degradation functions, and use the same division ratio for the element, that is, the mesh size is equal to 0.1 mm. See the figure below for the specific numerical simulation results.

**FIGURE 9**. **(A)** Boundary conditions and geometric parameters **(B)** actual experimental results **(C)** Meshing diagram (Tight zone dimensions=0.1, other=1) **(D)** Crack propogation path (S=5) **(E)** Crack propogation path (S=200) **(F)** Crack propogation path (traditional degradation).

In notched rock slabs with eccentric holes, obtaining accurate crack information is crucial, as multiple factors affect the crack propagation process. In this simulation, we employed the traditional phase field model of Drucker-Prager fracture criterion and the traditional degradation function for numerical simulation, and the simulation results are depicted in Figures 9, 10. Comparing Figures 10 B, C, if we exaggerate the phase field length while maintaining the original set mesh density, our crack path will shift. At this time, we can alleviate this problem by keeping the proportion of the enlarged grid density and the phase field length the same, that is, the grid density is less than half of the phase field scale. Even so, the mesh density after densification is still much smaller than the initially given mesh density, which is more conducive to obtaining the correct crack path and more efficient numerical simulation. However, we can overcome this issue by densifying the mesh. Through multiple sets of experimental simulations, we determined that a phase field scale with a grid density less than half can accurately track the crack propagation path.

**FIGURE 10**. **(A)** Crack propagation path S = 1.107 (element size 0.4). **(B)** Crack propagation path S = 1.107 (element size 0.4). **(C)** Crack propagation path S = 1.107 (element size 1.0).

## 5 Conclusion

This study effectively accelerates the fracture phase field modeling for compressive failure of rocks materials. By splitting the strain energy using the approach introduced by [18, 20], the asymmetric tension-compression behaviors in rock fracture processes can be characterized. To decouple the phase-field length with the physical length scale, the degradation function proposed by [18, 20] were adopted and thus the mesh density is reduced for the structures much larger than the physical process zone. In the future, we will extend the present approach to the applications of hydraulic fracturing in quasi-brittle materials [43]. In addition, the combination of isogeometric analysis will be considered to solve high-order problems [44].

## 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

BL: optical experiment, ZL and LY: guidance and revision.

## Funding

The authors are grateful for the support of National Natural Science Foundation of China (Nos. 52174213, 51574174).

## Conflict of interest

The authors declare 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

2. León Baldelli AA, Maurini C. Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. *J Mech Phys Sol* (2021) 152:104424. doi:10.1016/j.jmps.2021.104424

3. Zhuang XY, Zhou SW, Sheng M, Li GS. On the hydraulic fracturing in naturally-layered porous media using the phase field method. *ENGINEERING GEOLOGY* (2020) 266:105306. doi:10.1016/j.enggeo.2019.105306

4. Zhang Q, Wang D, Zeng F, Guo Z, Wei N. Pressure transient behaviors of vertical fractured wells with asymmetric fracture patterns. *JOURNAL ENERGY RESOURCES TECHNOLOGY-TRANSACTIONS ASME* (2020) 142(4). doi:10.1115/1.4045226

5. Navidtehrani Y, Betegon C, Martinez-Paneda E. A unified abaqus implementation of the phase field fracture method using only a user material subroutine. *MATERIALS* (2021) 14(8):1913. doi:10.3390/ma14081913

6. Navidtehrani Y, Betegon C, Martinez-Paneda E. A general framework for decomposing the phase field fracture driving force, particularised to a Drucker-Prager failure surface. *THEORETICAL APPLIED FRACTURE MECHANICS* (2022) 121. doi:10.1016/j.tafmec.2022.103555

7. Navidtehrani Y, Betegon C, Zimmerman RW, Martinez-Paneda E. Griffith-based analysis of crack initiation location in a Brazilian test. *INTERNATIONAL JOURNAL ROCK MECHANICS MINING SCIENCES* (2022) 159:105227. doi:10.1016/j.ijrmms.2022.105227

8. Cundall PA. A computer model for simulating progressive large-scale movements in blocky rock systems. In: Proceedings of the Symposium of the International Society for Rock Mechanics, Society for Rock Mechanics (ISRM); 1971; France, Paris (1971). p. 11–8.

9. Li H, Huang Y, Yang Z, Yu K, Li QM. 3D meso-scale fracture modelling of concrete with random aggregates using a phase-field regularized cohesive zone model. *Int J Sol Structures* (2022) 256:111960. doi:10.1016/j.ijsolstr.2022.111960

10. Carlsson J, Isaksson P. A statistical geometry approach to length scales in phase field modelling of fracture and strength of porous microstructures. *Int J Sol Structures* (2020) 200-201:83–93. doi:10.1016/j.ijsolstr.2020.05.003

11. Borden MJ, Hughes TJR, Landis CM, Anvari A, Lee IJ. Corrigendum to “A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. *Comp Methods Appl Mech Eng* (2017) 324:712–3. doi:10.1016/j.cma.2017.06.023

12. Jiang D, Kyriakides S, Landis CM. Propagation of phase transformation fronts in pseudoelastic NiTi tubes under uniaxial tension. *Extreme Mech Lett* (2017) 15:113–21. doi:10.1016/j.eml.2017.06.006

13. Jiang D, Kyriakides S, Landis CM, Kazinakis K. Modeling of propagation of phase transformation fronts in NiTi under uniaxial tension. *Eur J Mech - A/Solids* (2017) 64:131–42. doi:10.1016/j.euromechsol.2017.02.004

14. Landis CM, Beyerlein IJ, McMeeking RM. Micromechanical simulation of the failure of fiber reinforced composites. *J Mech Phys Sol* (2000) 48(3):621–48. doi:10.1016/S0022-5096(99)00051-4

15. Landis CM, Huang R, Hutchinson JW. Formation of surface wrinkles and creases in constrained dielectric elastomers subject to electromechanical loading. *J Mech Phys Sol* (2022) 167:105023. doi:10.1016/j.jmps.2022.105023

16. Li W, Landis CM. Nucleation and growth of domains near crack tips in single crystal ferroelectrics. *Eng Fracture Mech* (2011) 78(7):1505–13. doi:10.1016/j.engfracmech.2011.01.002

17. Li W, McMeeking RM, Landis CM. On the crack face boundary conditions in electromechanical fracture and an experimental protocol for determining energy release rates. *Eur J Mech - A/Solids* (2008) 27(3):285–301. doi:10.1016/j.euromechsol.2007.08.007

18. Lo Y-S, Borden MJ, Ravi-Chandar K, Landis CM. A phase-field model for fatigue crack growth. *J Mech Phys Sol* (2019) 132:103684. doi:10.1016/j.jmps.2019.103684

19. Woldman AY, Landis CM. Thermo-electro-mechanical phase-field modeling of paraelectric to ferroelectric transitions. *Int J Sol Structures* (2019) 178-179:19–35. doi:10.1016/j.ijsolstr.2019.06.012

20. Lo Y-S, Hughes TJR, Landis CM. Phase-field fracture modeling for large structures. *J Mech Phys Sol* (2023) 171:105118. doi:10.1016/j.jmps.2022.105118

21. Ai W, Wu B, Martínez-Pañeda E. A coupled phase field formulation for modelling fatigue cracking in lithium-ion battery electrode particles. *J Power Sourc* (2022) 544:231805. doi:10.1016/j.jpowsour.2022.231805

22. Cui C, Ma R, Martínez-Pañeda E. A generalised, multi-phase-field theory for dissolution-driven stress corrosion cracking and hydrogen embrittlement. *J Mech Phys Sol* (2022) 166:104951. doi:10.1016/j.jmps.2022.104951

23. Isfandbod M, Martínez-Pañeda E. A mechanism-based multi-trap phase field model for hydrogen assisted fracture. *Int J Plasticity* (2021) 144:103044. doi:10.1016/j.ijplas.2021.103044

24. Kristensen PK, Niordson CF, Martínez-Pañeda E. A phase field model for elastic-gradient-plastic solids undergoing hydrogen embrittlement. *J Mech Phys Sol* (2020) 143:104093. doi:10.1016/j.jmps.2020.104093

25. Martínez-Pañeda E, Fuentes-Alonso S, Betegón C. Gradient-enhanced statistical analysis of cleavage fracture. *Eur J Mech - A/Solids* (2019) 77:103785. doi:10.1016/j.euromechsol.2019.05.002

26. Navidtehrani Y, Betegón C, Martínez-Pañeda E. A simple and robust Abaqus implementation of the phase field fracture method. *Appl Eng Sci* (2021) 6:100050. doi:10.1016/j.apples.2021.100050

27. Navidtehrani Y, Betegón C, Martínez-Pañeda E. A general framework for decomposing the phase field fracture driving force, particularised to a Drucker–Prager failure surface. *Theor Appl Fracture Mech* (2022) 121:103555. doi:10.1016/j.tafmec.2022.103555

28. Tan W, Martínez-Pañeda E. Phase field fracture predictions of microscopic bridging behaviour of composite materials. *Compos Structures* (2022) 286:115242. doi:10.1016/j.compstruct.2022.115242

29. Espadas-Escalante JJ, Van Dijk NP, Isaksson Pjcs T. A phase-field model for strength and fracture analyses of fiber-reinforced composites. *Compos Sci Technol* (2019) 174(APR.12):58–67. doi:10.1016/j.compscitech.2018.10.031

30. Bourdin B, Francfort GA, Marigo JJ. Numerical experiments in revisited brittle fracture. *J Mech Phys Sol* (2000) 48(4):797–826. doi:10.1016/S0022-5096(99)00028-9

31. Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. *Int J Numer Methods Eng* (2010) 83(10):1273–311. doi:10.1002/nme.2861

32. Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. *Comp Methods Appl Mech Eng* (2010) 199(45):2765–78. doi:10.1016/j.cma.2010.04.011

33. Bernard PE, Moës N, Chevaugeon N. Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings. *Comp Methods Appl Mech Eng* (2012) 233-236:11–27. doi:10.1016/j.cma.2012.02.020

34. Gurtin ME. Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance. *Physica D: Nonlinear Phenomena* (1996) 92(3):178–92. doi:10.1016/0167-2789(95)00173-5

35. Coleman BD, Noll W. The thermodynamics of elastic materials with heat conduction and viscosity. In: W Noll, editor. *The foundations of mechanics and thermodynamics: Selected papers*. Berlin, Heidelberg: Springer Berlin Heidelberg (1974). p. 145–56.

36. Feng D-C, Wu J-Y. Phase-field regularized cohesive zone model (CZM) and size effect of concrete. *ENGINEERING FRACTURE MECHANICS* (2018) 197:66–79. doi:10.1016/j.engfracmech.2018.04.038

37. Vinh Phu N, Wu J-Y. Modeling dynamic fracture of solids with a phase-field regularized cohesive zone model. *COMPUTER METHODS APPLIED MECHANICS ENGINEERING* (2018) 340:1000–22. doi:10.1016/j.cma.2018.06.015

38. Wu J-Y. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. *JOURNAL MECHANICS PHYSICS SOLIDS* (2017) 103:72–99. doi:10.1016/j.jmps.2017.03.015

39. Wu J-Y. A geometrically regularized gradient-damage model with energetic equivalence. *COMPUTER METHODS APPLIED MECHANICS ENGINEERING* (2018) 328:612–37. doi:10.1016/j.cma.2017.09.027

40. Wu J-Y. Robust numerical implementation of non-standard phase-field damage models for failure in solids. *COMPUTER METHODS APPLIED MECHANICS ENGINEERING* (2018) 340:767–97. doi:10.1016/j.cma.2018.06.007

41. Wu J-Y, Vinh Phu N. A length scale insensitive phase-field damage model for brittle fracture. *JOURNAL MECHANICS PHYSICS SOLIDS* (2018) 119:20–42. doi:10.1016/j.jmps.2018.06.006

42. Wu J-Y, Vinh Phu N, Nguyen CT, Sutula D, Sinaie S, Bordas SPA. Phase-field modeling of fracture. In: SPA Bordas, and DS Balint, editors, 53 (2020). p. 1–183.*ADVANCES APPLIED MECHANICS*

43. Yang J, Lian H, Nguyen VP. Study of mixed mode I/II cohesive zone models of different rank coals. *Eng Fracture Mech* (2021) 246:107611. doi:10.1016/j.engfracmech.2021.107611

Keywords: phase field method, fracture, Drucker-Pluger, degradation function, large scale

Citation: Liu B, Liu Z and Yang L (2023) Accelerating fracture simulation with phase field methods based on Drucker-Prager criterion. *Front. Phys.* 11:1159566. doi: 10.3389/fphy.2023.1159566

Received: 06 February 2023; Accepted: 08 March 2023;

Published: 24 March 2023.

Edited by:

Leilei Chen, Huanghuai University, ChinaReviewed by:

Qingwei Bu, Anhui University of Science and Technology, ChinaZhijun Liu, Heilongjiang University of Science and Technology, China

Copyright © 2023 Liu, Liu and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhenghe Liu, liuzhenghe@tyut.edu.cn