Analysis of Tunnel Reliability Based on Limit Strain Theory

The construction of tunnel projects entails significant costs and long service periods, making the study of their stability a subject of widespread interest among scholars both domestically and internationally. A multitude of uncertainties exist in the construction process of tunnel engineering, for which reliability methods offer an effective means of quantitative analysis. Currently, the principal methods for reliability analysis of tunnel engineering include the response surface method, stochastic finite element method, and Monte Carlo method. For instance, Park and Park [1] proposed a back-analysis method based on the continuous response surface method for identifying geomechanical parameters of rock masses near underground chambers, which boasts high accuracy and computational efficiency in parameter identification. Li and Yang [2] applied the response surface method to study the failure mechanisms of subway tunnels, identifying the impact patterns of both non-random and random variables on the failure probability. Li and Yang [3] introduced an improved response surface method for the reliability analysis of the working face in fractured soft rock tunnels, demonstrating that the stability of the tunnel face is more sensitive to the internal friction angle than to cohesion. Wang and Li [4] developed a stochastic response surface method that considers non-Gaussian random variables under incomplete probability information in tunnel excavation reliability analysis, extending the application scope of the stochastic response surface method. Lyu et al. [5] evaluated the failure probability of each failure mode in tunnel systems using the Monte Carlo method and response surface method, finding that reliability primarily depends on the timing of support installation. Liu et al. [6] investigated the application of combining multiple response surface methods with the Monte Carlo method in analyzing the reliability of tunnel structures across different failure modes. Li et al. [7] explored the application of the non-intrusive stochastic finite element method in the reliability analysis of deformations in underground chambers. Liu and Low [8] employed the neural network response surface method for the reliability analysis of circular tunnels, which proved to be more efficient and accurate compared to the second-order response surface method.

Currently, scholars commonly employ deformation response values of tunnels, such as crown settlement, arch shoulder settlement, and surface settlement, as indicators for tunnel reliability analysis studies. It is determined that tunnel structures become unstable when deformation values exceed permissible limits [6], [7], [8], [9]. During the operational phase of tunnels, stress concentration occurs around the tunnel [10]. When the surrounding rock stress gradually surpasses the tunnel's load-bearing capacity, the lining may undergo partial or total failure, weakening the tunnel's load-bearing capacity and consequently leading to instability. Therefore, tunnel failure generally lags behind the failure of the lining, and using surface or crown settlement values exceeding permissible limits as the criterion for instability does not effectively represent the tunnel instability phenomena induced by lining failure. Based on the limit strain theory, this article adopts shear strain penetration at the tunnel haunches as the criterion for tunnel structural instability. A new method for tunnel reliability study, integrating neural network response surface methodology and Monte Carlo simulation, is proposed and its feasibility is validated through an engineering case study.

Figure 1 illustrates the stress-strain relationship curve for ideal plastic materials. As stress increases, the material initially enters the elastic stage within its load-bearing limit, where strain grows linearly, and the majority of energy absorbed by the material is transformed into elastic energy stored internally [11], [12]. With further increase in stress, the material exits the elastic stage and enters the initial phase of the plastic stage, at which point the material's stress reaches its peak strength. This state retains certain elastic characteristics, known as the elastic limit state, and the strain at this point, $\gamma_y$, is referred to as the elastic limit strain. As plasticity gradually progresses and the material's strain accumulates to $\gamma_y$, the accumulated elastic energy is instantaneously released, leading to material failure. This strain, $\gamma_y$, is termed the plastic limit strain [13].

For a given material, when the strain value of any of its elements reaches the plastic limit strain, it signifies that the element has met the conditions for failure, resulting in local damage and the formation of cracks. However, due to the restraining effect of adjacent elements, the failure of a single element does not lead to macroscopic damage of the material. It is only when the strain values of all elements within a certain area reach the limit strain, thus forming a zone of plastic limit strain penetration, that the material is truly considered to have failed [13], [14]. The limit strain values for concrete materials of various strength grades are detailed in Table 1.

The rationale behind the response surface method involves constructing an easily manageable response surface to substitute the actual complex true function, thereby enhancing computational efficiency. Presently, the response surface method primarily comprises polynomial response surface methodology and neural network response surface methodology. The polynomial response surface methodology calculates the true function values for a certain number of regularly selected input variable samples, then fits a mathematical polynomial, using the polynomial function to progressively approximate and substitute the true function as the response surface [15]. Neural network response surface methodology calculates the true function values for a large number of input variable samples, utilizing a back propagation (BP) neural network to construct the relationship between variables and the true function values, with the neural network serving as the response surface in place of the true function. When the true function is complex and involves numerous input variables, the neural network response surface offers higher fitting precision and more accurately simulates the true function compared to the polynomial response surface.

The Monte Carlo method obtains a large number of samples, N, through random sampling and calculates the true function for each sample, tallying the number of unstable samples, M. The ratio M/N represents the probability of instability, P. Compared with traditional Monte Carlo reliability methods, substituting the true function calculation with a neural network response surface allows for the computation of a large volume of samples in a short time, significantly improving Monte Carlo operation efficiency. The research approach involves: numerically simulating a large number of samples to obtain the actual values of the characteristic parameters of the function, then constructing hidden layer neurons, using the computed samples to train the established neural network, obtaining the neural network mapping relationship between input and output parameters. This neural network response surface substitutes the real function's response surface, and finally, reliability analysis is conducted using the Monte Carlo method [15], with the specific process as follows:

(a) The number of random variables for tunnel surrounding rock and their statistical characteristics are determined, obtaining the parameter samples required for constructing the response surface.

(b) A numerical model of the tunnel is established. For each sample parameter, calculations are conducted individually within the numerical model to obtain the corresponding true response values of characteristic parameters. The training samples for constructing the response surface consist of the random variable values of each sample and the true response values of the characteristic parameters.

(c) Based on the training samples, a neural network response surface between random variables and characteristic parameters is constructed.

(d) According to the statistical characteristics of random variables, N random samplings of the random variable values are performed.

(e) Each set of sampled values is substituted into the constructed neural network response surface to calculate the response values of characteristic parameters.

(f) Based on the permissible values of characteristic parameters, the number of samples among the N random samplings whose response values exceed permissible values, i.e., the number of unstable samples M, is tallied, with the instability probability being P=M/N.

The Shilong Tunnel in Chongqing City traverses a ridge mountain range, primarily through Grade IV surrounding rock. The tunnel has a width of 12 m and a height of 10 m, with a burial depth of 150 m [16]. FLAC3D software is employed for the reliability analysis of the tunnel structure, simplifying the tunnel excavation issue into a two-dimensional plane strain model, as illustrated in Figure 2. The model dimensions are 140 m in length, 80 m in height, and 2 m in width, with a set burial depth of 40 m. The weight of the overlying 110 m rock mass is applied as an equivalent load on the model's top. C20 concrete is used for the lining, with a thickness of 0.3 m. Specific parameters are provided in Table 2.

Following the limit strain theory introduced in Section 2.1, a criterion for limit strain failure was proposed, stating, "The failure of certain units in geotechnical engineering indicates local failure of the structure; when limit strain units are interconnected, total failure is achieved." This criterion was applied to the failure of concrete cubic specimens, validating its feasibility for practical application [13].

This study extends the application of this limit strain failure criterion to the failure of tunnel-lining structures. It is considered that a tunnel becomes unstable when the overall shear strain value across the thickness of a certain area of the lining exceeds the limit strain value. The lining in the arch waist area of the tunnel is most prone to damage. Hence, the criterion of limit strain penetration in the tunnel lining's arch waist area is adopted as the indicator of tunnel instability.

The lining material used for the tunnel is C20 concrete, for which the corresponding limit strain value is identified as 1.06‰ in Table 1. According to the established numerical analysis model, strain monitoring points are uniformly placed at the lining's arch waist position. If the shear strain values at any point across the thickness exceed the limit strain of 1.06‰, it is determined that the lining has failed, leading to the tunnel becoming unstable.

Figure 3 presents a comparison between the shear strain of the tunnel lining and the vertical (Z-direction) displacement for a subset of the calculation samples. In the shear strain contour map, the red areas indicate parts where the shear strain value exceeds the ultimate shear strain of 1.06‰ for C20 concrete, signifying the shear strain failure zone; the blue areas represent sections where the shear strain does not exceed 1.06‰, indicating the non-failure zones. The vertical displacement contour map reveals that the crown settlement of the four samples gradually increases from 0.67 cm to 1.52 cm. Concurrently, the shear strain failure zone begins to develop downwards from the arch waist, eventually extending through the entire area between the arch waist and the arch foot, indicating total shear strain penetration failure. Table 3 lists the corresponding lining failure areas and crown settlement values for each sample.

According to the standard, a control value of 1.5 cm for crown settlement is stipulated for tunnel support structures [17]. In Sample 1, where the tunnel lining did not undergo shear strain penetration failure, the crown settlement was 0.67 cm, not reaching the deformation control value of 1.5 cm. For Sample 2, shear strain penetration failure occurred in the arch waist area at a height of 2-4 m, indicating the onset of tunnel instability, with a crown settlement of 1.05 cm, which did not reach the 1.5 cm deformation control value. In Sample 3, the range of lining damage expanded further, with shear strain penetration failure occurring in the arch waist area at a height of 1.5-5.1 m, indicating a more severe state of tunnel instability, where the crown settlement was 1.24 cm, still below the 1.5 cm deformation control value. For Sample 4, extensive shear strain penetration failure occurred in the arch waist area at a height of 0.9-5.8 m, resulting in a significant loss of the original load-bearing capacity of the lining and severe tunnel instability, with the crown settlement reaching 1.52 cm, exceeding the 1.5 cm tunnel deformation control value. Based on the "1.5 cm deformation control value" stipulated by the standard, only Sample 4 of the tunnel was determined to be unstable, while the other three tunnel samples were considered stable [17]. However, according to the limit strain criterion, the tunnel linings of Samples 2, 3, and 4 all underwent shear strain penetration failure, indicating that Tunnels 2, 3, and 4 were all unstable. Clearly, using the crown settlement control value defined by the standard as the criterion is potentially hazardous, as it incorrectly judges the unstable Tunnel Samples 2 and 3 as stable, a misjudgment that could lead to engineering settlement disasters [17]. Therefore, it is considered that the limit strain penetration instability criterion is superior to the instability criterion based on the crown settlement control value stipulated by the standard [17]. The limit strain penetration instability criterion not only can predict the safety state of the tunnel in advance but also identifies potential tunnel damage areas, thereby providing a reference for the support design of tunnel engineering.

In this study, the modulus of elasticity ($E$), cohesion ($c$), and the angle of internal friction ($\varphi$) of the rock mass were selected as the subjects of study. They were treated as random variables for constructing the response surface and calculating the probability of instability. Specific parameters are listed in Table 4.

Following the response surface construction method described in Section 1.2, a neural network response surface for the shear strain values at the monitoring points of the tunnel lining's arch waist was developed, with the construction carried out in Matlab. To validate the fitting accuracy of the response surface, 30 verification points were randomly selected, and a scatter plot was drawn with FLAC3D calculation results as the $x$-axis and response surface calculation results as the $y$-axis, as illustrated in Figure 4. The coefficient of determination, $R^2$, which measures the response surface's fitting accuracy, approaches 1 as the accuracy increases. For the shear strain response surface of the arch waist lining, the fitting coefficient $R^2$ is 0.99, indicating high accuracy and allowing it to substitute the calculations by the FLAC3D numerical software.

Upon completion of the response surface construction, the Monte Carlo method was employed to analyze the tunnel's reliability, calculating the probability of instability for different numbers of samplings. The number of samplings was plotted on the x-axis, and the calculated probability of instability on the $y$-axis, with the resulting curve depicted in Figure 5.

When the number of samplings was less than 20,000, the calculated probability of tunnel instability showed significant fluctuations. As the number of samplings increased from 20,000 to 50,000, the fluctuations decreased and gradually stabilized. When the number of samplings exceeded 50,000, the calculated values essentially remained unchanged, indicating convergence. To ensure accuracy of the results, the number of samplings was set to 100,000, ultimately determining the probability of tunnel instability to be 10.90%.

In actual tunnel engineering, the spatial variability of rock and soil parameters, such as modulus of elasticity, cohesion, and angle of internal friction, is relatively significant [18]. To study the sensitivity of the shear strain values at the tunnel arch waist to these three parameters, 15 working conditions as listed in Table 5 were designed. The instability probability for each was calculated, and the distribution of shear strain values was analyzed.

Figure 6 illustrates the trend of instability probability as a function of the coefficients of variation for parameters $E$, $c$, and $\varphi$. With the increase in the coefficient of variation for these three parameters, the calculated probability of instability shows a rising trend.

Among them, when the coefficient of variation ranges from 0.1 to 0.3, the change in the probability of instability for the modulus of elasticity $E$ is the most significant, ranging from 2.31% to 20.08%. The change in the probability of instability for the angle of internal friction $\varphi$ follows, ranging from 8.81% to 15.28%. The change in the probability of instability for cohesion $c$ is the least, ranging from 10.85% to 11.49%. Based on these calculation results, it is concluded that the variation in the modulus of elasticity $E$ has the greatest impact on the overall tunnel instability probability, followed by the angle of internal friction $\varphi$, with cohesion c having the least impact.

Figure 7, Figure 8, Figure 9 depict histograms of the frequency distribution of shear strain values in the arch waist lining, based on 100,000 Monte Carlo simulations under different working conditions.

Figure 7 shows the changes in the distribution of shear strain values in the arch waist lining of the tunnel when the coefficient of variation of the rock mass's modulus of elasticity E changes. Under working condition E_1 (COV=0.1), the distribution range of shear strain values is from 0.4‰ to 1.5‰, with median and 95th percentile values of 0.692‰ and 0.962‰, respectively. Under working condition E_3 (COV=0.2), the distribution range expands to 0.4‰ to 2.3‰, with median and 95th percentile values of 0.698‰ and 1.301‰, respectively. For working condition E_5 (COV=0.3), the range further extends to 0.4‰ to 2.8‰, with median and 95th percentile values of 0.700‰ and 1.832‰, respectively. As the coefficient of variation of the modulus of elasticity E increases, the changes in the size of the shear strain values in the tunnel’s arch waist lining become more pronounced, with the distribution range shifting from the initial 0.4‰ to 1.5‰ to 0.4‰ to 2.8‰. This significant overall change indicates that the modulus of elasticity E has a very pronounced impact on the size of the shear strain values in the tunnel’s arch waist lining.

Figure 8 shows the changes in the distribution of shear strain values in the arch waist lining of the tunnel as the coefficient of variation for the rock mass's cohesion c changes. As the coefficient of variation for cohesion c increases, the distribution of shear strain values across the three conditions does not show significant changes, with the overall distribution range remaining between 0.4‰ and 2.3‰. The median values are respectively 0.691‰, 0.698‰, and 0.706‰, while the 95th percentile values are respectively 1.292‰, 1.313‰, and 1.341‰. This indicates that the size of shear strain values in the tunnel’s arch waist lining is minimally affected by changes in the rock mass's cohesion c. Variations in the rock mass's cohesion c do not significantly alter the size of the shear strain values in the arch waist lining.

Figure 9 shows the changes in the distribution of shear strain values in the arch waist lining of the tunnel as the coefficient of variation for the rock mass's internal friction angle $\varphi$ changes. Under condition $\varphi$_1 (COV=0.1), the distribution range of shear strain values is from 0.4‰ to 2.1‰, with median and 95th percentile values of 0.691‰ and 1.223‰, respectively. Under condition $\varphi$_3 (COV=0.2), the range is from 0.4‰ to 2.2‰, with median and 95th percentile values of 0.697‰ and 1.301‰, respectively. For condition $\varphi$_5 (COV=0.3), the range extends from 0.4‰ to 2.7‰, with median and 95th percentile values of 0.710‰ and 1.537‰, respectively. The distribution of shear strain values becomes significantly more pronounced as the coefficient of variation for the internal friction angle $\varphi$ changes, with the overall values noticeably increasing, indicating that the size of shear strain values in the tunnel arch waist lining is significantly affected by the internal friction angle $\varphi$.

The distribution changes of shear strain values reveal that when the coefficients of variation for the modulus of elasticity E, cohesion c, and internal friction angle $\varphi$ are altered, the impact of the modulus of elasticity E on the distribution of shear strain values is the most significant, followed by the internal friction angle $\varphi$, with cohesion c having the least impact. This demonstrates that the sensitivity of shear strain values to these three parameters decreases in the order of modulus of elasticity E, internal friction angle $\varphi$, and cohesion c.

As the distribution range of shear strain values widens and the values become more dispersed, the number of samples exceeding the permissible shear strain values (limit strain) for concrete lining increases, thereby leading to an increase in the probability of instability. The changes in the coefficient of variation for the modulus of elasticity E cause the most significant changes in the distribution of shear strain values, resulting in the greatest impact of the modulus of elasticity E on the probability of instability, as shown in Figure 6.

The following conclusions have been drawn from this study:

(a) Based on limit strain theory, with the shear strain penetration of the lining structure as the criterion for tunnel instability, a novel method for tunnel structure reliability analysis combining neural network response surface methodology and Monte Carlo simulation is proposed.

(b) Compared to using tunnel deformation values as the control indicator for tunnel instability, employing shear strain penetration failure of the tunnel lining as the criterion for instability more effectively determines the unstable state of tunnels, enhancing the accuracy of tunnel reliability analysis.

(c) The probability of tunnel structure instability increases with the rise in the coefficient of variation of rock mass parameters. Among these, the modulus of elasticity E has the most significant impact on the probability of tunnel structure instability.

(d) The criterion of shear strain penetration failure for tunnel lining instability more accurately represents phenomena of tunnel structure instability induced by events such as earthquakes and blasting. This offers a new perspective for the reliability analysis of tunnel structures.