# Research on Crack Propagation of Deep Geologic Mass Disturbed by Excavation Based on Phase Field Method

Authors: Ningdong Chang; Jinan Wang; Fei Li
Journal: Geofluids
 (2022)

---

## Abstract

In recent years, the phase field fracture model has been widely studied and applied. It has good convergence in crack propagation simulation. Comparing with other methods, the phase field method has advantages in simulating crack intersection, bifurcation, and three-dimensional propagation. Based on the phase field method, the influence of excavation disturbance on crack initiation of rock mass is realized in this paper. The phase field fracture variational model is built by using user-defined element interface (UEL) and user material subroutine (UMAT) in ABAQUS. Firstly, the prefabricated crack propagation simulation is carried out to verify the algorithm. The fracture initiates in a butterfly shape and then expands along the horizontal direction. The results show that the maximum support reaction decreases with the gradual increase ofl, which is compared with the results obtained by Miehe et al. The result proved the correctness and reliability of the algorithm. In this paper, the phase field fracture model of a flat plate with a reserved small hole under the upper tension is established. The results show that the crack finally produces a crack in the lower left and upper right directions of the square hole and continues to extend to the model boundary, which proves the feasibility of crack independent initiation and propagation by the phase field method. The stress formed a butterfly region until the fracture occurs. And the butterfly stress distribution was still present at the end of crack propagation. The maximum vertical stress was 1.7×103 MPa. Based on the South-to-North Water Transfer Project, the simulation of tunnel crack propagation under excavation disturbance is realized for the first time, which is based on the phase field method. The results show that the influence area of excavation disturbance will increase after considering crack development. Comparing the simulation results without considering crack propagation with the simulation results considering crack propagation, it is found that the stress level in the excavation disturbance area around the tunnel is greatly affected by cracks. When the crack is not considered, the maximum vertical stress is2.16×105 Pa, and the maximum horizontal stress is 9.35×105 Pa, which occurs at the waist of the tunnel on the horizontal axis. When the crack is considered, the maximum vertical stress is 2.53×105 Pa, and the maximum horizontal stress is 1.10×106 Pa. It shows that the stress at the dome increases greatly. The vertical stress reaches 3.68×105 Pa, and the horizontal stress is up to 3.07×103 Pa. For the rock mass far away from the excavation disturbance area, because part of the elastic strain energy is absorbed by the surface crack, the stress level considering the crack is lower than that without the crack. But it is basically similar, indicating the accuracy of the phase field fracture model. This paper realizes the simulation of crack propagation under excavation disturbance and provides a way for the application of phase field fracture model in rock mechanics. This paper proves that phase field method has broad prospects in simulating rock crack propagation and provides the possibility for the popularization of phase field method.


---

## Body

## 1. Introduction

In rock engineering, crack initiation and development are the main failure modes of rock. The crack development eventually leads to rock fracture [1, 2]. The crack development also leads to local stress disturbance. Therefore, numerical simulation of fracture is of great significance in the field of geotechnical engineering.At present, finite element methods for crack propagation can be divided into two categories: One is geometric crack description method based on grid expansion, such as element deletion method [3] and interface element method. The element deletion method is a relatively simple method to describe the crack by discontinuous means. In this method, the stress of the element becomes 0 when it reaches the fracture critical [4]. The interface element method [5, 6] inserts cohesive interface elements between elements. When reached the critical condition of fracture, the interface element will be disconnected. The interface element method only allows the crack to expand on the boundary of the element [7]. This type of approach allows Raven to extend only at the boundaries of the cell. The other is the nongeometric description method, which is the diffusion crack model, such as the smooth particle hydrodynamics (SPH), the diffuse element method (DFM), the extended finite element method (XFEM), and phase field methods. Lucy and Gingold first applied the smooth particle hydrodynamics in the field of astrophysics [8, 9]. Nayrole et al. proposed the scatterer method in 1992 [10]. This method is collectively known as the meshless method. Although the meshless method does not need element structure, its crack propagation needs manual setting, and its computational efficiency is very low. Extended finite element method (XFEM) [11, 12] combined displacement field related parameters with variational function and test function by using the shape function. The fracture has nothing to do with the grid by using XFEM. The expansion of cracks can be achieved within the grid, but its calculation is limited to a few relatively simple bifurcation problems of crack. It is difficult to solve the three-dimensional model of crack propagation simulation [13, 14]. The phase field method [15] used the dispersed phase boundary to approximate the actual sharp boundary, which implemented a model describing fracture with continuous functions. By using phase field variables, the phase field method can explicitly track crack, instead of tracking the surface of crack. The propagation path of crack can be obtained by the evolution of order parameters. So, the phase field method of crack propagation simulation is not affected by meshing. Moreover, it has good convergence [16] and can also realize crack propagation in the grid.In recent years, phase field model has been applied to fracture field. It has attracted extensive attention in this field. Emilio [17] realized the operation of phase field fracture model in user-defined unit interface (UEL) and user material subroutine (UMAT) of ABAQUS, but the model is simple. Hofacker [18, 19] established the evolution of complex crack mode for dynamic problems and applied interlacing algorithm to phase field simulation. He provided the basis for solving dynamic problems. Park [20] realized the refinement model of adaptive grid and carried out the simulation of cohesive cracks under dynamic loading. Liu Guowei [21] used ABAQUS platform to realize the step algorithm of phase field fracture model and analyzed the problem of two parallel air foil crack intersection. Cao Yakuo [22] studied the fracture process of elastic-plastic plate with holes at different spatial distances based on the elastic-plastic fracture theory of metal. It provided a thought for solving the elastic-plastic problem. Liu Jia [23] established a phase field model of hydraulic coupling based on porous elasticity theory and energy minimization theory. It solves the problem of phase field seepage.The phase field fracture method describes the physical process of fracture through a series of differential equations, so as to avoid the tedious crack surface tracking. It has great advantages in simulating crack initiation, propagation, and bifurcation. At present, most of the phase field method simulation is a two-dimensional simple model of single homogeneous material, and there is no correlation analysis of excavation disturbance. Based on the phase field method, this paper established a composite stratum model and realized the excavation simulation by element deletion.In this paper, the user-defined element (UEL) of finite element software ABAQUS is used to calculate the propagation of prefabricated crack and the crack propagation of perforated plate during tensile, which verified the correctness of the code. The deep rock mass disturbed by excavation will crack on the excavation boundary, which will lead to local stress concentration of rock mass. Finally, rock mass will break and be instability. Therefore, it is of great significance to study the stress distribution of deep rock mass disturbed by excavation.


## 2. Phase Field Method Fracture Model

### 2.1. Fracture Variational Theory

The phase field method is based on Ginzburg-Laudau theory. Hakim and Karma proposed the phase field model [15], which defines the free energy of the fracture system as follows:
(1)Fu,s=∫gϕψ0εu−ψc+Vϕ+12Dϕ∇ϕ2dx,where ϕ is the order parameter; gϕ is the interpolation function, gϕ=ϕ2+α; ψ0 is strain energy density; ψc is the critical strain energy density; Vϕ is the two potential well function, Vϕ=ϕ21−ϕ2/4; and Dϕ is surface parameter. The first two terms are the body free energy term and the last term is the surface energy term.In 1920, Griffith proposed Griffith energy theory based on energy [24], but this theory is suitable for simple model. It is difficult to analyze the complex crack model. In 1998, G. A. Francfor and J. J. Marigo proposed a variational method for fracture problems, which was based on the Griffith energy theory [25]. The total potential energy Π of elastomer Π can be divided into elastic energy of elastomer and surface energy of fracture:
(2)∏u,Γ=∫ΩψeεudΩ+∫ΓGcdΓ,where ψe is the elastic energy density; εu is the strain tensor; and Gc is the critical energy release rate.In the fracture variational model proposed by Francfor and Marigo, the crack propagation is controlled by the principle of minimum free energy. However, the boundary integral of surface energy in this method is not easy to deal with in the case of unknown crack boundary. In 2000, Bourdin and Francfort introduced the order parameter [26] and defined one in the interval [0, 1]. When ϕ=0, it means that the material is intact. When ϕ=1, it means that the material is completely fractured, as shown in Figure 1. The crack surface density can be expressed as
(3)γϕ,∇ϕ=12lϕ2+l2∇ϕ2,where lis the characteristic length parameter and the time length parameter determines the diffusion degree of the crack. When l⟶0, the characterization is closer to the tip crack. However, the size of L does not represent the actual width of crack diffusion.Figure 1
Crack described with a phase field mode.Then, the total surface energy of cracks in the elastic body can be expressed as(4)∫ΓGcdΓ≈∫ΩGcϕ−122l+l2∇ϕ2dΩ,where Gc is the critical value of Griffith energy release rate of material per unit area.To analyze the fracture phase field, the total energy function of the material can be expressed as(5)Ψu,ϕ=Ψbu,ϕ+Ψsϕ,where Ψbu,ϕ is the elastic strain energy stored in the material and Ψsϕ is the surface strain energy related to the crack.The elastic strain energy can be expressed by the elastic strain as(6)Ψbu,ϕ=∫Ωψεu,ϕdV,where u is the displacement of material; ϕ is the fracture phase field variable; energy storage function ψof material per unit volume can be expressed as ψε,ϕ=gϕψ0ε; ψ0 is the elastic strain energy density; gϕ can be expressed as gϕ=1−ϕ2+k, and k is a very small constant.Based on the above formula, the total energy function of materials can be expressed as(7)Ψu,ϕ=∫Ω1−ϕ2+kψ0ε+Gc12lϕ2+l2∇ϕ2dV,


### 2.2. Fracture Governing Equation by Phase Field Method

The energy of the system can be divided into external work generated by the application of load and internal work generated by deformation. The external energy incrementδWext can be expressed as
(8)δWext=∫VbδudV+∫δVhδu∂V,where b is the body force per unit volume unit and h is the boundary surface force per unit area.The change of internal potential energy increment is(9)∂Wint=∂ψ=∂ψ∂εδε+∂ψ∂ϕδϕ.Combined with Equation (7), the internal potential energy increment can be expressed as
(10)∂Wint=∫Ωσδε−21−ϕδϕψε+Gc1lϕδϕ+l∇ϕ⋅∇δϕdV,where σ is Cauchy stress tensor,
(11)σ=∂ψ∂ε=gϕσ0=1−ϕ2+kσ0.Assuming that the model is a quasi-static process, the virtual work done by the external force of the structure is equal to the virtual work done by the internal force of the structure.(12)∂Wint−∂Wext=0.Combining Equations (8) with Equation (10), Dirichlet boundary conditions are considered, and Gauss theorem is used to obtain
(13)∫∂ΩGcl∇ϕ⋅n⋅δϕdA+∫Ω−Divσ+b⋅δu−21−ϕψε−Gc1lϕ+Divl∇ϕδϕdV+∫∂Ωhσ⋅n−h⋅δudA=0,where n is the unit vector perpendicular to the plane.Considering the Neumann boundary conditionsσ⋅n=h (on ∂Ωh) and ∇ϕ⋅n=0 (on ∂Ω), the strong form of the governing equation of the phase field method fracture model is
(14)Divσ+b=0,Gc1lϕ+l∇ϕ−21−ϕψε=0.


## 2.1. Fracture Variational Theory

The phase field method is based on Ginzburg-Laudau theory. Hakim and Karma proposed the phase field model [15], which defines the free energy of the fracture system as follows:
(1)Fu,s=∫gϕψ0εu−ψc+Vϕ+12Dϕ∇ϕ2dx,where ϕ is the order parameter; gϕ is the interpolation function, gϕ=ϕ2+α; ψ0 is strain energy density; ψc is the critical strain energy density; Vϕ is the two potential well function, Vϕ=ϕ21−ϕ2/4; and Dϕ is surface parameter. The first two terms are the body free energy term and the last term is the surface energy term.In 1920, Griffith proposed Griffith energy theory based on energy [24], but this theory is suitable for simple model. It is difficult to analyze the complex crack model. In 1998, G. A. Francfor and J. J. Marigo proposed a variational method for fracture problems, which was based on the Griffith energy theory [25]. The total potential energy Π of elastomer Π can be divided into elastic energy of elastomer and surface energy of fracture:
(2)∏u,Γ=∫ΩψeεudΩ+∫ΓGcdΓ,where ψe is the elastic energy density; εu is the strain tensor; and Gc is the critical energy release rate.In the fracture variational model proposed by Francfor and Marigo, the crack propagation is controlled by the principle of minimum free energy. However, the boundary integral of surface energy in this method is not easy to deal with in the case of unknown crack boundary. In 2000, Bourdin and Francfort introduced the order parameter [26] and defined one in the interval [0, 1]. When ϕ=0, it means that the material is intact. When ϕ=1, it means that the material is completely fractured, as shown in Figure 1. The crack surface density can be expressed as
(3)γϕ,∇ϕ=12lϕ2+l2∇ϕ2,where lis the characteristic length parameter and the time length parameter determines the diffusion degree of the crack. When l⟶0, the characterization is closer to the tip crack. However, the size of L does not represent the actual width of crack diffusion.Figure 1
Crack described with a phase field mode.Then, the total surface energy of cracks in the elastic body can be expressed as(4)∫ΓGcdΓ≈∫ΩGcϕ−122l+l2∇ϕ2dΩ,where Gc is the critical value of Griffith energy release rate of material per unit area.To analyze the fracture phase field, the total energy function of the material can be expressed as(5)Ψu,ϕ=Ψbu,ϕ+Ψsϕ,where Ψbu,ϕ is the elastic strain energy stored in the material and Ψsϕ is the surface strain energy related to the crack.The elastic strain energy can be expressed by the elastic strain as(6)Ψbu,ϕ=∫Ωψεu,ϕdV,where u is the displacement of material; ϕ is the fracture phase field variable; energy storage function ψof material per unit volume can be expressed as ψε,ϕ=gϕψ0ε; ψ0 is the elastic strain energy density; gϕ can be expressed as gϕ=1−ϕ2+k, and k is a very small constant.Based on the above formula, the total energy function of materials can be expressed as(7)Ψu,ϕ=∫Ω1−ϕ2+kψ0ε+Gc12lϕ2+l2∇ϕ2dV,


## 2.2. Fracture Governing Equation by Phase Field Method

The energy of the system can be divided into external work generated by the application of load and internal work generated by deformation. The external energy incrementδWext can be expressed as
(8)δWext=∫VbδudV+∫δVhδu∂V,where b is the body force per unit volume unit and h is the boundary surface force per unit area.The change of internal potential energy increment is(9)∂Wint=∂ψ=∂ψ∂εδε+∂ψ∂ϕδϕ.Combined with Equation (7), the internal potential energy increment can be expressed as
(10)∂Wint=∫Ωσδε−21−ϕδϕψε+Gc1lϕδϕ+l∇ϕ⋅∇δϕdV,where σ is Cauchy stress tensor,
(11)σ=∂ψ∂ε=gϕσ0=1−ϕ2+kσ0.Assuming that the model is a quasi-static process, the virtual work done by the external force of the structure is equal to the virtual work done by the internal force of the structure.(12)∂Wint−∂Wext=0.Combining Equations (8) with Equation (10), Dirichlet boundary conditions are considered, and Gauss theorem is used to obtain
(13)∫∂ΩGcl∇ϕ⋅n⋅δϕdA+∫Ω−Divσ+b⋅δu−21−ϕψε−Gc1lϕ+Divl∇ϕδϕdV+∫∂Ωhσ⋅n−h⋅δudA=0,where n is the unit vector perpendicular to the plane.Considering the Neumann boundary conditionsσ⋅n=h (on ∂Ωh) and ∇ϕ⋅n=0 (on ∂Ω), the strong form of the governing equation of the phase field method fracture model is
(14)Divσ+b=0,Gc1lϕ+l∇ϕ−21−ϕψε=0.


## 3. Realize Fracture Finite Element by Using Phase Field Method

In order to solve the numerical solution of partial differential Equation (14), the weak form of the governing equation is solved by finite element method.
(15)∫Ωσδε−b⋅δudV+∫∂Ωhh⋅δudA=0,∫Ω−21−ϕδϕψε+Gc1lϕδϕ+l∇ϕ⋅∇δϕdV=0.Using the Voigt-notation method for discretization, displacement field variablesu and phase field variables ϕon cell nodes I can be expressed as
(16)u=∑i=1nNiuui,ϕ=∑i=1nNiuϕi,where n is the number of nodes on a cell; Ni is the corresponding shape function on the node I, and the shape function matrix in two-dimensional case is
(17)Niu=Ni00Ni.Accordingly, the derivative of the form function can be discretized as(18)ε=∑i=1nBiuui,∇ϕ=∑i=1nBiϕϕi,where ε=εxx,εyy,εxyT.Strain-displacement matrix and phase field matrix are(19)Biu=Ni,x00Ni,yNi,yNi,x,Biϕ=Ni,xNi,y.Then, the gradient of displacement variable and phase field variable can be expressed as(20)δε=∑i=1nBiδuui,∇ϕ=∑i=1nBiϕδϕi.To ensure that∂Wint−∂Wext=0 for any u and ϕ, and, therefore, to ensure balance, the residual of its equation can be expressed as
(21)Riu=∫Ω1−ϕ2+kBiuTσdV−∫ΩNiuTbdV−∫∂ΩhNiuThdA,Riϕ=∫Ω−21−ϕNiψε+Gc1lNiϕ+lBiϕT∇ϕdV.In order to get the residual close to 0, Newton-Raphson method is used for incremental iteration.(22)uϕt+Δt=uϕt−KuuKuϕKϕuKϕϕt−1RuRϕt.The governing equation is a partial differential equation system composed of quasi-static equilibrium equation and phase field equation. However, in order to prevent crack closure when the elastic body is under pressure or unloading springback and ensure the irreversible process of crack evolution (ϕt+Δt≥ϕt), a historical state variable H is introduced to realize the irreversible process.
(23)H=ψε,ψε>HtHt,ψε≤Ht.Therefore, after considering the historical state variables, the residual corresponding to phase field crack evolution is(24)Riϕ=∫Ω−21−ϕNiH+Gc1lNiϕ+lBiϕT∇ϕdV.The residual of the displacement field is(25)Riu=∫Ω1−ϕ2+kBiuTσdV−∫ΩNiuTbdV−∫∂ΩhNiuThdA.The stiffness matrix gradually degenerates in the process of crack initiation and propagation, which will result in constant rearrangement of the stress field. The implicit solver cannot be used to obtain a stable equilibrium solution. Therefore, the displacement field and phase field in Equation (14) are considered as coupling fields and solved separately. Considering the historical state variables, Equation (14) can be expressed as
(26)Gc1lϕ+l∇ϕ−21−ϕH=0.By using the Newton-Raphson method for iteration, the Equation (23) can be expressed as
(27)uϕt+Δt=uϕt−Kuu00Kϕϕt−1RuRϕt,where each stiffness matrix is
(28)Kijuu=∂Riu∂uj=∫Ω1−ϕ2+kBiuTC0BjudV,Kijϕϕ=∂Riϕ∂ϕj=∫Ω2H+GclNiNj+GclBiϕTBjϕdV.In this paper, ABAQUS is used to realize the phase field method fracture simulation. This research realized the phase field method model by using UEL of ABAQUS. The relevant stiffness matrix is updated in UEL. The specific process is shown in Figure2. Due to the use of user element subroutine, the integral points defined by it cannot be visualized independently in ABAQUS postprocessing. However, the element stiffness matrix and stress-strain results on each node can be obtained. In order to realize the visualization of the results, this paper creates an auxiliary grid by using UMAT to define the relevant material parameters at each integral point of the auxiliary grid. Then, transmit the data of user unit subroutine to user material subroutine, and output SDV variables for final visualization. The stress component and stiffness matrix of the auxiliary virtual grid are zero, which have no influence on the solution result.Figure 2
The process for solving phase field model.


## 4. Numerical Simulation Examples and Verification

### 4.1. A Case of Prefabricated Crack Propagation

To verify the correctness of the code, a square plate with prefabricated cracks on the left side was created. The specific size of the plate is shown in Figure3. The lower boundary of the model was fixed. The upward displacement plate of 0.01 mm was applied to the upper boundary. The Young’s modulus E=210GPa. Poisson’s ratio is μ=0.3. The critical energy release rate is Gc=2.7×10−3kN/mm.Figure 3
Geometric model and the location of crack (mm).In this paper, the crack length parameters were, respectively, selected asl=0.015mm and l=0.0075mm. The crack propagation results are shown in Figures 4 and 5. In the figure, the blue area represented no failure, while the red area represented complete failure. The fracture initiates in a butterfly shape and then expands along the horizontal and vertical direction. Finally, the crack reaches the right side of the model.Figure 4
The propagation of crack whenl=0.015mm. (a) Crack initiation stage. (b) Crack development stage. (c) Crack penetration stage.
(a)(b)(c)Figure 5
The propagation of crack whenl=0.0075mm. (a) Crack initiation stage. (b) Crack development stage. (c) Crack penetration stage.
(a)(b)(c)Figure6 shows the comparison between the reaction force and displacement curve calculated in this paper and Miehe [27]. At that time, the maximum support reaction force calculated was 0.71 kN, and the maximum support reaction force calculated was 0.74 kN. The calculation results were basically consistent with the crack growth and support reaction displacement curve obtained by Miehe. Although there is a small deviation on the curve, the maximum reaction force is basically similar, which is consistent with the results of Miehe. It is found that with the increase of l, while the maximum reaction force decreases. The result verifies the correctness of the phase field code.Figure 6
Comparison of load-displacement curves of differentl.


### 4.2. A Case of Tensile Crack Propagation of Perforated Plate

A square plate with side length of 40 m was established. And a square hole with one side length of 4 m was established in the center of the plate, as shown in Figure7. Young’s modulus E=210GPa, Poisson’s ratio μ=0.3, and critical energy release rate Gc=2.7×10−3kN/mm were used to fix the lower boundary of the CPE4 unit model, and upward displacement was applied to the upper boundary of the model. The lower boundary of the CPE4 element model is fixed, and the upper boundary of the model is subjected to upward displacement.Figure 7
Prefabricated cavity plate model (m).Characteristic lengthl=0.04m. The crack extension results are shown in Figure 8. The results showed that symmetrical butterfly-shaped microcracks appear around the hole in the initial stage. At this time, the structure damage was small. As the tension continued to be applied, the butterfly area expanded, eventually creating cracks in the upper right and lower left corners of the square hole. And the crack continued to expand and develops to the edge of the model, which proved the feasibility of crack initiation and self-propagation.Figure 8
The result of crack propagation. (a) Crack initiation stage. (b) Crack initiation stage. (c) Crack development stage. (d) Crack penetration stage.
(a)(b)(c)(d)The variation of vertical stress with crack growth in this example is shown in Figure9. The results showed that there was a great stress concentration at the four corners of the hole at the initial stage of crack development. And the stress distribution likes the butterfly. The maximum vertical stress was 4.7×102 MPa. The stress on both sides of the hole gradually increased and gradually formed a butterfly region until the fracture occurs. The stress on both sides of the hole is positive, and the stress above and below the hole is negative. The maximum vertical stress was 1.23×103 MPa. Crack propagation begins at this time. After that, the crack propagated gradually. And the butterfly stress distribution was still present at the end of crack propagation. At this time, the maximum vertical stress was 1.7×103 MPa.Figure 9
The vertical stress when the crack grows. (a) Crack initiation stage. (b) Crack initiation stage. (c) Crack development stage.
(a)(b)(c)


### 4.3. A Case of Crack Propagation in Composite Stratum Disturbed by Excavation

#### 4.3.1. Engineering Background

As a strategic project in China, the South-to-North Water Diversion Project is an extremely large infrastructure project to solve the serious shortage of water resources in the north and implement the optimal allocation of water resources. It plays a very important role in alleviating the bearing pressure of water resources and improving the water supply security and guarantee rate in Beijing.Under the influence of excavation disturbance, some stress concentration occurs in the geological body of Eastern Canal in Beijing, which leads to lining damage and pipeline leakage. Therefore, the pile no. 21+979 and its surrounding area of Eastern Canal of the South-to-North Water Diversion Project were selected for the analysis. Crack propagation and stress disturbance were analyzed in this case. The left and right width of the numerical model near pile no. 21+979 is 60 m. The tunnel buried depth is 19.62 m. The depth below the axis of the tunnel is 22.38 m and the tunnel radius is 3 m.According to the site geological survey data, the strata were simplified, and a composite stratum tunnel model was established. The tunnel size and stratum distribution are shown in Figure10. The material parameters are shown in Table 1. Material critical energy release rate is Gc=2.7×10−3kN/mm. Controlling the parameters of diffuse crack width is l=0.04m.Figure 10
Initial model around pile No. 21+979.Table 1
Parameters of geological model around pile No. 21+979.
NumberStratumDensity (103kg/m3)Elasticity modulus (MPa)Poisson’s ratio1Silty clay1.965.940.32Fine medium sand2.0013.50.2


#### 4.3.2. Realization of Phase Field Method

In order to simulate the initial stress state of rock and soil mass, the ground stress balance was carried out first. Then, the excavation was carried out. The specific process of the excavation of rock and soil mass simulated based on the phase field method is as follows:(1)
Establish the initial composite formation model (without tunnel excavation). Fix the bottom boundary for the model,ux=0 and uy=0. The left and right boundaries of the model are fixed with normal displacement, ux=0. The initial geostress field of the model was obtained by the ground stress balance analysis step under the dead weight stress(2)
Delete some excavated units, renumber nodes and grids, and establish tunnel excavation model(3)
In the .inp file of ABAQUS, the fracture parameters related to the layer phase field method were defined and the virtual visualization grid was generated. Apply the initial in situ stress field calculated in the first step to the tunnel excavation model, and define the output variables to simulate the excavation of rock and soil mass, which considered the phase field methodConsidering the crack growth, the crack growth is shown in Figure11. After tunnel excavation, the cracks of the model are mainly distributed around the tunnel, and the weak links are on the horizontal sides of the tunnel. The crack initiation occurs firstly in horizontal direction and then gradually develops. After excavation, the stress of the tunnel is released, and the shape of tunnel is ellipse. Due to the maximum stress on both sides of the middle line of the tunnel, there is a large damage in this part. The crack begins to expand at the horizontal direction.Figure 11
Crack propagation under the disturbance of excavation.The horizontal stress cloud diagram after excavation is shown in Figure12. Figure 12(a) is the horizontal stress result obtained without considering the phase field method, that is, without crack development; Figure 12(b) is the horizontal stress result obtained with crack development by using the phase field method. Both cases indicate that stress disturbance occurs in the surrounding area after tunnel excavation. When the crack is not considered, the maximum horizontal stress is 2.16×105 Pa, which occurs near the horizontal axis of the tunnel waist. When considering the crack, the maximum stress on the horizontal axis is 2.53×105 Pa, and it is found that the stress at the tunnel dome increases significantly, up to 3.68×105 Pa. It indicates that the crack development has a great influence on the stress.Figure 12
Horizontal stress distribution after excavation. (a) No cracks. (b) Consider cracks.
(a)(b)The vertical stress cloud is shown in Figure13. In situ stress is disturbed after excavation. When do not consider the crack, the maximum stress can reach 9.35×105 Pa. It occurs at the horizontal axis of the tunnel edge. When considering the crack development, the stress level around the tunnel increases due to the influence of the crack development. The maximum stress can reach 1.10×106 Pa. At the horizontal axis of the tunnel edge, the stress at the dome changes under the influence of the crack. It produces a vertical upward stress of 3.07×103 Pa.Figure 13
Vertical stress distribution after excavation. (a) No cracks. (b) Consider cracks.
(a)(b)After the phase field method is adopted to simulate the excavation, the stress around the tunnel changes greatly from considering the cracks, which is due to the development of cracks around the tunnel. In the area far away from the tunnel, it is less affected by the crack, as shown in Figures14 and 15. However, its total stress level is reduced, which is due to part of the elastic energy into the surface energy of the crack. As the cracks developed along the horizontal axis, it can be found in Figure 14(b) and Figure 15(b) that when the distance from the center of the tunnel is 3 m (tunnel boundary), the stress level of the cracks is significantly different from that when the crack is not considered. The vertical stress reaches −7.90×105 MPa when the crack is considered. The vertical stress reaches −8.69×105 MPa when the crack is not considered. The horizontal stress is 5.76×104 MPa when the crack is considered and −1.24×105 MPa when the crack is not considered.Figure 14
Comparison of vertical stress after excavation. (a) Above the vertical axis of the tunnel. (b) Right of horizontal axis of tunnel.
(a)(b)Figure 15
Comparison of horizontal stress after excavation. (a) Above the vertical axis of the tunnel. (b) Right of horizontal axis of tunnel.
(a)(b)


## 4.1. A Case of Prefabricated Crack Propagation

To verify the correctness of the code, a square plate with prefabricated cracks on the left side was created. The specific size of the plate is shown in Figure3. The lower boundary of the model was fixed. The upward displacement plate of 0.01 mm was applied to the upper boundary. The Young’s modulus E=210GPa. Poisson’s ratio is μ=0.3. The critical energy release rate is Gc=2.7×10−3kN/mm.Figure 3
Geometric model and the location of crack (mm).In this paper, the crack length parameters were, respectively, selected asl=0.015mm and l=0.0075mm. The crack propagation results are shown in Figures 4 and 5. In the figure, the blue area represented no failure, while the red area represented complete failure. The fracture initiates in a butterfly shape and then expands along the horizontal and vertical direction. Finally, the crack reaches the right side of the model.Figure 4
The propagation of crack whenl=0.015mm. (a) Crack initiation stage. (b) Crack development stage. (c) Crack penetration stage.
(a)(b)(c)Figure 5
The propagation of crack whenl=0.0075mm. (a) Crack initiation stage. (b) Crack development stage. (c) Crack penetration stage.
(a)(b)(c)Figure6 shows the comparison between the reaction force and displacement curve calculated in this paper and Miehe [27]. At that time, the maximum support reaction force calculated was 0.71 kN, and the maximum support reaction force calculated was 0.74 kN. The calculation results were basically consistent with the crack growth and support reaction displacement curve obtained by Miehe. Although there is a small deviation on the curve, the maximum reaction force is basically similar, which is consistent with the results of Miehe. It is found that with the increase of l, while the maximum reaction force decreases. The result verifies the correctness of the phase field code.Figure 6
Comparison of load-displacement curves of differentl.


## 4.2. A Case of Tensile Crack Propagation of Perforated Plate

A square plate with side length of 40 m was established. And a square hole with one side length of 4 m was established in the center of the plate, as shown in Figure7. Young’s modulus E=210GPa, Poisson’s ratio μ=0.3, and critical energy release rate Gc=2.7×10−3kN/mm were used to fix the lower boundary of the CPE4 unit model, and upward displacement was applied to the upper boundary of the model. The lower boundary of the CPE4 element model is fixed, and the upper boundary of the model is subjected to upward displacement.Figure 7
Prefabricated cavity plate model (m).Characteristic lengthl=0.04m. The crack extension results are shown in Figure 8. The results showed that symmetrical butterfly-shaped microcracks appear around the hole in the initial stage. At this time, the structure damage was small. As the tension continued to be applied, the butterfly area expanded, eventually creating cracks in the upper right and lower left corners of the square hole. And the crack continued to expand and develops to the edge of the model, which proved the feasibility of crack initiation and self-propagation.Figure 8
The result of crack propagation. (a) Crack initiation stage. (b) Crack initiation stage. (c) Crack development stage. (d) Crack penetration stage.
(a)(b)(c)(d)The variation of vertical stress with crack growth in this example is shown in Figure9. The results showed that there was a great stress concentration at the four corners of the hole at the initial stage of crack development. And the stress distribution likes the butterfly. The maximum vertical stress was 4.7×102 MPa. The stress on both sides of the hole gradually increased and gradually formed a butterfly region until the fracture occurs. The stress on both sides of the hole is positive, and the stress above and below the hole is negative. The maximum vertical stress was 1.23×103 MPa. Crack propagation begins at this time. After that, the crack propagated gradually. And the butterfly stress distribution was still present at the end of crack propagation. At this time, the maximum vertical stress was 1.7×103 MPa.Figure 9
The vertical stress when the crack grows. (a) Crack initiation stage. (b) Crack initiation stage. (c) Crack development stage.
(a)(b)(c)


## 4.3. A Case of Crack Propagation in Composite Stratum Disturbed by Excavation

### 4.3.1. Engineering Background

As a strategic project in China, the South-to-North Water Diversion Project is an extremely large infrastructure project to solve the serious shortage of water resources in the north and implement the optimal allocation of water resources. It plays a very important role in alleviating the bearing pressure of water resources and improving the water supply security and guarantee rate in Beijing.Under the influence of excavation disturbance, some stress concentration occurs in the geological body of Eastern Canal in Beijing, which leads to lining damage and pipeline leakage. Therefore, the pile no. 21+979 and its surrounding area of Eastern Canal of the South-to-North Water Diversion Project were selected for the analysis. Crack propagation and stress disturbance were analyzed in this case. The left and right width of the numerical model near pile no. 21+979 is 60 m. The tunnel buried depth is 19.62 m. The depth below the axis of the tunnel is 22.38 m and the tunnel radius is 3 m.According to the site geological survey data, the strata were simplified, and a composite stratum tunnel model was established. The tunnel size and stratum distribution are shown in Figure10. The material parameters are shown in Table 1. Material critical energy release rate is Gc=2.7×10−3kN/mm. Controlling the parameters of diffuse crack width is l=0.04m.Figure 10
Initial model around pile No. 21+979.Table 1
Parameters of geological model around pile No. 21+979.
NumberStratumDensity (103kg/m3)Elasticity modulus (MPa)Poisson’s ratio1Silty clay1.965.940.32Fine medium sand2.0013.50.2


### 4.3.2. Realization of Phase Field Method

In order to simulate the initial stress state of rock and soil mass, the ground stress balance was carried out first. Then, the excavation was carried out. The specific process of the excavation of rock and soil mass simulated based on the phase field method is as follows:(1)
Establish the initial composite formation model (without tunnel excavation). Fix the bottom boundary for the model,ux=0 and uy=0. The left and right boundaries of the model are fixed with normal displacement, ux=0. The initial geostress field of the model was obtained by the ground stress balance analysis step under the dead weight stress(2)
Delete some excavated units, renumber nodes and grids, and establish tunnel excavation model(3)
In the .inp file of ABAQUS, the fracture parameters related to the layer phase field method were defined and the virtual visualization grid was generated. Apply the initial in situ stress field calculated in the first step to the tunnel excavation model, and define the output variables to simulate the excavation of rock and soil mass, which considered the phase field methodConsidering the crack growth, the crack growth is shown in Figure11. After tunnel excavation, the cracks of the model are mainly distributed around the tunnel, and the weak links are on the horizontal sides of the tunnel. The crack initiation occurs firstly in horizontal direction and then gradually develops. After excavation, the stress of the tunnel is released, and the shape of tunnel is ellipse. Due to the maximum stress on both sides of the middle line of the tunnel, there is a large damage in this part. The crack begins to expand at the horizontal direction.Figure 11
Crack propagation under the disturbance of excavation.The horizontal stress cloud diagram after excavation is shown in Figure12. Figure 12(a) is the horizontal stress result obtained without considering the phase field method, that is, without crack development; Figure 12(b) is the horizontal stress result obtained with crack development by using the phase field method. Both cases indicate that stress disturbance occurs in the surrounding area after tunnel excavation. When the crack is not considered, the maximum horizontal stress is 2.16×105 Pa, which occurs near the horizontal axis of the tunnel waist. When considering the crack, the maximum stress on the horizontal axis is 2.53×105 Pa, and it is found that the stress at the tunnel dome increases significantly, up to 3.68×105 Pa. It indicates that the crack development has a great influence on the stress.Figure 12
Horizontal stress distribution after excavation. (a) No cracks. (b) Consider cracks.
(a)(b)The vertical stress cloud is shown in Figure13. In situ stress is disturbed after excavation. When do not consider the crack, the maximum stress can reach 9.35×105 Pa. It occurs at the horizontal axis of the tunnel edge. When considering the crack development, the stress level around the tunnel increases due to the influence of the crack development. The maximum stress can reach 1.10×106 Pa. At the horizontal axis of the tunnel edge, the stress at the dome changes under the influence of the crack. It produces a vertical upward stress of 3.07×103 Pa.Figure 13
Vertical stress distribution after excavation. (a) No cracks. (b) Consider cracks.
(a)(b)After the phase field method is adopted to simulate the excavation, the stress around the tunnel changes greatly from considering the cracks, which is due to the development of cracks around the tunnel. In the area far away from the tunnel, it is less affected by the crack, as shown in Figures14 and 15. However, its total stress level is reduced, which is due to part of the elastic energy into the surface energy of the crack. As the cracks developed along the horizontal axis, it can be found in Figure 14(b) and Figure 15(b) that when the distance from the center of the tunnel is 3 m (tunnel boundary), the stress level of the cracks is significantly different from that when the crack is not considered. The vertical stress reaches −7.90×105 MPa when the crack is considered. The vertical stress reaches −8.69×105 MPa when the crack is not considered. The horizontal stress is 5.76×104 MPa when the crack is considered and −1.24×105 MPa when the crack is not considered.Figure 14
Comparison of vertical stress after excavation. (a) Above the vertical axis of the tunnel. (b) Right of horizontal axis of tunnel.
(a)(b)Figure 15
Comparison of horizontal stress after excavation. (a) Above the vertical axis of the tunnel. (b) Right of horizontal axis of tunnel.
(a)(b)


## 4.3.1. Engineering Background

As a strategic project in China, the South-to-North Water Diversion Project is an extremely large infrastructure project to solve the serious shortage of water resources in the north and implement the optimal allocation of water resources. It plays a very important role in alleviating the bearing pressure of water resources and improving the water supply security and guarantee rate in Beijing.Under the influence of excavation disturbance, some stress concentration occurs in the geological body of Eastern Canal in Beijing, which leads to lining damage and pipeline leakage. Therefore, the pile no. 21+979 and its surrounding area of Eastern Canal of the South-to-North Water Diversion Project were selected for the analysis. Crack propagation and stress disturbance were analyzed in this case. The left and right width of the numerical model near pile no. 21+979 is 60 m. The tunnel buried depth is 19.62 m. The depth below the axis of the tunnel is 22.38 m and the tunnel radius is 3 m.According to the site geological survey data, the strata were simplified, and a composite stratum tunnel model was established. The tunnel size and stratum distribution are shown in Figure10. The material parameters are shown in Table 1. Material critical energy release rate is Gc=2.7×10−3kN/mm. Controlling the parameters of diffuse crack width is l=0.04m.Figure 10
Initial model around pile No. 21+979.Table 1
Parameters of geological model around pile No. 21+979.
NumberStratumDensity (103kg/m3)Elasticity modulus (MPa)Poisson’s ratio1Silty clay1.965.940.32Fine medium sand2.0013.50.2


## 4.3.2. Realization of Phase Field Method

In order to simulate the initial stress state of rock and soil mass, the ground stress balance was carried out first. Then, the excavation was carried out. The specific process of the excavation of rock and soil mass simulated based on the phase field method is as follows:(1)
Establish the initial composite formation model (without tunnel excavation). Fix the bottom boundary for the model,ux=0 and uy=0. The left and right boundaries of the model are fixed with normal displacement, ux=0. The initial geostress field of the model was obtained by the ground stress balance analysis step under the dead weight stress(2)
Delete some excavated units, renumber nodes and grids, and establish tunnel excavation model(3)
In the .inp file of ABAQUS, the fracture parameters related to the layer phase field method were defined and the virtual visualization grid was generated. Apply the initial in situ stress field calculated in the first step to the tunnel excavation model, and define the output variables to simulate the excavation of rock and soil mass, which considered the phase field methodConsidering the crack growth, the crack growth is shown in Figure11. After tunnel excavation, the cracks of the model are mainly distributed around the tunnel, and the weak links are on the horizontal sides of the tunnel. The crack initiation occurs firstly in horizontal direction and then gradually develops. After excavation, the stress of the tunnel is released, and the shape of tunnel is ellipse. Due to the maximum stress on both sides of the middle line of the tunnel, there is a large damage in this part. The crack begins to expand at the horizontal direction.Figure 11
Crack propagation under the disturbance of excavation.The horizontal stress cloud diagram after excavation is shown in Figure12. Figure 12(a) is the horizontal stress result obtained without considering the phase field method, that is, without crack development; Figure 12(b) is the horizontal stress result obtained with crack development by using the phase field method. Both cases indicate that stress disturbance occurs in the surrounding area after tunnel excavation. When the crack is not considered, the maximum horizontal stress is 2.16×105 Pa, which occurs near the horizontal axis of the tunnel waist. When considering the crack, the maximum stress on the horizontal axis is 2.53×105 Pa, and it is found that the stress at the tunnel dome increases significantly, up to 3.68×105 Pa. It indicates that the crack development has a great influence on the stress.Figure 12
Horizontal stress distribution after excavation. (a) No cracks. (b) Consider cracks.
(a)(b)The vertical stress cloud is shown in Figure13. In situ stress is disturbed after excavation. When do not consider the crack, the maximum stress can reach 9.35×105 Pa. It occurs at the horizontal axis of the tunnel edge. When considering the crack development, the stress level around the tunnel increases due to the influence of the crack development. The maximum stress can reach 1.10×106 Pa. At the horizontal axis of the tunnel edge, the stress at the dome changes under the influence of the crack. It produces a vertical upward stress of 3.07×103 Pa.Figure 13
Vertical stress distribution after excavation. (a) No cracks. (b) Consider cracks.
(a)(b)After the phase field method is adopted to simulate the excavation, the stress around the tunnel changes greatly from considering the cracks, which is due to the development of cracks around the tunnel. In the area far away from the tunnel, it is less affected by the crack, as shown in Figures14 and 15. However, its total stress level is reduced, which is due to part of the elastic energy into the surface energy of the crack. As the cracks developed along the horizontal axis, it can be found in Figure 14(b) and Figure 15(b) that when the distance from the center of the tunnel is 3 m (tunnel boundary), the stress level of the cracks is significantly different from that when the crack is not considered. The vertical stress reaches −7.90×105 MPa when the crack is considered. The vertical stress reaches −8.69×105 MPa when the crack is not considered. The horizontal stress is 5.76×104 MPa when the crack is considered and −1.24×105 MPa when the crack is not considered.Figure 14
Comparison of vertical stress after excavation. (a) Above the vertical axis of the tunnel. (b) Right of horizontal axis of tunnel.
(a)(b)Figure 15
Comparison of horizontal stress after excavation. (a) Above the vertical axis of the tunnel. (b) Right of horizontal axis of tunnel.
(a)(b)


## 5. Conclusion

Based on ABAQUS platform, this paper realized the phase field method fracture model and simulated the crack growth of composite stratum under excavation disturbance. This paper realized the excavation simulation of composite stratum based on the phase field method. The main conclusions are as follows:(1)
By using UMAT/UEL of ABAQUS secondary development interface, prefabricated crack growth simulation, prefabricated hole crack growth simulation, and composite stratum crack growth simulation under excavation disturbance were carried out. Based on the previous phase field model, the simulation of various materials and the disturbance analysis of excavation to the model are realized in this paper(2)
The correctness and reliability of the phase field fracture model code were verified by comparing the results of prefabricated crack propagation simulation with those in other literature. It was found that the crack propagated in the upper right corner and the lower left corner of the hole through the tensile example of the perforated plate, which demonstrated the feasibility of the phase field method for the self-initiation and propagation of cracks. By comparing the crack growth model under excavation disturbance with the model without crack growth, it was found that the stress disturbance was great at the crack development. When considering the crack, the maximum stress on the horizontal axis was2.53×105 Pa. It was found that the stress at the tunnel dome increased greatly, up to 3.68×105 Pa(3)
Based on the crack propagation model of composite strata established by phase field method, the crack distribution of rock mass after excavation disturbance and the fine in situ stress distribution around the tunnel are obtained. It provides an idea for simulating crack propagation in geotechnical engineering and lays a foundation for simulating in situ stress analysis in geotechnical excavation engineering



---

*Source: 5791006-2022-05-09.xml*