Phase Field Modeling of Anisotropic Tension Failure of Rock-Like Materials

The tensile fracture is a widespread feature in rock excavation engineering, such as spalling around an opened tunnel. The phase field method (PFD) is a non-local theory to effectively simulate the quasi-brittle fracture of materials, especially for the propagation of a tensile crack. This work is dedicated to study the tensile failure characteristics of rock-like materials by the PFD simulation of the Brazilian test of the intact and fissure disk samples. The numerical results indicate that the tensile strength of the disk sample is anisotropic due to the influence of pre-existing cracks. The peak load decreases at first and then increases with the increase of the inclination angle, following the U-shaped trend. The simulation results also indicate that the wing crack growth is the main failure characteristic. Moreover, the crack propagation path initiates at the tip of the pre-existing crack when the inclination angle is less than 60°. Crack propagation initiates near the tip of the pre-existing crack when the angle is 75°, and it initiates at the middle of the pre-existing crack when the angle is 90°. Finally, all cracks extend to the loading position and approximately parallel to the loading direction. This process is in agreement with the Brazilian test of pre-existing cracks in the laboratory, which can validate the effectiveness of the PFD in simulating the tensile fracture of rock-like materials. This study can provide a reference for the fracture mechanism of the surrounding rock in the underground excavation.


INTRODUCTION
Due to the geological tectonic movement or artificial disturbance, the fracture of engineering rock mass exhibits a rich variety of crack patterns under loading or unloading conditions [1][2][3]. Tensile and compression-shear fractures often occur in practical rock engineering. Because rock-like materials have the unilateral effect (the tensile strength is much lower than the compression strength), the tensile fracture is more likely to occur in practical rock engineering. Moreover, the tensile crack is also found in rock mass under compressive condition due to the effect of its heterogeneity or pre-existing cracks [4,5]. For example, rock block spalling appears near a cave wall after rock excavation [6]. To avoid the occurrence of engineering disasters, corresponding solution should be developed by numerical simulation in advance. It is therefore very important to develop a numerical method to accurately capture the crack propagation of the surrounding rock.
Rock fracture failure behaviors can be understood by means of on-site monitoring, laboratory test, theoretical analysis, and numerical simulation. Up to now, the mechanics of rock materials can be summarized as strain hardening/softening [7], brittle-ductile transition [8], time-dependent [9], unilateral effect [10], anisotropic behaviors [11], etc. The fracture mechanism of rock has been understood qualitatively and quantitatively by experimental and theoretical methods. The analysis results are limited to a simple configuration or simple stress paths. However, the failure pattern and crack propagation are also affected by the geometry of the engineering problem. The failure characteristics will become more complex. For example, the localization damage and discontinuous deformation failure are often found in rock engineering [12]. The description of these non-linear behaviors needs to improve the existing models, which result in the difficulty in the theoretical analysis and the complexity of model formulation. Moreover, there is a difference between the analytical results and test data due to a simple assumption of uniform stress around cracks [13]. Therefore, it is difficult to satisfy the practical demand only from the aspects of experimental or theoretical analysis. As with the development of numerical methods, the fracture problem of rock-like materials may be described by combining a non-local numerical method (such as, phase field (PF)) and a simple constitutive model.
The phase field method (PFD) is a non-local theory and can capture the crack initiation, propagation, and coalescence of quasi-brittle materials [14,15], especially for the tension brittle fracture. The discontinuous deformation near cracks can be represented by a phase field variable. Propagation and branching of cracks can be reflected directly in the phase field modeling. The hybrid PFD is further proposed by decomposing the contributions of the driving energy of phase field evolution into tensile and compressive parts [16][17][18]. So far, PFD has been successfully applied in the simulation of quasi-brittle fracture of materials [15, 19, and 20] and multi-field coupling problems [21,22]. Compared with the direct approaches for modeling crack propagation, the PFD has the following advantages: without an enrichment function, crack tracking algorithm, and ad hoc criterion for crack initiation.
The purpose of this work was to further promote the application of the PFD in rock mechanics. First, a brief literature overview of mechanical characteristics of rock-like materials and the background of phase field development is given in the introduction. Second, the theoretical framework about the PFD and its anisotropic formulation are introduced to describe the brittle fracture. Then, the finite element discretization of the PFD is derived in Section 3. The Brazilian disk test samples with intact or pre-existing cracks are numerically simulated by the PFD, which is further validated by comparing it with the laboratory Brazilian splitting test results. Finally, some conclusions and further research on the PFD are given at the end of this work.

PHASE FIELD METHOD (PFD) Theoretical Foundation Description
The rock-like material failure is the result of crack initiation and propagation. How to accurately predict the crack growth path of engineering materials has always been the focus of scholars. The PFD is an attractive theory to represent discontinuous characteristics near the crack tip by a phase field variable. An elastic body Ω ⊂ R n dim (n dim ∈ {1, 2, 3}) is shown in Figure 1. The body boundaries contain an external boundary zΩ and an internal discontinuous boundary Γ within the material. The displacement (Dirichlet) boundary zΩ u and the traction (Neumann) boundary zΩ t constitute the external boundary, that is, zΩ u ∩ zΩ t zΩ. n is the outward unit vector normal to boundary zΩ. The crack width is abstractly represented by the regularization parameter ℓ 0 , which is used as a width scale of the crack in the numerical simulation. The body force b is applied throughout the body, and the traction force t acts on zΩ t . The sharp discontinuous field is represented by a phase field variable ϕ, which changes from 0 to 1. The material element is fully fractured when ϕ equals to 1, whereas the material element is undamaged when ϕ equals to 0.
According to the first law of thermodynamics, the brittle fracture of the material is the result of the system from a nonequilibrium state to equilibrium state, with the energy conversion from the elastic storage energy to the dissipation energy [23]. Moreover, the energy release always drives the total energy to minimize based on the theorem of minimum potential energy. The crack begins to initiate and propagates when the storage energy exceeds the material resistance. The fracture process can be described by adding an auxiliary field according to the variational approach [24]. The total potential energy Ψ(ϕ, u) of the system consists of two parts: internal potential energy Ψ int (ϕ, u) and external energy Ψ ext (u), namely, where Ψ int can be expressed as the sum of the free energy and fracture energy in Eq.2.a, and Ψ ext is formulated by the body force b and the boundary force t, presented in Eq.2.b. Frontiers in Physics | www.frontiersin.org December 2021 | Volume 9 | Article 809417 where ψ(ε) is the strain energy density, c(ϕ, ∇ϕ) represents the crack surface density per unit volume, and G c is the critical fracture energy density. ψ(ε) needs to further consider the degradation of phase field in fracture modeling. And it can be defined by ψ(ε) g(d)ψ 0 (ε). Similar to the previous published literature [25], the most used form of g(d) is taken as the following formulation.
where parameter κ is selected as a small value to ensure the numerical stability. The elastic strain energy ψ 0 (ε) of an undamaged material is given by the volumetric-deviatoric decomposition.
where λ and μ are Lamé coefficients, expressed as λ and μ E 2(1+υ) . And the stress tensor of the bulk matrix can be deduced by As same as the published literatures [19,26], the crack surface density per unit volume c(ϕ, ∇ϕ) is defined as The crack propagation is a process to minimize the total energy functional. The extreme condition can be determined by the variational principle. The differential form of the total energy with respect to displacement u, phase field ϕ, and its gradient ∇ϕ can be expressed as Eq. (7) is satisfied for any increments such as δϕ, δu, and δ∇ϕ. The governing equations and boundary conditions are derived by substituting Eq. (1), (2), and (4) into Eq. (7).
The governing equations are as follows: And the boundary conditions are as follows where the first equation in Eq. (9) is also called the Dirichlet (displacement) boundary condition and the last two equations are called the Neumann boundary condition. Given the irreversibility of cracks growth, a history-field H(x, t) is adopted for driving the crack propagation and ensuring the monotonical accumulation of the phase field variable.
Anisotropic Formulation of PFM Due to the unilateral effect of rock-like materials, the mechanical properties of material deterioration are only restricted in the tensile condition in this work. The formulation of the strain energy function ψ(ε) requires to be written as the tension and compression parts, namely.
where the tensile part and the compression part of ψ(ε) are further expressed as where the bracket (· + ) represents the positive value and (· − ) represents the negative value. Therefore, the stress tensor can be divided into two parts as seen in.
where the constitutive relation of the bulk matrix is further expressed by the volumetric-deviatoric decomposition in the tensile and compression states.
Eq. (10) can prevent crack propagation healing. However, the same formulation of the strain energy density is not obviously suitable for various crack modes. In order to capture the tension and shear crack modes, a modified phase field model is proposed by distinguishing the critical release rates for these two crack modes [17]. The phase field governing equation can be rearranged as

FE IMPLEMENTATION OF THE PHASE FIELD MODEL
The PFD has been successfully solved by the finite element method (FEM) [15], material point method (MPM) [27], and numerical manifold method (NMM) [28]. Because the FEM is widely adopted to simulate in geotechnical engineering, the numerical implementation of the PFM is briefly introduced in the standard FEM in this work. The small deformation assumption is assumed. The strain tensor ε can be defined in terms of the gradient of displacement u.
The displacement field and the phase field are approximately discretized by the test function and nodal variables, as follows: where N u I and N ϕ I are, respectively, the shape functions for displacement field and phase field, andû andφ are the displacement vector and the phase field vector of element nodes, respectively. The strain tensor and the phase field gradient ∇ϕ can be expressed as The matrix forms of N u I , N where x and y are coordinate variables, and n is the number of element nodes. The brittle fracture modeling of the PFD is transformed to a multi-field problem (displacement field and phase field). By the extremum condition in Eq. (7), the residual vectors corresponding to u and ϕ can be derived in Eq. (22) and (23), respectively.
The above equations can also be given by the weak form of the governing equations. Then the staggered algorithm is widely applied for an alternate updating the displacement increment δu and of phase field increment δϕ during the phase field modeling. That is, the following equation system is iteratively solved by the Newton-Raphson method. where

NUMERICAL SIMULATION AND ANALYSIS
In order to avoid the disturbance damage caused by sample preparation, the tensile fracture is usually studied by an indirect test (such as Brazilian splitting test) in the laboratory. The tensile fracture of rock-like materials is first simulated by the PFD modeling of the Brazilian test of the intact sample. And then, the numerical modeling is further applied in simulating the Brazilian test of the disk sample with a single inclination crack. The anisotropic tensile fracture of rock-like materials is analyzed and compared with the test results. This section will study the tensile fracture of rock-like materials by PHD modeling. The plane strain condition is taken for all numerical tests.

PFD Simulation of Intact Disk Sample
In order to study tensile behaviors of rock-like materials, disk samples were prepared from Portland pozzolana cement (PPC), fine sands, and water [29]. The test results indicate that the tension strength of the intact rock-like material sample is 3.81 MPa, Young's modulus is 15 GPa, and the Poisson's ratio is 0.21. The above mechanical indices will be regarded as the basis of determination of model parameters in the following PFD simulation. The test structure consists of two rigid jaws and a disk sample with a radius of 50 mm. In the process of numerical simulation, the sample thickness is taken as 1 mm. The critical fracture energy can be approximately estimated as 1.1 J/m 2 to match the tensile strength of the material in the test. Because the shear energy is much more than the tensile energy, the relationship between the two is set as G II c 10G I c , as the same with the previous literature [17]. Therefore, G I c is 0.1 J/m 2 and G II c is 1 J/m 2 . The displacement-controlled mode is applied with the edge of the top jaw in Figure 2. And the loading rate is 1.0 × 10 −3 mm/step. To facilitate the numerical convergence of this problem, the top and bottom nodes of the disk sample are constrained in the horizontal direction. The structure is discretized into 81,960 mixed triangle and quadrilateral elements. The minimum mesh size is 0.25 mm at axisymmetric lines of the sample. The length scale parameter ℓ 0 affects the structure bearing capacity and crack width. In order to ensure that the simulation strength and test tensile strength are consistent, ℓ 0 is taken as 0.5 mm about twice the minimum size by repeated calculation in this work. In particular, ℓ 0 does not represent the actual crack width because of the smooth processing of the crack domain. The parameter κ 1.0 × 10 −7 is taken to avoid numerical singularity.
The PFD simulation results of the Brazilian test of intact disk samples are presented in Figure 3. The obvious force drop is found in the force-displacement curves of fissure rocks with various inclined cracks. The peak load is 0.598 KN, and the tensile strength can be calculated as 3.807 MPa, which is consistent with the tensile results in the laboratory. The crack path is middle in the sample and is parallel to the loading direction. This phenomenon is also the same with the failure mode in the test. Therefore, the effectiveness of the PFD simulation of Brazilian splitting is verified.

PFD Simulation of Fissure Disk Sample
The anisotropic tensile strength is further studied by the PFD modeling of Brazilian disk samples with various inclination cracks. The fissure samples are performed by inserting one crack in the intact disk. The inclination angle θ is between the length direction of the crack and the vertical direction of the   The force-displacement curves are plotted in Figure 5. The results show that the peak load decreases at first and then increases with the increasing inclination angle. The maximum peak load appears at 0°, and the minimum value is at 45°. Therefore, the pre-existing crack causes the anisotropy of tensile strength. In order to represent the anisotropic degree of the tensile strength, a strength ratio is defined by the ratio of peak load between the Brazilian test of the fissure sample and intact    sample. The strength ratio changing with the inclination angle is plotted in Figure 6. The changing trend follows the U-shaped distribution, which is also found in the compression test of jointed or bedded rock mass [11,30].
Based on the Brazilian tests of the disk sample with various inclination angles for a single crack [29], the failure modes between the PFD simulation and test are compared in Figure 7. The result indicates that the final crack path in the numerical simulation process is consistent with the test results. The wing crack is only found in the Brazilian test and the phase field simulation of the pre-existing crack sample. Crack propagation initiates at the crack tip when the inclination angle changes from 0°to 60°. Crack propagation initiates near the crack tip when the angle is 75°, and it initiates at the middle of the crack when the angle is 90°. And then all cracks grow to the loading position and coalesce to form a penetrating crack. Its growth path is curvilinear when the crack angle is at 15°, 30°, 45°, 60°, and 75°, while the path is a straight line, parallel to the loading direction at 0°and 90°. It should be noted that the final crack path is parallel to the loading direction regardless of the crack angle. This is also the main feature of the wing crack development.

CONCLUSION
The tensile failure behaviors of rock-like materials are studied by the PFD simulation of the Brazilian disk test in this work. The anisotropy tensile fracture is considered by a single crack with various inclination angles in the anisotropic PFD simulation. The simulation results indicate that the peak load is weakened due to the influence of pre-existing cracks. And the changing of peak load with the inclination angle follows a U-shaped trend. The wing crack in the Brazilian test can be reproduced by the PFD simulation. The crack propagation path is affected by the pre-existing crack. Crack propagation initiates at the crack tip when the inclination angle changes from 0 to 60. Crack propagation initiates near the crack tip when the angle is 75, and crack propagation initiates at the middle of the crack when the angle is 90°. Finally, all cracks grow to the loading position and coalesce to form a penetrating crack. The PFD can better simulate the tensile failure of rock-like materials. More attention will be performed on the PFD simulation on compression-shear failure of rock-like materials in the future.

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.