Numerical Modeling of Stress Disturbance Characteristics during Tunnel Excavation

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China School of Engineering Science, University of Chinese Academy of Sciences, Beijing 10049, China School of Resources and Environment, North China University of Water Resources and Electric Power, Zhengzhou 450046, China College of Water Conservancy and Hydropower, Hebei University of Engineering, Handan, Hebei 056001, China


Introduction
It is widely known that tunnel excavation can lead to stress changes and possible damage ahead of a tunnel face. Increased attention has been given to this issue during rock mechanics research and underground excavation projects [1][2][3][4][5][6]. Excavation-induced stress is commonly described by analyzing the progressive development and evolution of the near-field stress distribution and stress paths during the advancement of a tunnel face. Many researchers have carried out studies on stress disturbance during tunnel excavation by means of various methods such as theoretical analysis and field tests. e main representative studies are the Underground Research Laboratory (URL) of Atomic Energy of Canada Limited (AECL) (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), which carried out a large number of field tests and monitoring work for Lac du Bonnet granite [7,8]; the Swedish Nuclear Fuel and Waste Management Company (SKB), which conducted a variety of tests and analysis work based on the Hard Rock Laboratory on the islandÄspö of southeast Sweden [9]. In the zone of excavation disturbance experiments (ZEDEX) conducted by SKB (Sweden), ANDRA (France) and UK Nirex (UK), the research suggests that the damage zone is in the proximity of a tunnel boundary (referred to as the near-field), where irreversible changes of rock mass properties occur. Moving away from the tunnel boundary, rock mass properties in the disturbance zone (relatively, in the far-field) undergo little change and most can be recovered [10]. In China, in regard to the deeply buried diversion tunnels of the Jinping II Hydropower Station, Jiang et al. [11] and Feng et al. [6] performed comprehensive in situ observation tests and dynamic feedback analysis to evaluate the cracking process and stability of surrounding rock; Liu et al. [12] carried out research on the tunnel face effect due to stress changes by monitoring the stress in deeply buried diversion tunnels during TBM tunneling and used the results to determine the range of the damage zone. e distribution and redistribution of excavation-induced stress involves not only the change in stress magnitudes but also the rotation of the principal stress axes as a tunnel excavation is approached. However, the rotation of principal stress axes is normally neglected as tunneling projects continue to progress into deeper and more complex geological environments. A wealth of geotechnical tests show that pure principal stress axis rotation (with no change in magnitudes of principal stresses) can lead to obvious plastic deformation [13]. In recent years, soil constitutive models that consider the principal stress axes rotation, as well as its effect on deformation, have been developed [14]. However, few studies on the principal stress axes rotation of rock mass, especially on rotation due to excavation disturbance, were found. Carter [15] and Germanovich and Dyskin [16] have demonstrated that the microfracturing process, related to damage and rock strength degradation/ failure, is largely controlled by the deviatoric stresses and the orientation of the principal stress axes by analyzing the excavation-disturbed near-field stress zone. Eberhardt [17] analyzed the change in principal stress values and the characteristics of principal stress axes rotation ahead of a tunnel face in circular tunnels during excavation under different geostress conditions. Generally, describing the three-dimensional space change of principal stress orientation following rock tunnel excavation is difficult using a two-dimensional analysis method or graphic representation method. In addition, conventional mechanical indexes such as the average normal stress, equivalent shear stress, and rod parameter are extensively used to analyze the excavationinduced stress ahead of the tunnel face, but these stress indexes do not consider the principal stress axes rotation and its effect, especially the stress state change of the excavationinduced stress.
is paper develops a numerical model of diversion tunnels for the Jinping II Hydropower Station to simulate the rotation of principal stress axes and the change in stress state ahead of the tunnel face by means of finite difference modeling. e pole diagrams prevailing in geology and the stress triaxiality commonly used to indicate the deformation and damage of metal materials are introduced to carry out the initial study on the three-dimensional characteristics and distribution pattern of stress disturbance during the advancement of the tunnel face.

Engineering Case and Numerical Model
e Jinping II Hydropower Station is located in the winding region of the Jinping Bend of theYalong River in Sichuan Province, China, with approximately a 310 m water head along a 150 km river section. As a diversion-type power station, four diversion tunnels were excavated to provide a 4800 MW installed capacity. e four tunnels are the key project simply because they are buried deep, with an average length of approximately 16.67 km, an excavation diameter D of approximately 13 m, a maximum buried depth of approximately 2525 m, and an overburden rock depth of generally 1500-2000 m. erefore, based on the literature and measured geostress and rock mass mechanical parameters [6], this paper takes the diversion tunnel excavation as an example to research stress disturbance near and ahead of the tunnel face. Furthermore, only one of the four diversion tunnels was selected to establish the numerical model for the sake of calculative simplification.
e Chainage 15 + 700 section of No. 2# tunnel is selected as a reference section, with a 100 m long tunnel segment near the reference section as a study case. is segment is buried by 1100 m of thick overburden rock, with a gravity stress of approximately 29.2 MPa. e major principal stress σ 1 is assumed to be parallel to the tunnel axis, and σ 1 is 1.3 times as large as the gravity stress, namely, approximately 38.0 MPa. e intermediate principal stress σ 2 is 1.2 times as large as the gravity stress, namely, approximately 34.8 MPa, forming an angle of 22.5°with the horizontal plane. e minor principal stress σ 3 is the same as the gravity stress, namely, approximately 29.2 MPa, forming an angle of 67.5°with the horizontal plane. e corresponding normal and shear stresses calculated by the elastic mechanical formula are applied on the boundary of the numerical model to meet the magnitudes and directions of the primary in situ stress field. e numerical model and the X, Y, and Z coordinate systems are established, as shown in Figure 1, with axis Z along the tunnel centerline, axis X normal to axis Z in the horizontal plane, and axis Y as vertical. e total model size is 130 m along axis X, 130 m along axis Y, and 100 m along axis Z, with 16,650 elements and 173,316 nodes. e reference section selected is at Z � 60 m, and the schematic advancing of the tunnel face is shown in Figure 2. e cyclical footage is set to be 2 m in the important study region (Z � 40-60 m) and 5 m in the other regions to better reveal the stress changes and principal stress axes rotation with the advancement of the tunnel face. e tunnel size is shown in Figure 3 and some monitoring points are set at typical positions such as the tunnel roof and tunnel wall, both on the tunnel boundary and along a parallel line 1m into the rock mass in the reference section (Z � 60 m). e Mohr-Coulomb model is used and Table 1 provides the material properties used in the elastoplastic model.

Principal Stress Axes Rotation ahead of an
Advancing Tunnel Face e typical principal stress orientation vectors of monitoring point A in the reference section are shown in Figure 4 during the advancement of the tunnel face. e direction of the principal stresses experiences rotations at different degrees when the tunnel face is far from, before, on, and behind the 2 Advances in Materials Science and Engineering reference section.
e stress tensor in the tunnel roof is generally unchanged in comparison to its initial in situ state when the tunnel face is far from the reference section ( Figure 4(a)). As the tunnel face approaches the reference section, the principal stress axes in the tunnel roof rotate significantly from the initial in situ condition (Figure 4(c)). en, the principal stress direction vectors are readjusted to be approximately parallel to the coordination of the numerical model after the tunnel face passes through the reference section, with the σ 1 stress axis rotated to point in the X-axis direction.
However, illustrating the rotation process of the principal stress axes during tunnel excavation in a precise and intuitive way using the direction vectors in Figure 4 is difficult. us, the pole diagram is introduced to represent the principal stress orientation and its rotational change. According to the geographical orientation of the actual project, the X, Y, and Z axes of the coordinate system in Figure 1 point to the north, vertical, and east, respectively. At the same time, η and ξ are assumed to represent the inclination and dip of a principal stress axis, as shown in Figure 5. e angles between the principal stress axis and X, Y, and Z are denoted by α, β, and c, respectively. In terms of the spatial relation, η and ξ can be expressed by α, β, and c as follows:  In the finite difference numerical calculation, three direction cosines cos α, cos β, and cos c are used to describe the direction of principal stress. rough equations (1) and (2), the direction of principal stress can be described in the twodimensional pole diagram by the two direction indexes of inclination and dip. Using this method, information about the principal stress axes at monitoring points is transformed into pole diagrams as the tunnel face moves from 0 m to 60 m, which can reflect the change process of principal stress axes in some detail. e relationship between geographic coordinates and computational coordinates in Figure 5 shows that the diversion tunnel axis extends from west to east. e pole diagrams of monitoring point A and B are shown in Figure 6 during the tunneling period from beginning to reference section. Let "WF" denote the position of the working face. For monitoring point A, the inclination of the σ 1 stress axis, which points to the tunnel axis at WF � 0, rotates gradually from 0°to approximately 60°at which the σ 1 stress axis cuts the tunnel axis when WF � 50 m and decreases rapidly to approximately 0°at WF � 60 m. From WF � 0 to WF � 60 m, the dip angle ξ of the σ 1 stress axis increases constantly from 0°to approximately 34°. e inclination of the σ 3 stress axis, which is perpendicular to the tunnel axis at WF � 0, rotates gradually from 90°to approximately 0°at WF � 60 m, while the dip angle ξ of the σ 3 stress axis changes slightly from 67.5°(WF � 0) to approximately 55°(WF � 60 m). Generally, for the tunnel roof at WF � 60 m, the inclination of the major and minor principal stress axes are approximately parallel to the tunnel axis, and ξ of the major and minor principal stress axes are approximately 34°and 55°, respectively. e inclination of the σ 2 stress axis is basically normal to the tunnel axis, with ξ equal to approximately 0°.

Advances in Materials Science and Engineering
For monitoring point B in the sidewall, the principal stress directions adjustment process is different from monitoring point A, but it shows some similar characteristics. e major and minor principal stress axes ahead of the tunnel face tend to intersect the excavation surface of the tunnel at angles ranging from 30°to 60°, with the intermediate principal stress axis almost parallel to the excavation surface (approximately normal to the tunnel axis). In addition, the principal stress directions will adjust themselves dramatically when the tunnel face is 10-12 m (approximately equal to the tunnel diameter) away from the reference section. When the tunnel face approaches WF � 60 m, the monitoring points and the working face are located at the same cross-section of the tunnel. e three-dimensional space relationship between the principal stress orientation and the excavation surface, which is beneficial to analysis of excavation stability near the working face, as shown in Figure 7, with schematic diagrams of the principal stress orientation and the tunnel boundary.

Stress State Change of Sequential
Tunnel Advancement e research above shows that the direction of the principal stress will change at different degrees during tunnel excavation. However, there are only a few reports on the effect of principal stress axis rotation on the stress state ahead of the tunnel face. Here, the stress triaxiality R σ and the Rod parameter μ σ are introduced, which are regularly used to feature the stress states of metal materials. In early fracture mechanics, Sih [18] presented the theory that the location with minimum strain energy density is always the location with smaller deformation energy density and larger dilation energy density where brittle rupture readily begins. Based on Sih's theory, Rice and Tracey [19] defined the expression of stress triaxiality R σ . Rice believed that the mixed-mode crack extends along the direction of maximum stress triaxiality, which agrees with the type I fracture toughness testing results. In the 1990s, Bao and Wierzbicki [20] and Bao [21] demonstrated that the value of stress triaxiality can reflect the stress state of the material in various force conditions. e stress triaxiality R σ and the Rod parameter μ σ are defined in equations (3) and (4), and Table 2 lists the values of R σ and μ σ that represent different stress states from triaxial tension to triaxial compression. R σ decreases gradually as the stress state changes from tension to compression. e material is under tension and subject to tensile failure if R σ is positive, and the degree of tension increases as R σ increases, and vice versa. However, there is no correspondent relationship between the rod parameter and stress state. erefore, the stress triaxiality R σ is relatively better for reflecting the stress state:

Advances in Materials Science and Engineering
e stress triaxiality R σ is the ratio of average normal stress I 1 /3 and the equivalent shear stress �� J 2 , where I 1 is the first stress invariant and J 2 is the second deviatoric stress invariant. Figure 8 shows the change curves of σ 1 , σ 2 , σ 3 , I 1 /3, �� J 2 , and R σ of the monitoring points as the tunnel face advances. e stress path of monitoring point A is shown generally as I 1 /3 decreases and �� J 2 increases, and the stress state of the tunnel roof changes from triaxial compression to biaxial or uniaxial compression.
is is the same for monitoring point B in the tunnel wall.
Careful observation of the 46-60 m (Figure 8(a)) or 48-60 m (Figure 8(b)) cyclical footage shows that the principal stress values of the monitoring points undergo little change and I 1 /3 is nearly unchanged at the beginning of the excavation process. e stress concentration begins to occur when the tunnel face is 4-6 m from the reference section and reaches a peak value as the tunnel face arrives at the reference section. However, during this process (46-60 m or 48-60 m cyclical footage), R σ and �� J 2 increase rapidly, as shown in Figure 8, and the principal stress axes experience dramatic rotation, as shown in Figure 6.
Obviously, the stress state influenced by the change in magnitudes of the principal stresses can be expressed using stress parameters (e.g.,I 1 /3, σ 1 , σ 2 , and σ 3 ). To determine how the stress state changes with the principal stress axis rotation, the stress increment as a result of only the principal stress axes rotation is analyzed in this paper. A two-dimensional stress state σis analyzed for simplification, as written in equation (5), where σ 1 and σ 2 denote the major and minor principal stress, respectively. It is assumed that σ 1 and σ 2 are constant and the principal stress axes rotate a certain angle θ. e new stress state σ ′ after the principal stress axes rotate is written as equation (6), and the relationship between σ and σ ′ as equation (7), where Q � cos θ −sin θ sin θ cos θ is the coordinate transformation matrix. e stress increment dσ ′ due to the principal stress axis rotation in the new principal stress space is written as equation (8), and the stress increment dσ (transformed from dσ ′ ) in the original principal stress space is written as equation (9): Equation (9) shows that the stress increments due to the principal stress axis rotation are mainly influenced by the deviatoric stresses and the rotation angle of the principal stress axes, and the shear stress component of (σ 1 − σ 2 )dθ explains that the principal stress axes rotation can result in an increase of �� J 2 . us, �� J 2 can be considered to be closely related to the stress state change of rock mass due to the principal stress axes rotation, and the reason R σ shares similar change laws with �� J 2 and with dramatic change between 46-60 m or 48-60 m cyclical footage in Figure 8 can be explained. It is clear from equation (3) that R σ denotes the stress state variation induced by the interaction of I 1 /3 and �� J 2 and can represent the stress state change due to both the variation of magnitude and orientation of the stress-field tensor.
According to the curve of R σ with the tunnel face advancing in Figure 8, the stress state change during tunnel excavation can be roughly divided into four stages. In Stage I, at 0-46 m or 0-48 m cyclical footage, R σ decreases gradually from the initial value of approximately -4.2, indicating a slight adjustment due to the tunnel excavation. In Stage II, at 46-60 m or 48-60 m cyclical footage, R σ increases quickly to approximately -1.67 mainly due to the principal stress axes rotation and the stress concentration. In Stage III, at 60-70 m cyclical footage, R σ increases to approximately -0.67 due to the stress release and the principal stress axis rotation as a result of the tunnel face passing through the reference section. en, R σ is approximately unchanged in Stage IV at 70-100 m cyclical footage. Figure 9 shows the nephogram of R σ when the tunnel face moves towards the reference section. According to the R σ curves in Figure 8 and the physical meaning of the stress state values in Table 2, some characteristic values are selected as thresholds for evaluating the degree and characteristics of stress disturbance. erefore, the stress disturbance zone of the excavated tunnel is divided into approximately two sections, as shown in Figure 9. One section is near the tunnel surface with R σ > − 1.67, similar to an even hollow cylinder and extended approximately half of a tunnel diameter (0.5D) from the tunnel boundary into the rock mass. is section suffers from quick stress release, strong stress direction    Advances in Materials Science and Engineering rotation, and the stress state variations from triaxial compression to biaxial or uniaxial compression due to the tunnel excavation. Another section is far from the tunnel surface, with −1.67 > R σ > − 4.2. e distribution range of this section is 1-4 times as large as the tunnel diameter (1-4D) along the radial direction and approximately as large as the tunnel diameter (1D) ahead of the tunnel face, and its space form may be affected by the primary in situ stress field. is zone is disturbed weakly, where the stress magnitude and stress directions after excavation have little change compared with the initial in situ condition, and the stress state changes from unequal triaxial compression to hydrostatic loading condition. e distribution of the strong stress disturbance zone nearly coincides with that of the damage zone of tunnel excavation in some literature [6,10], which indicates that the stress triaxiality can reasonably be used to describe the change of stress state and the stress disturbance zone of the excavated tunnel.

Conclusions
e study on the principal stress orientation and the stress state changes during the advancement of the tunnel face in the diversion tunnels of the Jinping II Hydropower Station presents some elementary conclusions on stress disturbance characteristics and mechanisms: (1) e pole diagram can describe the rotation of the principal stress axes as the tunnel face advances and the analysis result shows that the orientation adjustments of principal stresses in different positions near the tunnel boundary have some common characteristics. e minor and major principal stresses ahead of the tunnel face will rotate to intersect the excavation surface of the tunnel at angles ranging from 30°to 60°, with the intermediate principal stress nearly parallel to the excavation surface (approximately normal to the tunnel axis). In addition, the principal stress directions will adjust themselves dramatically at approximately 10-12 m (approximately equal to the tunnel diameter) ahead of the tunnel face, as the tunnel face moves forward.
(2) e stress triaxiality, which is commonly used to indicate the deformation and damage of metal materials, is introduced to describe the stress state changes of the primary in situ stress field. e stress triaxiality is found to represent the stress state change due to the variation of both magnitude and orientation of the stress-field tensor. According to the physical meaning and the change law of stress triaxiality, the stress disturbance during tunnel excavation can be divided into four stages. (3) In terms of the disturbance mechanism and spatial distribution of the stress triaxiality, the stress disturbance zone is divided into the strong disturbance zone and the weak disturbance zone with different disturbance characteristics. e strong disturbance zone extends approximately half the tunnel diameter from the tunnel boundary into the rock mass, which suffers from quick stress release, strong stress direction rotation, and stress state variations from triaxial compression to biaxial or uniaxial compression due to the tunnel excavation.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.