Mechanical Parameter Inversion in Sandstone Diversion Tunnel and Stability Analysis during Operation Period

A large number of experimental studies show that the mechanical parameters of deep buried surrounding rock show significant attenuation characteristics with the increase of strain from the rheological acceleration stage to the attenuation stage. However, the existing numerical models all take mechanical parameters as constants when describing the rheological behavior of surrounding rocks, which can only be applied to the stability analysis of the shallowly buried tunnel. Therefore, this work proceeding from the actual project, improved the sandstone rheological constitutive model and optimized the algorithm of parameter inversion, and put forward a long-term stability analysis model that can accurately reflect the rheological characteristics of surrounding rocks under the complex geological condition including high stress induced by great depth and high seepage pressure. In the process, a three-dimensional nonlinear rheological damage model was established based on Burgers rheological model by introducing damage factors into the derivation of the sandstone rheological constitutive model to accurately describe the rheological behaviors of the deep buried tunnel. And BP (Back Propagation) neural network optimized by the multi-descendant genetic algorithm is used to invert the mechanical parameters in the model, which improves the efficiency and precision of parameter inversion. Finally, the rheological equation was written by using parametric programming language and incorporated into the general finite element software ANSYS to simulate the rheological behavior of the tunnel rock mass at runtime. The results of the model analysis are in good agreement with the monitoring data in the later stage. The research results can provide a reference for the stability analysis of similar projects.


Introduction
At present, the cross-basin water diversion project, cross-river and cross-sea passageway, and large-scale hydropower project under construction and planning in China often involve the problem of deep-buried tunnels, which has become the main project, or control project of the whole project construction. For example, the Wawushan Hydropower Station requires the construction of a high dam and deep excavation of a diversion tunnel, with the maximum buried depth of about 510m, and abundant underground water. In such a high stress and high external water pressure environment, the rheological characteristics of the surrounding rock will become very significant and complex, especially the mechanical parameters of the surrounding rock which shows remarkable attenuation with the increase of creep strain. The long-term stability analysis of the surrounding rock is bound to be a prominent engineering problem. Many scholars have made special studies on these problems and technical difficulties. Chen Guoqing et al. [1] used the intelligent optimization method to obtain the rheological model of the combination of viscoelastic model and CWFS (bond strength weakening -friction strength strengthening) model, which was used for the long-term stability analysis of the diversion tunnel of Jinping II Hydropower Station. Hu Quanguang et al. [2] established a power exponential creep model applicable to argillaceous siltstone through the indoor triaxial creep test and analyzed the long-term stability of the deep diversion tunnel of the Pakistan N-J hydropower station. Liu Zhiyong et al. [3] used the power-law rheological equation to describe the rheological behavior of mudstone, schist and other soft rock mass to predict the long-term stability of surrounding rock of the diversion tunnel. Based on the fluid-solid coupling theory, Jinpeng et al. [4] analyzed the long-term stability of tunnel surrounding rock by using the finite element software ABAQUS. Feng Shiguo et al. [5] used the power function creep mechanical model to carry out a three-dimensional simulation study on the long-term stability of tunnel surrounding rock. Based on the theory of unloading rock mass mechanics, Ru Long et al. [6] used the FLAC3D finite difference program to carry out three-dimensional simulation research on the long-term stability of tunnel surrounding rock.
A large number of experimental studies show that the mechanical parameters of the deep buried surrounding rock show significant attenuation characteristics with the increase of strain from the rheological attenuation stage to the acceleration stage [7][8][9][10][11]. However, the existing models for predicting the long-term stability of surrounding rocks all take mechanical parameters as constants for calculations and their models can only be applied to environments where the depth of surrounding rocks is relatively shallow. Direct application of existing research results and theories to the long-term stability analysis of deep buried surrounding rock has difficulty meeting the needs of practical engineering. Therefore, taking the complex environment of high in-situ stress and the characteristics of the sandstone rock mass of the deep buried diversion tunnel of Wawushan Hydropower Station as the research background, damage factors were introduced into the deduction of the sandstone rheological constitutive model based on Burgers rheological mode, and a three-dimensional nonlinear rheological damage model was established to describe the rheological characteristics of sandstone accurately. And BP [12] neural network optimized by multi-generation genetic algorithm [13][14][15] is used to invert the mechanical parameters in the model, which ensures the efficiency and precision of parameter inversion. Finally, the rheological equation was written by using parametric programming language and incorporated into the general finite element software ANSYS [16][17][18][19][20] to simulate the rheological, mechanical behavior of the tunnel rock mass under the condition of water filling, and the surrounding rock displacement, stress, plastic zone and support stress of the tunnel were obtained, so as to comprehensively evaluate the long-term stability of the project.

Project Overview and Geological Conditions
Wawushan Hydropower Station is located in Hongya County, Meishan City, Sichuan Province. The exact location of the study area is shown in Figure 1. The hub consists of a dam, a spillway tunnel on the left and right bank, a water diversion system, and a factory building. The water diversion tunnel is circular, with an inner design diameter of 6.2m, a lining thickness of 40 ~ 80cm, and a total length of 4813m. The overlying rock mass is generally buried at a depth of 150 ~ 450m, and the maximum buried depth of the tunnel is about 510m. It has the characteristics of large buried depth and long tunnel line and has become a key control project of Wawushan Hydropower Station.

Figure 1. Location map of case study
According to geological survey data, surrounding rock types include: Class III surrounding rock, length 3260 m, station number: 0+239.97~3+499.97, mainly composed of thick layer to medium-thick layer sandstone, some thin sandstones sandwiching argillaceous rocks; Class IV surrounding rock, length 1425 m, station number: 0+079.97~0+239.97, 3+499.97~4+439.97, 4+499.97~4+825.292, mainly composed of medium-thick argillaceous siltstone or medium-thick silty mudstone; V-type surrounding rock, length 140m, station number: 0+000~0+079.97; 4+439.97~4+499.97, mainly composed of weathered and unloaded medium-thick layer of silty siltstone with thin layer silty mudstone. Classes Ⅲ and class Ⅳ accounted for about 97% of the total length of the surrounding rock. The rock is dense, hard, and complete. There are still some V-type surrounding rocks in the III and IV surrounding rocks. The main mechanical parameters of rock mass are shown in Table 1.  Figure 2 shows the layout of monitoring points and lines of the tunnel. Combined with engineering geology and construction conditions, the main monitoring scheme is determined to be hole convergence and anchor stress monitoring. The deformation monitoring cross-section, according to the principle of subsection control, selects six deformation monitoring sections in the surrounding rock of Class III and IV [21][22][23][24]. The monitoring points are embedded in the design section near the tunnel face, and the monitoring points are arranged according to the five measuring points and eight measuring lines of the circular section. Based on the fixed point arbitrary triangle method, the displacement of the surrounding rock of the diversion tunnel is calculated. The fixed point arbitrary triangle method assumes that the monitoring point 5 is at the origin of the coordinate and does not deform. So monitoring point 1 and monitoring point 5 are in the same horizontal position, and the deformation of monitoring points 2, 3, and 4 are decomposed into horizontal and vertical displacement.

Figure 2. The layout of monitoring points and monitoring lines of the tunnel
In this paper, the downstream surface (K3+950) of 3# branch tunnel was selected as the research object. The section rock is a type III surrounding rock, accounting for 69% in three main types surrounding rock, with certain representativeness; only temporary support will be provided after the excavation of the section is stopped. During the period, the stress state of the surrounding rock does not change much, and it is easy to simulate and analyze; After the excavation activity is stopped, the convergence value of each line is still increasing, and the rheological properties of the sandstone are fully expressed, which is beneficial to the inversion analysis of the rheological parameters of the surrounding rock.

Constitutive Model and Mechanical Parameters of Sandstone
Rock samples show instantaneous elastic deformation under external load and show obvious rheological properties. Among the existing constitutive models, the Burgers model can best reflect the rheological and mechanical properties of sandstone. However, a large number of experimental studies showed that the mechanical parameters of sandstone decreased with the increase of strain from the rheological attenuation stage to the acceleration stage. It was not practical to calculate the mechanical parameters as a constant in the Burgers model. Therefore, based on the Burgers model, a three-dimensional nonlinear rheological damage model was established to accurately describe the rheological characteristics of sandstone by introducing damage factors into the derivation of the rheological constitutive model of sandstone.
The rock is mainly in the state of compression. Therefore, the damage mainly affects the shear modulus, and its influence on the volume modulus can be negligible. Because of the generation and propagation of cracks, the viscosity coefficient of rock changes. So the damage also affects the viscosity coefficient. The component composition of the improved Burgers model is shown in Figure 4.

Figure 4. Creep damage model of sandstone
According to the improved Burges model in Figure 4, the rheological equation of sandstone under three-dimensional conditions can be deduced as follows: Where GK=viscoelastic shear modulus; ηK=viscoelastic shear modulus; GM =shear modulus; ηM=shear modulus; K= bulk modulus; K and GM can be expressed as: Where EM = elasticity modulus; μ=poisson ratio.
In Equation 2, damage evolution equation D is defined as: Where, A, B, and C are damage parameters, which are determined by fitting the test curve; εthr is the strain damage threshold, which is 8.38×10 -3 in this paper; εequ is the equivalent strain, calculated as follows: By combining (1) and (2), the constitutive equation of rheological damage under the three-dimensional condition of sandstone can be obtained:  Where, η1, E1, η1, G2 and A need to be determined by parameter inversion. Based on geological data, triaxial test results, and engineering experience, the value ranges of four parameters are determined, as shown in Table 2.

BP Neural Network Model Optimized by Multi-generation Genetic Algorithm
The BP neural network model optimized by the multi-generation genetic algorithm is designed and programmed through MATLAB. The establishment of the network model is mainly divided into two parts.

1) BP neural network optimized by multi-generation genetic algorithm.
The main ideas of BP neural network optimization by the multi-generation genetic algorithm are as follows. To improve the probability of offspring producing more fitness individuals, the parent population composed of weights and thresholds crossed to produce multiple offspring populations with more than the number of parents. The search interval of network weights and thresholds is adjusted adaptively according to the change of fitness of current multi-offspring individuals. Compared with the traditional genetic algorithm optimization neural network, the optimization BP neural network algorithm based on multi-generation genetic algorithm has higher learning accuracy and faster convergence speed.
Interval adaptive adjustment process. Set the initial interval of the weight was [a,b], and randomly generate the initial weight population ( w1, w2,..., wn ) within the interval, with the scale of n. The network error of the weight w at point b is better than that at point a. Suppose the optimal weight corresponding to the minimum network error is w, and w'>b. Because the optimal weight w is beyond the range of the interval, the traditional genetic algorithm cannot search for the optimal weight. However, after several iterations, as the network error of point b is better than point a, the individuals in the population will gradually move closer to point b, and the distance between w1, the worst individual in the population, and point a becomes larger and larger. Therefore, a new search interval [a1,b1]needs to be formed through Equation 7 to complete an interval movement. For w'<a is similar to the above method.

2) Formulate the neural network structure
The structure of BP neural network includes an input layer, an output layer, and a hidden layer. The number of nodes in the input layer and the output layer is determined according to the actual problem, while the number of nodes in the hidden layer needs to be determined by Equation 8. a n m N    Where, is a constant between 1 and 10, m and n are the number of input layer nodes and the number of output nodes, respectively.
According to the layout requirements of BP neural network structure, the mechanical parameters are taken as the target vector of the neural network, and the number of output nodes is 5. The convergence value of displacement monitoring is taken as the input vector of the neural network, and the number of input nodes is 350. According to Equation 8, the initial value of hidden layer nodes is set as 20 layers.
Based on the uniform design method, the uniform design scheme U11(11 7 ) was selected to carry out the experimental design combination of the mechanical parameters to obtain the parameters of 11 groups of different combinations. The displacement values of the monitoring points corresponding to the 11 groups of mechanical parameters were calculated by numerical analysis. The specific simulation process can be referred to in section 6.
The 11 groups of mechanical parameters and the corresponding displacement values of monitoring points jointly constitute the training samples of the artificial neural network, providing initial samples for the inversion of mechanical parameters. Due to limited space, only the training samples for combination η1=1.03×10 18 Pa·s, E1=39Pa, η2=1.21×10 16 Pa·s, G2=9.80×10 10 Pa, and A=0.02 are listed for reference, as shown in Figure 6.

Mechanical Parameter Inversion and Rationality Verification
The BP neural network method optimized by the multi-generation genetic algorithm was designed and programmed by MATLAB. The 70-day monitoring values of the four monitoring points of K3+950 monitoring section were imported into the initial training sample as an input vector and four mechanical parameters were output. The results were as follows: η1 =7.11×10 17 Pa·s, E1 =39 Pa, η2 =1.27×10 16 Pa·s and G2 =1.80×10 11 Pa. Other mechanical parameters required for the establishment of the rheological model are shown in Table 3.  The inversion parameters were imported into the calculation model, and the rheological displacement of each monitoring point in the section for 70 days was obtained through numerical simulation. Figure 7 shows the comparison between the displacements calculated value and the monitored value. The calculation results of each monitoring point have a good fitting degree with the monitoring data, and the relative error is controlled within 10%, which meets the needs of the project. Therefore, the rheological damage model of sandstone can be used to describe the deformation prediction of the surrounding rock at the downstream of no.3 branch tunnel of the diversion tunnel of the hydropower station during the operation period.

Numerical Model
ANSYS software based on the MOTIF graphical user interface, intelligent menu guidance and help, etc. apply parametric design language APDL to write parametric user programs, providing users with powerful solid modeling and grid division tools, and can analyze various physical field quantities.
The rheological equation was written by using parametric programming language and incorporated into the general finite element software ANSYS. The APDL is called once every iteration of the program, and the stress and strain are updated to achieve the introduction of the sandstone rheological model.  Entity modeling and mesh generation. Since the surrounding rock lithology of the diversion tunnel is single in a wide range along the axis, the problem can be simplified as a plane strain problem. The upper and lower boundaries of the mesh model are set to be greater than 10 times the hole diameter to reduce the boundary effect. According to the design data, the section of the downstream side of the 1# branch tunnel of the diversion tunnel is round with an inner diameter of 6.2 m and a lining thickness of 0.6m. The calculation range of the model is set as 80m×80m, and the simulation unit of surrounding rock and lining is planar element Plane42. The finite element mesh of the model is shown in Figure 8.
 Normal constraints were applied on both sides and bottom.
 Establish the initial stress field before excavation. The initial stress field is dominated by the gravity stress field.
 Establish the initial stress field during operation. After the simulation of the given initial stress field and external water pressure is completed, the full section excavation simulation and tunnel support simulation are carried out successively. Since the interval between excavation and lining construction is long and the surrounding rock deformation is large, creep simulation should be carried out during the excavation and lining construction and after completion.
Through the above construction simulation, the initial stress field during the operation period is established, and on this basis, the internal water pressure is exerted on the tunnel. Finally, the rheological conditions of the tunnel in different periods (T=0d, T=10d, T=30d, T=60d, T=120d, T=240d, T=360d) are simulated.  To systematically analyze the rheological law of surrounding rock, the displacement cloud charts at different moments are analyzed to obtain the rheological displacement of characteristic points of surrounding rock at the corresponding moment, as shown in Table 4. From the quantitative value analysis, the overall deformation of the surrounding rock is small. After 360 days of rheological calculation, the displacement of the characteristic point of the surrounding rock is about 4.5-6.9 mm. From the perspective of the rheological trend, the rheological trend of each characteristic point of surrounding rock of the tunnel is relatively consistent, and the initial deformation is large. The deformation amount in the first ten days is up to 2.6-4.2 mm, and then the deformation rate gradually decreases.

Analysis of Calculation Results
The deformation enters the rheological attenuation period from 10d to 30d. The deformation tends to be stable after the 30th day and enters into the stable creep stage. The average deformation rate is 0.005 mm/d ~ 0.006 mm/d. The analysis results show that the deformation of the surrounding rock of this section can tend to be stable within 30 days, the stable deformation rate is very low, and there is no new plastic zone around the tunnel with the passage of rheological time. Therefore, the rheology deformation of surrounding rock during operation period meets the requirements of longterm stability of surrounding rock of diversion tunnel.  Table 5 gives the statistical results of the maximum stress around the cave at different rheological times. In terms of quantitative value, the surrounding rock and lining structure during the operation period are mainly under pressure. The tensile stress is inversely proportional to the rheological time, and the compressive stress is proportional to the rheological time. The maximum compressive stress of 3rd principle stress at each moment appeared at the right side of the bottom and the left side of the roof, and there was no tensile stress area around the rock. The maximum compressive stress of 1st principle stress appeared in the 3 m ~ 4 m area on both sides of the surrounding rock at each moment. The tensile stress mainly appeared on the right side of the lining structure floor, and the stress range decreases with the passage of rheological time. Therefore, the surrounding rock stress during the operation period meets the requirements of the long-term stability of the surrounding rock of the diversion tunnel.

Conclusion
Given the complex environment of high in-situ stress and the characteristics of sandstone rock mass in the deep water diversion tunnel of Wawushan Hydropower Station, the rheological damage model is used to describe the rheological characteristics of sandstone. BP neural network optimized by the multi-generation genetic algorithm was used to invert the mechanical parameters in the model. Finally, the rheological equation was written, and ANSYS was connected to simulate the rheological, mechanical behavior of the tunnel during its operation. The calculation results show that the diversion tunnel has good overall stability. The displacement value changes little, the deformation convergence speed is fast, and the stable deformation rate is low. The stress value tends to be uniform, the stress concentration is gradually reduced, the plastic area is small, and no new plastic area appears after entering the rheological stability period. The stress of the support structure is also safe and reliable.

Funding
This work was supported by the Fund of National Key R&D Program of China "high-efficient development and utilization of water resource." (No. 2018YFC0406803).