Analysis of Crack Formation and Growth in Tunnel Linings Using Double-K Fracture Criterion

: Empirical criteria and fracture/damage mechanics are used to evaluate the safety of lining cracks in the conventional methods. However, the former lacks a scientiﬁc basis, and the latter requires complicated mechanical calculations. To overcome the above shortcomings, this paper proposes a new method to perform crack analysis of plain concrete linings, based on the double-K fracture criterion. The proposed method uses two crack width indices, i.e., initiation and unstable fracture widths, to divide the fracture process of lining into three stages: initiation stage, stable propagation stage, and instability propagation stage. These two crack width indices are calculated by the equivalent transformation of fracture toughness. Using the proposed criterion, the safety state of the concrete lining can be determined by comparing the ﬁeld measurement width and crack width indices. A speciﬁc code based on the extended ﬁnite element method (XFEM) is developed to simulate the fracture process of concrete lining. Several numerical experiments are carried out to evaluate the proposed fracture criterion. The results show that the two fracture indices of the proposed criterion can accurately identify two demarcation points of the three stages of the lining fracture process, including the nonlinear starting point and the unstable fracture point of the load–displacement curve. Compared with conventional methods, the proposed method uses the geometric parameter to estimate the mechanical state of cracks, so the complicated mechanical calculation can be avoided.


Introduction
Cracks, leakage, spalling, and displacements are the main manifestations of tunnel anomalies, which affect the stability of the tunnel [1,2]. Among them, the lining crack is the most common and adverse tunnel anomaly, which is frequently adopted as the key indicator of tunnel safety. Moreover, the lining crack is also one of the main causes of other tunnel anomalies [3]. In order to avoid the potentially unfavorable impact on traffic caused by the lining crack, the relationship between the existing crack and tunnel safety must be understood.
One of the most commonly used methods for tunnel safety evaluation is based on field measurement results combined with empirical criteria [4], such as the maintenance manual of road tunnels suggested by Japan Road Association [5] and the maintenance technical specifications proposed by the professional standards compilation group of China [6]. However, these safety evaluation methods are semi-quantitative and more often qualitative because they neglect the mechanical essence of crack generation and evolution. Another method used for tunnel safety evaluation adopts the complex variable method, fracture mechanics, or damage mechanics to describe the fracture behavior of lining cracks [7][8][9][10][11]. For instance, Khoteev et al. used fracture mechanics to analyze the crack resistance of a fiber-reinforced concrete tunnel lining [10]. Amorim et al. developed a simplified numerical method based on lumped damage mechanics to model the inelastic behavior of lining cracks [11]. Although these mechanical methods are usually quantitative, they cannot take full advantage of field measurement results, such as the crack width and depth. Furthermore, a large number of experiments show that direct use of linear elastic and elastoplastic fracture theory to describe the fracture behavior of concrete has an obvious size effect [12][13][14][15][16]. For details of the size effect, interested readers are referred to the research work of Bazant [13][14][15][16][17][18][19].
For all of these reasons, analyzing the evolution of cracks in concrete linings is extremely important for tunnel safety evaluation, which should be carried out using concrete fracture mechanics due to the special fracture behavior of concrete. On the other hand, the field measurement results should also be fully utilized in the mechanical analysis of tunnel safety evaluation, because they are the most direct parameters indicating the mechanical state of cracks. Accordingly, this paper proposes a modified double-K fracture criterion to analyze the fracture process of concrete linings. The proposed criterion effectively combines the fracture mechanics of concrete with the field-measured data without complicated mechanical calculations. Moreover, a specific code based on the XFEM-FCM algorithm is developed to simulate the fracture process of concrete lining. The proposed theory of the modified double-K fracture criterion is verified by several numerical experiments using the XFEM-FCM algorithm.

Fracture Characteristics of Concrete
As the main ingredient of tunnel lining, concrete has significantly different fracture characteristics compared with other homogeneous materials. At the material level, concrete is a porous material composed of fine aggregate, coarse aggregate, and cement paste. Due to the poor strength of the bonding interface between aggregate and cement paste, a large number of microcracks are generated in the bonding interface before concrete cracking. At the macro level, the microcracks are clustered at the front end of the crack (physical crack tip) and form the so-called fracture process zone, which plays an important role in the fracture process of concrete [20][21][22].
The composition of concrete cracks is plotted in Figure 1, where the fracture process zone is located at the front of the physical crack tip. Due to a large number of micro-cracks clustered in the fracture process zone, the stiffness and strength of concrete are greatly reduced. This phenomenon is usually called stress softening. The stress softening behavior of concrete can be expressed by different stress-crack width softening curves, such as linear softening curves, bilinear softening curves, and exponential softening curves. When the stress at the tip of the fracture zone (fictitious crack tip) reaches tensile strength, the fracture process zone will further expand. When the open width of the physical tip increases to the critical value, the fracture process zone near the physical tip cannot transmit stress and fracture completely into free fractures [21,23]. A typical load-crack width curve (P-δ curve) is plotted in Figure 2. Generally, the fracture process of concrete consists of three stages, i.e., crack initiation, stable propagation, and unstable fracture, which correspond to the linear segment (OB), nonlinear segment (BC), and softening segment (CE) in Figure 2, respectively. In the crack initiation stage, the size of the fracture process zone can be neglected, and the crack width increases linearly with the increase of the load P. In the stable propagation stage, the fracture process zone starts to extend steadily and the nonlinear characteristics of the P-δ curve begin to appear. The approaches of linear elastic fracture theory are no longer available at this stage. In the unstable fracture stage, the free crack starts to cause rapid extension, and the load-bearing capacity of concrete decays drastically [24].

Basic Theory of the Double-K Fracture Criterion
A large number of experiments show that direct use of linear elastic or elastoplastic fracture theory to describe the fracture behavior of concrete has an obvious size effect [12,26], due to the existence of the fracture process zone. Over the past few years, many fracture theories have been developed to describe the fracture behavior of concrete, including the fictitious crack model by Hillerborg et al. [21], crack band model by Bazant et al. [17], two-parameter fracture criterion by Jenq et al. [27], and double-K fracture criterion by Xu et al. [24,25,28]. Since the double-K fracture criterion has physically meaningful parameters and simple closed analytic solutions, it is employed to analyze the fracture behavior of tunnel lining in this paper. A brief introduction to the double-K fracture criterion is provided as follows: The double-K fracture criterion introduces the following two basic hypotheses [25]: (1) The nonlinear characteristic of the P-δ curve is caused by the fracture process zone; (2) The fracture analysis of concrete is calculated by the elastic equivalent method, where the equivalent crack is composed of a free crack and an equivalent fictitious crack.
In the double K fracture criterion, the fracture process zone is equivalent to an elastic crack and the stress intensity factor is calculated by linear elastic fracture mechanics with an equivalent crack. Using the three-point bending beam in Figure 3 as an example, the detailed calculation process is as follows [25]: (1) Measure the crack width δ and the load P; (2) Calculate the equivalent crack depth a e by solving the following linear elastic fracture equation: where l, h, and b are the length, height, and width of the specimen, respectively; E is the modulus of elasticity; and V(·) is the shape function.
(3) Calculate the stress intensity factor by the following linear elastic fracture formula: where F(·) is the shape function.
Since the fracture process of concrete has three distinct stages, the double-K fracture criterion introduces two fracture toughness parameters, i.e., initiation toughness K ini and unstable fracture toughness K un , to evaluate the mechanical state of cracks. The corresponding fracture criterion is as follows [24]: (1) When K < K ini , crack initiation stage: the crack width increases linearly with the load; (2) When K ini ≤ K < K un , stable propagation stage: the fracture process zone extends steadily; (3) When K ≥ K un , unstable fracture stage: the crack develops drastically.
The two fracture toughness can be calculated by the experiment of the three-point bending beam, as shown in Figure 3. The specific calculation formulas are given by: where P ini is the load of the nonlinear starting point of the P-δ curve, a 0 is the depth of the free crack, P un is the load of the unstable fracture point of the P-δ curve, and a e is the depth of the equivalent crack [24,25]. From the above analysis, it can be seen that the external load and crack width should be obtained before applying the double K fracture criterion to evaluate the mechanical state of the lining cracks. The crack width can be measured by specific instruments, but the external loads acting on the lining are almost impossible to obtain. On the other hand, as an underground structure, the tunnel has a complicated interaction with the surrounding rock, and its structural form is extremely complex, so it is quite difficult to calculate the stress intensity factor using the double-K fracture theory. In view of these problems, the double-K fracture criterion should be modified to adapt it for applications in tunnel engineering.

Modified Double-K Fracture Criterion for Cracks in Tunnel Linings
If the whole lining is taken as the analysis object, it is difficult to evaluate the mechanical state of cracks by the analytical method. To overcome this, a small segment of lining containing a crack is separated from the whole model, as shown in Figure 4. This lining segment is subjected to the bending moment M, axial force N, shear force S, and external load F. Ignoring the effect of the shear force and external load on the crack, the crack width δ and the stress intensity factor K can be given by the following equations [29]: The shape functions f M (·), f N (·), v M (·), and v N (·) are given as follows: where b is the longitudinal length of the lining (b = 1), h is the lining thickness, E is the elastic modulus, a e is the equivalent crack depth, and s = a e /h. The double-K fracture criterion adopts the stress intensity factor calculated by the elastic equivalent method to evaluate the mechanical state of cracks. Although use of the above formulas to calculate the stress intensity factor is convenient, the axial force N and bending moment M still cannot be measured accurately. On the other hand, as a statically indeterminate structure, the extension of cracks will change the stiffness of the lining, which results in redistribution of the internal force. Therefore, the stress intensity factor is not suitable for the mechanical analysis of lining cracks.
As field monitoring of lining cracks usually only includes the width and depth, it is more appropriate to take the crack width as a mechanical state parameter to determine the fracture stage of the lining. To this end, the initiation toughness K ini and unstable fracture toughness K un presented in the double-K fracture criterion should be transformed into two crack width indices, i.e., initiation width δ ini and unstable width δ un . The transformation method is described in detail below.
Plug Equation (6) into Equation (5), the relationship between crack width and stress intensity factor can be expressed as: where e is the eccentricity, e = M/Nh. The initiation and unstable widths can be obtained respectively by substituting the initiation toughness and unstable fracture toughness with the corresponding equivalent crack depths into the above formula.
For the initiation width, the size of the fracture process zone can be neglected, so the equivalent crack depth is equal to the free crack depth a 0 . It takes the form: For the unstable width, the equivalent crack depth can be obtained by solving the following equation: where the shape functions F(·) and stress distribution σ(·) in the fracture process zone are given as follows, respectively: where u = x/a e , s = a e /h, t = a 0 /a e , f t is the tensile strength of concrete; σ s is the cohesive force at the physical crack tip; and w is the crack width at the physical crack tip. K c is the additional fracture toughness caused by the cohesive force in the fracture process zone, namely, cohesion toughness. It has the following relations with K ini and K un [24,25]: Since the integral Equation (13) is somewhat complicated, the equivalent crack depth a e for an unstable width may only be calculated using the numerical method, such as dichotomy. After obtaining the equivalent crack depth, the unstable width takes the form: Eventually, the detailed calculation steps of the modified double-K fracture criterion for lining cracks are as follows: (1) Measure the crack width δ and the free crack depth a 0 ; (2) Obtain two fracture toughness values and then transform them into the crack width indices; (3) Evaluate the mechanical state of the lining cracks by the following criteria: When δ < δ ini , crack initiation stage: the crack width increases linearly with the load; (b) When δ ini ≤ δ < δ un , stable propagation stage: the fracture process zone extends steadily; (c) When δ ≥ δ un , unstable fracture stage: the crack develops drastically.

Remark 1.
Although the above theory is proposed for the I mode crack, it is still applicable for the I-II mixed mode crack if the crack width δ adopts the vector sum of crack opening displacement δ I and slip displacement δ II , as shown in Figure 5. This method was first proposed by Jenq and Shah for the application of two-parameter fracture criterion in the I-II mixed mode crack [30,31].

Remark 2.
Although the above derivation is based on the three-point bending beam, this method is also applicable to other shapes of concrete specimens. In practical application, the expression of crack width and stress intensity factor in Equations (5)-(10) needs to be replaced by those of other specimens.

Numerical Experiments of the Fracture Behavior of Tunnel Linings
Owing to the existence of the fracture process zone, the fracture behavior of concrete linings become extremely complicated. In this section, the extended finite element method (XFEM) and fictitious crack model (FCM) are adopted to simulate the complicated fracture behavior of concrete lining. A specific code is developed using Matlab Language to perform calculations with the XFEM and FCM. Furthermore, the proposed theory of modified double-K fracture criterion is verified by analyzing the variation of the crack width and depth during the whole fracture process.
The XFEM is a partition of the unity-based method, where the enrichment functions describing the discontinuity on both sides of a crack are incorporated locally into the standard finite element method. The geometry of internal cracks is separated from the finite element mesh, so there is no need to regenerate the mesh as the crack expands. This feature makes XFEM a powerful tool for analyzing crack propagation [32][33][34][35].
The FCM divides the concrete crack into two parts: physical crack and fictitious crack. The physical crack represents the traction-free crack, and the fictitious crack simplifies the fracture process zone into a separate line crack. Applying the stress softening relation of concrete to the fictitious crack can effectively calculate the complicated fracture behavior of concrete linings caused by the fracture process zone [21]. For more details on XFEM and FCM, the reader is referred to the following references [21,[36][37][38].

Basic Parameters of the Numerical Experiments
In this section, a plain concrete lining of a double-track railway tunnel in China is considered as the research object as shown in Figure 6. The material parameters of rgw concrete lining used in the numerical experiment are listed in Table 1, in which the elastic modulus, Poisson's ratio, tensile strength, and maximum crack opening width are derived directly from the experiments in [39]. In order to obtain the initiation toughness and unstable fracture toughness, a numerical experiment of the three-point bending beam with the above material parameters is carried out using the XFEM-FCM algorithm, and the resulting P-δ curve is plotted in Figure 7. According to the double-K fracture criterion introduced in Section 2.2, two fracture toughness values are calculated and listed in Table 1.   Two numerical experiments are performed to study the crack propagation mechanism of the lining and demonstrate the feasibility of the modified double K fracture model. The interaction between the concrete lining and geotechnical material is simulated by a Winkler foundation spring, where the spring stiffness is 150 MPa/m. In the first numerical experiment, the arch crown is subjected to an additional load and a 100 mm initial crack is reserved on the arch crown. In the second experiment, the arch waist is subjected to an additional load and a 100 mm initial crack is reserved on the arch waist. In the whole numerical experiment, the additional load will gradually increase until the crack completely penetrates the lining. The maximum value of the additional load is about 1.5 MPa. The quadrilateral elements are used in both numerical experiments, where the global mesh size in the radial direction is 20 mm and that in the circumferential direction is about 180 mm. In order to improve the calculation accuracy of crack propagation, a fine finite element mesh is used near the crack, where the local mesh size is about 10 mm both in the radial and circumferential directions. The schematic diagrams of the calculation model are plotted in Figure 8.  Figure 9 shows the propagation paths of the cracks at the arch crown and arch waist under the additional load. It is seen that the crack at crown is the I mode crack while the crack at waist is the I-II mixed mode crack. The relationship between the crack depth and eccentricity at the cross-section of the lining crack is shown in Figure 10. It is illustrated that the arch crown is subjected to an internal force of eccentric tension while the arch waist is subjected to an internal force of eccentric compression. The two groups of the resulting eccentricity increase with the additional load at first, and then decrease with the additional load when the additional load reaches a certain value, which is due to the decrease of the bending stiffness caused by crack propagation. In the initial stage of crack propagation, the decrease of the bending stiffness is negligible, so the eccentricity increases with the additional load. With the increasing crack depth, the bending stiffness decreases sharply, so the eccentricity begins to decrease with the additional load. Since the lining is a statically indeterminate structure, the external force acting on it is not equal to that acting on the crack. For convenience, let us introduce an equivalent load to characterize the force acting on the lining crack:

Analysis of the Fracture Process of Tunnel Linings
This equivalent load is the weighted sum of the maximum stress caused by the bending moment and the axial force. The curves of the equivalent load-crack width and equivalent load-crack depth are respectively shown in Figures 11 and 12. The vertical dashed lines in Figure 11 represent two crack width indices calculated by the proposed method, and the horizontal dashed lines in Figure 12 represent the corresponding equivalent loads. The solid-line curves in Figures 11 and 12 represent the results obtained by the finite element method with the XFEM-FCM algorithm. It is seen that the proposed width indices accurately divide the fracture process of the lining into three stages. In the crack initiation stage, the crack depth is almost not developed, and the crack width increases linearly with the equivalent load. The crack initiation equivalent loads of the arch crown crack and arch waist crack are 0.68 and 0.51 Mpa, respectively. In the stable propagation stage, the fracture process zones of the arch crown crack and arch waist crack extend 59 and 54.6 mm, respectively, which are on the same order of magnitude as the free crack depth. Meanwhile, both load-width curves present the intense nonlinear feature. In the unstable fracture stage, the equivalent loads acting on the lining cracks begin to decrease, but the crack width and depth increase dramatically. It implies that the lining reaches the damage stage. The unstable fracture equivalent loads of the arch crown crack and arch waist crack are 1.30 and 1.17 Mpa, respectively.  Overall, the proposed theory of the modified double-K fracture criterion can accurately identify two demarcation points of the three stages of the lining fracture process, including the nonlinear starting point and the unstable fracture point. Therefore, use of the proposed criterion is feasible when analyzing the safety of tunnel cracks.

Conclusions
Due to the complexity of geological formations, highly indeterminate characteristics of the surrounding soils/rocks, and the complicated interaction between the lining and surrounding soils/rock, the external load acting on the tunnel is usually difficult to determine, which results in great difficulty regarding the fracture analysis of lining. This paper presents a novel analysis method, namely, the modified double-K fracture criterion, for analysis of the fracture mechanism of plain concrete linings. The proposed criterion uses two crack width indices to divide the fracture process into three different stages, i.e., crack initiation, stable propagation, and unstable fracture. The crack width indices including initiation and unstable widths are derived from the equivalent transformation of the fracture toughness. The current fracture stage of the lining can be determined by comparing the field measurement width and two crack width indices. Since this method uses the geometric parameter of the crack width to estimate the mechanical state of cracks, it can avoid complicated mechanical calculations.
A specific code based on the XFEM-FCM algorithm was developed to simulate the fracture process of concrete lining. The existence of three stages of the lining fracture process was confirmed by numerical experiments and the feasibility of the proposed theory was also verified.

Conflicts of Interest:
The authors declare no conflict of interest.