The Baihetan Hydropower Station is located in Qiaojia County, Yunnan Province, China. It is a hydropower station of over ten million kilowatts similar to the Three Gorges Hydropower station or Xiluodu Hydropower Station. The bedrock of the dam site is dominated by cryptocrystalline basalt, almond basalt, and breccia lava, with hard lithology; the columnar joints are developed in some lithologic sections. The rock masses of the abutment and foundation of the dam are poorly developed, and the joints, interlayer, and intralayer dislocation zones are widely distributed. The faults with steep inclination and structural planes cut each other to form multiple wedges, which makes the rock masses exhibit high anisotropy characteristics (as shown in Figure 1), and affects the stability of the arch dam and safety of the underground cavern group.

As shown in Figure 2, the Baihetan hydropower station has five diversion tunnels; three diversion tunnels are set on the left bank and two diversion tunnels are set on the right bank, with elevations ranging from 574 to 605 m. The red part in Figure 2 represents the section of the diversion tunnel where the columnar jointed basalt is exposed, and its range is relatively large. Figure 3 shows the typical columnar jointed basalt in the diversion tunnel. The irregular columns have obvious contours, and their diameters range from 13 to 25 m. The inclination angle of columns are 70°–80°. The cross section of the columns are mainly irregular pentagonal and quadrilateral. Parallel cylindrical longitudinal micro-cracks are developed in the columns.

In the excavation process of the diversion tunnels, unloading phenomena such as structural plane opening, tensile fracture, and rock mass relaxation appear in the surrounding rock near the excavation face. After the adjustment of the stress field near the surrounding rock, the original structural planes open and new fractures form. The widths of unloading fractures gradually decrease from the entrance to the end of the tunnel. Figure 4 shows the comparison of relaxation characteristics of the columnar jointed basalt in different sections of the diversion tunnel. The unloading fractures in the entrance section are generally developed, and are filled with gravel, rock debris, and secondary mud (Figure 4A). The relaxed columnar jointed basalt shown in Figure 4B maintains a certain integrity. Although most of the original columnar joints are open and the columns are relaxed, the rock blocks are still interlocking.

Taking the No.4 diversion tunnel as an example, this paper studies the unloading relaxation of columnar jointed basalt. The diversion tunnel is excavated in three steps from top to bottom. The upper layer is 9 m high, the middle layer is 9 m high, and the lower layer is 6.2 m high. After excavation, the three-point displacement gauges are used to monitor the deformation of multiple sections of the tunnel. The deformation displacement was recorded at the depth of the hole at 0, 2, and 9 m to obtain excavation response characteristics of the columnar jointed rock mass. The layout and deformation of the multi-point displacement gauges of the 1 + 075 section are shown in Figure 5. After excavation, the top arch has the largest deformation, the shallow deformation reaches 40.71 mm, and the deep deformation also reaches 29.89 mm, which is due to the stress concentration at the arch top. The deformations of the side walls are smaller than that of the arch top. The displacement gauge in Myd6-2 measuring hole is damaged and has no recording. The displacements recorded by gauges in the remaining three monitoring holes show that the deformations of the side walls gradually decrease from shallow to deep. And the deformation of the lower part of the side wall is slightly higher than that of the upper part of the side wall, and the maximum deformation appears in the lower part of the left side wall, which is 22.38 mm.

Figure 6 shows the displacement sequence curves at the lower part of the left-side wall and the upper part of the right-side wall. At the initial stage after excavation, the deformation rate of the surrounding rock continues to increase. After a period of time, the deformation rate slows down and gradually becomes stable, and the deformation rate of the deep part is lower than that of the shallow part.

Considering the columnar jointed rock mass is a transverse isotropic material, the above rock mass has 10 mechanical parameters including deformation modulus E 1 , E 3 , Poisson’s ratio γ 12 , γ 13 , shear modulus G 12 , G 13 , cohesion C , internal friction angle φ , and anisotropic azimuths d d , d i p , where the anisotropic azimuths are constant values. Therefore, there are still eight unknown mechanical parameters. For the convenience of engineering application, the parameters in the model need to be selected before back analysis.

Through parameter sensitivity analysis, we know that deformation modulus E 1 , E 3 and shear modulus G 13 are more sensitive to displacement, so these three mechanical parameters are selected as the parameters of back analysis.

The displacements of monitoring points MYD6-1 (0 m), MYD6-1 (2 m), and MYD6-5 (2 m) were selected as the target values of the back analysis. The three parameter levels of 85, 60, and 35%, and the orthogonal sample scheme shown in Table 1 are considered. The anisotropic elastoplastic model is used for calculation, and the simulation results of the sample are shown in Table 2.

The mapping relationship between deformation modulus E 1 , E 3 , shear modulus G 13 and displacement is established by support vector machines (Zhao et al., 2003). On this basis, a differential evolution algorithm is used as an optimization tool to search and monitor the deformation parameters of columnar jointed rock mass corresponding to displacement (Storn and Price, 1997). Table 3 is the target value of back analysis, and Table 4 is the result of back analysis.

The results of back analysis are used as mechanical parameters to simulate the diversion tunnel deformation, and the calculated displacements of the monitor points are recorded and compared with the measured displacements (Figure 7). The comparison results show that the displacements calculated by the results of back analysis are very close to the measured displacements, therefore, the mechanical parameter obtained by back analysis can be used to determine the relaxation scale of the diversion tunnel.

Pietruszczak and Morz extended the classical isotropy criterion to anisotropy by combining the spatial distribution of strength parameters and established anisotropy yield criterion based on the microstructure tensor (Pietruszczak 1999). In the model, each scalar of anisotropy parameter is composed of mixed variables of stress and microstructure tensor. This means the physical and geometric properties of materials and the loading direction can be used as a set of variables.

The microstructure tensor a i j is used to mark the material structure, and its deviatoric tensor is: A i j = ( a i j − 1 3 δ i j a k k ) 1 3 a k k (1) where a k k is the principal value of the microstructure tensor a i j . The loading direction can be calculated by the stress component: { L 1 = σ 11 2 + σ 12 2 + σ 13 2 L 2 = σ 21 2 + σ 22 2 + σ 23 2 L 3 = σ 31 2 + σ 32 2 + σ 33 2 (2)

The generalized loading vector l i can be expressed as: l i = L i / ( L k L k ) 1 / 2 (3)

The anisotropy parameter η = a i j l i l j is introduced to represent the projection of the microstructure tensor on the loading direction.

Ning (2008) proposed an anisotropic yield criterion considering tensile strength: f s = σ 1 − σ 3 N φ + 2 c N φ   ( σ 1 ≥ σ 2 ≥ σ 3 ) (4) where, N φ = 1 + sin ⁡ φ 1 − sin ⁡ φ

Anisotropic Tensile Stress Yields Function: f t = σ t − σ 1 (5)

The shear potential function of columnar jointed rock mass corresponds to the non-associated flow rule, which can be expressed as: g s = σ 3 − σ 1 N ψ (6) where ψ is the dilatancy angle, N ψ = 1 + sin ⁡ ψ 1 − sin ⁡ ψ

The tensile potential function g t corresponds to the associated flow law, which can be expressed as follows: g t = − σ 1 (7)

The strain increment can be expressed as the superposition of elastic part and plastic part: d ε i j = d ε i j e + d ε i j p (8)

The incremental stress-strain relationship of the elastic part is: d σ i j = K i j k l d ε k l e (9) where K i j k l is the elastic stiffness matrix.

The strain increment of the plastic part is determined by the plastic flow rule: d ε i j p = d λ ∂ g ∂ σ i j (10) d λ = 1 H ∂ f ∂ σ i j K i j k l d ε k l (11) H = ∂ f ∂ σ i j K i j k l ∂ g ∂ σ k l (12)

Therefore, the incremental stress-strain relationship can be expressed as: d σ i j = K i j k l e p d ε k l = ( K i j k l − 1 H K i j m n ∂ g ∂ σ m n ∂ f ∂ σ p q K p q k l ) d ε k l (13)

Relaxation factor is introduced to describe the deterioration of elastic modulus and shear modulus of rock after excavation. When the rock has no relaxation, the relaxation factor is 0; when the rock mass is completely relaxed, the relaxation factor is 1. For the rock in engineering, the relaxation factor is a value in the range of 0–1.

The expression of the degradation of elastic modulus used in this paper is E = ( 1 − S ) E 0 (Eberhardt et al., 1999); S is the relaxation factor, and E 0 is the initial elastic modulus. The relaxation factor describing the relaxation degree of rock mass can be defined as: S = 1 − E E 0 .

According to the aforementioned results of back analysis, the deformation modulus E 1 , E 3 and shear modulus G 13 are affected by relaxation. Therefore, three relaxation factors are needed to describe the different deterioration rates of these three mechanical parameters: S 1 = 1 − E 1 E 1 0 (14) S 2 = 1 − E 3 E 3 0 (15) S 3 = 1 − G 13 G 13 0 (16) where S i ( i =1, 2, 3) are relaxation factors, E i 0 and E i ( i =1, 2) are the initial and relaxation values of the two deformation moduli of rock mass, and G 13 0 and G 13 are the initial and relaxation values of the shear moduli of rock mass, respectively.

Considering the influence of the three principal stresses comprehensively, the relaxation threshold is defined from the respect of equivalent strain. When the equivalent strain of rock is greater than the relaxation threshold, the columnar jointed basalt enters the relaxation stage. The equivalent strain is ε ¯ = ε 1 2 + ε 2 2 + ε 3 2 .Here, ε i ( i =1, 2 and 3) are the values of the three principal strains of rock, and the relaxation criterion is: S { = 0   ,   ε ¯ ≤ ε ¯ ∗ >0 ,  ε ¯ > ε ¯ ∗ (17) where ε ¯ ∗ is the relaxation threshold. Referring to the research results of Zhou et al. (Zhou et al., 2009), a relaxation evolution equation directly related to strain is obtained, S = 1 − e − m ( ε ¯ − ε ¯ ∗ ) , where m is the parameters of columnar jointed basalt related to relaxation.

Three relaxation evolution equations can be obtained as: S 1 = 1 − e − m 1 ( ε ¯ − ε ¯ ∗ ) (18) S 2 = 1 − e − m 2 ( ε ¯ − ε ¯ ∗ ) (19) S 3 = 1 − e − m 3 ( ε ¯ − ε ¯ ∗ ) (20)

There are four parameters ( m 1 、 m 2 、 m 3 、 ε ¯ ∗ ) that need to be determined in the relaxation evolution equation.

The corresponding relaxation factor values can be calculated according to the results of back analysis in Table 4. Table 5 is the calculated relaxation factor, which can be used as the target value of back analysis of relaxation evolution equation.

Considering the limitation of the length of this paper, the specific method of back analysis is not given a more detailed description, although one can be found in the literature (Shen and Chun, 2014). The optimal results of back analysis are m 1 = 14 , m 2 = 20 , m 3 = 14 , ε ¯ ∗ = 1 E − 05 .

The average value of three relaxation factors S ¯ = 1 3 ( S 1 + S 2 + S 3 ) is used to determine the scale of unloading relaxation zone, and the value is 0.49.

Figure 8A shows the calculation model of the section of the diversion tunnel in software Flac3d; the model has a total of 14,012 elements. The model size is 60 m*60 m, and the ground stress was applied to the boundaries. The parameters of the model are shown in Table.6.

Figures 8B–E shows the scale of the relaxation zone of the surrounding rock after each step of excavation. After the first step of excavation, the relaxation zone is mainly distributed on both sides of the top arch (Figure 8B). After the second step of excavation, a small amount of relaxation zone appears on the upper part of the top arch, and the relaxation zone on the side walls expands rapidly (Figure 8C). After the third step of excavation, the relaxation zone of the top arch expands into a through relaxation ring with a maximum thickness of 0.65 m. The relaxation zone on the side walls is thicker, and the thickest part of the relaxation zone appears in the middle of the side walls. The maximum thickness of the relaxation zone on the left-side wall is 6.17 m, and the maximum thickness of the relaxation zone on the right-side wall is 6.20 m (Figure 8D). Figure 8E shows the expansion process of the relaxation zone after each step of excavation. It can be seen that the next step of excavation will cause the expansion of the relaxation zone to induce by the previous step of excavation.

Table 7 shows the relaxation depth of the diversion tunnel obtained by acoustic detection. It can be found that the thickness and distribution of the relaxation zone calculated by the customized anisotropic elastic-plastic model are well consistent with the results of the in-situ acoustic detection, especially at the side walls. The simulation results are in good agreement with the measured values. At the top arch, the relaxation depth calculated by the model is slightly different from the measured values, but the error is still acceptable, and fully meets the engineering requirements.