Skip to main content

BRIEF RESEARCH REPORT article

Front. Phys., 24 March 2023
Sec. Statistical and Computational Physics
This article is part of the Research Topic Moving Boundary Problems in Multi-physics Coupling Processes View all 16 articles

Accelerating fracture simulation with phase field methods based on Drucker-Prager criterion

Bin LiuBin Liu1Zhenghe Liu
Zhenghe Liu1*Lusheng YangLusheng Yang2
  • 1Key Laboratory of In-situ Property Improving Mining of Ministry of Education, Taiyuan University of Technology, Taiyuan, China
  • 2Shanxi 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 [57] 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 [810]. 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 [1120]. 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 [2128].

In this paper, we combine the work of Navidtehrani et al. [57] 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

Πu,Γ=ΩψeεudΩ+ΓGcdΓ(1)

where: ψe is the elastic energy density of the elastomer, εu Is the strain tensor, u is the displacement, the range of u is ∈ Rd ({1,2,3}), and Gc is the fracture energy release rate of the material. We can obtain from the variational method that at time t ∈ [0, t], crack I (x, t) has any behavior of x ∈ Rd at any position, Once the crack is formed, it cannot be recovered, so an irreversible condition needs to be imposed. For energy, minimizing the total potential energy will only increase but not decrease. This condition is Γ Γ (x, s) (x, t), (s < t).

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 ϕ to represent the topology of the crack, when ϕ = 1 to represent the crack, and when ϕ = 0 to represent the material is intact, and then use ϕ to represent the crack surface density in unit volume γ [31]:

γϕ,ϕ=ϕ22l0+l02|ϕ|2(2)

L0∈R+ is an important model parameter to control the range of crack diffusion fracture transition zone (0< ϕ <1), as shown in Figure 1. l0 characterizes the crack diffusion range. In short, when l0 is small, the crack is thinner, while the larger l0 is, the fatter the crack is. In this way, the crack surface Γ acts as a time-varying scalar field, or phase field, in space.

FIGURE 1
www.frontiersin.org

FIGURE 1. Fissures described with a phase-field model.

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

ΓGcdΓΩGcϕ122l0+l02|ϕ|2dΩ(3)

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]:

ε±=i=1dεi±nini(4)

In the formula, ε+ and ε are tensile strain tensors and compression strain tensors respectively. εi and ni are the main strain values and their corresponding directions. d is an independent parameter, namely the spatial dimension. When d = 2, this is a two-dimensional problem. When d = 3, this is a three-dimensional problem. Macaulay brackets in the formula are defined as: < • >+ = (• + | • |)/2,< • >−= (•− | •|)/2.

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

ψe±ε=λ2trε±2+μtrε±2(5)

Here λ and μ is the Lame constant, tr (•) represents the trace of the matrix. In the process of material failure, the stiffness will also decrease. At this time, the material stiffness weakening is related to the phase field variable. If we assume that only the tensile strain energy density in the elastic body receives the weakening of the phase field variable, the elastic strain energy at this time is

ψeε=1ϕ2+κψe+ε+ψeε(6)

Where: κ is a parameter model parameter, 0 ≤ κ<< 1. Avoid generating numerical singularity when ϕ = 1.

2.3 Characterization of Drucker-Prager fracture surface

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

J2σ=A+BI1σ(7)

I1 and J2 are the first tensor and the second partial derivative of the invariant, respectively.

A is the equation about uniaxial tension (σt). B is the equation about uniaxial compression (σc). The specific form is as follows

A=23σcσtσc+σt;B=13σtσcσc+σt(8)

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 gϕ in the numerical calculation:

Aϕ=23gϕσcgϕσtgϕσc+gϕσt=gϕ23σcσtσc+σt=gϕAϕ=0Bϕ=13gϕσtgϕσcgϕσc+gϕσt=13σtσcσc+σt=Bϕ=0(9)

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

Aϕ=1=0;Bϕ=1=Bϕ=0(10)

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:

Tu˙=Ω12ρu˙u˙dV(11)

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

L=TΠ=Ω12ρu˙u˙dVΩ1ϕ2+κψe+ε+ψeεdVΩGcϕ22l0+l02|ϕ|2dV(12)

According to Hamilton’s principle and ignoring the physical force, the Lagrangian function L takes the variation of {u˙, ϕ}, and the control equation of the phase field fracture model can be obtained:

Divσ=ρu¨Gcl0+2ψe+ϕGcl0Δϕ=2ψe+(13)

In the formula: σ is Cauchy stress tensor. We can obtain the stress from the partial derivative of the elastic strain energy corresponding to the strain

σ=1ϕ2+κεψe+ε+εψeε=1ϕ2+κλtrε+I+2με++λtrεI+2με(14)

Where: I is the second order unit tensor. The control Eq. 13 is a set of partial differential equations composed of the dynamic balance equation and the phase field evolution equation. In order to make the control equation equal to the aforementioned variational method, the irreversible condition should be added, that is, when the elastomer is under pressure or unloading, the crack healing should be prevented. The simplest and most effective way to deal with the irreversibility of variables is to introduce historical state variables:

Hx,t=maxs0,tψe+εx,s(15)

The state variable H represents the maximum tensile elastic strain energy. The maximum value of the H variable from loading to the current time is also related to the relative position and time, that is, Hx,t. The irreversibility condition can be satisfied by replacing the value of ∈ ψe+ in Eq. 13 with Hx,t in the control equation. It is also because of the previous irreversible and only increase monotonically. Because of the irreversibility of the state variable H, ϕ also has the same irreversibility and will only increase monotonically. At this time, the control equation should be rewritten as:

Divσ+b=ρu¨Gcl0+2HϕGcl0Δϕ=2H(16)

Where: b is the body force. The phase field fracture model can describe any behavior of the fracture only by a set of partial differential equations. The separated fracture is equivalent to the fracture field coupled with the displacement field, and no additional tracking of the geometric shape of the tracking crack is required. Therefore, any propagation path of the fracture can be calculated.

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 σ, It is related to the strain ε Work conjugation. In addition, the traction force T is defined as partial Ω on the solid boundary. We introduce scalar stress into the fracture model ω. This stress ω And phase field ϕ and phase field micro-stress vector ξ, which conjugate with the phase field gradient ϕ. It is assumed that the phase field ϕ is only driven by the solution of the displacement problem. Therefore, there is no external traction and ϕ relevant. In the absence of external force, the principle of virtual work is as follows:

Ω{σ:δε+ωδϕ+ξδϕ}dV=ΩTδudS(17)

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:

σ=0ξω=0in(18)

boundary conditions:

σn=Tξn=0on(19)

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:

σji,j+bi=0inVandti=σjinjonS(20)

V is the volume of the object, S is its boundary surface, ni is the unit component (perpendicular to the surface), bi is the physical strength per unit volume, ti is the traction vector, σji is the component of Cauchy stress tensor. Assuming small deformation and deformation gradient, make infinitesimal strain tensor εij can be determined by the deformation gradient of the displacement vector ui is expressed as,

εij=12ui,j+uj,i(21)

Phase field parameters usually introduced by phase field method ϕ .It is used to describe the degradation of materials, and is generally set ϕ When the boundary changes from 1 to 0, 1 indicates that the material is intact without cracks ϕ = 0 indicates that the material has been completely destroyed. According to Gurtin’s [34] research, the micro-forces cited at the same time are: a group of micro-forces on the outer surface λ and body micro-force γ, the micro force on the inner surface is π and ξi, The following is the point direction balance equation of these micro forces:

ξi,i+γ+π=0inVandξini=λonS(22)

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

Vψ˙dVVb1u˙i+γϕ˙dV+Stiu˙i+λϕ˙dS(23)

Assume that for any volume ψ depends on ε˙ij, ϕ, and ϕ˙. Under the condition of applying divergence theory, we can determine that Eq. 23 is bound to hold.

ψεijε˙ij+ψϕϕ˙+ψϕ,iϕ˙,i+ψϕ˙ϕ¨σjiε˙ijπϕ˙+ξiϕ˙,i(24)

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

σji=ψεij,ξi=ψϕ,i,andψϕ=0(25)

The remaining items lead to unequal reduced dissipation,

π+ψϕϕ˙0.(26)

The attenuation of this dissipative inequality always satisfies the following conditions

π=ψϕβϕϕ˙,ϕ,i,εijϕ˙(27)

In the fracture phase field method of brittle materials similar to rock, the value of the β is zero.

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

ψ=gϕψe++ψe+Gc38l01ϕ+38l0ϕ,iϕ,i(28)

Among ψe+ represents the effect of tensile stress on elasticity, ψe represents the effect of compression on elastic strain energy. The study of strain energy decomposition by Borden et al is very profound. For the numerical simulation of brittle fracture, we usually introduce the physical quantity representing the tensile or compressive strength of the material within the elastic limit, namely the elastic modulus E, into the linear elastic isotropic material, as well as the elastic constant Poisson’s ratio ν reflecting the transverse deformation of the material. For our phase-field simulation, we also need to introduce physical quantities such as the effective fracture surface energy Gc per unit area, and the length scale parameter l0. For phase-field fracture simulation, we have mentioned that the variation range of phase-field parameters is between 0 and 1, and the length scale affects the influence range of phase-field parameters.

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 [3642], and its manifestation in the phase-field fracture model is 1ϕ2. With the failure of the material, the physical properties of the material itself will also change. At this time, we need to use the degradation function, namely g, to determine the weakening of the physical properties of the material. Obviously, our degradation function is related to the phase field parameters. Just as the phase field parameters ϕ only change from 0 to 1, so is our degradation function g. When the material is completely destroyed, the degradation function gϕ = 0, we can easily realize the weakening of the elastic property by multiplying it by the elastic energy. The following is the function of the degradation function and the elastic energy equation:

gϕ=1ϕ2(29)

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

σc=3GcE8l0(30)

The final governing equation is as follows:

ϕ2ψe+εij+ψeεij+bi=βϕ˙inVϕ2ψe+εij+ψeεijnj=tionS(31)

And

38Gcl0ϕ,ii+3Gc8l02ϕψe++γ=0inV,38Gcl0ϕ,ini=λonS(32)

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 σc. Breaking energy Gc is related to Young’s modulus E, and also to the phase field length scale l0. Because of our need for numerical simulation of large-scale rock phase field fracture, the phase field length scale is usually a very small physical quantity. This has brought a great burden to our actual calculation. We can easily associate whether we can decouple the peak stress and the phase field length scale if we adopt a new degradation function to ensure more efficient calculation results. Therefore, we adopt the following degradation function:

gϕ=s1s1sϕ2(33)

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.

σc*=3GcE8l0s1lnss1(34)

FIGURE 2
www.frontiersin.org

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 l0 and the physical length scale lp. In the traditional phase field simulation, l0=lp, and we choose different s values, the corresponding phase field length scale l0 and the physical length scale lp ratio are also different, and the corresponding physical length scale lp formula is

lp=s1lnss1l0(35)

At this time, the ratio of peak stress is also different. We use the new peak stress symbol σc* to distinguish,

σc*=3GcE8lp.(36)

For the value of s, we can find that when s ≈ 1.0148, l0 increase by 16 times with the peak stress unchanged. This has increased by an order of magnitude, and if we want to increase l0 100 times while maintaining the same peak stress, when s ≈ 1.00155. With different values of s, we can easily simulate a larger scale model, which has many advantages for the phase field modeling of our large-scale rock mass model.

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/mm2, a Poisson’s ratio of 0.2, and a fracture toughness of Gc = 2.5 × 10−3 kN/mm2. Using displacement loading method, t = 1s, ΔU = 0.15 mm until complete failure. During the calculation, l0 = 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
www.frontiersin.org

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
www.frontiersin.org

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
www.frontiersin.org

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
www.frontiersin.org

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/m2 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
www.frontiersin.org

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
www.frontiersin.org

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 l0 = 0.25 mm and energy release rate Gc = 0.15 kJ/m2. 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
www.frontiersin.org

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
www.frontiersin.org

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

1. Jaeger JC, Cook NJ. Fundamental of rock mechanics (1979). doi:10.1016/j.mineng.2009.02.015

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

44. Chen L, Wang Z, Peng X, Yang J, Wu P, Lian H. Modeling pressurized fracture propagation with the isogeometric BEM. Geomechanics Geophys Geo-Energy Geo-Resources (2021) 7(3):51. doi:10.1007/s40948-021-00248-3

CrossRef Full Text | Google Scholar

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, China

Reviewed by:

Qingwei Bu, Anhui University of Science and Technology, China
Zhijun 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, bGl1emhlbmdoZUB0eXV0LmVkdS5jbg==

Disclaimer: 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.