Numerical Study on the Fracture Properties of Concrete Shield Tunnel Lining Segments

Shield tunnel lining structure is usually under very complex loading conditions in the underground space. As a kind of the common concrete structures, any defect in the tunnel lining segment may deteriorate its bearing capacity and even cause severe disasters. Three-dimensional numerical models of shield tunnel lining segments with initial cracks are built using the Symmetric Galerkin Boundary Element Method(SGBEM-) Finite Element Method (FEM) Alternating Method. The cracking load and ultimate load of the tunnel segments are obtained, and crack propagation under fatigue load is also simulated by employing the Paris Fatigue Law. Results show that loading eccentricity has a very large influence on the bearing capacity of the cracked lining segment; the larger the loading eccentricity, the smaller the bearing capacity. Deformation and damage of the lining segment show obvious phases, which consist of the initial crack, cracking stage, steady crack propagation, unsteady crack propagation, and eventual failure.


Introduction
As the construction technology in underground engineering develops rapidly nowadays, the shield tunnel has gradually become the main form of urban metro tunnel. The supporting structure of the shield tunnel is a precast reinforced concrete lining segment, which may contain various flaws such as initial cracks during the pouring of concrete or the assembling of the lining ring. These flaws obviously reduce the effective force-bearing area of the lining segment, thus deteriorating the safety level of the tunnel supporting system. In the worst situation, it may even cause entire failure of the shield tunnel.
Over the past several decades, plenty of numerical and experimental studies have been performed to investigate the mechanical behavior of shield tunnel lining structures. A simplified comparative static computation of an underground metro concrete tunnel lining was shown by Hudoba [1], and different means of assembly and coupling were used. Three lining models were investigated and compared in the paper by Kavvadas et al. [2], which included contin-uous shell without joints, shell with aligned joints (2D joint configuration), and shell with staggered (rotated) joints, which is the most realistic 3D lining model. Yan et al. [3] conducted extended finite element analysis to investigate the cracking and failure characteristics of the segmental lining structure of an underwater shield tunnel upon a derailed high-speed train impact, allowing for the optimization of tunnel design under accidental conditions. Wu et al. [4] investigated the thermal behavior of tunnel segment joints when exposed to fire, and the effectiveness of using concrete-filled steel tubes (CFSTs) to restore joint strength after a fire was also investigated. Wang et al. [5] proposed a novel multiscale modeling method in which potentially damaged and undamaged zones were recognized according to preanalysis results and then simulated at different scales. In their model, mesoscopic features within the potentially damaged zones were modeled using a mesostructure cohesive zone method, while the other zones were simulated as macroscopic homogeneous materials. To investigate the mechanical behavior of the segment joint in a shield TBM tunnel, full-scale joint tests were conducted, and a mechanical model of segment joints has been proposed by Ding et al. [6]. Li et al. studied the development of longitudinal joint opening with bending moment under different axial stress levels, and investigated the longitudinal joint opening in the Ultimate Limit State (ULS), by both experimental and analytical methods [7]. A series of customized model tests were conducted by Zheng et al. to investigate the mechanism underlying tunnel collapse due to the construction of a connecting passage [8]. Liu et al. designed and performed full-scale experiments relating to the structural bearing capacity of stagger joint assembled shield tunnels based on an unloading situation [9], which combined engineering practices relevant to the current urban metro with circumferential disturbances. Based on the field investigation and analysis, Wang et al. compared and analyzed the mechanical characteristics of the segment structure in shield tunnels under the circumstances of different concrete spalling regions by the method of the similarity model experiment [10]. A novel loading setup was developed by Huang et al. [11] for performing full-scale loading tests on "standing" segmental lining of a subrectangular shield tunnel, which for the first time allowed the mechanical behavior of the segmental lining subject to self-weight to be assessed in a full-scale loading test. To study the influence of the dynamic load on underground tunnel stability, a theoretical method to predict underground tunnel behavior considering the peak particle velocity (PPV) and the stress distribution was presented, and the influence of explosion-induced wave from an underground charge tunnel on an adjacent tunnel was explicitly considered in order to evaluate tunnel stability [12]. Tian et al. conducted a prototype load test on a cracked shield tunnel lining segment to study its fracture properties and found that the lining segment with an initial crack finally showed the I-type fracture property, and the failure mode was dominated by the initial crack [13]. Besides, many other researches [14][15][16][17][18][19] were conducted by numerical simulation or in situ test to study the fracture evolution and failure process of underground engineering.
Almost all of these researches are mainly focused on intact tunnel lining structures, while little research has been done on the mechanical behavior and fracture properties of the shield tunnel lining segment with initial flaws, and researches on the propagation law of existing cracks and the failure mode of a cracked tunnel lining structure are rare. Considering the existing body of the literature, there exists the possibility and need for further study on the fracture and fatigue properties of concrete shield tunnel lining structures with various flaws. In this paper, cracked tunnel lining segments under different loading conditions are studied by numerical methods. Three-dimensional numerical models of lining segments with initial cracks are built. The Stress Intensity Factors of the crack front are computed during each step of the crack increment by employing the SGBEM-FEM Alternating Method. The crack growth rates are determined by the Paris Fatigue Law. The crack propagation paths and number of loading cycles are also predicted to estimate the fatigue life of the tunnel lining segment.  [20,21], the Symmetric Galerkin Boundary Element Method (SGBEM) has many advantages. SGBEM has a symmetrical coefficient matrix of the equation system, and there is no need to specially treat sharp corners for numerical simulation. The early work on SGBEM mainly focused on the regularization of hypersingular integrals [22][23][24][25]. The nonhypersingular integral equations for tractions were developed in [26,27], then a systematic way was proposed to obtain the weakly singular symmetric Galerkin boundary integral equations in [28,29]. This method was employed to make analyses on cracked three-dimensional solids with various surface flaws [30,31].

Theoretical Background
As shown in Figure 1, for a domain of interest with source point x and target point ξ, three-dimensional weakly singular symmetric Galerkin boundary integral equations for displacements and tractions are developed as the following.

Geofluids
The displacement boundary integral equation is And the corresponding traction boundary integral equation is In Equations (1) and (2), D a is a surface tangential operator: and H * abq are kernel functions and can be found in reference [31]. These kernel functions are all weakly singular, so the implementation of the current boundary integral equations is very simple and convenient.
As shown in Figure 2, a symmetric system of equations can be finally obtained by applying Equation (1) In the above equations, p, q, and r denote the unknown tractions at S u , unknown displacements at S t ′, and unknown displacement discontinuities at S c , respectively.
Since the coefficient matrix of SGBEM is fully populated, it is unwise to simulate large scale structures using the pure SGBEM. A better way is to make use of the advantages of both SGBEM and FEM. As shown in Figure 3, the uncracked global structure is modeled with FEM, and the cracked local subdomain is modeled with SGBEM. An alternating procedure was conducted to impose the residual stresses at the global and local boundaries. By superposing the solution of each individual subproblem, the solution of the original problem was finally obtained. This method is named as the SGBEM-FEM Alternating Method. By employing the SGBEM-FEM Alternating Method, Tian et al. conducted three-dimensional fracture and fatigue analyses of many typical structure components in civil and mechanical engineering and obtained a series of numerical results with reference values for practical engineering [32][33][34].

Stress Intensity Factors.
In Linear Elastic Fracture Mechanics, Stress Intensity Factors (SIFs) are the most important coefficients to describe the fracture properties of the structure. By employing the SGBEM-FEM Alternating Method, the SIFs of three-dimensional complex cracks can be acquired very conveniently. As shown in Figure 4, the quarter-point singular element is used to capture the singularity at the crack front, which is a specially designed 8noded quadrangular element, with its middle nodes moved 1/4 length of its side to the crack front.
The corresponding SIFs can be captured by the displacements near the crack front, using the following equations: where E denotes the elastic modulus, ν denotes Poisson's ratio, r denotes the distance to the crack front, and u 1 , u 2 , and u 3 denote the discontinuities of displacement at the crack surface.

Simulation of Fatigue-Crack
Propagation. SGBEM is very efficient for modeling cracks and their fatigue propagation. For fatigue-crack propagation, it is unnecessary to use any other special technique such as the Level Sets to depict the crack surface since the crack surface is already efficiently modeled by boundary elements. In every fatigue analysis step, the Stress Intensity Factors along the crack front are obtained by the SGBEM-FEM Alternating Method. Then, the crack propagation is simulated by adding a layer of additional elements at the crack front, in the direction determined by the Eshelby-force vector [35], with the size determined by the Paris Fatigue Law.
As shown in Figure 5, the direction of crack growth is determined by the direction of vector ΔJ, which can be acquired by the following equations:

Geofluids
Based on the Paris Fatigue Law, the fatigue life of the cracked structure is estimated as Thus, at each fatigue-growth step, the maximum value of ΔK along the crack front is computed by the SGBEM-FEM Alternating Method, and Equation (7) is numerically evaluated using the trapezoidal rule.
Indeed, plasticity-induced crack closure can be very useful to account for the effects of plate thickness on crack growth rates. Here, the simple Paris Law is used to predict the fatigue growth rate. However, other models which account for the effect of the thickness can also be incorporated in the framework of the current SGBEM-FEM Alternating Method, which will be our future study. Moreover, some 3D effects of fracture mechanics are also neglected in this study, such as 3D corner singularity and mode II and III coupling, which are expected to have an insignificant effect for the damage tolerance of the cracked structures.

Numerical Simulation
3.1. Geometry of the Specimen and Numerical Model. The present numerical models of shield tunnel lining segments are built based on the prototype shield tunnel of Shanghai Metro Line 12 in China. The supporting structure of the tunnel is a concrete lining ring, which has an outer radius of 3.1 m, inner radius of 2.75 m, and width of 1.2 m. The concrete strength grade is C55, elastic modulus of the concrete is E = 3:55 × 10 10 Pa, and Poisson's ratio ν = 0:2. The initial crack is generated by cutting a straight crack on the inner surface of the lining segment, which is 50 mm in depth and runs through the entire lining segment. The schematic diagram of the loaded specimen is shown in Figure 6.
Since the specimen is symmetrically loaded, only half of it is used to make the force analysis. As shown in Figure 7, the moment at the centroid of the cross-section is where M and N are the bending moment and the axial force on the cross-section of the lining segment, respectively; P is the vertical load; and W is the self-weight of the segment.
The loading eccentricity is defined as

Geofluids
where e denotes the loading eccentricity and M and N are the corresponding bending moment and axial force, respectively.
The numerical models of the tunnel lining segment with initial cracks are built by employing the SGBEM-FEM Alternating Method. As shown in Figures 8 and 9, the global structure with no crack is modeled using finite elements (10noded tetrahedron elements), and the initial crack surface is modeled independently using boundary elements (8-noded quadrilateral elements). Two different numerical models with different loading eccentricities are built in this study. When loading eccentricity e is 0.24 m, the numerical model is named Z1, while the numerical model with e of 0.12 m is named Z2.

Fracture
Properties of the Lining Segment. Double-K fracture criterion of concrete [36] is employed to analyze the fracture properties of the cracked tunnel lining segment. The SIFs along the crack front are captured by the SGBEM-FEM Alternating Method and are constantly compared with the fracture toughness of concrete C55. According to the double-K fracture criterion of concrete, the fracture toughness of concrete C55 is set as K ini Ic = 0:53 × 10 6 Pa ffiffiffiffi m p and K un Ic = 1:42 × 10 6 Pa ffiffiffiffi m p . By simultaneously adjusting the load N and P, to make the computed SIFs reach the concrete's fracture toughness, the initial cracking load N ini and P ini and the ultimate load N un and P un are finally obtained. Here, the cracking load of the specimen is defined as the load that makes the initial crack start to grow. The ultimate load is defined as the largest load that the specimen can bear.
Since it is a mixed-mode fracture problem, the efficient SIF is defined as   For the shield tunnel lining structure, its surrounding pressure usually changes seasonally, such as the water pressure and temperature pressure. Besides, the lining structure will bear the vibration load of the train in the long term. All of these loads will bring fatigue problems. A microcrack may initiate inside the concrete structure under all these fatigue loads, and the existing cracks will propagate. Thus, it is of importance to study the fracture mechanism and fatiguecrack propagation of the shield tunnel lining structure under fatigue loads. Paris Fatigue Law can be used for the fatigue analyses of concrete structures [37]. For the concrete structure with ini-tial cracks, its Paris Fatigue Law is In the above equation, da/dN is the crack growth during the cyclic load step, ΔK is the change of SIF, a 0 is the initial depth of the crack, and D is the thickness of the concrete structure.
The amplitude of fatigue load is set to 0:7N max ; N max is the ultimate load of the specimen for monotonic loading. The stress ratio, which is the ratio of the minimum stress to the maximum stress, is set to 0.1. The numerical model Z1 7 Geofluids is firstly studied under such a fatigue load. Figure 10 illustrates the crack growth along with the load cycles, and the process of fatigue-crack growth is shown in Figure 11.
As shown in Figure 10, after about 34,184 load cycles, the SIF along the crack front of specimen Z1 reaches the fracture toughness of the concrete, and the structure is broken. The increment of the fatigue crack is about 0.102 m. For numerical model Z2, its fatigue-crack growth is shown in Figure 12.
As can be seen from Figures 11 and 12, for the case of small loading eccentricity, the fatigue crack deflects more seriously during the propagation process compared with the case of large loading eccentricity, and it finally grows to a local crack, not as the cut-through crack in Figure 11. So, a simple conclusion can be drawn: for a cracked shield tunnel lining segment, large loading eccentricity is more dangerous than small loading eccentricity for the structure.

Results and Discussion
In Figure 13, the numerical results in this study and the experimental results by Tian et al. [13] on crack propagation in the lining segment under large loading eccentricity are compared, which show good agreement with each other. It can be seen that the initial crack does not propagate strictly along the original direction; in fact, it deflects and extends to the midpoint of the vertical load.
The cracking loads of the tunnel lining segment obtained by experiment [13] and this numerical simulation are listed 8 Geofluids and compared in Table 1; quite good agreement is reached with each other.

Conclusions
In the present paper, three-dimensional numerical models of shield tunnel lining segments with initial surface cracks are built, for which the Stress Intensity Factor analyses and fatigue-crack growth simulations under different load conditions are carried out by employing the three-dimensional SGBEM-FEM Alternating Method. The following conclusions can be drawn: (1) The SGBEM-FEM Alternating Method requires independent and very coarse meshes for both the uncracked tunnel lining segment and the crack surface; it only requires very minimal human-labor cost as well as minimal computational burden to make fracture analyses and model the fatigue-crack propagation in shield tunnel lining segments (2) By employing the SGBEM-FEM Alternating Method, the whole crack growth path up to failure of the lining segment can be easily simulated (3) Loading eccentricity has a very large influence on the bearing capacity of the lining segment; the larger the loading eccentricity, the smaller the bearing capacity of the lining segment. For the case of the small loading eccentricity, the fatigue crack deflects more seriously during the propagation process compared with the case of the large loading eccentricity, and it finally grows to a local crack, not as the cut-through crack for the large loading eccentricity. So, for a cracked shield tunnel lining segment, the large loading eccentricity is more dangerous than the small loading eccentricity of the structure (4) The lining segment with an initial crack mainly shows the I-type fracture property, and its failure mode is dominated by the initial crack

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.