Numerical Simulation on the Progressive Failure Processes of Foundation Pit Excavation Based on a New Particle Failure Method

The predictions of failure zone during the foundation excavations will provide important guidance for the safety constructions of engineering structures. Based on this background, the smoothing kernel function in the traditional SPH method has been improved. The failure mark η is introduced into the program to realize the failure characteristics of particles at meso – scale. The “ Killing Particle Method ” has also been proposed, which can realize the simulations of complex excavation processes. The whole progressive failure processes of the excavation of a foundation pit are numerically simulated and the results show that (1) the failure zone of the excavated foundation pit without retaining walls appears at the corner and then gradually develops into the deep. However, the failure zone of the excavated foundation pit with retaining walls only develops longitudes along the retaining wall. (2) The sti ﬀ ness of retaining wall has a great impact on the failure zone of foundation pit excavation. The greater the sti ﬀ ness of retaining wall, the greater the damage degree. (3) The rationality of the proposed method is veri ﬁ ed by the comparisons of the simulation results of the proposed method with the ABAQUS numerical examples and the engineering practices. Future research directions should focus on developing the 3D parallel IKSPH programs. The research results can provide some references for the applications of SPH method into predicting the failure zone of foundation pit excavations and ensuring the safety of engineering constructions.


Introduction
With the accelerated development of China's economy, the constructions of urban underground space, large-scale water conservancy, and civil engineering projects are becoming the current focus of China, where the safety and stability of foundation pit during the excavation processes have become the main problems [1,2]. However, due to the complex geological conditions of foundation pit, it is easy to collapse during the processes of excavations. For example, the foundation pit excavation accident that occurred in Suzhou, China, in 2008 caused the destructions of the supporting structure and the damage of construction equipment, as shown in Figure 1(a). The collapse of the foundation pit in the Xiaoshan, Hangzhou in 2008 caused 21 people dead and 24 injured, as shown in Figure 1(b). Therefore, understanding and grasping the failure mechanisms of foundation pit engineering will undoubtedly provide an important guidance for ensuring the safety of engineering constructions and people's lives.
Previous works on the deformation and instability of foundation pit excavations mainly focused on three aspects: (1) experimental studies, (2) theoretical research, and (3) numerical simulation. Experimental studies are regarded as the most important and direct means to obtain the deformation and failure rules of foundation pit excavations, which can be divided into the field test and model test. Field test can directly obtain the deformation laws of actual engineering, but the cost is much. For example, Yang et al. [3] obtained the 3D deformation laws of the foundation pit in Hangzhou based on the monitoring data; Li [4] summarized the general rules of deep foundation pit deformations based on the measured data of the Shanghai subway foundation pit. Model test simplifies the complex and changeable factors of the engineering practice to a certain extent, so that the mechanics law of the whole model can be analyzed comprehensively. For example, Zhang and Qian [5] studied the laws of surface subsidence during the processes of foundation pit excavations under different displacement modes of rigid retaining wall through model test; Xia [6] carried out the model test on the penetration depth of underground diaphragm wall in cohesive soil and conducted a single study on many factors affecting the rebound of pit bottom. However, experimental studies only obtain the macroscopic deformation characteristics, but cannot quantitatively describe the internal mechanisms of excavation failure. Theoretical researches are based on the experimental results and summarize and refine the quantitative mathematical expressions of foundation pit deformation characteristics and failure rules. For example, Jiang [7] derived the analytical formula of pipeline deformation and internal force in the processes of foundation pit excavation by using the elastic foundation beam method; Gao et al. [8] established a lateral displacement prediction model for deep foundation pit excavation based on the combination of weighted first-order local method and trust domain; Zhang et al. [9] proposed a two-stage simplified analysis method for the longitudinal deformation of the adjacent existing tunnel caused by the excavation of double foundation pits based on the Pernak foundation model and analyzed the influence of the excavation of double foundation pits in soft soil on the vertical settlement of the tunnel. However, theoretical research can only obtain the analytical solution of the excavation deformation under simple boundaries and geometric shapes, and complex excavation steps as well as the complex foundation pit shapes will lead to extremely complex mathematical expressions. Meanwhile, previous experimental and theoretical studies on foundation pit excavation rarely paid attention to the progressive failure processes.
Numerical simulation can not only verify the correctness of theoretical studies but also can quantitatively reflect the inherent mechanisms of experimental research, which has been regarded as the "third method" of scientific researches [10,11]. The finite element method (FEM) was the first method used to study the foundation pit excavation [12][13][14]; however, FEM has limitations in dealing with excavation failures. Foundation excavation is a progressive failure process which contains the treatments of discontinuous properties such as crack propagation [15]. Therefore, the mesh refinements should be applied to crack tips, and mesh redivisions should also be applied to every step of progressive failure process, which costs huge amounts of computational resources. Meanwhile, for complex crack propagation paths (such as crack intersecting, etc.), the mesh grids will be extremely distorted, leading to the low calculation accuracy or even calculation failure. Different from FEM, the discrete element method gets free from the mesh grids [16,17], which discretizes the whole computational domain into particles. The interactions of different particles are characterized by the establishments of the contact model between particles, which can be well applied to modeling the progressive failure processes of foundation pit excavation. However, the DEM has many mesoscopic parameters with no actual physical meanings and requires complex parameter calibrations before numerical simulation, which is inconvenient to apply to engineering practice. The newly proposed discontinuous numerical method such as numerical manifold method (NMM) [18,19], peridynamics (PD) [20,21], and material point method (MPM) [22,23], which all have certain applications in the foundation pit excavation, but also have their limitations: the crack tips of NMM must be on the mesh nodes; the bond-based PD method has some theoretical defects which leads to the Poisson's ratio being constant; the MPM still needs background grids.
In this paper, based on the existing researches, the smoothing kernel function in the traditional SPH method has been improved to realize the failure characteristics of particles, which can reflect the progressive failure processes of the foundation pit during its excavation. Therefore, this method can also be called the Improved Kernel of Smoothed Particle Hydrodynamics (IKSPH). The "Killing Particle Method" has also been put forward to realize the simulations of the complex processes of foundation pit excavation. Based on the engineering practice of foundation excavation in Niulanjiang pumping station, the numerical simulation of progressive failure processes during foundation pit excavation is carried out, and influences of the stiffness of the retaining walls on the failure zone are also discussed. The research results can provide some references for the understandings of internal mechanisms of foundation excavation failure and ensuring the safety of the project.

Basic Principles of IKSPH
2.1. Solid Elastic Equations. The total stress tensor σ αβ in IKSPH can be expressed as the combinations of shear stress tensor τ αβ and the isotropic stress p, which can be written as: where the isotropic stress p can be expressed as:

Geofluids
where p H is the Hugoniot function, and Γ is the Gruneisen parameter.
The shear stress τ αβ is calculated by updating the shear stress rate of every step, which can be written as: where τ _ is the shear stress rate, and the shear stress tensor σ αβ can be calculated by multiplying the shear stress rate τ _ and the time step t. B is the shear modulus. ε αβ is the strain tensor. δ is the Kronecker delta. R αβ is the torsion tensor, which can be written as: 2.2. Governing Equations. IKSPH assigns to each particle its corresponding physical properties, including density, velocity, energy, and position coordinates. These particles should satisfy the following governing equations: where ρ is the density of the particle; m is the mass of the particle. v is the velocity tensor of the particle; e is the energy of the particle; x is the position of the particle; T is the artificial viscous term; W is the smoothing kernel function, which can be written as: where h is the smoothing length and R is the radio of average distances between particles and the smoothing length h.
2.3. Particle Pairing Method. Before IKSPH calculation, the first step to be carried out is the pairing between different particles. For example, the smoothing kernel function W in the governing equation (5) is calculated by particle pairing. This is also why IKSPH is different from FEM: the mesh grids in FEM have already been divided in advance, but the particles can move freely in the IKSPH method. While pairing particles in the IKSPH method, the number of paired particles inside the influencing domain should be determined first, and then the position relationship between two different particles can be then calculated.
The direct searching method is the simplest and most direct particle search method, which performs a full pair search by traversing all particles in each time step. We can find that the complexity order of this method is O ðN 2 Þ, which costs huge amounts of calculation resources when the particle number is relatively large.
The Linked-cell list method has the advantages of high efficiency and low memory saving compared with fullpaired searching method and is suitable for parallel computing. The details are as follows: the temporary searching grids are firstly laid onto the computing domain, as shown in Figure 2. The length of the temporary searching grid is defined as the searching radius of the particle, which is 2 h in our paper. For 1D, 2D, and 3D problems, the searching grids of any given particle (blue particle in Figure 2 as an example) are 3, 9, and 27 grids adjacent to it. By cycling through each particle, all pairs of particles can be found, whose complexity order is O ðNÞ.

Time Integration.
The Leap-frog integrating method is adopted in IKSPH, which has the advantage of low storage required for calculation. Meanwhile, only one optimization estimation is required for each calculation step. Therefore, the density, energy, velocity, and position of each particle can be obtained by cyclic iterations of the following formula: Search radius Base particle Critical base particle Searching grid 2 h Figure 2: The linked-cell list method.

Geofluids
In IKSPH, the determination of time step is related to the of material state change process, which can be estimated by the following equation: where f is the average force applied to the particle.

Failure Treatments of Particles
3.1. Failure Criteria. There is no unified criterion for soil failure in foundation pit excavation at present. Therefore, the improved Mohr-Coulomb criterion is selected in this section, which has two advantages: (1) the formula form is simple and does not need complex derivations; (2) the parameters are less and easy to access, which can be well applied to the engineering practice. The formula can be written as: where σ f and τ f are the tensile and shear stress on the failure surface. σ t is the tensile strength of the particle. c is the cohesion of the particle. φ is the internal friction angle of the particle. While judging whether the particle failure happens, equation (9) is firstly determined, which means that the tensile failure of the particle is easier to happen. When equation (9) is not satisfied, then equation (10) is determined whether the shear failure happens.
3.2. Treatments of Particle Failure. As can be seen from governing equation (5), the derivative of the smoothing kernel function ∂W ij,β /∂x β i governs the transfer of physical properties between different particles. Therefore, in order to reflect the failure characteristics of particles, the failure mark η is defined. When the particle failure occurs, η = 0, otherwise, η = 1, which can be clearly shown in Figure 3. The relationship between the improved smoothing kernel function D and the original smoothing kernel function W can be written as: Therefore, the final IKSPH governing equations considering the particle failure can be expressed as:

Killing Particle Method
To model the excavation processes of the foundation pit, different excavation parts should be grouped. In this section, a particle searching algorithm suitable for the IKSPH method is proposed, which can realize the particle grouping. The details are as follows: (1) The searching area of the target particles should be determined firstly, as shown in the purple area of Figure 4 (2) The searching points are generated uniformly on the target area (yellow points in Figure 4). The average spacing between searching points should be less than that of real particles, which is set to be 1/2 of the spacing between real particles (3) For every searching point, a searching radius r s is assigned. What should be noticed is that r s should 4 Geofluids be less than the average spacing of real particles, which is set to be 1/2 of the spacing between real particles (4) For every real particle covered by the radius of searching points, it is moved to the group of excavation particles. While operating the excavation process, the failure mark η of the excavation part is set to be 0 according to equation (11) to "kill" the particles

Stress Boundaries
The stress boundaries adopt the method of stress mapping, and more than 5 layers of "stress particles" are laid outside the solid particles. The "stress particles" shall have the following characteristics, as shown in Figure 5.
(1) "Stress particles" participates in the calculation of internal forces in IKSPH, and the particle density, mass, energy, and positions are updated according to equation (12) in each time step (2) In every time step, stress is reassigned to the "stress particles," that is, although the "stress particle" participates in the parameter updating of solid particles, its stress changes conform to the preset stress boundary requirements (3) A layer of "type I virtual particles" with velocity v_inf set to 0 should be laid outside the stress particles

Verification of IKSPH Method
In order to verify the proposed IKSPH method, a simple 2D cube model is established. The model size is 1 m × 1 m, and a crack with a length of 0.1 m is prefabricated in the center of the model, the dip angle of which is 45°. The model boundaries are subjected to 1 MPa confining pressure. Figure 6 shows the comparisons between the IKSPH results and the Abaqus results, which shows that compressive stress concentrates at the crack tip. Meanwhile, the maximum principal stress distributions calculated by the IKSPH program are consistent with the Abaqus results, which verify the proposed method.

Numerical Models
Based on the engineering practice, the corresponding 2D foundation model is established. The model size is 180 m × 70 m, which is greater than 3 times of the foundation pit depth, as shown in Figure 7. The whole model is divided into 360 × 140 = 50400 particles. Two retaining walls are set on the two sides of the model, and 4 excavation parts are arranged, the depth of which is 3 m, 4.5 m, 4 m, and 3 m, respectively. The detailed excavation simulation steps are as follows: first, 5000 steps of in situ stress balance are carried out. Then, the excavation steps are operated. The 5000-7000 th steps are the excavation of excavation part 1, the 7000-9000 th steps are the excavation of excavation part 2, the 9000-11000 th steps are the excavation of excavation part 3, and the 11000-13000 th steps are the excavation of excavation part 4.   Figure 8 shows the progressive failure process of the foundation pit excavation without retaining walls (Figure 8(a)) and with retaining walls (Figure 8(b)). As can be seen, for the condition without retaining walls, the failure zone firstly appears at the corner of the excavation part 2 and then develops to the deep of the foundation. Meanwhile, after the excavation of the last part, the failure zone appears on the slope surface. For the condition with retaining walls, the failure part of the foundation pit is similar to that of the condition without retaining walls, which appears at the corner of the excavation part 2. What is different is that the failure of the foundation only develops longitudes along the retaining wall, which is due to the fact that the existence of the retaining wall restricts the development trend of the failure zone to the deep. At the same time, after the final excavation, there is no failure particle on the exposed surface of the foundation pit, which is due to the fact that the retaining walls limit the large deformation of the soil, so the disturbance of the soil inside the foundation pit is relatively small, and the foundation pit is stable.

Influence of Retaining Wall Stiffness on the Failure Zone.
To quantitatively characterize the influence of different retaining wall stiffness on the failure zone of foundation pit excavation, the radio of the retaining wall elasticity modulus and the foundation elasticity modulus E R /E F is set to be 2, 5, 10, and 20, and corresponding numerical models are established for simulating. Figure 9 shows the failure zone under different conditions. As can be seen, the stiffness of retaining wall does have great impacts on the failure modes of the foundation pit. When the retaining wall stiffness is relatively small, the failure zone is limited to the small range in front of the retaining walls. This is because coordination deformation of retaining wall and foundation is good. However, when the retaining wall stiffness is large, the degree of incongruous deformation increases, and the failure range of soil becomes larger gradually. What should be noticed is that when E R / E F = 20, a large failure zone appears behind the retaining wall, which indicates that risk of foundation pit instability is increasing under this circumstance. Figure 10 shows the damage counts under different retaining wall stiffness. As can be seen, with the increasing of E R /E F , the damage counts increase accordingly, which means that high stiffness retaining wall will have negative effect on the foundation pit stability.    Figure 11 shows the comparisons between the IKSPH results and the engineering practice. As can be seen, the failure positions are all at the corner of the foundation pit, which verifies the proposed method. Figure 12 shows the maximum principal stress distributions and the particle velocity distributions during excavation. We can find that the tensile stress concentrates at the foundation pit corner, meanwhile, the velocity vector deviates at the corners of the foundation pit, indicating that the failure of the corner is the tensile failure. Therefore, in practical engineering, support should be set at the corner points of the foundation pit to prevent tensile failure.

Application Prospects of IKSPH in Failure Prediction of
Foundation Pit Excavation. In this paper, by improving the smoothing kernel function in the traditional SPH method, we can realize the failure modeling of particles at microscale. Compared with traditional FEM, IKSPH gets free from mesh grids, which can therefore well reflect the large deformation, failure, damage, and other discontinuous characteristics of rock and soil. Meanwhile, compared with DEM, its parameters have definite physical meanings. Therefore, IKSPH has a wide application prospect in rock and soil failure simulation and prediction. What should be stressed is that this paper only considers the numerical realization of foundation pit excavation under simple circumstances. In practical foundation pit engineering, there are many kinds of supports and the soil properties are different. Therefore, the real construction processes of actual engineering should be taken into consideration in the subsequent research. Meanwhile, the actual foundation pit engineering is a complex 3D problem, and simplifying into a 2D problem will miss lots of useful information. However, the computational efficiency of the 3D IKSPH program is one of the difficult problems in the field of geotechnical simulation. So future research should focus on the development of 3D parallel IKSPH program and its applications in 3D foundation pit excavation simulation.

Conclusions
(1) The failure mark η is introduced in this paper, and the smoothing kernel function in the traditional SPH method has been improved, which can realize the simulation of progressive failure process of rock and soil (2) The "Killing Particle Method" has been put forward to realize the complex excavation processes of foundation pit (3) The progressive failure processes of the foundation pit are simulated. The failure zone of the excavated foundation pit without retaining walls appears at the corner and then gradually develops into the deep. However, the failure zone of the excavated foundation pit with retaining walls only develops longitudes along the retaining wall (4) The stiffness of the retaining wall has a great influence on the failure zone of foundation pit excavation. The greater the stiffness of the retaining wall, the greater the damage degree (5) The numerical simulation results of IKSPH are consistent with the calculation results of commercial software ABAQUS and engineering practice, which verifies the correctness of the proposed method.
Meanwhile, future research directions should focus on the development of a 3D parallel IKSPH program

Data Availability
The program data used to support the findings of this study are restricted by Hohai University in order to protect the privacy. Data are available from yushuyang_hhu@163.com for researchers who meet the criteria for access to confidential data.

Conflicts of Interest
None of the authors have any conflicts of interest.