A New Reliability Analysis Model of the Chegongzhuang Heat-Supplying Tunnel Structure Considering the Coupling of Pipeline Thrust and Thermal Effect

Based on the operating Chegongzhuang heat-supplying tunnel in Beijing, the reliability of its lining structure under the action of large thrust and thermal effect is studied. According to the characteristics of a heat-supplying tunnel service, a three-dimensional numerical analysis model was established based on the mechanical tests on the in-situ specimens. The stress and strain of the tunnel structure were obtained before and after the operation. Compared with the field monitoring data, the rationality of the model was verified. After extracting the internal force of the lining structure, the improved method of subset simulation was proposed as the performance function to calculate the reliability of the main control section of the tunnel. In contrast to the traditional calculation method, the analytic relationship between the sample numbers in the subset simulation method and Monte Carlo method was given. The results indicate that the lining structure is greatly influenced by coupling in the range of six meters from the fixed brackets, especially the tunnel floor. The improved subset simulation method can greatly save computation time and improve computational efficiency under the premise of ensuring the accuracy of calculation. It is suitable for the reliability calculation of tunnel engineering, because “the lower the probability, the more efficient the calculation.”


Introduction
Heat-supplying tunnels are built for district heating pipelines in cold regions of the world, which is what distinguishes it from traffic and other municipal tunnels in its function and environment. High-temperature fluid transmission throughout the pipe leads to the temperature rising in heat-supplying tunnels. In addition, the pressure pipeline also produces a huge longitudinal thrust on the tunnel structure through the fixed brackets in the tunnel. The coupling of large thrust by the operating pipeline and thermal effect is always accompanied by cold and hot cycles and dry and wet alternation, which leads to the decrease of structural reliability and residual life. However, due to the fact that there is no direct contact between the public and the heat-supplying tunnel, the reliability and durability of these tunnels are seldom addressed by the public and scholars. In Beijing, for example, there were a large number of heat-supplying tunnels built in the 1990s. Heat-supplying tunnels in these regions were subject to structural degradation over time, causing latent risks to the above-grade traffic safety. Therefore, it is necessary to study the reliability of the tunnel lining structure under the coupling of large thrust by operating pipeline and thermal effect.
As far as reliability is concerned, the failure probability of the tunnel structure is often below 10 −4 as a statically indeterminate system. A large number of sample points are needed in the document. The properties of the three kinds of overlaying soil are listed in Table 1. The cross section is shown in Figure 1. During operation, the fluid transfers heat to the pipe and its outer insulating layer through heat conduction. Because of the convection and radiation, the temperature of the inner wall of the tunnel rises and transfers heat to the outer wall and then to the surrounding rock. On the other hand, highpressure fluid acts on the pressure in the pipe. Due to the limitation of the fixed bracket near the tunnel portal, the energy will act on the cross arm and the column of the fixed bracket in the form of a huge thrust and then the fixed bracket will transmit the force to the lining structure of the tunnel. Taking the Chegongzhuang tunnel as an example, the clearance size of the standard section of the main tunnel is 5.0 m (width) × 3.0 m (height). A set of fixed brackets in the tunnel is at the 1.5 m position of the north side of Manhole 15. The two fixed brackets are staggered 1 m front and back and the lining is thickened with 20 cm in 3 m range before and after the brackets' zone as shown in Figure 2. The longitude thrust on the brackets is 2200 kN and 2150 kN respectively, pointing to the tunnel portal. The action diagram is shown in Figure 3. During operation, the fluid transfers heat to the pipe and its outer insulating layer through heat conduction. Because of the convection and radiation, the temperature of the inner wall of the tunnel rises and transfers heat to the outer wall and then to the surrounding rock. On the other hand, high-pressure fluid acts on the pressure in the pipe. Due to the limitation of the fixed bracket near the tunnel portal, the energy will act on the cross arm and the column of the fixed bracket in the form of a huge thrust and then the fixed bracket will transmit the force to the lining structure of the tunnel. Taking the Chegongzhuang tunnel as an example, the clearance size of the standard section of the main tunnel is 5.0 m (width) × 3.0 m (height). A set of fixed brackets in the tunnel is at the 1.5 m position of the north side of Manhole 15. The two fixed brackets are staggered 1 m front and back and the lining is thickened with 20 cm in 3 m range before and after the brackets' zone as shown in Figure 2. The longitude thrust on the brackets is 2200 kN and 2150 kN respectively, pointing to the tunnel portal. The action diagram is shown in Figure 3.

Performance Functions and Variable Parameters
Based on the Chinese current design specification of concrete structure GB50010-2010 [18], the performance functions were set up. The general section of the lining of the tunnel was treated as the eccentric compression member. The large and small eccentricity was judged by Equation (1) at first and then Equations (2) and (3) were used respectively. In particular, in view of the longitudinal thrust response of Section 1, the longitudinal local bearing capacity of concrete Equation (4) was added. ' ' y s y s 0 1 y cm When F1 < 0, it belongs to the large eccentricity and uses Equation (2): y y 0 y 0 s 0 a s cm (N-f As +f As ) h F =(N-f As +f As )h -+f As (h -a )-N(e +e + -a ) 2bf 2 (2) When F1 > 0, it belongs to the small eccentricity and uses Equation (3):

Performance Functions and Variable Parameters
Based on the Chinese current design specification of concrete structure GB50010-2010 [18], the performance functions were set up. The general section of the lining of the tunnel was treated as the eccentric compression member. The large and small eccentricity was judged by Equation (1) at first and then Equations (2) and (3) were used respectively. In particular, in view of the longitudinal thrust response of Section 1, the longitudinal local bearing capacity of concrete Equation (4) was added.
When F1 < 0, it belongs to the large eccentricity and uses Equation (2): y y 0 y 0 s 0 a s cm (N-f As +f As ) h F =(N-f As +f As )h -+f As (h -a )-N(e +e + -a ) 2bf 2 When F1 > 0, it belongs to the small eccentricity and uses Equation (3):

Performance Functions and Variable Parameters
Based on the Chinese current design specification of concrete structure GB50010-2010 [18], the performance functions were set up. The general section of the lining of the tunnel was treated as the eccentric compression member. The large and small eccentricity was judged by Equation (1) at first and then Equations (2) and (3) were used respectively. In particular, in view of the longitudinal thrust response of Section 1, the longitudinal local bearing capacity of concrete Equation (4) was added.
When F 1 < 0, it belongs to the large eccentricity and uses Equation (2): Particularly, the longitudinal local pressure performance function of the concrete in Section 1: where: M yy , N xx = the bending moment and axial force of lining structure during the operation period are derived from numerical simulation results; N yy = the longitudinal axial force of lining in Section 1, data from the said numerical simulation; b = the width of lining cross section. b = 1 m; c = the width of bracket cross section. c = 0.2 m; e 0 = the eccentricity of axial pressure on the center of gravity of the lining section, e 0 = M yy N xx ; and e a = the additional eccentricity, 20 mm.
The variable parameters showed in Table 2 were all from the mechanical properties test of the in-situ specimens. Figure 4 showed compressive and tensile tests of concrete and steel bar specimens respectively. By testing the physical and mechanical properties of the concrete and steel bars of the lining structure, the real denaturation of the structural reliability can be determined.

Analysis Model and Basic Assumption
(1) The establishment of the analysis mechanism. In order to obtain the mechanical response of a heat-supplying tunnel under the coupling of temperature and large thrust during the operation period, the indirect analysis method was used in this paper. Firstly, the thermal analysis should be carried out to calculate the distribution of the temperature field. Secondly, the results of the thermal analysis are then applied to the model as the temperature load and the large thrust of the pipe. Finally, the internal force distribution of the lining structure in the operation period is obtained to provide internal force parameters for the performance function of reliability analysis.
Three-dimensional heat conduction partial differential equation [19]: ∂U k ∂U ∂U ∂U = ( + + ) ∂t ρc ∂ x ∂ y ∂ z (5) where: U = temperature at any point; t = time; ρ = density; c = specific heat capacity; and k = thermal conductivity. Xu Z.S. [17] gives the relationship between k, c and U: ( ) 0.6 k U =1.6-U 850 (6) 420U c(U) = 840 + 850 (7) Furthermore, this paper assumes that: a. As far as the constitutive model was concerned, the concrete lining was a uniform linear elastic body and the soil was an elastoplastic body. The concrete lining and fixed bracket obeyed Hooke's law and Table 3 offers the parameters of structure.

Analysis Model and Basic Assumption
(1) The establishment of the analysis mechanism. In order to obtain the mechanical response of a heat-supplying tunnel under the coupling of temperature and large thrust during the operation period, the indirect analysis method was used in this paper. Firstly, the thermal analysis should be carried out to calculate the distribution of the temperature field. Secondly, the results of the thermal analysis are then applied to the model as the temperature load and the large thrust of the pipe. Finally, the internal force distribution of the lining structure in the operation period is obtained to provide internal force parameters for the performance function of reliability analysis.
Three-dimensional heat conduction partial differential equation [19]: where: U = temperature at any point; t = time; ρ = density; c = specific heat capacity; and k = thermal conductivity.
Xu Z.S. [17] gives the relationship between k, c and U: Furthermore, this paper assumes that: a.
As far as the constitutive model was concerned, the concrete lining was a uniform linear elastic body and the soil was an elastoplastic body. The concrete lining and fixed bracket obeyed Hooke's law and Table 3 offers the parameters of structure. The soil mass obeyed the Mohr-Coulomb criterion. Mohr-Coulomb criterion can be expressed as [20,21]: where: τ f = the shear strength of soil; σ = the vertical compressive stress on soil; c = the cohesion of soil and ϕ = internal friction angle of soil, as shown in Table 1. b.
The gradual heating process in the early stage of heating was not considered. It was assumed that the inner wall of the lining was constant and the same. There was no heat loss between the outer wall of a tunnel and the soil. The temperature transition is uniform on the contact surface of the lining and the soil mass. c.
The saturation, porosity, soil moisture, thermal conductivity and the volume-specific heat of the soil never change with the change of the temperature in the calculation. d. A finite soil model takes the place of an infinite soil model.   The soil mass obeyed the Mohr-Coulomb criterion. Mohr-Coulomb criterion can be expressed as [20,21]: f τ = c+ σtanφ (8) where: τ = the shear strength of soil; σ = the vertical compressive stress on soil; c = the cohesion of soil and φ = internal friction angle of soil, as shown in Table 1.
b. The gradual heating process in the early stage of heating was not considered. It was assumed that the inner wall of the lining was constant and the same. There was no heat loss between the outer wall of a tunnel and the soil. The temperature transition is uniform on the contact surface of the lining and the soil mass. c. The saturation, porosity, soil moisture, thermal conductivity and the volume-specific heat of the soil never change with the change of the temperature in the calculation. d. A finite soil model takes the place of an infinite soil model.  According to the field monitoring, Figure 6 showed the distribution of control sections like the central section of the fixed brackets, Y = 0 m, Y = 12 m and the expansion section Y = 23 m. They were named Sections 1, 2 and 3 by order. According to the field monitoring, Figure 6 showed the distribution of control sections like the central section of the fixed brackets, Y = 0 m, Y = 12 m and the expansion section Y = 23 m. They were named Sections 1-3 by order. (3) The selection of boundary conditions. 1. The mechanical boundary. Horizontal constraints were imposed on the vertical boundary. The edge sections of the tunnel were also limited by the horizontal constraints. Fixed constraints were imposed on the bottom boundary and the upper boundary was free. 2. The thermal boundary. It was measured that the temperature of the inner wall of the tunnel lining was 50 °C. According to the statistical research results in the literature [22], the temperature of land surface and initial boundary temperature of surrounding soil in the Beijing area were assumed 0 °C and 20 °C respectively. The distribution variation of the temperature value in the vertical direction and the influence of the heat-supplying operation of the branch tunnel were discarded.

Field Monitoring
Due to the need for excavation and dynamic monitoring, the team reserved survey points in the construction period of 2006. The longitudinal and lateral strain of the concrete were monitored by using concrete strain gauges, which were reserved for the lining within five sections of the Y = 0 m (Section 1), Y = 6 m, Y = 12 m (Section 2), Y = 18 m and Y = 23 m (Section 3). The layout of the monitoring points is shown in Figure 7.

Subset Simulation Process
If G is the failure domain of performance function F(x): G = {x|F(x) ≤ 0}. By introducing the threshold value q1 > q2 > q3 > … > qn = 0, a series of failure events Gk are formed. Gk needs to obey these rules: (1) Gk = {x|F(x) ≤ qk} (k = 1,2,3,…,n); (3) The selection of boundary conditions. 1. The mechanical boundary. Horizontal constraints were imposed on the vertical boundary. The edge sections of the tunnel were also limited by the horizontal constraints. Fixed constraints were imposed on the bottom boundary and the upper boundary was free. 2. The thermal boundary. It was measured that the temperature of the inner wall of the tunnel lining was 50 • C. According to the statistical research results in the literature [22], the temperature of land surface and initial boundary temperature of surrounding soil in the Beijing area were assumed 0 • C and 20 • C respectively. The distribution variation of the temperature value in the vertical direction and the influence of the heat-supplying operation of the branch tunnel were discarded.

Field Monitoring
Due to the need for excavation and dynamic monitoring, the team reserved survey points in the construction period of 2006. The longitudinal and lateral strain of the concrete were monitored by using concrete strain gauges, which were reserved for the lining within five sections of the Y = 0 m (Section 1), Y = 6 m, Y = 12 m (Section 2), Y = 18 m and Y = 23 m (Section 3). The layout of the monitoring points is shown in Figure 7. (3) The selection of boundary conditions. 1. The mechanical boundary. Horizontal constraints were imposed on the vertical boundary. The edge sections of the tunnel were also limited by the horizontal constraints. Fixed constraints were imposed on the bottom boundary and the upper boundary was free. 2. The thermal boundary. It was measured that the temperature of the inner wall of the tunnel lining was 50 °C. According to the statistical research results in the literature [22], the temperature of land surface and initial boundary temperature of surrounding soil in the Beijing area were assumed 0 °C and 20 °C respectively. The distribution variation of the temperature value in the vertical direction and the influence of the heat-supplying operation of the branch tunnel were discarded.

Field Monitoring
Due to the need for excavation and dynamic monitoring, the team reserved survey points in the construction period of 2006. The longitudinal and lateral strain of the concrete were monitored by using concrete strain gauges, which were reserved for the lining within five sections of the Y = 0 m (Section 1), Y = 6 m, Y = 12 m (Section 2), Y = 18 m and Y = 23 m (Section 3). The layout of the monitoring points is shown in Figure 7.

Subset Simulation Process
If G is the failure domain of performance function F(x): G = {x|F(x) ≤ 0}. By introducing the threshold value q 1 > q 2 > q 3 > . . . > q n = 0, a series of failure events G k are formed. G k needs to obey these rules: (1) G k = {x|F(x) ≤ q k } (k = 1,2,3, . . . ,n); . . ,n). The failure probability can be expressed as a form of Markov chain as follows: If P 1 = P(G 1 ), P i = P(G i |G i−1 ) , P f can also be expressed as Equation (10): From Equation (10) we can know that when the order of magnitude of P i is 10 −1 and n = 5, the order of magnitude of P f will be 10 −5 . As a result, the small probability problem can be transformed into the product of the larger conditional probability by using the subset simulation method, which can improve the efficiency of estimation of a small failure probability event.

Select the Intermediate Failure Event
From the theory of the subset simulation method, failure event can be expressed as a set of a series of intermediate failure events like G 1 ,G 2 ,G 3 , . . . ,G n . Therefore, it was very important to select appropriate intermediate failure events when using the subset simulation method to calculate the reliability problems. Au put forward a preset value of conditional probability P 0 and the method of automatic segmentation, in which the number of sampling points of each layer need to be the same [1,2]. Like N i = N 0 , both P 0 N 0 and (1-P 0 )N 0 must be integers.

Algorithm Implementation
Based on the principle of subset simulation method, the automatic hierarchical process of subset simulation was shown as follows: (1) The Carlo Monte method was used to generate N sample points {X (1) j : j = 1, 2, 3, . . . , N 0 }. These sample points were independent and obeyed the basic random variables of the probability density function f x (X).
(3) Starting from the number of P 0 N 0 sample points that fell within the G i−1 (i = 2, 3, 4, . . . , n), each sample point can be simulated by a Markov chain. Therefore, the number of 1/P 0 sample points were generated by each Markov chain. Therefore, the number of samples generated per layer was maintained at N 0 . The specific process was shown in Figure 8, where transfer acceptance rate of Markov chain Acc = min{1, F(X(i)) }. (4) Analog the second step and use the N sample points captured by the third step to calculate the value F(x (i) j ) of each function. In descending order, take the number (1-P 0 )N 0 value as the critical values q i of failure event G i . Here, q i = F(x (i) (1−p 0 )N 0 ) and P i = P 0 .

(5) Repeat
Step (3) and Step (4), until the critical value q n in Layer l was less than 0. Ultimately, the number of failure events was N f in N 0 sample points. Therefore, the conditional probability estimated: P n = N f N 0 . (6) Final failure probability estimated:P f = P 0 n−1 P n = P 0 n−1 · N f N 0 .  ) and Pi = P0.
(5) Repeat Step (3) and Step (4), until the critical value qn in Layer l was less than 0. Ultimately, the number of failure events was Nf in N0 sample points. Therefore, the conditional probability estimated:

Problems of the Existing Method
The N0 sample points are generated by the Monte Carlo method in the general subset simulation method. These sample points are pseudorandom numbers, which has a certain correlation. The distribution of the sample points is not uniform and the deviation can reach -0.5 ο(n ) [23]. For the calculation of small probability events represented by the structural failure probability, the convergence speed and accuracy are affected by PRN sampling points. On the one hand, the low convergence speed determines a huge time consuming in calculation. Especially in the reliability analysis of an underground structure engineering, the operation time cost may reach ten hours or

Problems of the Existing Method
The N 0 sample points are generated by the Monte Carlo method in the general subset simulation method. These sample points are pseudorandom numbers, which has a certain correlation. The distribution of the sample points is not uniform and the deviation can reach o(n −0.5 ) [23]. For the calculation of small probability events represented by the structural failure probability, the convergence speed and accuracy are affected by PRN sampling points. On the one hand, the low convergence speed determines a huge time consuming in calculation. Especially in the reliability analysis of an underground structure engineering, the operation time cost may reach ten hours or even days without the subset simulation process. On the other hand, the reliability of the calculation results can be doubtful because of the accuracy problems.
where p should be a prime number and p ≥ 2m + 3. The point set P n has a deviation o(n −0.5− 1 2(m−1) +ξ ). One hundred sample points generated randomly in two-dimensional unit space by using the Hua-Wang point set and the PRN method respectively were shown in Figure 9. It is found that the point distribution generated by the Hua-Wang point set method is more uniform. Therefore, the improvement of the subset simulation algorithm is embodied in Step (1) of Section 3.1.2 by using the Hua-Wang point set instead of Monte Carlo's PRN to generate more independent and uniformly distributed sample points.
items have a deviation ⋅ -1+ξ n D(n,P ) ≤ c(γ,ξ) n , then Pn can be called a good point set and γ is a good point.
where p should be a prime number and p ≥ 2m+3 . The point set Pn has a deviation 1 0.5 2(m 1) One hundred sample points generated randomly in two-dimensional unit space by using the Hua-Wang point set and the PRN method respectively were shown in Figure 9. It is found that the point distribution generated by the Hua-Wang point set method is more uniform. Therefore, the improvement of the subset simulation algorithm is embodied in Step (1) of Section 3.1.2 by using the Hua-Wang point set instead of Monte Carlo's PRN to generate more independent and uniformly distributed sample points.

Numerical Simulation
(1) Before the operation period, according to the construction steps, the excavation and support of the main tunnel, the expansion section and then the branch tunnel were simulated. After the installation of the fixed brackets, the state of the main tunnel before heating was obtained. Because of the thickening of the lining, the maximum principal stress σ of Section 1 is smaller in Table 4. The longitudinal strain of Section 1 and Section 2 is very small, which is in line with the plane strain characteristics of tunnel mechanics. The maximum tensile stress of the spandrel of Section 3 is up to 2.1 MPa, which is greatest affected by the excavation of the expansion section and branch tunnel among three sections.
(2) During the operation period, based on the indirect analysis method in said material, first consider the effect of temperature field separately. We can get the temperature field distribution by numerical calculation. The temperature field distribution of the main tunnel standard section (Section 2) is shown in Figure 10. In Figure 10, the crown and spandrel of the tunnel structure are dense areas of equal value, which indicates that the temperature gradient in crown and spandrel is large. In other words, the temperature difference between inner and outer lining at crown and spandrel is high. Because the shallow tunnel is greatly influenced by the land surface temperature, the temperature

Numerical Simulation
(1) Before the operation period, according to the construction steps, the excavation and support of the main tunnel, the expansion section and then the branch tunnel were simulated. After the installation of the fixed brackets, the state of the main tunnel before heating was obtained. Because of the thickening of the lining, the maximum principal stress σ 1 of Section 1 is smaller in Table 4. The longitudinal strain of Sections 1 and 2 is very small, which is in line with the plane strain characteristics of tunnel mechanics. The maximum tensile stress of the spandrel of Section 3 is up to 2.1 MPa, which is greatest affected by the excavation of the expansion section and branch tunnel among three sections.
(2) During the operation period, based on the indirect analysis method in said material, first consider the effect of temperature field separately. We can get the temperature field distribution by numerical calculation. The temperature field distribution of the main tunnel standard section (Section 2) is shown in Figure 10. In Figure 10, the crown and spandrel of the tunnel structure are dense areas of equal value, which indicates that the temperature gradient in crown and spandrel is large. In other words, the temperature difference between inner and outer lining at crown and spandrel is high. Because the shallow tunnel is greatly influenced by the land surface temperature, the temperature difference can reach 4.9 • C and spandrel is 4.6 • C. Against this backdrop, the temperature difference between the inner and outer walls is the cause of the temperature stress in the tunnel structure and the development of the relevant temperature crack should be prevented. difference can reach 4.9 °C and spandrel is 4.6 °C. Against this backdrop, the temperature difference between the inner and outer walls is the cause of the temperature stress in the tunnel structure and the development of the relevant temperature crack should be prevented. (3) Then, after applying thermal load and large thrust of pipeline, the mechanical response of lining structure of heat-supplying tunnel under the coupling condition of large thrust and temperature field during operation is obtained. The difference between before and after heating in Table 4 can be seen as the characterization of the coupling effect. Table 4 shows that because the huge pipe thrust is directly on the fixed bracket, it has a great influence on Section 1. The minimum principal stress of the floor in Section 1 is increased from −14.3 MPa to −17.8 MPa. The longitudinal effect is particularly prominent and the longitudinal compressive strain εyy reaches the −56.7 με. Furthermore, the stress factor applies the greatest influence on the floor, the crown and on the side wall in the respective order. Comparing the results before and after operation, the variation of longitudinal stress and strain of Section 2 is less and within 5.5 με at each place. The lateral strain changes in the range of 3 με. It shows that the coupling effect has been markedly attenuated at the distance 12 m from the bracket. For Section 3, there is little change in the longitudinal strain before and after operation. It can be seen that for the lining structure away from the fixed bracket, the influence of the coupling effect is mainly reflected by the single temperature field in the form of circumferential stress and strain. Results show that its influence is very limited and its maximum lateral strain can reach −5.6 με at the foot of sidewall. According to the results of numerical simulation, the internal force distribution of the lining structure can be obtained as Table 5 shows. It can be plugged into the performance functions in reliability calculation. The results indicate that due to the multi-field coupling effect in the operation period, the design of the heat-supplying tunnel at the present stage is blind only by enlarging the whole section thickness. The floor of the tunnel under the fixed bracket should be the key part. (3) Then, after applying thermal load and large thrust of pipeline, the mechanical response of lining structure of heat-supplying tunnel under the coupling condition of large thrust and temperature field during operation is obtained. The difference between before and after heating in Table 4 can be seen as the characterization of the coupling effect. Table 4 shows that because the huge pipe thrust is directly on the fixed bracket, it has a great influence on Section 1. The minimum principal stress of the floor in Section 1 is increased from −14.3 MPa to −17.8 MPa. The longitudinal effect is particularly prominent and the longitudinal compressive strain ε yy reaches the −56.7 µε. Furthermore, the stress factor applies the greatest influence on the floor, the crown and on the side wall in the respective order. Comparing the results before and after operation, the variation of longitudinal stress and strain of Section 2 is less and within 5.5 µε at each place. The lateral strain changes in the range of 3 µε. It shows that the coupling effect has been markedly attenuated at the distance 12 m from the bracket. For Section 3, there is little change in the longitudinal strain before and after operation. It can be seen that for the lining structure away from the fixed bracket, the influence of the coupling effect is mainly reflected by the single temperature field in the form of circumferential stress and strain. Results show that its influence is very limited and its maximum lateral strain can reach −5.6 µε at the foot of sidewall. According to the results of numerical simulation, the internal force distribution of the lining structure can be obtained as Table 5 shows. It can be plugged into the performance functions in reliability calculation. The results indicate that due to the multi-field coupling effect in the operation period, the design of the heat-supplying tunnel at the present stage is blind only by enlarging the whole section thickness. The floor of the tunnel under the fixed bracket should be the key part. In this paper, the tensile stress and strain are positive and the compressive stress and strain are negative.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Section 1, Section 2 and Section 3 is −65 με, −6.4 με and −0.7 με respectively and the results of numerical simulation are −56.7 με, −5.5 με and −0.8 με. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 με, 4.9 με and 1.9 με respectively and the results of numerical simulation are 25.5 με, 6.1 με and 3.1 με. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Section 1, Section 2 and Section 3 is −65 με, −6.4 με and −0.7 με respectively and the results of numerical simulation are −56.7 με, −5.5 με and −0.8 με. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 με, 4.9 με and 1.9 με respectively and the results of numerical simulation are 25.5 με, 6.1 με and 3.1 με. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Section 1, Section 2 and Section 3 is −65 με, −6.4 με and −0.7 με respectively and the results of numerical simulation are −56.7 με, −5.5 με and −0.8 με. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 με, 4.9 με and 1.9 με respectively and the results of numerical simulation are 25.5 με, 6.1 με and 3.1 με. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Section 1, Section 2 and Section 3 is −65 με, −6.4 με and −0.7 με respectively and the results of numerical simulation are −56.7 με, −5.5 με and −0.8 με. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 με, 4.9 με and 1.9 με respectively and the results of numerical simulation are 25.5 με, 6.1 με and 3.1 με. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Section 1, Section 2 and Section 3 is −65 με, −6.4 με and −0.7 με respectively and the results of numerical simulation are −56.7 με, −5.5 με and −0.8 με. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 με, 4.9 με and 1.9 με respectively and the results of numerical simulation are 25.5 με, 6.1 με and 3.1 με. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Section 1, Section 2 and Section 3 is −65 με, −6.4 με and −0.7 με respectively and the results of numerical simulation are −56.7 με, −5.5 με and −0.8 με. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 με, 4.9 με and 1.9 με respectively and the results of numerical simulation are 25.5 με, 6.1 με and 3.1 με. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Field Monitoring
Corresponding to the changes in the longitudinal and lateral strain of concrete before and after the operation respectively, Figure 11 shows the measured value of the coupling effect on the concrete strain. It can be found that: 1. the variation of longitudinal and lateral strain on concrete is roughly consistent with the increase of distance from the center of the fixed brackets and the decrease of strain. In the range of 6 m, the longitudinal strain on the floor attenuates by about 60% and the lateral strain attenuates by 50%. Therefore, the 6 m range of the fixed bracket can be considered as the focus of monitoring and protection in the operation period of the heat-supplying tunnel and the coupling effects can be less beyond 6 m. 2. For longitudinal strain, taking the floor as an example, the average strain of the measured floor in Sections 1-3 is −65 µε, −6.4 µε and −0.7 µε respectively and the results of numerical simulation are −56.7 µε, −5.5 µε and −0.8 µε. 3. Similarly, for lateral strain, the average strain on the floor is 21.2 µε, 4.9 µε and 1.9 µε respectively and the results of numerical simulation are 25.5 µε, 6.1 µε and 3.1 µε. 4. Comparing the computed results with the monitored results, it is not difficult to find that the results of the numerical simulation are consistent with the measured laws in general. Therefore, the reasonableness of the numerical simulation in said materials has been verified and its result can be used for reliability calculation and analysis.

Failure Probability and Reliability Index
In the study of reliability, it is generally believed that when the number of samples is large enough, the accuracy of the Monte Carlo method can be high enough to be a contrast standard. When the general subset simulation method and its improved method are used, the conditional failure probability P0 = 0.1 and 1000 sample points were pumped into each layer. The failure probability of each section and the corresponding reliability index beta were obtained. The results were shown in Table 6, compared with the traditional Monte Carlo method.

Failure Probability and Reliability Index
In the study of reliability, it is generally believed that when the number of samples is large enough, the accuracy of the Monte Carlo method can be high enough to be a contrast standard. When the general subset simulation method and its improved method are used, the conditional failure probability P 0 = 0.1 and 1000 sample points were pumped into each layer. The failure probability of each section and the corresponding reliability index beta were obtained. The results were shown in Table 6, compared with the traditional Monte Carlo method.
The results of Sections 1 and 2 in Table 6 show that the floor near the fixed bracket of the heat-supplying tunnel is the section with the lowest reliability in the structure. The reliability index of the floor in Section 1 calculated by Monte Carlo method, improved subset simulation method and general subset simulation method are 3.04, 2.93 and 2.91 respectively. For the general heat-supplying underground structure, the target design reliability index β t = 3.2. In engineering, 0.85β t = 2.72 is generally used as a control reliability index of the tunnel [24]. Therefore, it is not difficult to find the great influence of the coupling effect on the reliability of the floor of the heat-supplying tunnel and it is necessary to reinforce the floor near the fixed bracket.

Accuracy and Dispersion
Accuracy and dispersion are the most important evaluating indicators of a reliability analysis method. The reliability indexes of all the sections calculated in Table 6 are extracted to integrate and contrast the calculation errors of the general subset simulation method and its improved method, as shown in Figure 12.
In Figure 12, the scatter points of the improved subset simulation approach are denser, approaching the diagonal line with a slope of 1. The slope of the fitting curve k reaches 0.98 over the general method's 0.88 and the correlation coefficient 0.98 is also greater than 0.9. It shows that the accuracy of the improved method is higher and the degree of dispersion is lower.
The key to the success of the improved method lies in the low deviation of the Hua-Wang point set. The existing methods, such as Monte Carlo and general subset simulation, are based on PRN. Table 6 indicates that the number of sample points required for the general subset simulation method and the improved method is in the range of 3000 to 6000 in this study. A fewer number of samples were derived from the reliability analysis. The six random variable parameters in Table 2 brings a high dimension problem for PRN. High dimension and small sample size can lead to the distribution of PRN points inconsistent with the true distribution. Against this backdrop, the general subset simulation method has a higher probability of producing erroneous results than that of the improved method as Table 6 shows.
These cases indicate that the subset simulation method improved by the introduction of the Hua-Wang point set is optimized in the calculation error. It could be a preferred method to ensure accuracy and dispersion. of the floor in Section 1 calculated by Monte Carlo method, improved subset simulation method and general subset simulation method are 3.04, 2.93 and 2.91 respectively. For the general heat-supplying underground structure, the target design reliability index βt = 3.2. In engineering, 0.85βt = 2.72 is generally used as a control reliability index of the tunnel [24]. Therefore, it is not difficult to find the great influence of the coupling effect on the reliability of the floor of the heat-supplying tunnel and it is necessary to reinforce the floor near the fixed bracket.

Accuracy and Dispersion
Accuracy and dispersion are the most important evaluating indicators of a reliability analysis method. The reliability indexes of all the sections calculated in Table 6 are extracted to integrate and contrast the calculation errors of the general subset simulation method and its improved method, as shown in Figure 12.
In Figure 12, the scatter points of the improved subset simulation approach are denser, approaching the diagonal line with a slope of 1. The slope of the fitting curve k reaches 0.98 over the general method's 0.88 and the correlation coefficient 0.98 is also greater than 0.9. It shows that the accuracy of the improved method is higher and the degree of dispersion is lower.
The key to the success of the improved method lies in the low deviation of the Hua-Wang point set. The existing methods, such as Monte Carlo and general subset simulation, are based on PRN. Table 6 indicates that the number of sample points required for the general subset simulation method and the improved method is in the range of 3000 to 6000 in this study. A fewer number of samples were derived from the reliability analysis. The six random variable parameters in Table 2 brings a high dimension problem for PRN. High dimension and small sample size can lead to the distribution of PRN points inconsistent with the true distribution. Against this backdrop, the general subset simulation method has a higher probability of producing erroneous results than that of the improved method as Table 6 shows.
These cases indicate that the subset simulation method improved by the introduction of the Hua-Wang point set is optimized in the calculation error. It could be a preferred method to ensure accuracy and dispersion.

Economic Efficiency
Economic efficiency is also an important evaluating indicator for the improvement of any reliability calculation methods. Reducing the sample quantity and improving the sampling efficiency

Economic Efficiency
Economic efficiency is also an important evaluating indicator for the improvement of any reliability calculation methods. Reducing the sample quantity and improving the sampling efficiency can be the basis for increasing economic efficiency. The CPU time calculated by each method of Section 1 is extracted as (1) Sampling quantity. The scientific counting method is generally used to record the failure probability P f , P f = a × 10 −b . The number of required sample points is usually between 10 b+1~1 0 b+2 in accordance with the experience of the existing Monte Carlo sampling times, N M−C ≥ 10 b+1 . Looking at the algorithm principle and Table 6, it can be concluded that the subset simulation method and its improved method need a sample number of N SS = N 0 × b. Therefore, the analytic relation between the number of sample points N SS and the number of Monte Carlo sample points N M-C for the subset simulation method is as follows: where, N 0 is the number of sampling points per layer in the subset simulation method; N 0 = 1000 in this paper. In Table 6, the failure probability of the sidewall is the lowest and that of the floor is the highest. In Figure 13, the maximum time utilization ratio (TUR) of the improved subset simulation method at the side wall can reach 158.8. The number of sampling points in the Monte Carlo method increases sharply when the failure probability is low, while the sample points needed by the subset simulation method increase little. Under these circumstances, the efficiency of the improved method is significantly reflected. However, owing to the required number of sample points by these three methods at the floor of tunnel is not large, the advantage of subset simulation is not obvious because of the complexity of per step logic operation.
(2) Sampling efficiency. From a global point of view, general subset simulation and the improved subset simulation have the same number of sampling points. Therefore, the efficiency of the sampling method determines the whole calculation efficiency. The results obtained in Figure 13 actually show the difference in efficiency of Hua-Wang point set and PRN sampling method. The analysis results suggest that in a six-dimensional space of 1000 sample points per layer, PRN method has a greater disadvantage in convergence speed than Hua-Wang point set. The improved method can save computing time to the reduction of 27.7%, 28.1%, 30.4%, 27.8% and 30.8% respectively at the locations in Section 1. In the analysis of the reliability of complex tunnel structures, the benefits "Lower deviation and correlation in small scale sample points high dimensional space" of Hua-Wang point set contribute to the economic efficiency more prominently.
In a complex underground engineering, a great many sections demand to be selected and the performance functions are repeatedly invoked. It often takes more than ten hours or even days to operate on a microcomputer when using the PRN sampling method. The above mentioned results prove that the subset simulation method modified by the Hua-Wang point set can be a preferred method to ensure the economic efficiency. Figure 13. can be the basis for increasing economic efficiency. The CPU time calculated by each method of Section 1 is extracted as Figure 13.
(1) Sampling quantity. The scientific counting method is generally used to record the failure probability Pf, Pf = a × 10 −b . The number of required sample points is usually between 10 b+1 ~ 10 b+2 in accordance with the experience of the existing Monte Carlo sampling times, NM−C ≥ 10 b+1 . Looking at the algorithm principle and Table 6, it can be concluded that the subset simulation method and its improved method need a sample number of NSS = N0 × b. Therefore, the analytic relation between the number of sample points NSS and the number of Monte Carlo sample points NM-C for the subset simulation method is as follows: NSS ≤ N0•(lgNM-C − 1) (12) where, N0 is the number of sampling points per layer in the subset simulation method; N0 = 1000 in this paper. In Table 6, the failure probability of the sidewall is the lowest and that of the floor is the highest. In Figure 13, the maximum time utilization ratio (TUR) of the improved subset simulation method at the side wall can reach 158.8. The number of sampling points in the Monte Carlo method increases sharply when the failure probability is low, while the sample points needed by the subset simulation method increase little. Under these circumstances, the efficiency of the improved method is significantly reflected. However, owing to the required number of sample points by these three methods at the floor of tunnel is not large, the advantage of subset simulation is not obvious because of the complexity of per step logic operation.
(2) Sampling efficiency. From a global point of view, general subset simulation and the improved subset simulation have the same number of sampling points. Therefore, the efficiency of the sampling method determines the whole calculation efficiency. The results obtained in Figure 13 actually show the difference in efficiency of Hua-Wang point set and PRN sampling method. The analysis results suggest that in a six-dimensional space of 1000 sample points per layer, PRN method has a greater disadvantage in convergence speed than Hua-Wang point set. The improved method can save computing time to the reduction of 27.7%, 28.1%, 30.4%, 27.8% and 30.8% respectively at the locations in Section 1. In the analysis of the reliability of complex tunnel structures, the benefits "Lower deviation and correlation in small scale sample points high dimensional space" of Hua-Wang point set contribute to the economic efficiency more prominently.
In a complex underground engineering, a great many sections demand to be selected and the performance functions are repeatedly invoked. It often takes more than ten hours or even days to operate on a microcomputer when using the PRN sampling method. The above mentioned results prove that the subset simulation method modified by the Hua-Wang point set can be a preferred method to ensure the economic efficiency.

Conclusions
In this paper, preliminary conclusions and recommendations are as follows: (1) The multi-field coupling effects in the operation period attenuate with the increase of the distance between the structure section and the fixed bracket. The tunnel lining structure in the range of the fixed bracket at 6 m is greatly influenced by the coupling and the floor has the lowest reliability. In the process of design, construction and maintenance of heat-supplying underground structures, special design and protective measures for the structural floor of the fixed bracket 6 m should be put forward according to the operation characteristics of heat-supplying engineering.
(2) This paper improves subset simulation based on Hua-Wang point set. The results show that this method can save computation time and improve computing efficiency on the premise of guaranteeing the computation error and its efficiency is more significant when the failure probability is lower. Considering that the failure probability of tunnel engineering is generally small, the improved subset simulation method has a certain application prospect in the calculation of tunnel engineering reliability.
(3) The reliability analysis system of the heat-supplying tunnel in the operation period established in this paper can provide some theoretical support for the hidden transformation of the general heat-supplying underground engineering.