^{1}

^{1}

^{2}

^{1}

^{2}

This article was submitted to Statistical and Computational Physics, a section of the journal Frontiers in Physics

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.

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.

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 [

The phase-field fracture model has limitations for simulating large-scale rock models [

In this paper, we combine the work of Navidtehrani et al. [

According to [

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 [

B. Bourdin et al. [

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

Fissures described with a phase-field model.

From

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. [

In the formula,

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

Here

Where:

According to Navidtehrani [

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

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

Where:

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:

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 [

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 [

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

The final governing equation is as follows:

And

For the sake of completeness, we introduce the degradation function proposed by Lo et al [

As shown in

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,

Consider the tension square plate containing initial cracks as shown in ^{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

Based on the numerical simulation results presented in

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

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

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

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 ^{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.

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

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 ^{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.

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

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 [

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

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

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

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.

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.