Numerical Analysis of the Hydraulic Fracturing of Pressure Tunnel Lining Based on the 2P-IKSPH Method

In order to investigate the fracture mechanisms of the pressure tunnel lining under water-stress coupling, based on the traditional smoothed particle hydrodynamics (SPH)method, the solid-water particle interactionmethod, and the particle damage conversion algorithm are proposed to realize the hydraulic fracturing process, which is called the 2P-IKSPH method. +e “particle domain searching method,” the “birth-and-death particle method,” and the “group discrimination searching method” have also been proposed to realize the simulations of complex processes of excavation, lining, and operation of the hydraulic tunnel. Taking the Guzeng hydraulic tunnel as an engineering example, the hydraulic fracturing of tunnel lining under different conditions is numerically simulated. Results show the following: (1) the 2P-IKSPHmethod can dynamically reflect the stress wave propagation processes of surrounding rock and the damage process of tunnel lining. (2) +e lining damage mainly occurs on the vault and the arch foot. (3)+e critical internal water pressure increases with the increase of the tunnel buried depth and the thickness of lining, but increases first and then decreases with the increase of the surrounding rock mass grade.+e research results can provide some references for the optimization designs of tunnel lining and reinforcement of similar projects. Meanwhile, developing 3D parallelization program based on 2P-IKSPH will be the future research directions.


Introduction
With the implementation of the strategy of energy development in western China, a number of high dams have been built. However, most of the projects have been built between high mountains, so safe and stable operations of long-distance hydraulic pressure tunnels under complex stress conditions have become a difficult problem [1,2]. As an important part of hydraulic tunnel, the lining structures can not only guarantee the stability and limit the deformations of surrounding rock [3] but also reduce the roughness, improve the flow conditions, and reduce the water head loss [4]. Nevertheless, under the actions of complex geological conditions, in situ stress and internal water pressure, the lining structures will inevitably crack, which will lead to serious inner water exosmosis. At the same time, the existence of cracks leads to the stress redistribution of lining structures. e stress concentrates at the crack tips, which will easily contribute to the disaster [5,6]. erefore, the understanding of the damage mechanisms of tunnel lining under water-stress coupling will provide guidance for the design of tunnel lining and reinforcement as well as ensure the safe and stable operations of pressure tunnels.
Field survey is the most direct way to evaluate the damage degrees of tunnel lining. However, it can only be carried out when the tunnel is emptied, and it has the disadvantages of long time period, high risk, and small observation ranges. Meanwhile, most of the observed results are after the lining failure, so the whole damage processes cannot be directly described. Although various means have been developed to detect the lining structure damage, such as underwater vehicles [7], which can be carried out without the limits of the tunnel status, its detection range is relatively small and the equipment is very expensive. A variety of observation techniques have also been proposed, such as acoustic wave method [8] and microseismic-acoustic emission method [9], which can quantitatively detect the damage extents, but the detection accuracy and ranges are still limited. e indoor geomechanical model test can make up for this disadvantage [10]; nevertheless, the similarity between small-scale test and engineering practice is still a difficult problem to be solved. At the same time, the experiment cost is too high and the dispersion of the test results is also large.
As "the third method" of scientific researches [11], numerical simulation can not only dynamically reflect the failure characteristics of tunnel lining but also clarify its internal mechanisms, which has been developed rapidly in recent years. e finite element method (FEM) has been one of the earliest measures to study the water-stress coupling damage of the lining structures [12], but special treatments should be applied to the crack tips, and the mesh of complex crack networks will be extremely distorted [13,14]. e Discrete Element Method (DEM) discretizes the analysis domain into particles, and the interactions of different particles are realized by various contact models [15], which avoids the mesh divisions. However, its mesoscopic parameters are numerous and have no practical meanings, which is inconvenient to apply to the engineering practice [16].
e recently new numerical measures such as the Numerical Manifold Method (NMM) [17], the Phase Field Method [18], the Peridynamics (PD) [19], and the Material Point Method (MPM) [20] all have their unique advantages in dealing with failure modeling of tunnel lining, but the limitations also exist. e crack tips cannot be inside the element in the NMM method [21]. e bond-based PD method has some theoretical defects, which leads to Poisson's ratio being constant [22]. e Phase Field Method still has some problems in dealing with multiple cracks propagation [23]. e MPM method needs to match the material points through background grids [24].
In this paper, an improved water-stress coupling numerical method based on SPH is proposed, which is called two-phase improved kernel of smoothed particle hydrodynamics (2P-IKSPH), to realize the water-stress coupling damage analysis of tunnel lining. Compared with previous studies, the water-stress coupling can be calculated in one set of solid equation, which reduces the calculation amount and can improve the calculation efficiency. e solid-water particle interaction method and the particle damage conversion algorithm have been put forward to realize the hydraulic fracturing processes.
e "particle domain searching method" has been proposed to realize the generations of water particles, lining particles, and excavation particles. e "group discrimination searching method" has been proposed to realize the applications of water pressure to the water particles and fissure particles.
e "birth-anddeath particle method" has been proposed to realize the complex processes of tunnel excavation, lining, and operating. Based on the engineering background of Guzeng pressure tunnel, numerical simulations of critical internal water pressure and the failure modes of tunnel lining under different conditions have been carried out. e research results can provide some references for the understanding of failure mechanisms of tunnel lining and ensure the safe operations of tunnels.

Basic Principles of 2P-IKSPH
2.1. Governing Equations. For each particle in 2P-IKSPH, the following equations need to be satisfied, namely, (1) continuity equation, (2) momentum equation, (3) energy equation, and (4) motion equation, which can be expressed as follows [25,26]: where ρ is the density of the base particle, t is the calculation time step, m is the mass of the base particle, v is the velocity of the base particle, x is the position of the base particle, σ αβ is the total stress tensor of the base particle, e is the energy of the base particle, τ αβ is the shear stress tensor of the base particle, ε αβ is the strain tensor of the base particle, T is the artificial viscous part, which can reduce nonphysical oscillations during calculations [27], and W is the kernel function of the IKSPH method, which can be expressed as follows [28]: where h is the smoothing length, which is usually taken as the distance between two particles, and R is the ratio of the base particle distance and the smoothing length h.

Solid Constitutive Equation.
σ αβ consists of two parts: (1) the hydrostatic pressure p and (2) the shear stress τ, which can be written as

Mathematical Problems in Engineering
where p can be obtained by calculating the solid state equation [26]: where p H is the Hugoniot curve function and Γ is the Grüneisen parameter. e shear stress rate τ can be updated by the following equation: where τ ⌢ is the stress rate tensor, B is the shear modulus, R is the torsion tensor, and ε αβ is the strain tensor, which can be expressed as R αβ in equation (6) is the torsion tensor, which can be written as

Solid-Water
Interaction. For the neighboring solid particle and water particle, the water particle applies a force p w to the solid particle. Different form the previous studies [27], the flow characteristics of water particles are no longer considered in this paper to simplify the calculation processes. e water-stress interaction is shown in Figure 1; as can be seen, the normal unit vector between water and solid particles can be expressed as where i represents the solid particle and j represents the water particle.
ij . e direction cosine between the normal unit vector and the x-axis n x,i as well as the y-axis n y,i can then be written as en, the stress components of the water particle to solid particle can be written as (1) control the transfer of information between different particles. When the solid particle damage occurs at the mesoscopic scale, setting the value of zW ij,β /zx β i to 0 will reflect the brittle fracture characteristics. erefore, a fracture mark variable ξ is introduced in the 2P-IKSPH program. When the particle is damaged, then ξ � 0; otherwise, ξ � 1 [29], and the brittle fracture of IKSPH particle is shown in Figure 2. e Mohr-Coulomb criterion with a tension cut-off is adopted here: where σ f and τ f are the maximum tensile stress and shear stress of the failure surface, σ t is the tensile strength of the particle, c is the cohesion of the particle, and φ is the internal friction angle of the particle. In our simulation, when the stress components satisfy equation (12), tensile failure first happens, if equation (12) is not satisfied, then the program will judge whether equation (13) satisfies, which means the shear failure. en, we can obtain the improved form of the derivations of smoothing kernel functions, which can be written as x x σ p y y Solid particle Water particle Solid-solid contact Solid-water contact where D is the improved form of the derivative of kernel function in 2P-IKSPH. erefore, the 2P-IKSPH governing equations can be rewritten as

Particle Domain Searching Method.
In numerical simulations of underground rock mechanics, the generations of the excavation part and the tunnel lining should be treated separately. In this section, the "particle domain searching method" is proposed, as shown in Figure 3, and the detailed steps are as follows: (1) Define the geometric locations of the excavation part and the tunnel lining. e circle shaped excavation part T i can be determined by the central point x ti and its radius r ti . e tunnel lining L i can be determined by the central point x li , the inner radius l si , and the outer radius l bi . (2) Generate a series of "searching particles" in T i and L i .
What should be noticed is that these "searching particles" should distribute uniformly, and the distance between two neighboring "searching particles" should be less than the minimum spacing between real particles. (3) A searching radius d is given to each "searching particle." Generally, d should be more than the average spacing between real particles. For those real particles inside the radius d, the excavation or lining attributes are given to the target particles.

Birth-and-Death Particle Method.
To realize the excavation-lining-operating processes of pressure tunnel, the particles need to be removed and then reactivated. In this section, the birth-and-death particle method is introduced here, as shown in Figure 4. For those particles that need to be "killed" (which is called the "death particles"), the particle domain searching method is utilized here for searching out the target particle groups, and equation (14) is used to "kill" these particles by setting ξ � 0. When the reactivation is required, ξ � 1, and the stress components of the "birth particles" need to be set to 0.

Group Discrimination Searching Method.
During the calculations of hydraulic fracturing processes, the lining damage is caused not only by hydraulic factors, but also the pure mechanical factors. erefore, the "group discrimination searching method" is put forward in this section to distinguish whether the particle damage is caused by the water pressure or the pure mechanical factors; only by this means can we decide to whether to add the water pressure. e steps are shown in Figure 5, which can be expressed as follows: (1) Define the arrays to store water particle group W i and fissure particle groups F i . (2) After the stress update of each particle, the newly damaged particles are transferred to the fissure particle groups F i temporarily. A searching radius d t is given to each water particle in W i , which should be more than 2 times the distance between real particles. (3) For those real particles inside the radius d t and belonging to the group F i , they are marked as water  particles and transferred to W i . For those fissure particles, no water pressure is applied to them.

Model Validation.
To verify the proposed numerical method, the numerical model of a cubic shape with one single crack is established. e model size is 1 m × 1 m, and one hydraulic fracture with a length of 1 m is prefabricated in the center, on which the internal water pressure of 1 MPa is acted on. e confining pressure of 1 MPa is acted on the model side. Figure 6 shows the maximum principal stress distribution of the model calculated by the 2P-IKSPH method proposed in this paper and the Abaqus software. As can be seen, the calculation results by 2P-IKSPH and Abaqus are consistent, which indicates that the proposed method is accurate and reasonable.

Engineering Background.
Guzeng Hydropower Station is located on the main stream of Muli River in Muli County, Liangshan Prefecture, Sichuan Province. e total installed capacity is 4 × 43 MW, the normal water level is 2215.00 m, and the total reservoir capacity is 484000 m 3 . e diversion tunnel of Guzeng Hydropower Station is located on the right bank of Muli River with a total length of 11060.119 m. e horseshoe-shaped tunnel and lining is adopted, and the cross section size is 6.1 m × 7.8 m. ere is a series of regional epimetamorphic rocks along the diversion tunnel of Guzeng Hydropower Station. From old to new, they are Ordovician Group (O1r), Ordovician Group (O1 W), Silurian Group (S1), Carboniferous Group (C1), and Triassic Group (T3q), which are shown in Figure 7.

Calculation Conditions.
e 2P-IKSPH model is set up according to the actual size of the Guzeng hydraulic tunnel, as shown in Figure 8. e model size is 60 m × 60 m, which is more than 5 times the diameter of the tunnel, as shown in Figure 8(a). e whole model is divided into 300 × 300 particles. e horseshoe-shaped tunnel is set in the middle of the model, as shown in Figure 8(c). e horseshoe-shaped lining is also shown in Figure 8(d).
According to the engineering practice of Guzeng tunnel, three calculation conditions are set in this section: A: different tunnel buried depths, B: different lining thickness, and C: different surrounding rock grades, as shown in Table 1.
e critical internal water pressure and the failure modes of tunnel lining are calculated in our simulation. What should be stressed is that the surrounding rock is treated as homogeneous and elastic media. e mechanical parameters of the lining are calibrated according to the C25 concrete. Considering the random distribution characteristics of concrete properties, the two-parameter Weibull functions are introduced here [30].
where x is the basic mechanical parameters of particles (e.g., modulus of elasticity, compressive strength, cohesion, etc.), x 0 is the mean value of basic mechanical parameters, and m is the heterogeneity extent of particles. By large numbers of trials, the simulated compressive strength of concrete is 24.6 MPa, which is similar to the compressive strength of C25 concrete 22.5 MPa. e calibrated parameters can be then used in the engineering practice.
e numerical simulation steps are as follows: firstly, the in situ stress balance of 4000 steps is carried out, where the boundary stress of the first 1000 steps is loaded to the predetermined stress and keeps unchanged in the later 3000 steps. en, the excavation operation of the tunnel is carried out in the 4000th step, the lining is applied in the 5000th step, and the internal water pressure is added in the 6000th step. e termination of the simulation is when the crack propagates through the lining. Figure 9 shows the minimum principal stress distributions of excavation-lining-operating processes of condition A2. As can be seen, 2P-IKSPH method can dynamically reflect the stress wave propagation processes of surrounding rock and the damage processes of tunnel lining. is is also the advantage of the proposed method. Figure 10 shows the failure modes of the tunnel lining under different buried depths, where the white color represents the tensile failure and the red color represents the Original kernel Improved kernel Water-solid contact Original particle Fissured particle Water particle Searching radius

Triassic group T3q
Silurian group S1 Carboniferous group C1 Ordovician group O1r Ordovician group O1w Hydraulic tunnel Joints  shear failure. As can be seen, under the combined action of surrounding rock pressure and internal water pressure, the damage mainly occurs on the vault and the arch foot of the lining. When the depth is less than 620 m, tensile and shear composite failure occur at the arch foot, but when the depth is more than 620 m, tensile failure mainly occurs both on the vault and the arch foot. Figure 11 shows the initiation and failure internal water pressure under different buried depths. As can be seen, the critical internal water pressure increases with the increase of buried depth, and the increase rate of the failure internal water pressure is much larger than the initiation internal water pressure. For the case of shallow buried depth, the lining will fail shortly after the crack initiation, but when the buried depth is deep, the process of crack propagation is relatively longer.

Effect of Different Lining ickness on the Hydraulic
Fracturing of Tunnel Lining. Figure 12 shows the failure modes of the tunnel lining under different lining thicknesses. As can be seen, tensile failure mainly occurs on the vault and the arch foot of the thin tunnel lining, for those lining thicknesses varying from 0.6 m to 0.8 m, tensile and shear composite failure happens on the arch foot. When the thickness is 1 m, only the tensile failure happens both on the vault and the arch foot. Figure 13 shows the initiation and failure internal water pressure under different lining thicknesses. As can be seen, the critical internal water pressure increases with the increase of the lining thickness. What should be noticed is that when the thickness is less than 0.8 m, the differences between the initiation and failure internal water pressure are relatively small, but when the thickness is more than 1 m, the differences increase rapidly. Figure 14 shows the failure modes of the tunnel lining under different surrounding rock mass grades. As can be seen, when the qualities of surrounding rock mass are relatively high (such as grade II), due to the fact that the surrounding rock bears a large part of the pressure, the stress applied to the tunnel lining is relatively small, and tensile failure happens on the vault and the arch foot. For the case that the surrounding rock mass grade is III or IV, the pressure applied to the tunnel lining is relatively larger, and the failure mode transforms from pure tensile failure to tensile and shear composite failure. When the qualities of surrounding rock mass are low, the stress applied to the tunnel lining is extremely high, and the stress concentrates on the arch foot, where the tensile failure mainly happens. Figure 15 shows the initiation and failure internal water pressure under different rock mass grades. As can be seen, when the quality of the surrounding rock mass becomes worse, the critical internal water pressure first increases and then decreases and reaches the largest when the rock grade is IV, which indicates that more attention should be paid to the lining support when the rock mass grade is high or low.

Verification of Numerical Simulation.
Previous studies focused less on the water-stress coupling damage of tunnel lining. Due to the limitations of field survey, the research results mostly focused on numerical simulation, and the tunnel shape is simplified into circles, but the horseshoeshaped tunnel is mostly used in practice. Figure 16 shows the comparisons between our simulation results and the simulation results in [31]. As can be seen, damage occurs on the arch foot, which is similar to the results in [31], and the rationality of the proposed method is verified.
However, previous studies only considered a few conditions, and the different factors affecting lining failure are not fully investigated. In this paper, the Guzeng hydraulic tunnel is taken as an engineering practice, and different factors such as buried depth, lining thickness, and surrounding rock mass grade are taken into consideration. e failure mode and the critical internal water pressure under different conditions are investigated, and the results can provide guidance for the design of the lining support and reinforcements. In this paper, the solid-water particle interaction method and the particle damage conversion algorithm are proposed to realize the hydraulic fracture processes. A series of numerical treatments such as "particle domain searching method," the "birth-and-death particle method," and the "group discrimination searching method" have also been proposed to realize the simulations of complex processes of excavation, lining, and operation processes of hydraulic tunnel. Compared with traditional FEM, the proposed 2P-IKSPH method can dynamically reflect the stress wave propagation processes of surrounding rock and the damage processes of tunnel lining, and no mesh redivision is required. erefore, applications of 2P-IKSPH method into rock mechanics are promising.
What should be stressed is that the external water pressure is not considered in this paper. Reasons are as follows: (1) the application of external water pressure is difficult in 2P-IKSPH method temporarily. (2) e buried depth is shallow, and the influence of the external water is small. So in order to illustrate the feasibility of numerical simulation, the calculation is a little simplified. Meanwhile, Step 4000     the practical rock mechanics engineering is usually a complex 3D problem. It has become a consensus that 2D model can not completely represent 3D. However, the computational efficiency of 3D model is recognized as a major difficulty in numerical simulation. e parallel programming of 2P-IKSPH method for 3D modeling of engineering practice will be the future directions.

Conclusions
(1) e solid-water particle interaction method and the particle damage conversion algorithm are proposed, and the hydraulic fracturing processes are realized. (2) e "particle domain searching method," the "birthand-death particle method," and the "group discrimination searching method" have been proposed to realize the complex processes of excavation, lining, and operation processes. Data Availability e programme 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
e authors declare that they have no conflicts of interest regarding the publication of this paper.