Abstract
The excavation damage of deep tunnels is one of the most important factors contributing to the failure of tunnel structures. In order to investigate the influence of tunnel shapes and fissure geometries, the kernel function in the traditional SPH method has been improved, which can realize the brittle fracture characteristics of particles and can be called the Improved Kernel of Smoothed Particle Hydrodynamics (IKSPH-2D). Meanwhile, the random fissure generation method in IKSPH has been put forward. Different tunnel shapes, fissure geometries, and locations are considered during the simulation of tunnel excavation, and results show that (1) the typical “V”-shaped shear damage zones appear after the tunnel excavation, which is consistent with engineering practice. Meanwhile, tunnel excavation also has an “activating” effect on the preexisting fissures. (2) The stability of circular-shaped tunnel is the best, while horseshoe shaped tunnel is worse, and the “U”-shaped tunnel is the worst. (3) Fissures with small and large dip angles have the greatest influence on the stability of tunnel excavation. With the increase of fissure numbers and lengths, the tunnel tends to be instable. (4) The IKSPH method gets free from traditional grids in FEM, which can dynamically reflect the fracture processes of tunnel excavation. Meanwhile, developing 3D IKSPH parallel program will be the future directions.
1. Introduction
A large number of large-scale hydropower stations have been developed in western China for its rich hydropower resources. However, these projects are mostly built between high mountains, and the excavations of tunnels will lead to redistributions of stress field, thus causing the excavation damage zone (EDZ) [1–3]. Meanwhile, rock is a typical heterogeneous material. Due to complex crustal movements, there exists large number of joints, fissures, and holes inside rock masses. Therefore, engineering disasters such as tunnel collapse, water inrush, and even rock burst are easy to occur during tunnel excavation, which will pose a great threat to the safety of tunneling equipment and human lives [4, 5]. Figure 1(a) shows the rock burst at Jinping drainage tunnel on 28 November 2009 [6], which resulted in 7 deaths and 1 injury, and caused the TBM to be buried. Figure 1(b) shows the collapse of Jinping No. 4 drainage tunnel [7]. Therefore, the evaluation of the damage range of tunnel excavation will undoubtedly provide guidance for the safe constructions of tunnel engineering.
(a)
(b)
Previous researches on the tunnel excavation damage mostly focused on the theory, experiment, and numerical simulation. Theoretical research can quantitatively analyze the influence of tunnel excavation disturbance on rock masses. For example, Castro et al. [8] utilized the deviatoric stress method to study the damage distributions after tunnel excavations in SNO underground laboratory in Canada; Zhou et al. [9] deduced the theoretical formula of tunnel excavation damage range according to the failure criterion obtained from engineering experience and the classical results of tunnel stress solutions. Cai et al. [10] put forward the stress threshold values of crack initiation in brittle rock masses after tunnel excavations, and the size of the excavation damage zone is judged according to the stress redistribution of surrounding rock masses. However, theoretical researches can only obtain the exact solution of homogeneous rock mass excavation under simple boundary conditions, while complex conditions will lead to complex mathematical expressions. Experimental studies can directly show the damage range of tunnel excavation. For example, Yan et al. [11] analyzed the characteristics and causes of surrounding rock damage zone under different excavation methods by in situ detection test; Zhang et al. [12] carried out in situ tests on Jinping deep tunnels, and the excavation damage of the tunnel was comprehensively evaluated by using stress monitoring, acoustic emission monitoring, displacement monitoring, and digital drilling camera technology; Chen et al. [13] carried out the acoustic emission monitoring tests during TBM construction to study the surrounding rock damage laws. However, previous experimental studies mostly used indirect methods to evaluate the damage degree of tunnel excavation, and the internal mechanisms cannot be intuitively exhibited.
As the “third method” of scientific research, numerical simulation has developed rapidly in recent years. Numerical simulation can not only verify experimental researches and theoretically studies, but also can directly reveal the variations of internal stress during rock failure, which breaks the “black box” existing in the rock mechanics. The tunnel excavation damage was firstly simulated by the elastoplastic finite element method [14, 15], and the damage range was characterized by plastic zone. The disadvantages are as follows: (1) quasistatic simulation of the elastoplastic finite element method cannot reflect the dynamic damage processes of surrounding rock during tunnel excavation, which is not consistent with the actual situations. (2) The mesh grids around the cracks or joints must be changed at every time step during the crack propagation processes [16, 17], and the computation load is relatively huge. Under the circumstances of complex joints or crossing cracks, it may lead to low calculation accuracy or calculation failure. Subsequently, the extended finite element method (XFEM) [18, 19] has been proposed to model the discontinuous features. However, mesh refinements are still needed. The discrete element method (DEM) gets free from the grids existing in traditional FEM, which can be well applied to the simulations of discontinuous problems [20, 21]. However, DEM has numerous microscopic parameters with no physical meanings, and complex parameter calibrations are needed before numerical simulation, which is inconvenient to apply. The other discontinuous numerical methods such as peridynamics [22], numerical manifold method [23], material point method [24, 25], and efficient remeshing techniques [26–29] all have their unique advantages in dealing with discontinuous problems, but still have some disadvantages. For example, the PD method has some theoretically defects which leads to Poisson’s ratio being constant [30]. The crack tips must be on the mesh nodes in NMM method [31]. The MPM method still need background grids [24].
Smoothed particle hydrodynamics (SPH) is a pure Lagrangian meshless method, which gets free from traditional grids in FEM and can deal with discontinuous problems. However, few applications are made in the field of rock mass excavations. Meanwhile, the traditional SPH method did not contain the treatments of particle damage. Thus, in this paper, an improved form of the SPH method has been proposed, which is named Improved Kernel of Smoothed Particle Hydrodynamics (IKSPH) [32]. By improving the traditional kernel function, the brittle fracture characteristics of particles have been realized, which can then be applied to modeling the rock fracturing processes. The “random fissure generation method” has also been put forward to generate the random distributed fissures in the model. The IKSPH model of Jinping tunnel is established, and the tunnel shapes, fissure geometries, and locations are considered during the excavation processes. The damage range, fracture type, and fracture counts are analyzed. The research results can provide some references for the understanding of tunnel excavation damage laws and also some guidance for its support measurements.
2. Basic Principles of IKSPH
2.1. Constitutive Equations
The total stress tensor consists of two different parts: (1) the isotropic pressure and (2) shear stress . The formula can be expressed as
where is the kronecker symbol; the isotropic pressure can be obtained by equation of state:
where is the Hugoniot function, and is the Gruneisen parameter.
IKSPH updates the stress components by updating the stress rates. By introducing the Jaumann rate, the shear stress tensor rate can be then expressed as
where −τ is the shear stress rate, is the shear modulus, and is the torsion rate, which can be written as
2.2. Governing Equations of IKSPH
For every particle in IKSPH, three equations should be satisfied, namely, (1) density equation, (2) momentum equation, and (3) motion equation, which can be written as
where is the density of the particle. is the calculation time step. is the mass of the particle. is the velocity of the particle. is the kernel function between particles. is the position of the particle. σ is the stress tensor of the particle. is the energy of the particle. is the artificial viscosity, which can be expressed below. What should be stressed is that upper Greek indices refer to tensor notation while the lower Latin letters refer to particle numbers.
where , , , , , and are the standard constants, whose values are usually taken as 1.0, is set here to prevent the base particles from approaching each other and can be expressed as , and is the sound speed.
3. Treatments of Particle Damage
3.1. Selection of Fracture Criterion
There has been no universally accepted fracture criterion applicable to all rock materials so far. Therefore, the improved form of maximum principal stress criterion has been adopted in our paper, and the reasons are as follows:(1)The form is simple and does not need complicated formula derivations(2)The formula has less parameters, and these parameters can be easily obtained, which can be well applied to engineering practices
Therefore, the improved form of maximum principal stress criterion can be written as
where σf and τf are the tensile and shear stress of failure surface. σt is the tensile strength of the rock masses. is the cohesive strength of the rock masses. φ is the internal friction angle of the rock masses. While using this criterion, equation (9) is firstly determined whether it is satisfied, which judges whether the tensile failure happens. If equation (9) is not satisfied, then equation (10) is determined whether shear failure happens.
3.2. Treatments of Particle Damage
As can be seen from IKSPH governing equations (5) ~(7), the kernel function in the traditional SPH method governs the transfer of information between particles. Therefore, to cut off the link between particles, a fracture mark ξ has been introduced in this paper, and the improved form of the kernel function can be expressed as
where is the improved kernel function considering the particle damage characteristics. When the particle is not damaged (the stress components do not satisfy equation (5) or equation (6)), the . When particle damage occurs (the stress components satisfy equation (5) or equation (6)), then the . Therefore, the final governing equation can then be expressed as below, and the fracture processes of IKSPH particles can then be illustrated in Figure 2.
4. Random Fissure Generation in IKSPH
Complex fissure geometries and locations exist in the surrounding rock masses in deep tunnel engineering. Therefore, it is necessary to utilize the statistical probability method to consider the influences of random fissure distributions on tunnel excavation damage. In this section, the random fissure geometry generation methods are introduced to provide preliminary preparations for the generations of random fissure particles in IKSPH below.
The first step in generating random fissures is to generate random numbers. The Monte Carlo Method is introduced, and the linear congruence method is used to generate uniformly distributed random numbers, which can be expressed as
where is the random numbers. is the multiplier. is the increasement. is the module. Mod represents the remainder of . is a random number between [0, 1]. Therefore, the normal distributions of random numbers can be written in the following form:
where is the random number that obeys normal distributions. μx is the mean value. σx is the standard deviation.
For 2D problems, each fissure can be characterized by its middle point coordinate (, ), fissure length , and dip angle θ. Therefore, we can obtain the two endpoints of the fissure (, ) and (, ):
Therefore, the coordinates of the crack endpoints are then determined by Equation (10).
5. Numerical Models and Parameters
5.1. Parameter Calibrations
In order to characterize the heterogeneity of rock materials, the heterogeneity coefficient is introduced, and the Weibull function can be expressed as
where represents the basic parameters of the particle, such as elastic modulus and compressive strength. is the heterogeneity coefficient. is the mean value of particle parameters.
For the convenience of calculations, only the heterogeneity of compressive strength is considered. The 2D model with a size of is established under the boundary of uniaxial compression, and the heterogeneity coefficient is set to be 5. Through trial-and-error tests, the stress-strain curve of the numerical simulation is obtained, which is consistent with previous experimental results [33], as shown in Figure 3. Meanwhile, the final failure modes of the numerical results are also consistent with experimental results [33], which mean that the calibrated parameters can be applied to engineering practice. The final compression strength σc is 155 MPa; Elasticity modulus is 50.59 MPa, Poisson’s ratio is μ 0.2, and the internal friction angle φ is 40°, as listed in Table 1.
5.2. Numerical Model and Calculation Conditions
The IKSPH numerical model is established according to the engineering practice of Jinping tunnels. The size of the model is . Four calculation conditions are considered, namely, (1) condition A: different hole shapes: A1: circular hole; A2: “U”-shaped hole, and A3: horseshoe-shaped hole. (2) Condition B: different fissure dip angles: B1: ; B2: ; B3: , and B4: . (3) Condition C: different fissure numbers: C1: , C2: , C3: , and C4: . (4) Condition D: different fissure lengths: , D2: , D3: , and D4: . The details are shown in Table 2. The whole model is divided into particles, which is shown in Figure 4.
The tunnel excavation steps in IKSPH program are as follows: Firstly, apply the in situ stress of Jinping tunnel to the model boundaries, the horizontal stress , the vertical stress , and the shear stress . Then, the stress balance is carried out in 7000 calculation steps. Finally, the excavation is carried out by setting the fracture mark ξ of tunnel particles to be 0.
6. Numerical Results
6.1. Influence of Tunnel Shapes on the Excavation Damage
Figure 5 shows the excavation damage ranges of different tunnel shapes. The white color represents the tensile failure, and the red color represents the shear failure. As can be seen, shear failure concentrates in the directions of 5 o’clock and 11 o’clock, which is consistent with the direction of the maximum principal stress. Meanwhile, the tunnel shapes have great impacts on the excavation damage. For the circular-shaped tunnel, the damage range reaches maximum in the directions of 5 o’clock and 11 o’clock, whose depths reach about 3.5 m. However, the damage range of other parts is relatively less, and the depths of which are within 1 m. For the horseshoe-shaped tunnel, apart from the directions of 5 o’clock and 11 o’clock, the damage range in the directions of 3 o’clock, 6 o’clock, and 9 o’clock are also large, the average depths of which is about 3 m. For the “U”-shaped tunnel, the damage ranges of side wall and bottom floor as well as the direction of 11 o’clock are relatively larger, and the average depths of which are about 3.5 m. What should be pointed out is that the excavation damage not only happens around the tunnel but also around the preexisting fissures, which is regarded as the “activating” effect and mostly the tensile damage.
(a)
(b)
(c)
Figure 6 shows the fracture counts of particles under different tunnel shapes. As can be seen, the ranks of damage degree under different tunnel shapes (judged by the fracture counts) are , indicating that the stability of circular-shaped tunnel is the best, the horseshoe-shaped tunnel is worse, and the “U”-shaped tunnel is the worst. In general, the tensile failure counts are more than shear failure counts. The tensile failure counts of horseshoe-shaped tunnel are more than that of circular-shaped tunnel, but the shear failure counts of which are less. The tensile failure counts of the “U”-shaped tunnel are more than that of horseshoe shaped tunnel, and the shear failure counts of which are equal.
(a)
(b)
(c)
6.2. Influence of Fissure Dip Angles on the Excavation Damage
Figure 7 shows the excavation damage ranges of different fissure dip angles. As can be seen, the cracks propagate through the tunnel and the preexisting fissures in the conditions of , , and . This is because distances between the fissures and tunnels are relatively small, which leads to the strong interactions. Therefore, in the actual engineering practice, special attentions should be paid to the hidden joints near the tunnel excavation, and necessary support measurements should be taken. Meanwhile, the damage degrees in the conditions of and are larger than that of and , which means that small and large dip angles have the greatest influence on the stability of tunnel excavation.
(a)
(b)
(c)
(d)
Figure 8 shows the fracture counts of particles under different fissure dip angles. As can be seen, the ranks of damage degree under different fissure dip angles (judged by the fracture counts) are , showing that small and large dip angles have the greatest influence on the stability of tunnel excavation. At the same time, the shear fracture counts are similar, but the tensile fracture counts are larger in the conditions and .
(a)
(b)
(c)
(d)
6.3. Influence of Fissure Numbers on the Excavation Damage
Figure 9 shows the excavation damage ranges of different fissure numbers. As can be seen, when the number of fissures is small, the excavation damage mainly happens around the tunnel, which is similar to the circumstances with no preexisting fissures. However, with the fissure number gradually increases, the excavated tunnel is more likely to interact with preexisting fissures, which contributes to the cracks propagating through the tunnel and the preexisting fissures. Notably, the interactions between tunnel and fissures are mainly the shear failure (the red marks in Figure 9).
(a)
(b)
(c)
(d)
Figure 10 shows the fracture counts of particles under different fissure numbers. As can be seen, the ranks of damage degrees under different fissure numbers are . In general, with the increase of fissure numbers, rock masses are more fragmented, which leads to the large damage degrees. What should be stressed is that the damage degree of is less than that of , and this is because the random fissures are closer to the excavated tunnel, which leads to the stronger interactions. Therefore, the locations of random fissures are also an important factor influencing tunnel stabilities.
(a)
(b)
(c)
(d)
6.4. Influence of Fissure Lengths on the Excavation Damage
Figure 11 shows the excavation damage range of different fissure lengths. As can be seen, when the fissure lengths are relatively small, it has limited influences on the tunnel excavation damage; meanwhile, the excavation process also has little impacts on the “activating” effect on the preexisting fissures. With the increase of fissure lengths, the probability of its interactions with the excavated tunnel increases, which will lead to the cracks propagating through the tunnel and the preexisting fissures, as shown in Figure 11(d). Meanwhile, the damage degrees under long fissure lengths are also large, which means that with the increase of fissure lengths, the stability of tunnel excavation gets worse.
(a)
(b)
(c)
(d)
Figure 12 shows the fracture counts of particles under different fissure lengths. As can be seen, the ranks of damage degree under different fissure numbers are . We can find that when the fissure length increases from 3 m to 6 m, the damage counts increase sharply. However, when the fissure length exceeds 6 m, the damage counts increase little, which indicates that the fissure length 6 m is a critical value. At the same time, when the fissure length is 3 m, the tensile fracture counts are almost equal to the shear fracture counts. With the increase of fissure lengths, the tensile fracture counts increase accordingly but the shear fracture counts remain unchanged.
(a)
(b)
(c)
(d)
7. Discussions
7.1. Influence of the Discretization
In order to study the influence of the discretization, the model A1 is selected for analysis, and different numbers of particles are set: A1-1: ; A1-2: , and A1-3: . Meanwhile, different time steps are also shown in Figure 13. We can see that the numerical results under different particle numbers are consistent, which verify the discretization method.
(a)
(b)
(c)
7.2. Validation of the IKSPH Method
In order to verify the rationality of the proposed method, the excavation damage of the tunnel with no fissures is simulated, and the numerical results are compared with the engineering practice [34], as shown in Figure 14. As can be seen, the excavation damage occurs in the direction of maximum principal stress, forming the typical “V”-shaped damage zone, which is consistent with the IKSPH numerical results.
(a)
(b)
Meanwhile, compared with the traditional FEM, the proposed IKSPH method can get free from the traditional grids, which do not need to apply special treatments to crack tips and can realize the random crack propagations. The IKSPH calculation properties can truly reflect the actual state of tunnel excavation and can be well applied to the rock mechanics engineering.
What should be stressed is that the actual rock engineering is the complex 3D problems. However, 3D IKSPH program has the disadvantages of low calculation efficiency. Therefore, future studies should focus on the development and application of 3D parallel IKSPH program.
8. Conclusions
(1)The kernel function in the traditional SPH method has been improved by introducing a fracture mark ξ, which can realize the brittle fracture characteristics of rock mechanics(2)The random fissure generating method has been put forward, which can realize the generations of random fissure particles in IKSPH and can better reflect the real characteristics of rock masses(3)The excavation of tunnels not only produce damage zone around tunnels but also have an “activating” effect on the preexisting fissures(4)The tunnel shapes, fissure lengths, dip angles, and fissure numbers all have great impacts on the tunnel excavation damage. The damage zone mainly occurs in the direction of maximum principal stress, and the counts of tensile failure are more than that of shear failure(5)The IKSPH method has been verified by comparing the numerical results with engineering practice. Meanwhile, developing 3D parallel IKSPH program will be the future research hotspot
Data Availability
The program data used to support the findings of this study have not been made available because it is applying for a patent.
Conflicts of Interest
Ren Xu-Hua, Yu Shu-Yang, Wang Hai-Jun, Zhang Ji-Xun, and Sun Zhao-Hua declare that they have no conflict of interest.
Authors’ Contributions
Ren Xu-Hua provided the ideas, Yu Shu-Yang did the numerical simulation and wrote the draft, and Wang Hai-Jun, Zhang Ji-Xun, and Sun Zhao-Hua edited the draft.
Acknowledgments
We acknowledge the financial supports of the National Natural Science Fund (Grant No. U1765204), “the Fundamental Research Funds for the Central Universities” (B210203078), and the National Natural Science Fund (51409170).