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 of , 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 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 is Pa, and the maximum horizontal stress is Pa, which occurs at the waist of the tunnel on the horizontal axis. When the crack is considered, the maximum vertical stress is Pa, and the maximum horizontal stress is Pa. It shows that the stress at the dome increases greatly. The vertical stress reaches Pa, and the horizontal stress is up to 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.
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: where is the order parameter; is the interpolation function, ; is strain energy density; is the critical strain energy density; is the two potential well function, ; and 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: where is the elastic energy density; is the strain tensor; and 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 , it means that the material is intact. When , it means that the material is completely fractured, as shown in Figure 1. The crack surface density can be expressed as where is the characteristic length parameter and the time length parameter determines the diffusion degree of the crack. When , the characterization is closer to the tip crack. However, the size of does not represent the actual width of crack diffusion.
Then, the total surface energy of cracks in the elastic body can be expressed as where 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 where is the elastic strain energy stored in the material and is the surface strain energy related to the crack.
The elastic strain energy can be expressed by the elastic strain as where is the displacement of material; is the fracture phase field variable; energy storage function of material per unit volume can be expressed as ; is the elastic strain energy density; can be expressed as , and is a very small constant.
Based on the above formula, the total energy function of materials can be expressed as
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 can be expressed as where is the body force per unit volume unit and is the boundary surface force per unit area.
The change of internal potential energy increment is
Combined with Equation (7), the internal potential energy increment can be expressed as where is Cauchy stress tensor,
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.
Combining Equations (8) with Equation (10), Dirichlet boundary conditions are considered, and Gauss theorem is used to obtain where n is the unit vector perpendicular to the plane.
Considering the Neumann boundary conditions (on ) and (on ), the strong form of the governing equation of the phase field method fracture model is
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.
Using the Voigt-notation method for discretization, displacement field variables and phase field variables on cell nodes can be expressed as where is the number of nodes on a cell; is the corresponding shape function on the node , and the shape function matrix in two-dimensional case is
Accordingly, the derivative of the form function can be discretized as where .
Strain-displacement matrix and phase field matrix are
Then, the gradient of displacement variable and phase field variable can be expressed as
To ensure that for any and , and, therefore, to ensure balance, the residual of its equation can be expressed as
In order to get the residual close to 0, Newton-Raphson method is used for incremental iteration.
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 (), a historical state variable is introduced to realize the irreversible process.
Therefore, after considering the historical state variables, the residual corresponding to phase field crack evolution is
The residual of the displacement field is
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
By using the Newton-Raphson method for iteration, the Equation (23) can be expressed as where each stiffness matrix is
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 Figure 2. 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.
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 Figure 3. 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 . Poisson’s ratio is . The critical energy release rate is .
In this paper, the crack length parameters were, respectively, selected as and . 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.
(a)
(b)
(c)
(a)
(b)
(c)
Figure 6 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 , while the maximum reaction force decreases. The result verifies the correctness of the phase field code.
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 Figure 7. Young’s modulus , Poisson’s ratio , and critical energy release rate 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.
Characteristic length . 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.
(a)
(b)
(c)
(d)
The variation of vertical stress with crack growth in this example is shown in Figure 9. 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 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 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 MPa.
(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 Figure 10. The material parameters are shown in Table 1. Material critical energy release rate is . Controlling the parameters of diffuse crack width is .
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, and . The left and right boundaries of the model are fixed with normal displacement, . 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 method
Considering the crack growth, the crack growth is shown in Figure 11. 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.
The horizontal stress cloud diagram after excavation is shown in Figure 12. 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 Pa, which occurs near the horizontal axis of the tunnel waist. When considering the crack, the maximum stress on the horizontal axis is Pa, and it is found that the stress at the tunnel dome increases significantly, up to Pa. It indicates that the crack development has a great influence on the stress.
(a)
(b)
The vertical stress cloud is shown in Figure 13. In situ stress is disturbed after excavation. When do not consider the crack, the maximum stress can reach 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 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 Pa.
(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 Figures 14 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 MPa when the crack is considered. The vertical stress reaches MPa when the crack is not considered. The horizontal stress is MPa when the crack is considered and MPa when the crack is not considered.
(a)
(b)
(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 was Pa. It was found that the stress at the tunnel dome increased greatly, up to 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
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
There is no conflict of interest regarding the publication of this paper.
Acknowledgments
The presented work has been financially supported by the National Key Research and Development Program of China (Grant No. 2017YFC1503104).