# Numerical Modeling of Stress Disturbance Characteristics during Tunnel Excavation

Authors: Zeqi Zhu; Qian Sheng; Yumin Zhang; Shiwei Liu
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

The stress state and principal stress axis changes of the stress-field tensor are analyzed during the advancement of a tunnel face on the basis of a given case study of the Jinping II Hydropower Station in China. First, the prevailing pole diagram in geology is used to illustrate the rotation of the stress axes as the tunnel face advances. The results show that the orientation adjustments of principal stresses in different positions near the tunnel boundary share common characteristics. The major and minor principal stress axes ahead of the tunnel face will rotate to intersect with the excavation surface at an angle, with the intermediate principal stress axis being almost parallel to the excavation surface. Furthermore, the stress triaxiality that is commonly used to indicate the deformation and damage of metal materials is introduced to describe the stress state change of the excavation-induced stress. The stress triaxiality is found to represent the stress state change due to the variation in both the magnitude and orientation of the stress-field tensor. According to the physical meaning and the change law of the stress triaxiality, stress disturbance during tunnel excavation can be divided into four stages, and the stress disturbance zone is divided into a strong disturbance zone and a weak disturbance zone. The disturbance characteristics of different stages and the distribution patterns of various zones are analyzed, which may be useful for practical application in the design and construction of rock tunnels.


---

## Body

## 1. 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–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. The main representative studies are the Underground Research Laboratory (URL) of Atomic Energy of Canada Limited (AECL) (1989–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.The 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 excavation-induced 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 excavation-induced stress.This 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. The 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.


## 2. Engineering Case and Numerical Model

The 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. The four tunnels are the key project simply because they are buried deep, with an average length of approximately 16.67 km, an excavation diameterD of approximately 13 m, a maximum buried depth of approximately 2525 m, and an overburden rock depth of generally 1500–2000 m. Therefore, 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.The 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. This segment is buried by 1100 m of thick overburden rock, with a gravity stress of approximately 29.2 MPa. The 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. The 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. The 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. The 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.The numerical model and theX, 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. The 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. The reference section selected is at Z = 60 m, and the schematic advancing of the tunnel face is shown in Figure 2. The 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.Figure 1
Three-dimensional finite difference mesh, model geometry, and in situ stress-field.Figure 2
Schematic advancing of the tunnel face.The tunnel size is shown in Figure3 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). The Mohr–Coulomb model is used and Table 1 provides the material properties used in the elastoplastic model.Figure 3
Tunnel geometry and layout of monitoring points in the reference section (Z = 60 m) (unit: m).Table 1
Finite difference mesh input parameters used in the three-dimensional tunnel model.
Model input parameter
Value
Density,ρ
2750 kg/m3
Young’s modulus, E
30.0 GPa
Poisson’s ratio,ν
0.20
Intact cohesion, c
2.0 MPa
Intact friction angle,ϕ
55.0∘
Tensile strength,σt
1.0 MPa


## 3. Principal Stress Axes Rotation ahead of an Advancing Tunnel Face

The typical principal stress orientation vectors of monitoring point A in the reference section are shown in Figure4 during the advancement of the tunnel face. The direction of the principal stresses experiences rotations at different degrees when the tunnel face is far from, before, on, and behind the reference section. The 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)). Then, 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.Figure 4
Principal stress orientation vectors of the tunnel roof at different advancing stages.(a) Tunnel face far from the reference section. (b) Tunnel face before the reference section. (c) Tunnel face on the reference section. (d) Tunnel face behind the reference section.
(a)
(b)
(c)
(d)However, illustrating the rotation process of the principal stress axes during tunnel excavation in a precise and intuitive way using the direction vectors in Figure4 is difficult. Thus, 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. The angles between the principal stress axis and X, Y, and Z are denoted by α, β, and γ, respectively. In terms of the spatial relation, η and ξ can be expressed by α, β, and γ as follows:(1)η=arctancosγcosα,if cosα>0,η=90∘,if cosα=0, cosγ>0,η=270∘,if cosα=0, cosγ<0,η=arctancosγcosα+180∘,if cosα<0,(2)ξ=π2−β.Figure 5
The principal stress axis represented in geographic coordinates and computational coordinates.In the finite difference numerical calculation, three direction cosinescosα, cosβ, and cosγ are used to describe the direction of principal stress. Through equations (1) and (2), the direction of principal stress can be described in the two-dimensional 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.The relationship between geographic coordinates and computational coordinates in Figure5 shows that the diversion tunnel axis extends from west to east. The 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°. The 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. The inclination of the σ2stress axis is basically normal to the tunnel axis, with ξ equal to approximately 0°.Figure 6
Pole diagrams of principal stress axes of monitoring points as the tunnel face approaches from WF = 0 to WF = 60 m. (a) Monitoring point A of the tunnel roof. (b) Monitoring point B of the tunnel wall.
(a)
(b)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. The 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. The 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 Figure7, with schematic diagrams of the principal stress orientation and the tunnel boundary.Figure 7
Schematic diagram of principal stresses orientation of monitoring points at WF = 60 m.


## 4. Stress State Change of Sequential Tunnel Advancement

The 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 triaxialityRσ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.The stress triaxialityRσ 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. The 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. Therefore, the stress triaxiality Rσ is relatively better for reflecting the stress state:(3)Rσ=2σ1+σ2+σ33σ1−σ22+σ2−σ32+σ3−σ12,(4)μσ=2σ2−σ1−σ3σ1−σ3.Table 2
R
σ and μσ under different stress states.
Stress state
Principal stress
R
σ
μ
σ
σ
1
σ
2
σ
3
Unequal triaxial tension
σ
σ
0.5σ
1.67
1
Biaxial tension
σ
σ
0
0.67
1
Simple tension
σ
0
0
0.33
−1
Pure shear
σ
0
−σ
0
0
Uniaxial compression
0
0
−σ
−0.33
1
Biaxial compression
0
−σ
−σ
−0.67
−1
Unequal triaxial compression
−0.5σ
−σ
−σ
−1.67
−1The stress triaxialityRσ is the ratio of average normal stress I1/3 and the equivalent shear stress J2, where I1 is the first stress invariant and J2 is the second deviatoric stress invariant. Figure 8 shows the change curves of σ1, σ2, σ3, I1/3, J2, and Rσ of the monitoring points as the tunnel face advances. The stress path of monitoring point A is shown generally as I1/3 decreases and J2 increases, and the stress state of the tunnel roof changes from triaxial compression to biaxial or uniaxial compression. This is the same for monitoring point B in the tunnel wall.Figure 8
The curves ofσ1, σ2, σ3, I1/3, J2, and Rσ with tunnel face advancing (Stage I: stress adjustment; Stage II: principal stress axis rotation and stress concentration; Stage III: stress release and principal stress axis rotation; Stage IV: stress redistribution).
(a)
(b)Careful observation of the 46–60 m (Figure8(a)) or 48–60 m (Figure 8(b)) cyclical footage shows that the principal stress values of the monitoring points undergo little change and I1/3 is nearly unchanged at the beginning of the excavation process. The 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 J2increase 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.,I1/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 σ2are constant and the principal stress axes rotate a certain angle θ. The 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. The 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):(5)σ=σ100σ2,(6)σ′=σ1′00σ2′,(7)σ′=QσQT,(8)dσ′=dQσQT=dQσQT+QdσQT+QσdQT,(9)dσ=QTdσ′Q=0σ1−σ2dθσ1−σ2dθ0.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−σ2dθ explains that the principal stress axes rotation can result in an increase of J2. Thus, J2can 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 J2 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 I1/3 and J2 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 ofRσ 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. Then, Rσ is approximately unchanged in Stage IV at 70–100 m cyclical footage.Figure9 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. Therefore, 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. This section suffers from quick stress release, strong stress direction 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. The 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. This 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. The 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.Figure 9
Nephogram of stress triaxiality as the tunnel face approaches the reference section.


## 5. Conclusions

The 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)
The 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. The 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)
The 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. The 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. The 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.



---

*Source: 4508905-2020-04-28.xml*