Damage Detection and Evaluation for an In-Service Shield Tunnel Based on the Monitored Increment of Neutral Axis Depth Using Long-Gauge Fiber Bragg Grating Sensors

It is difficult to detect and evaluate the structural damage in a shield tunnel during operation because many traditional techniques based on the observation of vibrations are limited in daily monitoring in tunnels. Thus, the curvature radius of a static longitudinal settlement curve is used to identify the residual health and safety of an in-service shield tunnel. However, there are still two problems. The curvature radius is suitable for a qualitative judgment rather than a quantitative evaluation for longitudinal damage detection. Moreover, the curvature radius, which is calculated from the measured settlements of three neighboring points, gives an average damage degree in a wide scope only and is difficult to use to identify the damage’s precise location. By means of the analysis of three kinds of longitudinal failure modes in a shield tunnel, this paper proposes: (1) a damage detection method based on the monitored increment of the neutral axis depth; and (2) an index to evaluate longitudinal damage. The index is composed of the residual ratios of the equivalent flexural stiffness (HFM1) and the equivalent shear stiffness (HFM3). The neutral axis position and the proposed damage index can be determined using long-gauge Fiber Bragg Grating sensors. Results from numerical simulations show that the deviation between the HFM1 and the true value residual ratio of the equivalent flexural stiffness is no more than 1.7%. The HFM3 is equal to its true value in the entire damage process. A loading experiment for a scaled-down model of a shield tunnel using long-gauge Fiber Bragg Grating sensors indicated that the errors in the HFM1 were no more than 5.0% in the case of early damage development (HFM1 ≥ 0.5). The maximum error did not exceed 9.0% even under severe damage conditions in the model. Meanwhile, the HFM3 also coincided with its true value in the entire testing process.


Introduction
Due to its safety, high efficiency in construction, and low impact on the urban environment, many metro tunnels have been built by the shield-driven method to satisfy the increasing demand for public transportation in many metropolises. However, some shield tunnels have suffered abrupt damage with severe consequences in recent years. Research points out that a shield tunnel may gradually deteriorate due to many potential risks [1,2], such as groundwater leakage [3,4], nearby excavation A shield tunnel in the longitudinal direction is typically simplified using two models, namely the beam-spring model [45] and the longitudinal continuous model [46]. In the beam-spring model, the segmental ring is viewed as a set of short beams, and the circumferential joint is replaced by axial, shear, and rotational spring elements. In the longitudinal continuous model, the tunnel is simplified to be a homogenous beam with reduced stiffness. The longitudinal continuous model has been widely applied to analyze the internal forces and deformation in shield tunnels based on the Euler-Bernoulli beam theory [47,48]. The homogenous continuous model is more practical to use in a theoretical analysis than the beam-spring model because the values of rigidities of springs in a joint that are required for the latter model are usually hard to determine.
The deformation of a shield tunnel in the longitudinal direction can be divided into two separate modes as shown in Figure 1 [9,49], and this can be done on the basis of the practically measured settlement data from long-term monitoring [8,9]. One mode is the bending mode, where the compression of the concrete segments occurs on one side, and the tension on the other side causes the joints to open and the bolts there to experience tension. The other mode of deformation is the dislocation mode. In this mode, the accumulation of dislocation between adjacent rings leads to differential settlement, and the circumferential joint mainly experiences a shear force. In the dislocation mode, uniform shear strain appears on the entire section of the joint. According to the mentioned analyses, Wu [50] proposed a Timoshenko beam simplified model (TBSM) to simplify the tunnel as a continuous Timoshenko beam with equivalent flexure stiffness (EI)eq and shear stiffness (κGA)eq. The TBSM can describe the dislocation of a shield tunnel in theory; so, this model is used as the basis for the subsequent derivation of the proposed damage index in this paper.
Sensors 2018, 18, x FOR PEER REVIEW 4 of 24 stiffness (κGA)eq. The TBSM can describe the dislocation of a shield tunnel in theory; so, this model is used as the basis for the subsequent derivation of the proposed damage index in this paper. Consideration of the two deformation modes illustrated in Figure 1 leads to the reasonable assumption of three different kinds of failure mode, denoted FM1, FM2, and FM3. FM1 is the bending failure mode, where yielding occurs from the outermost bolt to the inner ones successively in the tensile region before concrete crushing occurs in the compressive region. FM2 is the shearing failure mode, where all of the bolts in the circumferential joint rupture abruptly due to a shearing action. FM3 represents the decrease in equivalent shear stiffness as a result of bolt failure caused by an excessive bending moment applied to the joint. The failures in FM1 and FM3 are ductile while the failure in FM2 is brittle. Thus, the failure of steel bolts in FM1 and FM3 occurs over a relatively long time, and the rupture of steel bolts in FM2 is abrupt.

Damage Index Derivation
The first step in structural damage identification is to divide the entire structure into several independent, similar elements. An element comprises two adjacent half rings and the circumferential joint between them as shown in the range of l s in Figure 2. This element can be further divided into two parts: the area influenced by the joint, and the uninfluenced area [51]. The area influenced by the joint has a tensile force that is resisted by the bolts, and the compressive force is withstood by the concrete segment. In the uninfluenced area, both the tensile force and the compressive force are borne by the concrete segment. Consideration of the two deformation modes illustrated in Figure 1 leads to the reasonable assumption of three different kinds of failure mode, denoted FM1, FM2, and FM3. FM1 is the bending failure mode, where yielding occurs from the outermost bolt to the inner ones successively in the tensile region before concrete crushing occurs in the compressive region. FM2 is the shearing failure mode, where all of the bolts in the circumferential joint rupture abruptly due to a shearing action. FM3 represents the decrease in equivalent shear stiffness as a result of bolt failure caused by an excessive bending moment applied to the joint. The failures in FM1 and FM3 are ductile while the failure in FM2 is brittle. Thus, the failure of steel bolts in FM1 and FM3 occurs over a relatively long time, and the rupture of steel bolts in FM2 is abrupt.

Damage Index Derivation
The first step in structural damage identification is to divide the entire structure into several independent, similar elements. An element comprises two adjacent half rings and the circumferential joint between them as shown in the range of l s in Figure 2. This element can be further divided into two parts: the area influenced by the joint, and the uninfluenced area [51]. The area influenced by the joint has a tensile force that is resisted by the bolts, and the compressive force is withstood by the concrete segment. In the uninfluenced area, both the tensile force and the compressive force are borne by the concrete segment. independent, similar elements. An element comprises two adjacent half rings and the circumferential joint between them as shown in the range of l s in Figure 2. This element can be further divided into two parts: the area influenced by the joint, and the uninfluenced area [51]. The area influenced by the joint has a tensile force that is resisted by the bolts, and the compressive force is withstood by the concrete segment. In the uninfluenced area, both the tensile force and the compressive force are borne by the concrete segment. The degree of residual fitness of a damaged shield tunnel is defined by the residual ratios of H FM1 -H FM3 that correspond to FM1-FM3, respectively. The superscripts ud and d denote that the variables used are in the undamaged state or the damaged state, respectively.
According to the TBSM, H FM1 is defined as follows: where EI is the flexural stiffness in the longitudinal direction and η is the reduction factor of the flexural stiffness. The value of η d ranges from 0.1 to 0.2 depending on the form of the segment assembly [51]. The value of H FM2 is expressed as: Note that H FM2 is 0 when the circumferential joint is completely damaged by the shearing action, while an H FM2 value of 1 means that the joint is undamaged. The discontinuous nature of Equation (2) reflects that damage to the circumferential joint occurs suddenly and that FM2 is brittle. The expression for H FM3 is as follows: where G b is the shear modulus of a bolt, G s is the shear modulus of the concrete, A b is the cross-sectional area of the bolt, A s is the cross-sectional area of the segment, l b is the length of a bolt, l s is the length of the segmental ring, n ud is the number of bolts in the circumferential joint, n d is the number of un-yielded bolts in the circumferential joint, κ b is the Timoshenko shear coefficient of the bolt and is 0.9 for a circular cross-section [50], and κ s is the Timoshenko shear coefficient of the segmental ring and is 0.5 for an annular cross-section [50]. The expression of α is as follows: where E b is the elastic modulus of the steel bolt, E s is the elastic modulus of the concrete segment, µ b is the Poisson's ratio of the steel bolt, and µ s is the Poisson's ratio of the concrete segment. It is noted that α is a constant for a given tunnel.
, and A s /A b ≥ 1200. Therefore, α is approximately 90-130, and it is much larger than n ud .
The following simplified approximation of Equation (3) can be used when the structural damage is not severe: The damage index D is defined as a combination of H FM1 -H FM3 , according to: The probability of FM2 occurring is generally small, so H FM2 is usually equal to 1. Substituting Equation (1) and Equation (5) into Equation (6) with the assumption that H FM2 is equal to 1 gives the simplified expression of D as follows: Two main conclusions can be drawn from the damage index D derived from Equations (1)-(7). The first conclusion is that D is mainly related to the state of the bolts and is not affected by the concrete segment. Thus, the primary goals when monitoring a shield tunnel should be to observe whether the bolts have yielded, and, if any bolts have yielded, to determine the number of yielded bolts. The second conclusion is that the positions of the neutral axis and yielding line relative to the cross-section of the circumferential joint are important to ascertain. The former is used for the calculation of η d /η ud , and the latter is applied to determine n d /n ud .

Damage Detection and Evaluation Using LFBG-Based Strain Measurements
The derivation process in this section is based on the following assumptions: (i) the strain distributions on the circumferential joint and concrete segment are approximately the same as those in the plane-section; (ii) the tensile force in the joint area is resisted entirely by the bolts, and the compressive force in the same area is withstood by the concrete segment; (iii) the pretension force in the bolts and the axial force in the longitudinal direction of the shield tunnel are not considered; (iv) the second failure mode (FM2) is not taken in account; and (v) the bolts are distributed uniformly in the circumferential joint and all bolts have the same mechanical and geometrical properties.

Relationship between Structural Damage Level and NAD
The neutral axis is a horizontal line within each cross-section of a beam where the normal stress and strain is zero. Figure 3 shows that the maximum tensile stress σ ud t and the maximum compressive stress σ ud c on the cross-section of a rectangular beam are less than the yield stress σ y in the elastic stage. The NAD h ud is constant while the bending moment applied to the cross-section increases. However, when the cross-section is in the plastic stage, σ d t reaches σ y , the NAD increases, and h d is greater than h ud . Thus, the NAD is sensitive to structural damage. The neutral axis moves towards the compressive region during the development of structural damage. Therefore, the damage level and NAD correspond in beam-like structures, and the damage level can be determined quantitatively by measuring the NAD in practice.
The NAD of the area influenced by the joint has been shown to be distinct from that of the uninfluenced area of the element [46,50]. However, the tensile strain is concentrated mainly at the circumferential joint rather than the segments, so it is suitable to choose the area that is influenced by the joint as the monitoring area for the NAD rather than the uninfluenced area. The influencing factor λ of the joints ranges from 0.45 to 0.65 for different forms of the segment assembly [51]. Consequently, if the l b is assumed to be about 60 cm, the length of the area influenced by the joint is greater than 20-30 cm and exceeds the gauge length of an electrical resistance strain gauge. A solution to this problem is to replace the electrical resistance strain gauge by an LFBG sensor. the elastic stage. The NAD h is constant while the bending moment applied to the cross-section increases. However, when the cross-section is in the plastic stage, σ d t reaches σy, the NAD increases, and h d is greater than h ud . Thus, the NAD is sensitive to structural damage. The neutral axis moves towards the compressive region during the development of structural damage. Therefore, the damage level and NAD correspond in beam-like structures, and the damage level can be determined quantitatively by measuring the NAD in practice. The NAD of the area influenced by the joint has been shown to be distinct from that of the uninfluenced area of the element [46,50]. However, the tensile strain is concentrated mainly at the circumferential joint rather than the segments, so it is suitable to choose the area that is influenced by the joint as the monitoring area for the NAD rather than the uninfluenced area. The influencing factor λ of the joints ranges from 0.45 to 0.65 for different forms of the segment assembly [51]. Consequently, if the lb is assumed to be about 60 cm, the length of the area influenced by the joint is greater than 20-30 cm and exceeds the gauge length of an electrical resistance strain gauge. A solution to this problem is to replace the electrical resistance strain gauge by an LFBG sensor.  [40] illustrates the structure of an LFBG sensor, the gauge length of which ranges from several centimeters to several meters. The most remarkable feature of this sensor is that the FBG is sleeved and fixed at both ends of a hollow tube that is protected by the surrounding fiber-reinforced plastic (FRP) package. Thus, the sensing length can be predetermined according to different requirements. Figure 5 shows a comparison of the strain measurement from a traditional short-gauge strain sensor and the LFBG sensor. The short-gauge stain sensor may be ruptured due to the large tensile strain caused by the joint opening. The LFBG sensor can measure the strain safely after the joint opens because the sensor has point fixation only at the two ends bonded to the

Introduction to LFBG Sensors
Figure 4 [40] illustrates the structure of an LFBG sensor, the gauge length of which ranges from several centimeters to several meters. The most remarkable feature of this sensor is that the FBG is sleeved and fixed at both ends of a hollow tube that is protected by the surrounding fiber-reinforced plastic (FRP) package. Thus, the sensing length can be predetermined according to different requirements. Figure 5 shows a comparison of the strain measurement from a traditional short-gauge strain sensor and the LFBG sensor. The short-gauge stain sensor may be ruptured due to the large tensile strain caused by the joint opening. The LFBG sensor can measure the strain safely after the joint opens because the sensor has point fixation only at the two ends bonded to the concrete surface. Thus, the LFBG sensor measures a uniform strain distribution within its gauge length and reduces the risk of rupture when there is excessive joint opening.

Damage Evaluation Using LFBG-Based Strain Measurements
The NAD, and therefore damage, can be determined using an array of LFBG sensors. A series of LFBG sensors, labeled as S1-Sn, can be fixed on the inner surface of a tunnel as shown in Figure 6a,b. Figure 6c shows that the method of least squares can be used to determine the best-fit of the straight line y = kε + c based on the height coordinate of S1-Sn and the corresponding strain measurements

Damage Evaluation Using LFBG-Based Strain Measurements
The NAD, and therefore damage, can be determined using an array of LFBG sensors. A series of LFBG sensors, labeled as S1-Sn, can be fixed on the inner surface of a tunnel as shown in Figure 6a,b. Figure 6c shows that the method of least squares can be used to determine the best-fit of the straight line y = kε + c based on the height coordinate of S1-Sn and the corresponding strain measurements ε1-εn. The parameter c in the fitted equation is equal to the NAD and is denoted as χ. The value εy is defined as the yield strain of the bolts. If ε is substituted by εy in the best-fit equation, the calculated 125 250 500

Damage Evaluation Using LFBG-Based Strain Measurements
The NAD, and therefore damage, can be determined using an array of LFBG sensors. A series of LFBG sensors, labeled as S1-Sn, can be fixed on the inner surface of a tunnel as shown in Figure 6a,b. Figure 6c shows that the method of least squares can be used to determine the best-fit of the straight line y = kε + c based on the height coordinate of S1-Sn and the corresponding strain measurements ε 1 -ε n . The parameter c in the fitted equation is equal to the NAD and is denoted as χ. The value ε y is defined as the yield strain of the bolts. If ε is substituted by ε y in the best-fit equation, the calculated y is the height coordinate γ of the yield line as shown in Figure 6c. It means that the bolts above the yield line have not yielded, but the bolts with positions lower than the yield line have yielded. Thus, the parameter n d can be obtained easily. In addition, it may be required that a few LFBG sensors, sometimes associated with some other distributed sensing technique [52,53], are placed in one monitoring position in practice to not only ensure accurate measurements for a long time but also eliminate the effect of some unpredicted factors, such as the shrinkage of concrete on the fiber and local compression on the Bragg sensor. The neutral axis is a horizontal line, so the central angle φ corresponding to the neutral axis is: where R and r are the outer radius and inner radius of the segment ring, respectively. Similarly, the central angle ω corresponding to the yield line is obtained as follows: These measurements allow for the determination of the important parameters used in Equation (7). First, the undamaged state is considered. According to the TBSM, η ud is calculated by Xu [51]: where Es is the elastic modulus of the concrete segment, and Is is the area moment of inertia of the cross-section of the concrete segment ring. The coefficient of rotational stiffness ζ ud for a circumferential joint is determined from Shiba [46]: The following theoretical method for the calculation of φ ud was proposed by Xu [51]: where Kj1=n ud kj1 and kj1 is the translational stiffness of the bolt at the joint. The value φ ud remains constant in the longitudinal direction of a non-destructive shield tunnel. The calculated φ ud is used to determine ζ ud and η ud . All of these values are independent of the strain measurements. Next, damage evaluation using the long-gauge strain measurements is examined. Shiba [46] proposed the expression of (EI) d eq as follows: where M is the moment applied on the segment ring, and Ny is the ultimate axial force of the segment ring. The parameter R1 in Equation (13) is given by: Therefore, η d can be obtained from: The equilibrium equation and compatibility equation at a joint are solved simultaneously, and the solution was determined by Shiba [46] to be: The neutral axis is a horizontal line, so the central angle ϕ corresponding to the neutral axis is: where R and r are the outer radius and inner radius of the segment ring, respectively. Similarly, the central angle ω corresponding to the yield line is obtained as follows: These measurements allow for the determination of the important parameters used in Equation (7). First, the undamaged state is considered. According to the TBSM, η ud is calculated by Xu [51]: where E s is the elastic modulus of the concrete segment, and I s is the area moment of inertia of the cross-section of the concrete segment ring. The coefficient of rotational stiffness ζ ud for a circumferential joint is determined from Shiba [46]: The following theoretical method for the calculation of ϕ ud was proposed by Xu [51]: where K j1 = n ud k j1 and k j1 is the translational stiffness of the bolt at the joint. The value ϕ ud remains constant in the longitudinal direction of a non-destructive shield tunnel. The calculated ϕ ud is used to determine ζ ud and η ud . All of these values are independent of the strain measurements. Next, damage evaluation using the long-gauge strain measurements is examined. Shiba [46] proposed the expression of (EI) d eq as follows: where M is the moment applied on the segment ring, and N y is the ultimate axial force of the segment ring. The parameter R 1 in Equation (13) is given by: Therefore, η d can be obtained from: The equilibrium equation and compatibility equation at a joint are solved simultaneously, and the solution was determined by Shiba [46] to be: where the parameter R 2 is calculated by: where K j2 = n ud k j2 . k j2 is the plastic stiffness of the bolts at a joint. Substituting Equation (15) into Equation (16), the following expression of η d is obtained: The value of η d in Equation (18) depends on the ϕ and ω as determined from Equation (8) and Equation (9) because R 1 and R 2 are related only to the mechanical and geometrical characteristics of the bolts. Thus, η d is determined by the long-gauge strain measurements shown in Figure 5.
The calculated η ud , η d , n ud , and n d are substituted into Equation (7) to determine the value of the damage index D.

Numerical Simulation Verification
A numerical model was used to simulate the shield tunnel in the longitudinal direction and verify the applicability of the damage index D for damage detection and evaluation. The decrease in the longitudinal stiffness is compared to the damage index D calculated from the long-gauge strain measurements in the model. Verification using this model is concerned with three main aspects. First, the plane-section assumption should be shown to be applicable before and after damage, and that the NAD is sensitive enough to indicate the damage. Second, the theoretical ϕ ud should be close to its true value at the undamaged state. Finally, the monitored damage level from the LFBG sensors should approximate the true damage level.

Numerical Model
The numerical segment ring built in the ANSYS software is based on a real segment ring in the Shanghai Metro Line No. 1. The mechanical properties of the segment and bolt are as follows:  Figure 6, the segment ring is composed of one Block K (22.5 • ) and five standard Block As (67.5 • ). As shown in Figure 7, there are 20 rings assembled by segments with straight joints. A ring contains 16 circumferential weld bolts (M24) and 12 longitudinal joint bolts (M27). The value of l b is 0.6 m. The two half-segment rings at both ends are not considered, and the model is uniformly divided into 19 elements denoted E1-E19.
The concrete segments are simulated using the Solid65 element. The shear transfer coefficient is set to 0.5 and 0.9 for crack opening and crack closing, respectively. The bolts are simulated by using the Solid185 element. In this model, the two ends of each bolt are embedded in the corresponding segments. Contact elements are used to simulate the contact and friction between the contacting surfaces. The end surfaces of each bolt are set as the target surface in the TARGE170 contact element, and the surfaces of the bolt holes are set as the contact surface in the CONTA173 contact element. The friction coefficient between the contacting surfaces is set to 0.3. Furthermore, the surface-to-surface contact is simulated between the contacting end surfaces of adjacent segments. One end surface is set as the target surface in the TARGE170 contact element, and the other is set as the contact surface in the CONTA173 contact element. The friction coefficient between the contacting surfaces is set to 0.7. The penetration coefficient in each contact element is set to 1. Further details of the contact settings in this model are given by Zhong [54].

Sensor Placement and Loading Mode
The LFBG sensors, denoted S1-S9, are fixed on the inner surface of the area influenced by the joint in E6 (one-third span) as shown in Figure 8. This figure also gives the height coordinates of S1-S9 by setting the horizontal axis at the center of the cross-section of the ring. The bolts located on the Block K Block A LFBG sensor

Sensor Placement and Loading Mode
The LFBG sensors, denoted S1-S9, are fixed on the inner surface of the area influenced by the joint in E6 (one-third span) as shown in Figure 8. This figure also gives the height coordinates of S1-S9 by setting the horizontal axis at the center of the cross-section of the ring. The bolts located on the outside of the sensors are named B1-B9. The sensor placement for E9 (mid-span) is the same. These sensors are simulated by the Beam188 element. All sensors have the uniform gauge length of 300 mm and the uniform width of 10 mm. The elastic modulus of the sensors is set to 8.0 MPa. The two ends of the sensors are glued to the inner surface of the segment, and the sensing region is separated from the segment.

Sensor Placement and Loading Mode
The LFBG sensors, denoted S1-S9, are fixed on the inner surface of the area influenced by the joint in E6 (one-third span) as shown in Figure 8. This figure also gives the height coordinates of S1-S9 by setting the horizontal axis at the center of the cross-section of the ring. The bolts located on the outside of the sensors are named B1-B9. The sensor placement for E9 (mid-span) is the same. These sensors are simulated by the Beam188 element. All sensors have the uniform gauge length of 300 mm and the uniform width of 10 mm. The elastic modulus of the sensors is set to 8.0 MPa. The two ends of the sensors are glued to the inner surface of the segment, and the sensing region is separated from the segment. Displacement constraints are applied in the z-direction at both ends of the model and in the x-direction on the outer surface of the ring line, as shown in Figure 9. The uniform load p is applied by seven successive loading steps from 0 to 30 kPa, 60 kPa, 90 kPa, 105 kPa, 120 kPa, 150 kPa, and 180 kPa.  Displacement constraints are applied in the z-direction at both ends of the model and in the x-direction on the outer surface of the ring line, as shown in Figure 9. The uniform load p is applied by seven successive loading steps from 0 to 30 kPa, 60 kPa, 90 kPa, 105 kPa, 120 kPa, 150 kPa, and 180 kPa.  Tables 1 and 2 show the tensile and compressive strain measurements for E9 and E6 at each loading step. The tensile strain is defined to be positive. Figure 10 shows that there are approximately linear relationships between the measured strains and the heights of the sensors. Thus, the method of least squares easily obtains the linear equations of y = kε+c for each linear relationship. These equations are used to construct Table 3 that shows the coefficient of determination R 2 and the parameter c that is equal to the monitored NAD χ for each equation. The NAD of E9 remains approximately 2.560 m when the applied load p is less than 105 kPa. The value of χ begins to increase when p reaches 105 kPa. For E6, the same trend of χ also occurs for p greater than 120 kPa. Thus, the proposed detecting method in Section 3 defines p = 105 kPa and p = 120 kPa as the threshold values of structural damage in E9 and E6, respectively. The stresses on bolts at different loading steps for E9 and E6 are listed in Tables 4 and 5, respectively. Table 4 shows that the stress of the outermost bolt B9 in the tensile region of E9 reaches 661 MPa for a p of 105 kPa. This stress value of this bolt slightly exceeds the yield stress of 640 MPa. Table 5 shows that the stress of B9 in the tensile region of E6 reached 651 MPa at a p of 120 kPa. This bolt stress is also larger than the 640 MPa yield stress. These stress values indicate that the NAD χ is a sensitive index that can indicate the occurrence of damage. Moreover, the applicability of the plane-section assumption is verified since  Tables 1 and 2 show the tensile and compressive strain measurements for E9 and E6 at each loading step. The tensile strain is defined to be positive. Figure 10 shows that there are approximately linear relationships between the measured strains and the heights of the sensors. Thus, the method of least squares easily obtains the linear equations of y = kε + c for each linear relationship. These equations are used to construct Table 3 that shows the coefficient of determination R 2 and the parameter c that is equal to the monitored NAD χ for each equation. The NAD of E9 remains approximately 2.560 m when the applied load p is less than 105 kPa. The value of χ begins to increase when p reaches 105 kPa. For E6, the same trend of χ also occurs for p greater than 120 kPa. Thus, the proposed detecting method in Section 3 defines p = 105 kPa and p = 120 kPa as the threshold values of structural damage in E9 and E6, respectively. The stresses on bolts at different loading steps for E9 and E6 are listed in Tables 4  and 5, respectively. Table 4 shows that the stress of the outermost bolt B9 in the tensile region of E9 reaches 661 MPa for a p of 105 kPa. This stress value of this bolt slightly exceeds the yield stress of  Table 5 shows that the stress of B9 in the tensile region of E6 reached 651 MPa at a p of 120 kPa. This bolt stress is also larger than the 640 MPa yield stress. These stress values indicate that the NAD χ is a sensitive index that can indicate the occurrence of damage. Moreover, the applicability of the plane-section assumption is verified since the coefficient of determination R 2 of each best-fitting straight line is greater than 0.997 for each loading step before and after damage.   S2  3  7  13  20  57  70  270  S3  30  63  130  183  347  800  1400  S4  73  147  310  427  933  1743  3210  S5  120  243  517  713  1523  2803  5277  S6  170  340  723  1000  2133  4067  7307  S7  213  423  903  1227  2590  5013  9033  S8  240  480  1017  1413  3010  5733  10,367  S9  250  500  1063  1467  3170

Accuracy Comparison of ϕ ud
The theoretical ϕ ud is obtained from the proposed models of Xu [51] and Shiba [46]. In the latter work, ϕ ud is determined from: where the parameter l in Equation (19) represents the whole width of the ring. Shiba [46] proposed that the length of the area influenced by the joint should be equal to l. However, Xu [51] believed that the length may be less than l. Therefore, the parameter l in Equation (19) is replaced by l b as shown previously in Equation (12). The two models produce different ϕ ud values, so it is necessary to ascertain which model has better accuracy. The mechanical and geometrical properties of the real tunnel are used for the comparison of the theoretical ϕ ud as calculated by Equation (19), Equation (12), and Equation (8). The results of the comparison are given in Table 6. The value of ϕ ud obtained from Xu [51] is closer to the result from the numerical model than that from Equation (19). Thus, Xu's model is more accurate than Shiba's model. According to Xu [51], λ is 0.6225 for the case of straight-jointed segmental tunnel lining. Substituting this value of λ into Equation (10) results in a η ud of 0.126.

Verification of the Damage Index's Accuracy
The damage index D is determined by H FM1 and H FM3 according to Equation (7). The verification for the accuracy of this D can be divided into two parts: verifying the accuracy of H FM1 and verifying the accuracy of H FM3 . The following equation is used to obtain the practical flexural stiffness at each loading step before or after damage: where M is the average moment of the monitored element, and θ is the average slope of the monitored element. The calculation for this case shows that the theoretical flexural stiffness of the cross-section of the segment is 7.544 × 10 8 kN·m 2 . Tables 7 and 8 give the comparisons between the (EI) d /(EI) ud and H FM1 for E9 and E6, respectively. The H FM1 value is very close to the value of (EI) d /(EI) ud before and after damage occurrence. The maximum error between the (EI) d /(EI) ud and the H FM1 in E9 and E6 is no more than 1.7% and 0.8%, respectively. The calculated γ is used for determining n d at each loading step. Table 9 gives the comparison of this determined n d in E9 to that established from the fitting equations based on the long-gauge strain measurements and the stress on the bolt shown in Table 4. Table 10 shows the same comparison for E6. It is found that the n d determined from the fitting equations is equal to the true n d in the entire damage process.

Experimental Verification
A scaled-down longitudinal concrete model was employed to simulate an actual shield tunnel. The damage levels are monitored and compared to the true damage levels. This comparison is important to the applicability of the proposed damage detection and evaluation methods for practical monitoring.

Scaled-Down Model
The Qingchun Road of Hangzhou shield tunnel was selected as the prototype for the scaled-down model's design. The similarity ratio of the geometric dimensions and the elastic modulus in the model is 7 and 1, respectively. The M40 bolts in the real tunnel were simulated by bent screws with a uniform diameter of 6 mm. The normal concrete in the real tunnel was modeled by a fine aggregate concrete. The mechanical properties of the segment and bolt were obtained by laboratory tests and are as follows: E s = 35 GPa, E b = 190 GPa, µ s = 0.2, µ b = 0.3, the compressive strength of the concrete = 41.4 MPa, the yield stress of a bolt = 710.08 MPa, the yield strain of a bolt = 3680 µε, and k j2 /k j1 = 0.001. The scaled-down model has 20 segment rings. The outside diameter is 860 mm and the inside diameter is 780 mm. The width and the thickness of a segment ring is 200 mm and 40 mm, respectively. The segment ring is composed of six standard Block As (60 • ) with staggered joints. The staggered angle between the two adjacent rings is 20 • . A joint of the lining ring contains 12 circumferential weld bolts and 18 longitudinal joint bolts. The value of l b is 0.1 m. Figure 11 shows the components of the model. 3680 με, and kj2/kj1 = 0.001. The scaled-down model has 20 segment rings. The outside diameter is 860 mm and the inside diameter is 780 mm. The width and the thickness of a segment ring is 200 mm and 40 mm, respectively. The segment ring is composed of six standard Block As (60°) with staggered joints. The staggered angle between the two adjacent rings is 20°. A joint of the lining ring contains 12 circumferential weld bolts and 18 longitudinal joint bolts. The value of lb is 0.1 m. Figure  11 shows the components of the model.

Sensor Placement and Loading Mode
The two half-segment rings at both ends were ignored, and the remaining model was uniformly divided into 19 elements denoted E1-E19. Figure 12a shows the placement locations of 10 LFBG sensors with a uniform length of 70 mm that were fixed on the inner surface of the area influenced by the joint in E8. These sensors are denoted S1-S10. The maximum strain that can be measured by these strain sensors is approximately 3600 με. The bolts on the outside of the sensors are labeled B1-B10. As shown in Figure 12b, each bolt of B1-B10 is attached to an electrical resistance strain gauge with a uniform size of 4 mm × 3 mm to measure the tensile or compressive strain in the bolt at each loading step. These sensors have a maximum strain measurement of approximately 4500 με. The sensor placement for E10 is the same as E8.

Sensor Placement and Loading Mode
The two half-segment rings at both ends were ignored, and the remaining model was uniformly divided into 19 elements denoted E1-E19. Figure 12a shows the placement locations of 10 LFBG sensors with a uniform length of 70 mm that were fixed on the inner surface of the area influenced by the joint in E8. These sensors are denoted S1-S10. The maximum strain that can be measured by these strain sensors is approximately 3600 µε. The bolts on the outside of the sensors are labeled B1-B10. As shown in Figure 12b, each bolt of B1-B10 is attached to an electrical resistance strain gauge with a uniform size of 4 mm × 3 mm to measure the tensile or compressive strain in the bolt at each loading step. These sensors have a maximum strain measurement of approximately 4500 µε. The sensor placement for E10 is the same as E8.
Two concrete supports with a curved groove were placed at both ends of the model as shown in Figure 13. Two curved steel sheets with a length of 650 mm were laid between the model and the groove to smooth the stress concentration and prevent shear failure. The concentrated load produced by an electro-hydraulic servo loading system was divided equally into two parts by a transfer steel board. Each part of the load passes through a steel plate designed to have the same curved surfaces as the steel sheets. Rubber pads were placed between the plates and the model to ensure a smooth load transfer process. The applied load P was increased by using displacement control from a linear variable differential transformer installed in the loading system. The applied displacement y of the loading system was increased by 1 mm for each step from y = 0 to y = 12 mm. The two midpoints of the uniform loads are located at the boundary point between E8 and E9 and the boundary point between E11 and E12. The two supporting points are located at the boundary point between E1 and E2 and the boundary point between E18 and E19. The increase in speed of the applied displacement is 0.5 mm per minute to prevent an accidental rupture that may occur at any circumferential joint. A schematic and the setup of the loading system are shown in Figure 13. The applied load P measured at each loading step is listed in Table 11. Two concrete supports with a curved groove were placed at both ends of the model as shown in Figure 13. Two curved steel sheets with a length of 650 mm were laid between the model and the groove to smooth the stress concentration and prevent shear failure. The concentrated load produced by an electro-hydraulic servo loading system was divided equally into two parts by a transfer steel board. Each part of the load passes through a steel plate designed to have the same curved surfaces as the steel sheets. Rubber pads were placed between the plates and the model to ensure a smooth load transfer process. The applied load P was increased by using displacement control from a linear variable differential transformer installed in the loading system. The applied displacement y of the loading system was increased by 1 mm for each step from y = 0 to y = 12 mm. The two midpoints of the uniform loads are located at the boundary point between E8 and E9 and the boundary point between E11 and E12. The two supporting points are located at the boundary point between E1 and E2 and the boundary point between E18 and E19. The increase in speed of the applied displacement is 0.5 mm per minute to prevent an accidental rupture that may occur at any circumferential joint. A schematic and the setup of the loading system are shown in Figure 13. The applied load P measured at each loading step is listed in Table 11.   Two concrete supports with a curved groove were placed at both ends of the model as shown in Figure 13. Two curved steel sheets with a length of 650 mm were laid between the model and the groove to smooth the stress concentration and prevent shear failure. The concentrated load produced by an electro-hydraulic servo loading system was divided equally into two parts by a transfer steel board. Each part of the load passes through a steel plate designed to have the same curved surfaces as the steel sheets. Rubber pads were placed between the plates and the model to ensure a smooth load transfer process. The applied load P was increased by using displacement control from a linear variable differential transformer installed in the loading system. The applied displacement y of the loading system was increased by 1 mm for each step from y = 0 to y = 12 mm. The two midpoints of the uniform loads are located at the boundary point between E8 and E9 and the boundary point between E11 and E12. The two supporting points are located at the boundary point between E1 and E2 and the boundary point between E18 and E19. The increase in speed of the applied displacement is 0.5 mm per minute to prevent an accidental rupture that may occur at any circumferential joint. A schematic and the setup of the loading system are shown in Figure 13. The applied load P measured at each loading step is listed in Table 11.

Verification of the Damage Index's Accuracy
The applied load P and the positions of loading points and supporting points are used to obtain the bending moment M acting on the test specimen. The curvatures κ of the testing specimens at each loading step are obtained by the fitting equations to the monitored long-gauge strains. The practical flexural stiffness of the monitored element before or after damage at each loading step is calculated according to: According to Xu [51], λ is 0.535 for the case of a stagger-jointed segmental tunnel lining. Substituting a λ of 0.535 into Equation (10), the value of η ud is determined to be 0.144. The comparison between the (EI) d /(EI) ud and the H FM1 for E10 and E8 from the scaled-down model is shown in Tables 17 and 18, respectively. It is significant that H FM1 is equal to (EI) d /(EI) ud in the non-damaged state, and the deviations between the two ratios are no more than 5.0% for the early-damaged state (H FM1 ≥ 0.5) in both E10 or E8. The maximum error would not exceed 9.0% during the full damage process.   Table 19 gives the comparison of n d for E10 from the fitting equations based on monitored long-gauge strains and from those listed in Table 14. Table 20 presents the same comparison of n d for E8. The results show that the n d determined from the fitting equations is equal to the n d as monitored from strains at each loading step. Thus, the proposed H FM3 in this paper is an accurate damage index in both the undamaged state and the damaged state.

Conclusions
This study proposes a new damaged detection and evaluation method based on NAD monitoring in the longitudinal direction for in-service shield tunnels. Verifications of the proposed method were performed using numerical simulations and experiments, and the following conclusions can be drawn:

1.
Based on the analysis of two deformation modes (i.e., bending mode and dislocation mode), three kinds of failure modes are assumed, including the bending failure mode, the shearing failure mode due to dislocation, and the shearing failure mode of the equivalent shear stiffness reduction, which resulted from the bolts being subjected to an excessive bending moment. These three failure modes and their development were successfully detected and evaluated by the proposed damage index D.

2.
The damage index D can be determined only by the monitored NAD, which is obtained from the long-gauge strain distribution of the LFBG sensor array. The calculation of the damage index value is independent of the loading mode and the mechanical properties of materials. 3.
The damage index D characterizes the development of damage in the tunnel accurately through calculating the coefficients H FM1 and H FM3 in most cases. The maximum deviation between the H FM1 and the true decrease of the equivalent flexural stiffness is less than 1.7% in the numerical simulation result and less than 9% in the experimental result. The calculated H FM3 matches well with its true value both in the undamaged state and the early-damaged state in the numerical simulation and the experimental results. The accuracy of the index D demonstrates the great applicability of the proposed damaged detection and evaluation method in the practical monitoring and maintenance of in-service shield tunnels.