# Phase-Field Modeling Fracture in Anisotropic Materials

Authors: Haifeng Li; Wei Wang; Yajun Cao; Shifan Liu
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The 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. Therefore, 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. The 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.


---

## Body

## 1. 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–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 Gc. 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. The cracking particles method (CPM) [7–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. The 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. The 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–24], quasi-brittle fracture [25, 26], ductile fracture [27–30], dynamic fracture [15, 31], and multi-physicals fracturing problem [32, 33] for various materials.The 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]. The crack propagation problem is sequentially recast as a standard multifield problem which can then be handled using the conventional finite-element method for both two-dimensional and three-dimensional cases. As a result, issues related to crack discontinuities are circumvented, and complex crack evolution can be treated naturally without difficulty. The phase-field model was developed within the framework of a variational principle of fracture [14], and this further enhances its attraction. The 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. This 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). The sequential appearance of wing and secondary cracks can be attributed to the considerable difference in the critical energy release rates for different crack modes. Thus, 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. This is achieved by partitioning the active energy density into distinct parts corresponding to different crack modes. The present paper is organized as follows: in Section2, 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.


## 2. Fundamentals of the Phase-Field Method

In this section, an arbitrary bounded computational domainΩ⊂ℝndimndim=2,3 is considered with external boundary ∂Ω⊂ℝndim−1 and internal discontinuity Γ⊂ℝndim−1 as illustrated in Figure 1. The external boundary ∂Ω is decomposed into two disjoint parts ∂Ωu and ∂Ωt, i.e., ∂Ωu∩∂Ωt=∅ and ∂Ωu∪∂Ωt¯=∂Ω. The domain Ω is subjected to the Dirichlet boundary conditions, ux for x∈∂Ωu, and the Neumann boundary ∂Ωt∈Ω, with the corresponding outward unit normal vectors nu and nt. The Neumann conditions impose the traction t∗x on ∂Ωt.Figure 1
Schematic depiction of a solid body Ω with a strong discontinuity Γ (a), modeled by the phase-field fracture method (b). The parameterℓ controls the width of the diffused fracture zone.
(a)(b)


### 2.1. 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 quasi-static conditions, the total energy functional Ψu,Γ of a solid is the sum of elastic energy ψε, fracture energy, and external work. Thus, the total power is written as(1)Ψu,Γ=∫Ωb⋅udΩ+∫∂Ωtt∗⋅udS−∫ΓGcdΓ−∫ΩψεεdΩ,with the linear strain tensor ε=εu given by(2)ε≡∇su=12∇u+∇uT,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(3)ψdevε=μεdev:εdev=12ε:ℂdev:ε,ψvolε=12Ktr2ε=12ε:ℂvol:ε,where εdev=ε−1/3trε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 ℂdev and ℂvol are the deviatoric and volumetric parts of the elasticity tensor ℂ, and defined as follows:(4)ℂvol=KI⊗I;ℂdev=ℂ−ℂdev.


### 2.2. 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]. Therefore, 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 dx,t∈0,1 as shown in Figure 1, which satisfies the following conditions:(5)d=0,if material is intact,1,if material is cracked.This variable indicates the damage in the material. The material is fully broken ford=1, and d=0 represents the intact state. A typical one-dimensional phase-field is approximated with the exponential function:(6)dx=e−x/ℓ.The 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](7)γd,∇d=d22ℓ+ℓ2∂d∂xi∂d∂xi.Thus, exploiting equation (6), the fracture energy in equation (1) can be rewritten as(8)∫ΓGcdS=∫ΩGcd22ℓ+ℓ2∂d∂xi∂d∂xidΩ.


### 2.3. 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. Therefore, based on the volumetric-deviatoric decomposition of the elastic energy, a more adequate choice for a split would be adopted [40]:(9)ψ+=12Ktrε+2+μεdev:εdev,ψ−=12Ktrε−2.We follow Ambati et al. [41] and assume that the phase-field affects the positive part of the elastic energy:(10)ψε=1−kgd+kψ++ψ−,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(11)Ψu,Γ=∫Ωb⋅udΩ+∫∂Ωtt∗⋅udS−∫ΓGcd22ℓ+ℓ2∂d∂xi∂d∂xidΩ−∫Ω1−kgd+kψ++ψ−dΩ.Following Miehe et al. [16], we obtain two coupled local equations:(12)divσ+b=0,1−kg˙dψ++Gcdℓ−ℓΔd=0,where σ is the stress tensor, defined as(13)σ=gdKtrε+I+2μεdev+Ktrε−I.The degradation functiongd characterizes the ratio of residual strain energy and total strain energy during the crack evolution. The degradation gd function is a monotonically g0=1 and g1=0. Actually, the selection of degradation function depends on the mechanical properties of materials. For this reason, inspired by the work of [42–44], a new degradation function is proposed in this study to describe a large range of failure processes:(14)gd=1−k1−d2M+M−11−d2+k.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:(15)ℋx,t=maxs∈0,tψε+,inΩ×0,T,where x is the material point in a reference body and t is the pseudotime. Thus the phase-field evolution equation is finally given as(16)lΔd=dℓ+1−kg˙dℋGc.Figure 2
The evolution of degradation function for differentM.The 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. Therefore, equation (16) can be rewritten as(17)lΔϕ=dℓ+1−kg˙dℋvol+Gt+ℋdev+Gt+ℋdev−Gs,where the parameters Gt and Gs are the critical energy rates for tensile and compressive-shear fracture, respectively. ℋvol+ and ℋdev± represent parts of the volumetric-deviatoric split of the elastic strain energy:(18)ℋvol+=12Ktrε+,ℋdev+=μεdev:εdevHtrε,ℋdev+=μεdev:εdev1−Htrε,where H⋅ is the Heaviside function.


## 2.1. 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 quasi-static conditions, the total energy functional Ψu,Γ of a solid is the sum of elastic energy ψε, fracture energy, and external work. Thus, the total power is written as(1)Ψu,Γ=∫Ωb⋅udΩ+∫∂Ωtt∗⋅udS−∫ΓGcdΓ−∫ΩψεεdΩ,with the linear strain tensor ε=εu given by(2)ε≡∇su=12∇u+∇uT,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(3)ψdevε=μεdev:εdev=12ε:ℂdev:ε,ψvolε=12Ktr2ε=12ε:ℂvol:ε,where εdev=ε−1/3trε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 ℂdev and ℂvol are the deviatoric and volumetric parts of the elasticity tensor ℂ, and defined as follows:(4)ℂvol=KI⊗I;ℂdev=ℂ−ℂdev.


## 2.2. 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]. Therefore, 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 dx,t∈0,1 as shown in Figure 1, which satisfies the following conditions:(5)d=0,if material is intact,1,if material is cracked.This variable indicates the damage in the material. The material is fully broken ford=1, and d=0 represents the intact state. A typical one-dimensional phase-field is approximated with the exponential function:(6)dx=e−x/ℓ.The 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](7)γd,∇d=d22ℓ+ℓ2∂d∂xi∂d∂xi.Thus, exploiting equation (6), the fracture energy in equation (1) can be rewritten as(8)∫ΓGcdS=∫ΩGcd22ℓ+ℓ2∂d∂xi∂d∂xidΩ.


## 2.3. 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. Therefore, based on the volumetric-deviatoric decomposition of the elastic energy, a more adequate choice for a split would be adopted [40]:(9)ψ+=12Ktrε+2+μεdev:εdev,ψ−=12Ktrε−2.We follow Ambati et al. [41] and assume that the phase-field affects the positive part of the elastic energy:(10)ψε=1−kgd+kψ++ψ−,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(11)Ψu,Γ=∫Ωb⋅udΩ+∫∂Ωtt∗⋅udS−∫ΓGcd22ℓ+ℓ2∂d∂xi∂d∂xidΩ−∫Ω1−kgd+kψ++ψ−dΩ.Following Miehe et al. [16], we obtain two coupled local equations:(12)divσ+b=0,1−kg˙dψ++Gcdℓ−ℓΔd=0,where σ is the stress tensor, defined as(13)σ=gdKtrε+I+2μεdev+Ktrε−I.The degradation functiongd characterizes the ratio of residual strain energy and total strain energy during the crack evolution. The degradation gd function is a monotonically g0=1 and g1=0. Actually, the selection of degradation function depends on the mechanical properties of materials. For this reason, inspired by the work of [42–44], a new degradation function is proposed in this study to describe a large range of failure processes:(14)gd=1−k1−d2M+M−11−d2+k.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:(15)ℋx,t=maxs∈0,tψε+,inΩ×0,T,where x is the material point in a reference body and t is the pseudotime. Thus the phase-field evolution equation is finally given as(16)lΔd=dℓ+1−kg˙dℋGc.Figure 2
The evolution of degradation function for differentM.The 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. Therefore, equation (16) can be rewritten as(17)lΔϕ=dℓ+1−kg˙dℋvol+Gt+ℋdev+Gt+ℋdev−Gs,where the parameters Gt and Gs are the critical energy rates for tensile and compressive-shear fracture, respectively. ℋvol+ and ℋdev± represent parts of the volumetric-deviatoric split of the elastic strain energy:(18)ℋvol+=12Ktrε+,ℋdev+=μεdev:εdevHtrε,ℋdev+=μεdev:εdev1−Htrε,where H⋅ is the Heaviside function.


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


### 3.1. Finite-Element Discretization

The weak forms of the governing equations are given by(19)∫Ω−σ:δεdV+∫Ωb⋅δudV+∫∂Ωtt∗⋅δudS=0,∫Ω1−kg′dℋGcδϕdΩ+∫Ωℓ∇d⋅∇δd+1ℓdδddΩ=0.The quadrilateral four-node element in 2D elements and the hexahedral eight-node in 3D elements are implemented to discretize the bulk domainΩ. The node values ui and di are discretized as follows:(20)u=∑i=1nNiui,d=∑i=1nNidi,where n is the total number of nodes per element. Ni denotes the shape function associated with node i. The corresponding matrices of the spatial derivatives can be expressed as(21)Biu=Ni,x00Ni,yNi,yNi,x,Bid=Ni,xNi,y.Then it is possible to express the gradients:(22)ε=∑i=1nBiuui,∇d=∑i=1nBiddi.According to21) and (22), the contributions of one-element at node i to the residual of the overall systems of equations are given as(23)Riu=Fiu,ext−Fiu,int=∫∂ΩtNit∗dS−∫ΩBiuσdΩ,Rid=−Fid,int=∫Ω1−kg′dℋGcNi+ℓBidT∇d+1ℓdNidΩ,where Fiu,ext and Fiu,int denote the external forces and inner forces, respectively, which correspond to the displacement, whereas Fid,int can be explained as inner forces to the phase-field. We used the Newton–Raphson procedure to obtain the solutions by making Ru=0 and Rd=0. The corresponding tangents on the element level can be obtained based on the inner forces:(24)Kijuu=∂Fiu,int∂uj=∫ΩBiuTDBjudΩ,Kijdd=∂Fid,int∂dj=∫ΩBidTℓBid+Ni1−kg2dℋvol+Gt+ℋdev+Gt+ℋdev−Gs+1ℓNjdΩ,where D is the fourth-order elasticity tensor given by(25)D=1−kgd+kℂvol+ℂdev,trε≥0,ℂvol+1−kgd+kℂdev,trε<0.


### 3.2. 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. The 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. Therefore, 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.Figure 3
Flowchart of the staggered solution procedure.


## 3.1. Finite-Element Discretization

The weak forms of the governing equations are given by(19)∫Ω−σ:δεdV+∫Ωb⋅δudV+∫∂Ωtt∗⋅δudS=0,∫Ω1−kg′dℋGcδϕdΩ+∫Ωℓ∇d⋅∇δd+1ℓdδddΩ=0.The quadrilateral four-node element in 2D elements and the hexahedral eight-node in 3D elements are implemented to discretize the bulk domainΩ. The node values ui and di are discretized as follows:(20)u=∑i=1nNiui,d=∑i=1nNidi,where n is the total number of nodes per element. Ni denotes the shape function associated with node i. The corresponding matrices of the spatial derivatives can be expressed as(21)Biu=Ni,x00Ni,yNi,yNi,x,Bid=Ni,xNi,y.Then it is possible to express the gradients:(22)ε=∑i=1nBiuui,∇d=∑i=1nBiddi.According to21) and (22), the contributions of one-element at node i to the residual of the overall systems of equations are given as(23)Riu=Fiu,ext−Fiu,int=∫∂ΩtNit∗dS−∫ΩBiuσdΩ,Rid=−Fid,int=∫Ω1−kg′dℋGcNi+ℓBidT∇d+1ℓdNidΩ,where Fiu,ext and Fiu,int denote the external forces and inner forces, respectively, which correspond to the displacement, whereas Fid,int can be explained as inner forces to the phase-field. We used the Newton–Raphson procedure to obtain the solutions by making Ru=0 and Rd=0. The corresponding tangents on the element level can be obtained based on the inner forces:(24)Kijuu=∂Fiu,int∂uj=∫ΩBiuTDBjudΩ,Kijdd=∂Fid,int∂dj=∫ΩBidTℓBid+Ni1−kg2dℋvol+Gt+ℋdev+Gt+ℋdev−Gs+1ℓNjdΩ,where D is the fourth-order elasticity tensor given by(25)D=1−kgd+kℂvol+ℂdev,trε≥0,ℂvol+1−kgd+kℂdev,trε<0.


## 3.2. 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. The 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. Therefore, 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.Figure 3
Flowchart of the staggered solution procedure.


## 4. 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. The 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.


### 4.1. One-Element

One 2D plane strain element is the simplest case, to verify the correctness of the proposed model. The geometry and the boundary conditions are illustrated in Figure4(a). The 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×105MPa, v = 0.3, Gt=5×10−3 kN/mm, Gs=5×10−2 kN/mm, and ℓ=0.1 mm as reported in [43, 46]. The loading history is divided into 1000 steps with a constant increment Δu=10−4 mm.Figure 4
(a) Studied one-element and (b) comparisons of strain-force curves and macroscopic damage evolution between analytical solution and numerical results.
(a)(b)Figure4(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. This shows that the proposed phase-field method can well describe the damage process in homogeneous materials.


### 4.2. Three-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. Therefore, some previous experiences and results are available for comparison. The geometry and loading conditions of the specimen are 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. The displacement imposed onto the notched beam is Δu=10−2 mm.Figure 5
Symmetric three-point bending test: geometry and boundary conditions.Figure 6
FE-meshes used for computation of three-point bending on notched beam: (a) coarse; (b) fine.
(a)(b)Young’s modulus is taken asE=4.83×104 MPa. The parameter Gt on the one hand has been provided experimentally as Gt=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 Gt=0.029 N/mm.Figure7 illustrates the crack evolution process predicted from the phase-field model. The 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 Gt=0.0290 N/mm is faster than the other one, and the corresponding compressive-shear critical energy release rates are Gs=0.290 N/mm and Gs=0.451 N/mm, respectively.Figure 7
Predicted crack paths of the three-point bending test. (a)Gt=0.0451 N/mm. (b) Gt=0.029 N/mm.
(a)(b)Figure8 compares the load-displacement response from this study to the experimental data. There is no detectable mesh dependency for the two meshes under consideration, and the beam is subject to noncyclic loading. In the case of Gt=0.0290 N/mm and Gs=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 Gt=0.0451 N/mm and Gs=0.451 N/mm, the peak load is about 13 kN which is 8.3% higher than the experimental result.Figure 8
Calculated load-displacement curves of the three-point bending test.The simulated forces both drop much faster than the experimental data after reaching the bearing capacity. These 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 load-displacement curves from both cyclic and noncyclic loadings.


### 4.3. L-Shaped Panel Test

In this section, we simulate crack propagation in an L-shaped slab. The geometry and boundary conditions of the problem are depicted in Figure9(a). The 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). The material parameters are chosen as follows [50]: Young’s modulus E=2.0×104 MPa, Poisson’s ratio v=0.18, ℓ=10 mm, tension fracture energy Gt=1.3×10−4 kN/mm, and compressive-shear fracture energy Gs=1.3×10−3 kN/mm [51]. The simulation is led with the step increment Δu=2×10−3 mm.Figure 9
L-shaped panel test. (a) Geometry and boundary conditions. (b) Experimentally observed crack pattern.
(a)(b)The computational domain is discretized using a total of 51,766 elements with fine meshes assigned to the critical zone. Figure10 illustrates the crack progression at several loading stages. The corresponding load-displacement curves are presented in Figure 11. For ℓ=10mm, the crack path from the simulation is located at the path observed in the real test. The corresponding load-displacement curves are shown in Figure 11(a); the global response matches the experimental peal loads. The 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 (ℓ=20mm). For ℓ=20mm, a diffusive crack path is obtained due to the relatively large length scale. Although the maximum force obtained by using ℓ=20mm is even lower as shown in Figure 11(b), the predicted crack path agrees well with the experimental data (Figure 10(b)). The proposed model is fairly in good agreement with the experimental results in terms of crack pattern.Figure 10
L-shaped panel test. Damage profiles for various length-scale parameters at displacementsu = 0.22 mm, 0.30 mm, 0.45 mm, and 1.0 mm. (a) ℓ=10mm. (b) ℓ=20mm.
(a)(b)Figure 11
L-shaped panel: curves of load versus displacement for various length-scale parametersℓ. (a) Length-scale parameter ℓ=10mm. (b) Length-scale parameter ℓ=20mm.
(a)(b)


### 4.4. 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. Figure12 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. The length and the width for the flaws are 12.7 mm and 0.1 mm, respectively. The 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.002mm is prescribed in line with laboratory tests [35].Figure 12
Geometries of the model rock specimen.The mechanical properties of the material areE=5×103 MPa, Poisson’s ratio v=0.24, ℓ=0.5942 mm, tension fracture energy Gt=5.0×10−6 kN/mm, and compressive-shear fracture energy Gs=1.0×10−2 kN/mm. The 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. The fracturing process is explained in more detail in Figure 13. Initially a stable growth of four wing cracks is seen (Figure 13(b)). The 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)). Then, 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].Figure 13
Crack evolution processes with displacement at (a)u = 0.312 mm; (b) u = 0.632 mm; (c) u = 0.752 mm; (d) u = 0.856 mm; (e) u = 0.912 mm.
(a)(b)(c)(d)(e)


## 4.1. One-Element

One 2D plane strain element is the simplest case, to verify the correctness of the proposed model. The geometry and the boundary conditions are illustrated in Figure4(a). The 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×105MPa, v = 0.3, Gt=5×10−3 kN/mm, Gs=5×10−2 kN/mm, and ℓ=0.1 mm as reported in [43, 46]. The loading history is divided into 1000 steps with a constant increment Δu=10−4 mm.Figure 4
(a) Studied one-element and (b) comparisons of strain-force curves and macroscopic damage evolution between analytical solution and numerical results.
(a)(b)Figure4(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. This shows that the proposed phase-field method can well describe the damage process in homogeneous materials.


## 4.2. Three-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. Therefore, some previous experiences and results are available for comparison. The geometry and loading conditions of the specimen are 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. The displacement imposed onto the notched beam is Δu=10−2 mm.Figure 5
Symmetric three-point bending test: geometry and boundary conditions.Figure 6
FE-meshes used for computation of three-point bending on notched beam: (a) coarse; (b) fine.
(a)(b)Young’s modulus is taken asE=4.83×104 MPa. The parameter Gt on the one hand has been provided experimentally as Gt=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 Gt=0.029 N/mm.Figure7 illustrates the crack evolution process predicted from the phase-field model. The 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 Gt=0.0290 N/mm is faster than the other one, and the corresponding compressive-shear critical energy release rates are Gs=0.290 N/mm and Gs=0.451 N/mm, respectively.Figure 7
Predicted crack paths of the three-point bending test. (a)Gt=0.0451 N/mm. (b) Gt=0.029 N/mm.
(a)(b)Figure8 compares the load-displacement response from this study to the experimental data. There is no detectable mesh dependency for the two meshes under consideration, and the beam is subject to noncyclic loading. In the case of Gt=0.0290 N/mm and Gs=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 Gt=0.0451 N/mm and Gs=0.451 N/mm, the peak load is about 13 kN which is 8.3% higher than the experimental result.Figure 8
Calculated load-displacement curves of the three-point bending test.The simulated forces both drop much faster than the experimental data after reaching the bearing capacity. These 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 load-displacement curves from both cyclic and noncyclic loadings.


## 4.3. L-Shaped Panel Test

In this section, we simulate crack propagation in an L-shaped slab. The geometry and boundary conditions of the problem are depicted in Figure9(a). The 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). The material parameters are chosen as follows [50]: Young’s modulus E=2.0×104 MPa, Poisson’s ratio v=0.18, ℓ=10 mm, tension fracture energy Gt=1.3×10−4 kN/mm, and compressive-shear fracture energy Gs=1.3×10−3 kN/mm [51]. The simulation is led with the step increment Δu=2×10−3 mm.Figure 9
L-shaped panel test. (a) Geometry and boundary conditions. (b) Experimentally observed crack pattern.
(a)(b)The computational domain is discretized using a total of 51,766 elements with fine meshes assigned to the critical zone. Figure10 illustrates the crack progression at several loading stages. The corresponding load-displacement curves are presented in Figure 11. For ℓ=10mm, the crack path from the simulation is located at the path observed in the real test. The corresponding load-displacement curves are shown in Figure 11(a); the global response matches the experimental peal loads. The 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 (ℓ=20mm). For ℓ=20mm, a diffusive crack path is obtained due to the relatively large length scale. Although the maximum force obtained by using ℓ=20mm is even lower as shown in Figure 11(b), the predicted crack path agrees well with the experimental data (Figure 10(b)). The proposed model is fairly in good agreement with the experimental results in terms of crack pattern.Figure 10
L-shaped panel test. Damage profiles for various length-scale parameters at displacementsu = 0.22 mm, 0.30 mm, 0.45 mm, and 1.0 mm. (a) ℓ=10mm. (b) ℓ=20mm.
(a)(b)Figure 11
L-shaped panel: curves of load versus displacement for various length-scale parametersℓ. (a) Length-scale parameter ℓ=10mm. (b) Length-scale parameter ℓ=20mm.
(a)(b)


## 4.4. 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. Figure12 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. The length and the width for the flaws are 12.7 mm and 0.1 mm, respectively. The 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.002mm is prescribed in line with laboratory tests [35].Figure 12
Geometries of the model rock specimen.The mechanical properties of the material areE=5×103 MPa, Poisson’s ratio v=0.24, ℓ=0.5942 mm, tension fracture energy Gt=5.0×10−6 kN/mm, and compressive-shear fracture energy Gs=1.0×10−2 kN/mm. The 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. The fracturing process is explained in more detail in Figure 13. Initially a stable growth of four wing cracks is seen (Figure 13(b)). The 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)). Then, 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].Figure 13
Crack evolution processes with displacement at (a)u = 0.312 mm; (b) u = 0.632 mm; (c) u = 0.752 mm; (d) u = 0.856 mm; (e) u = 0.912 mm.
(a)(b)(c)(d)(e)


## 5. Conclusions

This 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. The 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. The 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.



---

*Source: 4313755-2021-07-30.xml*