Phase-Field Modeling Fracture in Anisotropic Materials

*e phase-field method is a widely used technique to simulate crack initiation, propagation, and coalescence without the need to trace the fracture surface. In the phase-field theory, the energy to create a fracture surface per unit area is equal to the critical energy release rate.*erefore, the precise definition of the crack-driving part is the key to simulate crack propagation. In this work, we propose a modified phase-field model to capture the complex crack propagation, in which the elastic strain energy is decomposed into volumetric-deviatoric energy parts. Because of the volumetric-deviatoric energy split, we introduce a novel form of the crack-driving energy to simulate mixed-mode fracture. Furthermore, a new degradation function is proposed to simulate crack processes in brittle materials with different degradation rates. *e proposed model is implemented by a staggered algorithm and to validate the performance of the phase-field modelling, and several numerical examples are constructed under plane strain condition. All the presented examples demonstrate the capability of the proposed approach in solving problems of brittle fracture propagation.


Introduction
Crack propagation is an active research topic in mechanical, energy, and environmental engineering, such as underground excavation, oil drilling, and nuclear waste storage [1][2][3], during the past decades. In particular, predictive investigations of crack-induced failure in rocks or rock-like materials are a complex problem due to the presence of preexisting fractures and voids which impact the strength and other mechanical properties. For geological materials like concrete and gypsum at no/low confinement, the failure mode is brittle fracturing which can be explained by Griffith's theory [4] and assuming that the energy to create a fracture surface per unit area is equal to the critical energy release rate G c . Based on the Griffith principle, many numerical computational methods have been developed. Numerical techniques for simulating crack propagation can be categorised into discrete and diffuse/smear methods depending upon how they handle the discontinuity.
Discrete methods attempt to capture the exact topology either in an explicit way or in an implicit manner. For instance, the extended finite-element method (XFEM) [5,6] has become a popular tool to consider the discontinuities. It enables the accurate approximation of solutions with jumps within elements through additional enrichment functions of discontinuous and asymptotic fields, thereby avoiding remeshing the cracked domain. Nevertheless, algorithmic tracking of the evolution of complex fracture surfaces is a tedious task in the numerical implementation. e cracking particles method (CPM) [7][8][9] is a pragmatic alternative to explicit modelling of crack surfaces in which a crack is represented by a set of cracking particles that can be easily updated when the crack propagates.
Diffuse methods for fracture modelling are based on the assumption that the discontinuity in the cracked material is not sharp but can be interpreted as smeared damage [10,11]. Recently, the phase-field method [12,13] has been attracting much attention because of its simplicity for numerical implementation. In the phase-field model, a smooth boundary of the phase-field is employed to approximate the internal discontinuity boundary of a crack. e use of the phase-field model for fracture can circumvent the complexity of tracking crack propagation that is typically required in discrete models. In this research, we adopt the phase-field method and propose novel modifications in order to elegantly simulate complicated fracture processes in geological materials. e phase-field method was proposed by Bourdin et al. [14], and further developed by Borden et al. [15] and Miehe et al. [16,17]. Due to its strong ability to simulate complex fracture processes such as nucleation, propagation, and branching, great efforts and extensions have also been done for brittle fracture [18][19][20][21][22][23][24], quasi-brittle fracture [25,26], ductile fracture [27][28][29][30], dynamic fracture [15,31], and multi-physicals fracturing problem [32,33] for various materials.
e fundamental concept behind this model is the introduction of a scalar damage field, which ranges from 0 (undamaged material) to 1 (fully damage material), to represent the degree of fracture or damage of the material [13]. e crack propagation problem is sequentially recast as a standard multifield problem which can then be handled using the conventional finite-element method for both twodimensional and three-dimensional cases. As a result, issues related to crack discontinuities are circumvented, and complex crack evolution can be treated naturally without difficulty. e phase-field model was developed within the framework of a variational principle of fracture [14], and this further enhances its attraction. e principle [12,34] can be regarded as a generalization of Griffith's theory which enables it to predict not only crack initiation but also the crack propagation path. Furthermore, the solution of the variational principle is globally rather than locally optimal; thus, any new crack nucleation can be detected naturally without being specified in advance [14,34].
Despite these contributions, these phase-field models assume that the critical energy release rates of different crack modes are the same, but in fact for many materials are not. In rock-like materials, such as concrete and gypsum, the critical energy release rate for Mode-I fracture is significantly lower than that for Mode-II fracture. is feature complicates the modelling of crack phenomena for rocks. In a rock specimen with a single inclined flaw under compression, wing cracks emerge first followed by secondary cracks. It has been reported [35,36] that the wing crack is a Mode-I crack (tensile crack), while the secondary crack is usually Mode-II crack (shear crack). e sequential appearance of wing and secondary cracks can be attributed to the considerable difference in the critical energy release rates for different crack modes.
us, it has been suggested [37] that this phenomenon cannot be captured by using the traditional critical energy release criterion, which does not account for the different critical energy release rates for mode-I and mode-II fractures.
In this paper, a modified phase-field model for brittle fracture is proposed to distinguish between the critical release rates for mode-I and mode-II cracks. is is achieved by partitioning the active energy density into distinct parts corresponding to different crack modes. e present paper is organized as follows: in Section 2, the fundamentals of the phase-field method are first briefly summarized. Section 3 presents the numerical implementation in detail. In Section 4, the accuracy of the numerical simulation is verified by using one-element example; then a number of classical experimental tests are simulated. Finally, Section 5 concludes the paper.

Fundamentals of the Phase-Field Method
In this section, an arbitrary bounded computational domain Ω ⊂ R n dim (n dim � 2, 3) is considered with external boundary zΩ ⊂ R n dim − 1 and internal discontinuity Γ ⊂ R n dim − 1 as illustrated in Figure 1.
e external boundary zΩ is decomposed into two disjoint parts zΩ u and zΩ t , i.e., zΩ u ∩ zΩ t � ∅ and zΩ u ∪ zΩ t � zΩ. e domain Ω is subjected to the Dirichlet boundary conditions, u(x) for x ∈zΩ u , and the Neumann boundary zΩ t ∈ Ω, with the corresponding outward unit normal vectors n u and n t . e Neumann conditions impose the traction t * (x) on zΩ t .

Energy Functional.
As already mentioned in Introduction, the phase-field approach to brittle fracture is based on the work of Bourdin et al. [13] and consists in the regularization of the variational formulation of Griffith's theory of brittle fracture, first proposed in 1998 by Francfort and Marigo [12]. Neglecting inertia effects and assuming quasistatic conditions, the total energy functional Ψ(u, Γ) of a solid is the sum of elastic energy ψ(ε), fracture energy, and external work. us, the total power is written as with the linear strain tensor ε � ε(u) given by where u denotes the displacement field of the body Ω and ∇(·) and ∇ s (·) denote the gradient and symmetric gradient operators, respectively, with (·) T being the transpose operator.
Regarding the choice of the split of the energy density function, we introduce the volumetric-deviatoric energy split proposed by Lancioni and Royer-Carfagni [38]. Assuming the solid is isotropic and linear elastic, the elastic energy density ψ ε (ε) � ψ dev (ε) + ψ vol (ε) is written as where ε dev � ε − (1/3)tr(ε)I is the deviatoric component of the strain tensor ε, K � λ + (2/3)μ is the bulk modulus of the material, λ and μ are the Lamé constants, I the second-order unit tensor, and C dev and C vol are the deviatoric and volumetric parts of the elasticity tensor C, and defined as follows:

Phase-Field Approximation for Fracture Energy.
In the phase-field model, a scalar field (phase-field) is used to diffuse the sharp crack topology [16,17] over a certain domain which avoids complex crack tracking procedures and an explicit representation of the crack surface as in discrete crack approaches [39]. erefore, a narrow transition band connects the fully fractured and intact domains with the displacement being still continuous. Within the context of quasi-static brittle fracture in elastic solids, the cracks are approximated as bands of finite thickness characterised by the phase-field d(x, t) ∈ [0, 1] as shown in Figure 1, which satisfies the following conditions: is variable indicates the damage in the material. e material is fully broken for d � 1, and d � 0 represents the intact state. A typical one-dimensional phase-field is approximated with the exponential function: e length-scale parameter ℓ plays an important role which controls the transition region between the fracture and intact material.
For 2D and 3D problems, the crack surface density per unit volume of the solid is given by [17] us, exploiting equation (6), the fracture energy in equation (1) can be rewritten as

Governing Equations for Evolution of the Phase-Field.
It is noticed that the phase-field formulation equation (1) does not distinguish fracture behaviour during tension and compression. erefore, based on the volumetric-deviatoric decomposition of the elastic energy, a more adequate choice for a split would be adopted [40]: We follow Ambati et al. [41] and assume that the phasefield affects the positive part of the elastic energy: where 0 < k ≪ 1 is the parameter that stabilizes the stiffness matrix to ensure the numerical convergence. Taking advantage of equations (8) and (10), equation (1) can be rewritten as Following Miehe et al. [16], we obtain two coupled local equations: where σ is the stress tensor, defined as e degradation function g(d) characterizes the ratio of residual strain energy and total strain energy during the crack evolution. e degradation g(d) function is a monotonically g(0) � 1 and g(1) � 0. Actually, the selection of degradation function depends on the Figure 1: Schematic depiction of a solid body Ω with a strong discontinuity Γ (a), modeled by the phase-field fracture method (b). e parameter ℓ controls the width of the diffused fracture zone.
Advances in Civil Engineering mechanical properties of materials. For this reason, inspired by the work of [42][43][44], a new degradation function is proposed in this study to describe a large range of failure processes: In this degradation function, M is a nondimensional material parameter, which characterizes the degradation rate of strain energy with the evolution of the phase-field, as shown in Figure 2. To prevent crack healing, we take advantage of the definition of local history field proposed by Miehe and Schänzel [45] to set up relationship between phase-field variable and maximum reference energy in history, whereby the following relation is given: where x is the material point in a reference body and t is the pseudotime. us the phase-field evolution equation is finally given as e critical energy release rate for brittle tension fracture is significantly lower than that for the compressive-shear fracture. To capture this feature, in this work, the model incorporates different contributions of energy components to crack growth, which is able to capture tensile and shear cracks according to the strain states. erefore, equation (16) can be rewritten as where the parameters G t and G s are the critical energy rates for tensile and compressive-shear fracture, respectively. H + vol and H ± dev represent parts of the volumetric-deviatoric split of the elastic strain energy: where H(·) is the Heaviside function.

Numerical Algorithm
In this section, we describe the numerical algorithm on implementation of the phase-field modelling of fracture propagation in isotropic medium. We use the finite-element method to discretize the spatial domain and a staggered scheme to solve the coupling equations and for pursuing a higher convergence rate.

Finite-Element Discretization.
e weak forms of the governing equations are given by e quadrilateral four-node element in 2D elements and the hexahedral eight-node in 3D elements are implemented to discretize the bulk domain Ω. e node values u i and d i are discretized as follows: where n is the total number of nodes per element. N i denotes the shape function associated with node i. e corresponding matrices of the spatial derivatives can be expressed as en it is possible to express the gradients:  Advances in Civil Engineering According to 21) and (22), the contributions of oneelement at node i to the residual of the overall systems of equations are given as where F u,ext i and F u,int i denote the external forces and inner forces, respectively, which correspond to the displacement, whereas F d,int i can be explained as inner forces to the phasefield. We used the Newton-Raphson procedure to obtain the solutions by making R u � 0 and R d � 0. e corresponding tangents on the element level can be obtained based on the inner forces: where D is the fourth-order elasticity tensor given by

Staggered Scheme.
A staggered scheme is adopted here; where at each increment step, the governing equation (23) is solved for u by freezing d. After updating the elastic strains, we solve the phase-field evolution for d. e step is repeated until a convergence criterion is reached. A strong coupling between the displacement field and phase-field introduces a relatively complicated numerical process. erefore, a N-R iteration procedure is needed for solving this nonlinear equation. For the sake of clarity, we show the flowchart of the N-R iterative scheme in Figure 3.

Verification of the Proposed Approach
In this section, starting with the simplest case where we compare different methods for one-element, more and more complex cases are introduced. In all cases, the relevant numerical parameters are summarized, then the results are shown and interpreted. In all cases (plane strain), the thickness of the element is 1 mm. e mesh is densified where the crack is expected to propagate, the size is specified in the text, and the mesh is shown in some of the figures. According to the results of Miehe et al. [17], the length-scale parameter is always taken two times larger than the smallest element around the crack path.

One-Element.
One 2D plane strain element is the simplest case, to verify the correctness of the proposed model. e geometry and the boundary conditions are illustrated in Figure 4(a). e bottom nodes are constrained in both directions, whereas we allow the top nodes to slide vertically. For the case of homogeneous materials, E � 2.1 × 10 5 MPa, v � 0.3, G t � 5 × 10 − 3 kN/mm, G s � 5 × 10 − 2 kN/mm, and ℓ � 0.1 mm as reported in [43,46]. e loading history is divided into 1000 steps with a constant increment Δu � 10 − 4 mm. Figure 4(b) shows the comparisons of macroscopic stress-strain relation and damage evolution between the analytical solution and the proposed method. Obviously, these two solutions well recover each other as illustrated in this figure. is shows that the proposed phase-field method can well describe the damage process in homogeneous materials.

ree-Point Bending Test.
In this section, an example of the problem of three-point bending test is presented performed by Perdikaris and Romeo [47] and widely taken as the benchmark model for numerical investigation on model-I fracture energy for plane concrete. erefore, some previous experiences and results are available for comparison. e geometry and loading conditions of the specimen are Advances in Civil Engineering shown in Figure 5. Two different finite-element meshes are considered and shown in Figure 6. It amounts to 6252 elements in coarse mesh and 18474 elements in the fine mesh. e displacement imposed onto the notched beam is Δu � 10 − 2 mm.
Young's modulus is taken as E � 4.83 × 10 4 MPa. e parameter G t on the one hand has been provided experimentally as G t � 0.0451 N/mm and, on the other hand, can be calculated by the linear elastic fracture mechanics in consideration of size effect [48]. According to the work elaborated in [48], one takes G t � 0.029 N/mm. Figure 7 illustrates the crack evolution process predicted from the phase-field model. e crack initiates from the notch tip and propagates towards the top surface of the beam. Comparing the crack evolution processes in Figure 6, the speed of crack growth in the case of G t � 0.0290 N/mm is faster than the other one, and the corresponding compressive-shear critical energy release rates are G s � 0.290 N/ mm and G s � 0.451 N/mm, respectively. Figure 8 compares the load-displacement response from this study to the experimental data. ere is no detectable mesh dependency for the two meshes under consideration,  Advances in Civil Engineering and the beam is subject to noncyclic loading. In the case of G t � 0.0290 N/mm and G s � 0.290 N/mm, the maximum force that the beam can sustain from the simulation is 12.3 kN which is very close to the experimental data (12.0 kN). However, taking G t � 0.0451 N/mm and G s � 0.451 N/mm, the peak load is about 13 kN which is 8.3% higher than the experimental result. e simulated forces both drop much faster than the experimental data after reaching the bearing capacity. ese deviations on one aspect stem from the fact that a linear fracture model is used in the simulation, whereas the real fracture is nonlinear; the model neglects any plastic deformation. It can be clarified by examining the loaddisplacement curves from both cyclic and noncyclic loadings.

L-Shaped Panel Test.
In this section, we simulate crack propagation in an L-shaped slab. e geometry and boundary conditions of the problem are depicted in Figure 9(a). e experimental results are taken from [49], and Figure 9(b) illustrates the crack path obtained from the phase-field modelling and superimposes the range of experimentally obtained crack paths (shaded region). e material parameters are chosen as follows [50]: Young's modulus E � 2.0 × 10 4 MPa, Poisson's ratio v � 0.18, ℓ � 10 mm, tension fracture energy G t � 1.3 × 10 − 4 kN/mm, and compressive-shear fracture energy G s � 1.3 × 10 − 3 kN/ mm [51]. e simulation is led with the step increment Δu � 2 × 10 − 3 mm. e computational domain is discretized using a total of 51,766 elements with fine meshes assigned to the critical zone. Figure 10 illustrates the crack progression at several loading stages. e corresponding load-displacement curves are presented in Figure 11. For ℓ � 10 mm, the crack path from the simulation is located at the path observed in the real test. e corresponding load-displacement curves are shown  Advances in Civil Engineering in Figure 11(a); the global response matches the experimental peal loads. e reaction force from the simulation increases steadily to a maximum value and then drops sharply. To further investigate the effect of the length scale, a second simulation is performed with a larger length scale (ℓ � 20 mm). For ℓ � 20 mm, a diffusive crack path is obtained due to the relatively large length scale. Although the maximum force obtained by using ℓ � 20 mm is even lower as shown in Figure 11(b), the predicted crack path agrees well with the experimental data ( Figure 10(b)). e proposed model is fairly in good agreement with the experimental results in terms of crack pattern.

Compression of a Rock Plate with Double Flaws.
To further verify the modified phase-field model proposed and highlight its capabilities, crack propagation involving fracture coalescence is considered. Figure 12 shows a schematic illustration of the problem where a specimen with double inclined open flaws is loaded under uniaxial compression. Such a test has been widely investigated experimentally [35] and numerically [52] by means of prefracture specimens of gypsum under uniaxial compression. In the test, the gypsum specimen is 76.2 mm long and 152.4 mm high. e length and the width for the flaws are 12.7 mm and 0.1 mm, respectively. e geometry of the flaws is represented by the terminology "β − s − c," as shown in Figure 12. Herein we consider the case of "45 -2a -2a" geometry (2a � 12.7 mm). Displacement increment Δu � 0.002 mm is prescribed in line with laboratory tests [35]. e mechanical properties of the material are E � 5 × 10 3 MPa, Poisson's ratio v � 0.24, ℓ � 0.5942 mm, tension fracture energy G t � 5.0 × 10 − 6 kN/mm, and compressive-shear fracture energy G s � 1.0 × 10 − 2 kN/mm. e numerical simulation gives a crack pattern similar to the experimental observation, that is, the growth of four wing cracks and the coalescence of a secondary and one wing crack. e fracturing process is explained in more detail in Figure 13. Initially a stable growth of four wing cracks is seen (Figure 13(b)). e outer wing cracks continue to grow for a while, and a shear crack is subsequently initiated close to the inner tip of the upper notch; it coalesces with the inner wing  crack of the lower notch ( Figure 13(d)). en, the inner wing crack of the upper notch is located in a strong shear stress state and approaches the preexisting flaw ( Figure 13(e)). Eventually, they connect with each other. Up to this point, the simulated crack pattern is quite similar to the experimental [35].

Conclusions
is paper presents a new framework of the phase-field method based on the split of the fracture energy release rate for simulation of crack propagation in geotechnical materials. e critical release rate for tensile cracks is significantly lower than the energy release rate for shear cracks. In the proposed approach, the crack-driving energy is identified, and a new degradation function is introduced, in which a nondimensional parameter is used to describe crack propagation of brittle materials with different weakening rates.
Several numerical examples are carried out. Firstly, the modified phase-field model is validated by the widely used benchmark example, and the simulation results well agree with the analytical results. Furthermore, to demonstrate the capability of the modified phase-field model in simulation crack propagation and bifurcation in brittle materials, several numerical examples are presented. e presented phenomenon of crack propagation shows that the modified phase-field fracture model gives results in good agreement with the experimental observations both with respect to crack patterns and critical stress loads.
In summary, the modified model is capable of describing the failure behaviour for brittle and presenting the failure processes. Moreover, it is noteworthy that the effects of length parameter and critical energy release rate on simulation results need to be investigated further [53]. In future work, the proposed approach can be extended to predict crack propagation in multi-physics problems, for example, hydraulic fracture propagation and heat transfer.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.