# A Slate Tunnel Stability Analysis considering the Influence of Anisotropic Bedding Properties

Authors: Jun Cai; Guangyin Du; Haiwang Ye; Tao Lei; Han Xia; Huangsong Pan
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

Bedding planes are the fundamental causes of anisotropic deformation and mechanical behaviors in slate, which will have great influence on tunnel stability. In order to analyse tunnel stability surrounded by slate, well-foliated slate in eastern Guizhou was taken as the specimen in tests. Microscopic analysis and test results show that slate can be regarded as a special continuous material. During the test, shear strength parameters and progressive failure varied when the direction of the bedding plane was changed, and two sets of reasonable shear strength were achieved by fitting. Numerical model verification is conducted before applying, and results indicate that the model can represent the anisotropic failure properties. So the model considering anisotropic shear strength simultaneously is utilized to analyse the tunnel stability in slate. When it is medium dip angle, the tunnel is significantly unstable especially for face and side walls, and at 45° (dip angle), the plastic zone depth ahead of the tunnel face can be the largest, being 1.7 times the tunnel height. The maximum deviator stress (σ1 − σ3) is centralized on the middle of the side wall, and also, the stress (σ1 − σ3) is the highest at 45° (dip angle), which will lead to shear failure.


---

## Body

## 1. Introduction

Layered rock masses with sedimentary structure can be found in many parts of the world, and they are often encountered in all types of underground engineering construction [1, 2]. During the process of deposition, layered rock masses possess different structures and coupling characteristics in different directions. One of the typical mechanical properties of a layered structure rock mass is anisotropy, which results from mineral particles being of different sizes and in compound mode [3]. These structures in the rock mass organized in alignment generate the bedding plane and foliation features. Because of the poor mechanical properties of bedding planes, the strength, failure patterns, and deformations are different in each direction [4], which becomes a great challenge for the construction of civil engineering. In order to understand this better, Tien et al. [5] discussed the failure of isotropic rocks to make the failure mechanism around a tunnel clear, which gives valuable insight into the collapse mode. Wang et al. [6] studied the failure mechanism for a circular tunnel in transversely isotropic rock by using RFPA, where the failure process is characterized by the initiation, propagation, and coalescence of cracks around the tunnel. Xu et al. [7] developed a transversely isotropic elastic-plastic model to describe the elastic response and post–peak failure behavior based on the Mohr–Coulomb and maximum tensile-stress criteria, which was utilized to evaluate rock mass excavation response in underground openings. Based on the evaluation of large scale of randomly distributed discontinuities, Yang et al. [8] studied the influence of joint plane angle.Numerical methods have been extensively utilized to give valuable insights into the evaluation of tunnel stability in many publications [9–12], and some calculation model achievements [13, 14] have been made. Besides a lot of fundamental anisotropic mechanical property, work about the layered rock has been made to make tunnel stability analysis proper. By performing a series of uniaxial and triaxial compression tests on the slate specimens with varying inclination angles, Chen et al. [15] examined the anisotropic behaviors and failure mechanisms of slate, based on which a micromechanical damage-friction model was proposed. Niandou et al. [16] probed into the failure pattern of Tournemire shale; together with some studies [17, 18], they have shown that the strength and failure pattern of layered rock are highly dependent on the bedding plane dip angle, θ.Even so many complicated excavation cases in anisotropic rock have been studied, but few work considered several sets of anisotropic shear strength simultaneously in one model. So in the current study, in order to analyse the tunnel stability in layered rock mass, a specific numerical model considering anisotropic shear strength that reflected in experiment tests was utilized. Before analysing tunnel stability, numerical model validation was conducted between experiment and numerical results. Then, shear loading process and the anisotropic failure properties shown in experiments were reproduced numerically. So the validated model was utilized to evaluate the stress and deformation of an anisotropic rock tunnel.


## 2. Shear Test

### 2.1. Specimen Preparation

A block of slate was collected from a slate quarry in Tongren city, eastern Guizhou Province, China. The slate is well developed in a laminated structure, and together with its sediment structure, the specimen was observed to be composed of macro- and microlayers (Figure1). Figure 1(b) shows a thin photomicrograph of the slate, and the foliation planes can be viewed. The array of layers is seen to be arranged directionally when magnified 100 times, and metamorphic minerals and shale particles between the layers are well cemented. The microscopic analysis indicates that the slate is mostly composed of calcite, quartz, and mica, among which the calcite is varied in size, and the fine grain size of the quartz and mica is about 25 μm.Figure 1
Macro- and microimages of slate. (a) Photomicrographs of bedding plane structure. (b). Thin section of the slate.
(a)
(b)


### 2.2. Test Program

The shear strength of the specimens can be tested by various methods, among which the confined shear test is extensively used in practice. Cylindrical or cubic specimens are usually used with separate corresponding loading methods.In the shear test, cylindrical specimens are mostly used with the device in Figure2(a), the normal and shear stress can be achieved by loading device independently. For cubic specimens, two components of normal and shear stress were achieved by changing the shear angle (α in Figure 2(b)); the specimen is placed between two beveled dies and tested along a predetermined shear plane. For the cylindrical device (Figure 2(a)), the progressive failure cannot be observed, and the failure pattern differences are not so significant when changing loading mode. Given all this, loading mode in Figure 2(b) was selected.Figure 2
Shear test loading modes comparison. (a). Shear test with cylindrical specimens. (b). Shear test with cube specimens.
(a)
(b)The cubic specimens (100 mm × 100 mm × 100 mm) were prepared (Figure3). Error of the length of the specimens is within ±0.1 mm, and the parallelism of the specimen ends is within ±0.5 mm after being polished. Before the test, the specimens were all stored in dry conditions at room temperature.Figure 3
Cube specimens of slate for shear tests.The strength of a rock mass not only is related to the rock properties but also can be affected by the properties of discontinuities and interfaces, and the properties can vary when the direction of discontinuities is changed [19–21]. A total of 54 cubic specimens, 100 mm in length, were classified into three groups according to the spatial location relationship (Figure 4). First group (Figure 4(a)) is the condition that the bedding plane is parallel to the shear surface, the second (Figure 4(b)) is that the bedding plane is perpendicular to the shear surface and the intersection line of these two planes is perpendicular to the shear stress, and the third (Figure 4(c)) is that the bedding plane is perpendicular to the shear surface and the intersection line of these two planes is parallel to the shear stress. Each of these three loading mode classifications (Figures. 4(a)–4(c)) were divided into different groups by the adjustment of the range of shear angle (α in Figure 2) as 20°, 30°, 40°, 50°, 60°, and 70°. By changing the shear angle (α), the test was devised to obtain different σ and τ. By changing the loading mode classifications (Figures. 4(a)–4(c)), the anisotropic shear strength and progressive failure characteristics of the slate can be achieved. The final results were optimized by choosing the intermediate in each group and eliminating the maximum deviation.Figure 4
Loading modes for shear tests.


## 2.1. Specimen Preparation

A block of slate was collected from a slate quarry in Tongren city, eastern Guizhou Province, China. The slate is well developed in a laminated structure, and together with its sediment structure, the specimen was observed to be composed of macro- and microlayers (Figure1). Figure 1(b) shows a thin photomicrograph of the slate, and the foliation planes can be viewed. The array of layers is seen to be arranged directionally when magnified 100 times, and metamorphic minerals and shale particles between the layers are well cemented. The microscopic analysis indicates that the slate is mostly composed of calcite, quartz, and mica, among which the calcite is varied in size, and the fine grain size of the quartz and mica is about 25 μm.Figure 1
Macro- and microimages of slate. (a) Photomicrographs of bedding plane structure. (b). Thin section of the slate.
(a)
(b)


## 2.2. Test Program

The shear strength of the specimens can be tested by various methods, among which the confined shear test is extensively used in practice. Cylindrical or cubic specimens are usually used with separate corresponding loading methods.In the shear test, cylindrical specimens are mostly used with the device in Figure2(a), the normal and shear stress can be achieved by loading device independently. For cubic specimens, two components of normal and shear stress were achieved by changing the shear angle (α in Figure 2(b)); the specimen is placed between two beveled dies and tested along a predetermined shear plane. For the cylindrical device (Figure 2(a)), the progressive failure cannot be observed, and the failure pattern differences are not so significant when changing loading mode. Given all this, loading mode in Figure 2(b) was selected.Figure 2
Shear test loading modes comparison. (a). Shear test with cylindrical specimens. (b). Shear test with cube specimens.
(a)
(b)The cubic specimens (100 mm × 100 mm × 100 mm) were prepared (Figure3). Error of the length of the specimens is within ±0.1 mm, and the parallelism of the specimen ends is within ±0.5 mm after being polished. Before the test, the specimens were all stored in dry conditions at room temperature.Figure 3
Cube specimens of slate for shear tests.The strength of a rock mass not only is related to the rock properties but also can be affected by the properties of discontinuities and interfaces, and the properties can vary when the direction of discontinuities is changed [19–21]. A total of 54 cubic specimens, 100 mm in length, were classified into three groups according to the spatial location relationship (Figure 4). First group (Figure 4(a)) is the condition that the bedding plane is parallel to the shear surface, the second (Figure 4(b)) is that the bedding plane is perpendicular to the shear surface and the intersection line of these two planes is perpendicular to the shear stress, and the third (Figure 4(c)) is that the bedding plane is perpendicular to the shear surface and the intersection line of these two planes is parallel to the shear stress. Each of these three loading mode classifications (Figures. 4(a)–4(c)) were divided into different groups by the adjustment of the range of shear angle (α in Figure 2) as 20°, 30°, 40°, 50°, 60°, and 70°. By changing the shear angle (α), the test was devised to obtain different σ and τ. By changing the loading mode classifications (Figures. 4(a)–4(c)), the anisotropic shear strength and progressive failure characteristics of the slate can be achieved. The final results were optimized by choosing the intermediate in each group and eliminating the maximum deviation.Figure 4
Loading modes for shear tests.


## 3. Anisotropic Properties Analysis

### 3.1. Experimental Results

The digital ultrasonic instrument (RSM-SY5) used for acoustic testing was in Figure5, and results of the acoustic wave test (Table 1) indicate that Vp¯ (average velocity of P-wave along the bedding plane) is 6058.4 m/s, and Vv¯ (average velocity of P-wave across the bedding plane) is 5274.7 m/s, Vv¯ is 12.9% less than Vp¯. The difference of the P-wave velocity clearly shows the anisotropic behavior of slate [15], and the reduction in velocity across the bedding plane means that the microcracks in discontinuities caused the diffraction, and these mostly exist in the bedding plane. Physical test results show that the average dry density is 2771 kg/m3, the elasticity modulus is 78.231 GPa, and the Poisson’s ratio is 0.274. Table 1 is the result of the shear test, and Figure 6 is the shear stress vs shear displacement, and loading modes (a), (b), (c) correspond to the three loading methods in Figure 4. During the loading process, the curve is an upward parabola. With the shear angle (α) increases, shear displacement increases in general except for condition a. For condition a, shear failure formed along bedding, and the test ended quickly; so shear displacement difference in condition a can be unapparent.Figure 5
Acoustic wave test.Table 1
Results of shear tests for slate.
Loading mode
Specimen
V
p (m/s)
V
v (m/s)
τ (MPa)
σ (MPa)
Loading mode
Shear angle (°)
(a)
20
L26−1
6208.485
5382.703
17.097
46.964
30
L4−2
5863.137
4881.301
12.066
20.899
40
L12−1
5976.078
4690.732
9.809
11.691
50
L4−1
6196.250
5425.926
10.668
8.952
60
L20−2
5961.176
5394.035
7.511
4.337
70
L18−1
6161.290
5405.614
7.468
2.719
(b)
20
L18−2
6196.250
5744.762
19.784
54.355
30
L24−1
6312.500
5387.368
14.380
24.908
40
L12−4
6052.157
5118.596
14.287
17.026
50
L32−1
5872.941
5485.926
13.408
11.250
60
L16−1
5995.686
5246.316
11.708
6.760
70
L17−2
6187.917
5324.675
11.362
4.135
(c)
20
L16−3
5842.745
5298.119
21.211
58.276
30
L31−2
6051.717
5227.719
33.621
58.234
40
L12−2
5998.462
5109.769
31.839
37.944
50
L20−1
6272.917
5422.222
36.694
30.790
60
L32−2
6011.667
5144.333
58.362
33.695
70
L32−3
5890.221
5254.019
55.212
20.096Figure 6
Shear stress vs displacement curves under loading modes: (a), (b), and (c).
(a)
(b)
(c)


### 3.2. Failure Pattern Analysis

The failure patterns (Figure7) clearly show the anisotropic behavior of slate. In condition a, when the stress peaked (Figure 7(a)), with the influence of the weak bedding plane, the failure formed along the direction of the bedding plane in the initial phase. Then, the test ended quickly without much bulk solids and fragments left on the test platform, leaving the shear plane of the specimen very smooth. In condition b, failure was on the two sides of the shear plane. As shear moving occurred between the upper and lower parts of the specimen, tensile failure appeared along the bedding plane.Figure 7
Failure patterns of slate with different loading modes in shear tests.
(a)
(b)
(c)In condition c, failure appeared along the bedding direction on the loading plane (Figure7(c)). With the influence of normal stress and shear stress, tensile strain formed in the direction perpendicular to the bedding plane, and tensile failure propagated along the bedding direction.Obviously, the cause of initial failure in condition c is different from the other two. Condition c is the transverse tensile failure under compression (Figure8), whereas condition a and condition b are the shear moving failure. Figure 8 shows the results of condition c and the uniaxial compression test. In uniaxial compression test, the bedding direction is parallel to the loading direction. The comparison illustrates that the failure pattern are the same, and they are all the transverse tensile failure along the bedding plane.Figure 8
Failure pattern comparison.
(a)
(b)
(c)The peak stress of the test results is shown in theτ−σ coordinate system, and Figures 9(a)–9(c) are the generalized sketch graph. Mohr–Coulomb equation is as follows:(1)τ=c+σ⋅tanφ,where c is the rock cohesion and φ is the rock internal friction angle. So in condition a, φ is 11.96° and c is 7.394 MPa, and in condition b, φ is 9.18° and c is 10.965 MPa. Reasonable small friction angle and cohesion were also achieved by Mao [22].Figure 9
Failure patterns and shear test results with fitting Mohr−Coulomb curves under loading modes: (a), (b) and (c).But for condition c, a reasonable Mohr–Coulomb curve cannot be fitted. During the shear loading, it is compression failure, not a shear failure [23]. So, conditions a, b, and c are all controlled by the bedding plane, but differently, in condition c, the transverse tensile strain under the compression leads to the damage of specimen.


### 3.3. Test Damage Failure Evolution Process

For anisotropic failure pattern, the crack propagation process is observed during the test. In condition a, specimen damaged rapidly, and failure evolution process is not obvious. Figure10 shows the failure evolution of condition b and c, and the observed zone is plotted.Figure 10
Failure evolution of loading modes: (b) and (c). (a). Failure evolution of condition b. (b). Failure evolution of condition c.
(a)
(b)For condition b, in the initial phase, a tensile crack formed on the two sides of the specimen. The crack across the matrix body propagated from the right side to the middle, and the left crack propagated along the bedding plane. The process is a mixed-mode crack propagation. When crack at the bottom reached the middle of the specimen, it arrested, as the upper half part was confined. Finally, a penetrating crack formed across the specimen along the shear plane, and the test ended.In condition c, initial microcracks formed on the shear stress loading plane, and the cracks propagated along the bedding plane. When cracks came to the middle of the specimen, influenced by preexisting microcracks, the cracks propagated to the bedding plane nearby. Then, more and more cracks developed on both side edges of the specimen.It can be concluded that weakened bedding planes are the fundamental causes of anisotropic mechanical behaviors of slate. The failure patterns are determined by anisotropic deformation in the bedding plane, and it is also illustrated by Zhou et al. [24].


## 3.1. Experimental Results

The digital ultrasonic instrument (RSM-SY5) used for acoustic testing was in Figure5, and results of the acoustic wave test (Table 1) indicate that Vp¯ (average velocity of P-wave along the bedding plane) is 6058.4 m/s, and Vv¯ (average velocity of P-wave across the bedding plane) is 5274.7 m/s, Vv¯ is 12.9% less than Vp¯. The difference of the P-wave velocity clearly shows the anisotropic behavior of slate [15], and the reduction in velocity across the bedding plane means that the microcracks in discontinuities caused the diffraction, and these mostly exist in the bedding plane. Physical test results show that the average dry density is 2771 kg/m3, the elasticity modulus is 78.231 GPa, and the Poisson’s ratio is 0.274. Table 1 is the result of the shear test, and Figure 6 is the shear stress vs shear displacement, and loading modes (a), (b), (c) correspond to the three loading methods in Figure 4. During the loading process, the curve is an upward parabola. With the shear angle (α) increases, shear displacement increases in general except for condition a. For condition a, shear failure formed along bedding, and the test ended quickly; so shear displacement difference in condition a can be unapparent.Figure 5
Acoustic wave test.Table 1
Results of shear tests for slate.
Loading mode
Specimen
V
p (m/s)
V
v (m/s)
τ (MPa)
σ (MPa)
Loading mode
Shear angle (°)
(a)
20
L26−1
6208.485
5382.703
17.097
46.964
30
L4−2
5863.137
4881.301
12.066
20.899
40
L12−1
5976.078
4690.732
9.809
11.691
50
L4−1
6196.250
5425.926
10.668
8.952
60
L20−2
5961.176
5394.035
7.511
4.337
70
L18−1
6161.290
5405.614
7.468
2.719
(b)
20
L18−2
6196.250
5744.762
19.784
54.355
30
L24−1
6312.500
5387.368
14.380
24.908
40
L12−4
6052.157
5118.596
14.287
17.026
50
L32−1
5872.941
5485.926
13.408
11.250
60
L16−1
5995.686
5246.316
11.708
6.760
70
L17−2
6187.917
5324.675
11.362
4.135
(c)
20
L16−3
5842.745
5298.119
21.211
58.276
30
L31−2
6051.717
5227.719
33.621
58.234
40
L12−2
5998.462
5109.769
31.839
37.944
50
L20−1
6272.917
5422.222
36.694
30.790
60
L32−2
6011.667
5144.333
58.362
33.695
70
L32−3
5890.221
5254.019
55.212
20.096Figure 6
Shear stress vs displacement curves under loading modes: (a), (b), and (c).
(a)
(b)
(c)


## 3.2. Failure Pattern Analysis

The failure patterns (Figure7) clearly show the anisotropic behavior of slate. In condition a, when the stress peaked (Figure 7(a)), with the influence of the weak bedding plane, the failure formed along the direction of the bedding plane in the initial phase. Then, the test ended quickly without much bulk solids and fragments left on the test platform, leaving the shear plane of the specimen very smooth. In condition b, failure was on the two sides of the shear plane. As shear moving occurred between the upper and lower parts of the specimen, tensile failure appeared along the bedding plane.Figure 7
Failure patterns of slate with different loading modes in shear tests.
(a)
(b)
(c)In condition c, failure appeared along the bedding direction on the loading plane (Figure7(c)). With the influence of normal stress and shear stress, tensile strain formed in the direction perpendicular to the bedding plane, and tensile failure propagated along the bedding direction.Obviously, the cause of initial failure in condition c is different from the other two. Condition c is the transverse tensile failure under compression (Figure8), whereas condition a and condition b are the shear moving failure. Figure 8 shows the results of condition c and the uniaxial compression test. In uniaxial compression test, the bedding direction is parallel to the loading direction. The comparison illustrates that the failure pattern are the same, and they are all the transverse tensile failure along the bedding plane.Figure 8
Failure pattern comparison.
(a)
(b)
(c)The peak stress of the test results is shown in theτ−σ coordinate system, and Figures 9(a)–9(c) are the generalized sketch graph. Mohr–Coulomb equation is as follows:(1)τ=c+σ⋅tanφ,where c is the rock cohesion and φ is the rock internal friction angle. So in condition a, φ is 11.96° and c is 7.394 MPa, and in condition b, φ is 9.18° and c is 10.965 MPa. Reasonable small friction angle and cohesion were also achieved by Mao [22].Figure 9
Failure patterns and shear test results with fitting Mohr−Coulomb curves under loading modes: (a), (b) and (c).But for condition c, a reasonable Mohr–Coulomb curve cannot be fitted. During the shear loading, it is compression failure, not a shear failure [23]. So, conditions a, b, and c are all controlled by the bedding plane, but differently, in condition c, the transverse tensile strain under the compression leads to the damage of specimen.


## 3.3. Test Damage Failure Evolution Process

For anisotropic failure pattern, the crack propagation process is observed during the test. In condition a, specimen damaged rapidly, and failure evolution process is not obvious. Figure10 shows the failure evolution of condition b and c, and the observed zone is plotted.Figure 10
Failure evolution of loading modes: (b) and (c). (a). Failure evolution of condition b. (b). Failure evolution of condition c.
(a)
(b)For condition b, in the initial phase, a tensile crack formed on the two sides of the specimen. The crack across the matrix body propagated from the right side to the middle, and the left crack propagated along the bedding plane. The process is a mixed-mode crack propagation. When crack at the bottom reached the middle of the specimen, it arrested, as the upper half part was confined. Finally, a penetrating crack formed across the specimen along the shear plane, and the test ended.In condition c, initial microcracks formed on the shear stress loading plane, and the cracks propagated along the bedding plane. When cracks came to the middle of the specimen, influenced by preexisting microcracks, the cracks propagated to the bedding plane nearby. Then, more and more cracks developed on both side edges of the specimen.It can be concluded that weakened bedding planes are the fundamental causes of anisotropic mechanical behaviors of slate. The failure patterns are determined by anisotropic deformation in the bedding plane, and it is also illustrated by Zhou et al. [24].


## 4. Modeling

### 4.1. Modeling of Slate with Bedding

Numerical methods can also be employed to evaluate the strength of anisotropically jointed media, and modeling can generally be accomplished in two ways: discontinuous and equivalent continuum approaches. Latha and Garaga [25] demonstrated that both approaches can appropriately model the conditions of triaxial compression tests on anisotropically jointed specimens.Slate in eastern Guizhou contains a high density of parallel bedding planes, and the spacing between the bedding planes can be assumed to be sufficiently close; so slate can be regarded as a special continuous material, and these agree with assumption of the model. The model is the achievement of Zienkiewicz and Pande in 1977 [26], and it has been embedded in ABAQUS [27].


#### 4.1.1. Stress and Strain of Bedding Structure

The stress state in the rock mass is assumed to beσ=σxτxyτxzτyxσyτyzτzxτzyσz. According to the research of Zienkiewicz and Pande in 1977 [26], a particular bedding plane is considered by the normal vector na, where na=l,m,nT, and l, m, n is the cosine of the angle between the normal vector and the three axes. The vector in the plane is defined as two units, taλ and λ=1,2, and similarly, taλ=laλ,maλ,naλT. So taλ and na are orthogonal vectors in the bedding plane, and they form the local coordinate system, as shown in Figure 11; dip angle (θ) is the angle between bedding plane and horizontal plane. So the normal pressure and shear stress of the bedding plane are given as(2)pa=na⋅σ⋅na,τaλ=na⋅σ⋅taλ.Figure 11
Coordinate systems used for transformation law of a layered rock mass. (a) Vectors of bedding plane. (b) Coordinate systems.
(a)
(b)The shear stress magnitude is defined as(3)τa=τaλ⋅τaλ.The local strain components are the normal strain across the bedding plane:(4)εan=na⋅ε⋅na,and the shear strain in direction λ (λ=1,2) of the bedding plane is(5)γaλ=na⋅ε⋅taλ+taλ⋅ε⋅na,where ε is the strain tensor.The linear strain rate in the bedding plane system can be divided into elastic and plastic parts:(6)dε=dεel+dεpl,where dε is the total strain rate, dεel is the elastic strain rate, and dεpl is the plastic strain rate. According to Hooke’s law,(7)εe=D−1σ,where D is the stiffness matrix of the elastic constant, and the plastic strain rate is linear:(8)dεpl=∑dεpl.


#### 4.1.2. Bedding Plane Opening

According to the ABAQUS documentation [27], when the pressure stress across the bedding plane becomes positive (tensile), exceeding the capacity of the bedding plane, the bedding plane is opened. Thus, opening the bedding plane creates an anisotropic elastic response at a point, and the joint system remains as(9)εanpsel≤εanel,where εanel is the component of direct elastic strain across the bedding plane and εanpsel is the component of hypothetical direct elastic strain across the bedding plane that is calculated as(10)εanpsel=−νEσa1+σa2,where E is the Young’s modulus of the material, ν is the Poisson’s ratio, and σa1,σa2 are the direct stresses in the bedding plane.


#### 4.1.3. Bedding Plane Sliding

Normally, shear failure forms in the bedding, so the sliding failure surface is defined by(11)fa=τa−σatanφa−ca=0,where φa is the friction angle of bedding system a, ca is the cohesion for system a, τa is the shear stress on the bedding plane, and σa is the pressure stress normal to the bedding plane; so it is the transformation of Mohr−Coulomb law.As long asfa<0, the bedding plane in system a does not slip. When fa=0, the bedding plane in system a slips. The plastic strain on the system is then given by(12)dγaλpl=dε¯aplτaλτacosψa,dεanpl=dε¯aplsinψa,where dγaλpl is the plastic shear strain rate in direction λ (λ=1,2 are the orthogonal directions) on the bedding plane, dε¯apl is the magnitude of the plastic strain rate, τaλ is the component of shear stress on the bedding plane, ψa is the dilation angle for the bedding system, and dεanpl is the plastic strain rate normal to the bedding plane.


### 4.2. Comparison between Experiment and Simulation

#### 4.2.1. Modeling

The size of model was 100 mm × 100 mm × 100 mm; pressure was applied on the upside and upper half side plane, respectively; the bottom and opposite lower half side plane is set with fixed boundary conditions (Figure12), and three-dimensional solid finite elements C3D8 were chosen as the model element algorithm.Figure 12
Three-dimensional numerical model mesh shape and boundary conditions.In order to make the shear effect obvious,σ and τ was the result of α = 70° under condition c. Pressure was applied on the upside and upper half plane, making σ = 20.096 MPa, τ = 55.212 MPa, bottom and lower half plane was fixed, as shown in Figure 12. According to Xie et al. [28], for slate, the reasonable friction angle can be 30–50°, and reasonable cohesion can be 2–20 MPa. In the simulation, cohesion (c = 7.394 MPa) and friction angle (φ=11.96°) were selected from the approximate result of condition a; even the friction angle may be small (φ=11.96°), and small friction angle is also acquired by Mao [22]. Table 2 lists the slate mechanical parameters for model calculation. Definition of the bedding plane direction relies on the three vectors of the plane (Figure 11), so simulation of the three conditions (Figure 4) in the shear tests could achieved.Table 2
Physicomechanical parameters of slate.
Material
Dry densityρd (kg/m3)
Poisson’s ratioυ
Elastic modulusEs (GPa)
Cohesionc (MPa)
Friction angleφ (°)
Slate
2771
0.27
78.23
7.394
11.96


#### 4.2.2. Simulation Results

Figure13 presents the simulation results, and the models in each group from left to right correspond to (a), (b), and (c) in Figure 4. In condition (a), the model computes 22 increment steps with nonconvergence; for condition (b), it is 24 increment steps with nonconvergence, and for condition (c), the computation is converged at 39 increment steps. In Figure 13, the shear stress are all centralized on the side edge of the shear plane, so this area are all prone to shear failure, which corresponds well to the final shear failure that occurs in (a), (b), and (c) in Figure 7.Figure 13
Results of numerical simulation. (a) Shear stress. (b) Logarithmic strain. (c) Logarithmic strain vector.
(a)
(b)
(c)Figures13(b) and 13(c) are the strain and its vector; (a) shows the strain concentrated on the side plane in the middle, and the orientation of the shear stress is in the middle; it fits well with the test results of condition a. Similarly, for (b), tensile strain occurs on the upper side plane, and it may results in the initial cracks, as shown in Figure 10. For (c), the strain is concentrated on the plane of Y0Z, and the direction of the strain vector is normal to the side planes, so it agrees with the initial failure mechanism in the tests (Figure 10); the bedding plane preferentially opens on both the front and back side planes.In Figure14, the curve is the point in the middle of the upper half loading plane (Figure 12). Even for condition (c), the maximum is different, and the value of the drop in the middle of the test curve is approximately the same. The good agreement between simulation and test shows that the model satisfactorily describes the anisotropic deformational and strength behaviors of slate.Figure 14
Comparison between test and simulation (α = 70°) under loading modes: (a), (b), and (c).


## 4.1. Modeling of Slate with Bedding

Numerical methods can also be employed to evaluate the strength of anisotropically jointed media, and modeling can generally be accomplished in two ways: discontinuous and equivalent continuum approaches. Latha and Garaga [25] demonstrated that both approaches can appropriately model the conditions of triaxial compression tests on anisotropically jointed specimens.Slate in eastern Guizhou contains a high density of parallel bedding planes, and the spacing between the bedding planes can be assumed to be sufficiently close; so slate can be regarded as a special continuous material, and these agree with assumption of the model. The model is the achievement of Zienkiewicz and Pande in 1977 [26], and it has been embedded in ABAQUS [27].


### 4.1.1. Stress and Strain of Bedding Structure

The stress state in the rock mass is assumed to beσ=σxτxyτxzτyxσyτyzτzxτzyσz. According to the research of Zienkiewicz and Pande in 1977 [26], a particular bedding plane is considered by the normal vector na, where na=l,m,nT, and l, m, n is the cosine of the angle between the normal vector and the three axes. The vector in the plane is defined as two units, taλ and λ=1,2, and similarly, taλ=laλ,maλ,naλT. So taλ and na are orthogonal vectors in the bedding plane, and they form the local coordinate system, as shown in Figure 11; dip angle (θ) is the angle between bedding plane and horizontal plane. So the normal pressure and shear stress of the bedding plane are given as(2)pa=na⋅σ⋅na,τaλ=na⋅σ⋅taλ.Figure 11
Coordinate systems used for transformation law of a layered rock mass. (a) Vectors of bedding plane. (b) Coordinate systems.
(a)
(b)The shear stress magnitude is defined as(3)τa=τaλ⋅τaλ.The local strain components are the normal strain across the bedding plane:(4)εan=na⋅ε⋅na,and the shear strain in direction λ (λ=1,2) of the bedding plane is(5)γaλ=na⋅ε⋅taλ+taλ⋅ε⋅na,where ε is the strain tensor.The linear strain rate in the bedding plane system can be divided into elastic and plastic parts:(6)dε=dεel+dεpl,where dε is the total strain rate, dεel is the elastic strain rate, and dεpl is the plastic strain rate. According to Hooke’s law,(7)εe=D−1σ,where D is the stiffness matrix of the elastic constant, and the plastic strain rate is linear:(8)dεpl=∑dεpl.


### 4.1.2. Bedding Plane Opening

According to the ABAQUS documentation [27], when the pressure stress across the bedding plane becomes positive (tensile), exceeding the capacity of the bedding plane, the bedding plane is opened. Thus, opening the bedding plane creates an anisotropic elastic response at a point, and the joint system remains as(9)εanpsel≤εanel,where εanel is the component of direct elastic strain across the bedding plane and εanpsel is the component of hypothetical direct elastic strain across the bedding plane that is calculated as(10)εanpsel=−νEσa1+σa2,where E is the Young’s modulus of the material, ν is the Poisson’s ratio, and σa1,σa2 are the direct stresses in the bedding plane.


### 4.1.3. Bedding Plane Sliding

Normally, shear failure forms in the bedding, so the sliding failure surface is defined by(11)fa=τa−σatanφa−ca=0,where φa is the friction angle of bedding system a, ca is the cohesion for system a, τa is the shear stress on the bedding plane, and σa is the pressure stress normal to the bedding plane; so it is the transformation of Mohr−Coulomb law.As long asfa<0, the bedding plane in system a does not slip. When fa=0, the bedding plane in system a slips. The plastic strain on the system is then given by(12)dγaλpl=dε¯aplτaλτacosψa,dεanpl=dε¯aplsinψa,where dγaλpl is the plastic shear strain rate in direction λ (λ=1,2 are the orthogonal directions) on the bedding plane, dε¯apl is the magnitude of the plastic strain rate, τaλ is the component of shear stress on the bedding plane, ψa is the dilation angle for the bedding system, and dεanpl is the plastic strain rate normal to the bedding plane.


## 4.1.1. Stress and Strain of Bedding Structure

The stress state in the rock mass is assumed to beσ=σxτxyτxzτyxσyτyzτzxτzyσz. According to the research of Zienkiewicz and Pande in 1977 [26], a particular bedding plane is considered by the normal vector na, where na=l,m,nT, and l, m, n is the cosine of the angle between the normal vector and the three axes. The vector in the plane is defined as two units, taλ and λ=1,2, and similarly, taλ=laλ,maλ,naλT. So taλ and na are orthogonal vectors in the bedding plane, and they form the local coordinate system, as shown in Figure 11; dip angle (θ) is the angle between bedding plane and horizontal plane. So the normal pressure and shear stress of the bedding plane are given as(2)pa=na⋅σ⋅na,τaλ=na⋅σ⋅taλ.Figure 11
Coordinate systems used for transformation law of a layered rock mass. (a) Vectors of bedding plane. (b) Coordinate systems.
(a)
(b)The shear stress magnitude is defined as(3)τa=τaλ⋅τaλ.The local strain components are the normal strain across the bedding plane:(4)εan=na⋅ε⋅na,and the shear strain in direction λ (λ=1,2) of the bedding plane is(5)γaλ=na⋅ε⋅taλ+taλ⋅ε⋅na,where ε is the strain tensor.The linear strain rate in the bedding plane system can be divided into elastic and plastic parts:(6)dε=dεel+dεpl,where dε is the total strain rate, dεel is the elastic strain rate, and dεpl is the plastic strain rate. According to Hooke’s law,(7)εe=D−1σ,where D is the stiffness matrix of the elastic constant, and the plastic strain rate is linear:(8)dεpl=∑dεpl.


## 4.1.2. Bedding Plane Opening

According to the ABAQUS documentation [27], when the pressure stress across the bedding plane becomes positive (tensile), exceeding the capacity of the bedding plane, the bedding plane is opened. Thus, opening the bedding plane creates an anisotropic elastic response at a point, and the joint system remains as(9)εanpsel≤εanel,where εanel is the component of direct elastic strain across the bedding plane and εanpsel is the component of hypothetical direct elastic strain across the bedding plane that is calculated as(10)εanpsel=−νEσa1+σa2,where E is the Young’s modulus of the material, ν is the Poisson’s ratio, and σa1,σa2 are the direct stresses in the bedding plane.


## 4.1.3. Bedding Plane Sliding

Normally, shear failure forms in the bedding, so the sliding failure surface is defined by(11)fa=τa−σatanφa−ca=0,where φa is the friction angle of bedding system a, ca is the cohesion for system a, τa is the shear stress on the bedding plane, and σa is the pressure stress normal to the bedding plane; so it is the transformation of Mohr−Coulomb law.As long asfa<0, the bedding plane in system a does not slip. When fa=0, the bedding plane in system a slips. The plastic strain on the system is then given by(12)dγaλpl=dε¯aplτaλτacosψa,dεanpl=dε¯aplsinψa,where dγaλpl is the plastic shear strain rate in direction λ (λ=1,2 are the orthogonal directions) on the bedding plane, dε¯apl is the magnitude of the plastic strain rate, τaλ is the component of shear stress on the bedding plane, ψa is the dilation angle for the bedding system, and dεanpl is the plastic strain rate normal to the bedding plane.


## 4.2. Comparison between Experiment and Simulation

### 4.2.1. Modeling

The size of model was 100 mm × 100 mm × 100 mm; pressure was applied on the upside and upper half side plane, respectively; the bottom and opposite lower half side plane is set with fixed boundary conditions (Figure12), and three-dimensional solid finite elements C3D8 were chosen as the model element algorithm.Figure 12
Three-dimensional numerical model mesh shape and boundary conditions.In order to make the shear effect obvious,σ and τ was the result of α = 70° under condition c. Pressure was applied on the upside and upper half plane, making σ = 20.096 MPa, τ = 55.212 MPa, bottom and lower half plane was fixed, as shown in Figure 12. According to Xie et al. [28], for slate, the reasonable friction angle can be 30–50°, and reasonable cohesion can be 2–20 MPa. In the simulation, cohesion (c = 7.394 MPa) and friction angle (φ=11.96°) were selected from the approximate result of condition a; even the friction angle may be small (φ=11.96°), and small friction angle is also acquired by Mao [22]. Table 2 lists the slate mechanical parameters for model calculation. Definition of the bedding plane direction relies on the three vectors of the plane (Figure 11), so simulation of the three conditions (Figure 4) in the shear tests could achieved.Table 2
Physicomechanical parameters of slate.
Material
Dry densityρd (kg/m3)
Poisson’s ratioυ
Elastic modulusEs (GPa)
Cohesionc (MPa)
Friction angleφ (°)
Slate
2771
0.27
78.23
7.394
11.96


### 4.2.2. Simulation Results

Figure13 presents the simulation results, and the models in each group from left to right correspond to (a), (b), and (c) in Figure 4. In condition (a), the model computes 22 increment steps with nonconvergence; for condition (b), it is 24 increment steps with nonconvergence, and for condition (c), the computation is converged at 39 increment steps. In Figure 13, the shear stress are all centralized on the side edge of the shear plane, so this area are all prone to shear failure, which corresponds well to the final shear failure that occurs in (a), (b), and (c) in Figure 7.Figure 13
Results of numerical simulation. (a) Shear stress. (b) Logarithmic strain. (c) Logarithmic strain vector.
(a)
(b)
(c)Figures13(b) and 13(c) are the strain and its vector; (a) shows the strain concentrated on the side plane in the middle, and the orientation of the shear stress is in the middle; it fits well with the test results of condition a. Similarly, for (b), tensile strain occurs on the upper side plane, and it may results in the initial cracks, as shown in Figure 10. For (c), the strain is concentrated on the plane of Y0Z, and the direction of the strain vector is normal to the side planes, so it agrees with the initial failure mechanism in the tests (Figure 10); the bedding plane preferentially opens on both the front and back side planes.In Figure14, the curve is the point in the middle of the upper half loading plane (Figure 12). Even for condition (c), the maximum is different, and the value of the drop in the middle of the test curve is approximately the same. The good agreement between simulation and test shows that the model satisfactorily describes the anisotropic deformational and strength behaviors of slate.Figure 14
Comparison between test and simulation (α = 70°) under loading modes: (a), (b), and (c).


## 4.2.1. Modeling

The size of model was 100 mm × 100 mm × 100 mm; pressure was applied on the upside and upper half side plane, respectively; the bottom and opposite lower half side plane is set with fixed boundary conditions (Figure12), and three-dimensional solid finite elements C3D8 were chosen as the model element algorithm.Figure 12
Three-dimensional numerical model mesh shape and boundary conditions.In order to make the shear effect obvious,σ and τ was the result of α = 70° under condition c. Pressure was applied on the upside and upper half plane, making σ = 20.096 MPa, τ = 55.212 MPa, bottom and lower half plane was fixed, as shown in Figure 12. According to Xie et al. [28], for slate, the reasonable friction angle can be 30–50°, and reasonable cohesion can be 2–20 MPa. In the simulation, cohesion (c = 7.394 MPa) and friction angle (φ=11.96°) were selected from the approximate result of condition a; even the friction angle may be small (φ=11.96°), and small friction angle is also acquired by Mao [22]. Table 2 lists the slate mechanical parameters for model calculation. Definition of the bedding plane direction relies on the three vectors of the plane (Figure 11), so simulation of the three conditions (Figure 4) in the shear tests could achieved.Table 2
Physicomechanical parameters of slate.
Material
Dry densityρd (kg/m3)
Poisson’s ratioυ
Elastic modulusEs (GPa)
Cohesionc (MPa)
Friction angleφ (°)
Slate
2771
0.27
78.23
7.394
11.96


## 4.2.2. Simulation Results

Figure13 presents the simulation results, and the models in each group from left to right correspond to (a), (b), and (c) in Figure 4. In condition (a), the model computes 22 increment steps with nonconvergence; for condition (b), it is 24 increment steps with nonconvergence, and for condition (c), the computation is converged at 39 increment steps. In Figure 13, the shear stress are all centralized on the side edge of the shear plane, so this area are all prone to shear failure, which corresponds well to the final shear failure that occurs in (a), (b), and (c) in Figure 7.Figure 13
Results of numerical simulation. (a) Shear stress. (b) Logarithmic strain. (c) Logarithmic strain vector.
(a)
(b)
(c)Figures13(b) and 13(c) are the strain and its vector; (a) shows the strain concentrated on the side plane in the middle, and the orientation of the shear stress is in the middle; it fits well with the test results of condition a. Similarly, for (b), tensile strain occurs on the upper side plane, and it may results in the initial cracks, as shown in Figure 10. For (c), the strain is concentrated on the plane of Y0Z, and the direction of the strain vector is normal to the side planes, so it agrees with the initial failure mechanism in the tests (Figure 10); the bedding plane preferentially opens on both the front and back side planes.In Figure14, the curve is the point in the middle of the upper half loading plane (Figure 12). Even for condition (c), the maximum is different, and the value of the drop in the middle of the test curve is approximately the same. The good agreement between simulation and test shows that the model satisfactorily describes the anisotropic deformational and strength behaviors of slate.Figure 14
Comparison between test and simulation (α = 70°) under loading modes: (a), (b), and (c).


## 5. Stability Analysis of Slate Tunnel

### 5.1. Modeling and Parameters

A tunnel construction is needed in Tongren city (Figure15). The stability analysis need to take all possible bedding plane dip angle (θ) into consideration, so the previous model was utilized for the numerical simulation. The size of model is set to be 55 m × 30 m × 50 m (Figure 16), the tunnel length is 30 m, and the depth of overlying rock above the tunnel is 30 m; the size was selected to ensure no interaction of the displacement field with boundaries. The surrounding rock properties and parameters were achieved in the shear test. It is widely believed that there is not much difference in friction angle between the rock and rock mass [29], so only cohesion c reduction is needed. Yang. et al. [30] concluded, for volcanic rocks and metamorphic rocks, the cohesion can be reduced by the following empirical formula:(13)cm=0.114e−0.48i−2+0.02ck,where i is the number of joint fractures per meter, ck is the rock cohesion (MPa), and cm is the weakened rock mass cohesion (MPa). As the average joint spacing of rock mass in the tunnel is about 0.4 m, so i is 2.5 and reduction coefficient is 0.1097. The reduction coefficient is often taken as 1/7 − 1/10 by experience [29]; for safety, the reduction coefficient of cohesion is set to be 0.1, and parameters of the two rock systems are demonstrated in Figure 17 and Table 3. In ABAQUS, the two systems are separated without interaction. The geometrical parameter of the tunnel considered is introduced by two intersecting circles (Figure 16). Considering the calculation amount in the simulations and the symmetry of the tunnel, only half of the total domain is considered.Figure 15
Construction scene of the tunnel.Figure 16
Schematic of shear strength parameters.Figure 17
Tunnel dimensions and 3D numerical model mesh shape with bedding plane direction.Table 3
Parameters of rock mass.
Material
Dry densityρd (kg/m3)
Poisson’s ratioυ
Elastic modulusEs (GPa)
Cohesionc (MPa)
Friction angleφ (°)
Slate rock Mass
Direction 1
2771
0.27
78.23
0.7394
11.96
Direction 2
1.0965
9.18Figure16 plots the relationship of the tunnel and the bedding plane. Dip angle (θ) between the bedding plane and horizontal plane is set to be 0°, 15°, 30°, 45°, 60°, 75°, and 90° to simulate the tunnel’s stability under all conditions. Also, three-dimensional solid elements (C3D8) are selected. For boundary conditions, gravity is applied to the whole model; the bottom of the model is fully constrained, the top of the model is unconstrained, and on the four side planes, only vertical deformations are allowed. For simplification, the effects of the excavation process, rock support, and reinforcement are not considered.


### 5.2. Modeling Results and Analysis

Figure18 presents the plastic zone under different dip angles; the failure region is centralized on the side wall. When it is 45°, the plastic area is the largest. Comparing the depth of the plastic zone on the side wall, the medium dip angle (45°) is the largest, followed by 60° and 30°. The sketch in the middle of Figure 18 shows the plastic zone ahead of the tunnel face, also the damage region of medium dip angle (45°) is the largest; the depth of the plastic area ahead of the tunnel is about 1.7 times the height, followed by 30° and 60° being 1.2 and 1.1 times the height, respectively.Figure 18
Tunnel plastic zone and plastic zone sketch ahead the tunnel under different dip angles.In order to investigate the stress and deformation of the tunnel, four monitoring points a, b, c, and d (Figure16) on each cross section along the tunnel axis are arranged.The displacements are plotted in Figure19, and the distance means the spacing apart from the tunnel face (Figure 16). Comparing the value, point d is the smallest, and it is the arch feet of the tunnel, so further comparison is needed for points a, b, and c. Generally, for tunnel excavation stability, when it is medium dip angle, the displacement is larger. Deformation is mainly concentrated on the vault and upside of the tunnel, especially for points a and b, displacement increases at an inconstant speed along the distance. 0–8 m is the rapid growth stage, so in the spacing of 0–8 m apart from tunnel face, upper side of tunnel excavation boundary need to control stability.Figure 19
Tunnel displacements variation of point a b, c, and d along the tunnel.
(a)
(b)
(c)
(d)Considering the characteristics of shear failure in slate, the maximum deviator stress is employed to analyse the tunnel mechanical stability. Figure20 plots the maximum deviator stress distribution along the tunnel, and the maximum points are marked with spots. As is shown, the red zone is concentrated on the side wall; it is the same as the plastic zone distribution (Figure 18), and the medium-dip angle (45°) is the largest. In this region (Figure 20), the maximum point position of 45° was selected, and the stress-strain curve at all dip angle (0°, 15°, 30°, 45°, 60°, 75°, and 90°) is plotted. In Figure 20, for medium dip angles (30°, 45°, and 60°), the curves yield at a lower stress. With strain increasing, the curve of the medium dip angles grows higher than the others.Figure 20
Stress (σ1−σ3) distribution of the tunnel and stress-strain curve of the selected point.In Figure21, the upper vertical axis is the stress (σ1−σ3) maximum of all values at all dip angles (0°, 15°, 30°, 45°, 60°, 75°, and 90°); the lower vertical axis is the dip angle that corresponds to the maximum value, and the horizontal axis is the distance from tunnel face. Indicated by the graph, the maximum is mainly centralized on the line of points c and d along the tunnel axis, and the stress (σ1−σ3) of the tunnel vault is minimum. The marked part (4–24 m) in Figure 21 illustrates that stress in the middle of tunnel length is higher, and the dip angle corresponds to marker part is medium dip angle similarly. So when medium dip angle is necessary, support is needed on the middle of sidewall.Figure 21
Maximum stress (σ1−σ3) and the corresponding dip angle variation along the tunnel.


## 5.1. Modeling and Parameters

A tunnel construction is needed in Tongren city (Figure15). The stability analysis need to take all possible bedding plane dip angle (θ) into consideration, so the previous model was utilized for the numerical simulation. The size of model is set to be 55 m × 30 m × 50 m (Figure 16), the tunnel length is 30 m, and the depth of overlying rock above the tunnel is 30 m; the size was selected to ensure no interaction of the displacement field with boundaries. The surrounding rock properties and parameters were achieved in the shear test. It is widely believed that there is not much difference in friction angle between the rock and rock mass [29], so only cohesion c reduction is needed. Yang. et al. [30] concluded, for volcanic rocks and metamorphic rocks, the cohesion can be reduced by the following empirical formula:(13)cm=0.114e−0.48i−2+0.02ck,where i is the number of joint fractures per meter, ck is the rock cohesion (MPa), and cm is the weakened rock mass cohesion (MPa). As the average joint spacing of rock mass in the tunnel is about 0.4 m, so i is 2.5 and reduction coefficient is 0.1097. The reduction coefficient is often taken as 1/7 − 1/10 by experience [29]; for safety, the reduction coefficient of cohesion is set to be 0.1, and parameters of the two rock systems are demonstrated in Figure 17 and Table 3. In ABAQUS, the two systems are separated without interaction. The geometrical parameter of the tunnel considered is introduced by two intersecting circles (Figure 16). Considering the calculation amount in the simulations and the symmetry of the tunnel, only half of the total domain is considered.Figure 15
Construction scene of the tunnel.Figure 16
Schematic of shear strength parameters.Figure 17
Tunnel dimensions and 3D numerical model mesh shape with bedding plane direction.Table 3
Parameters of rock mass.
Material
Dry densityρd (kg/m3)
Poisson’s ratioυ
Elastic modulusEs (GPa)
Cohesionc (MPa)
Friction angleφ (°)
Slate rock Mass
Direction 1
2771
0.27
78.23
0.7394
11.96
Direction 2
1.0965
9.18Figure16 plots the relationship of the tunnel and the bedding plane. Dip angle (θ) between the bedding plane and horizontal plane is set to be 0°, 15°, 30°, 45°, 60°, 75°, and 90° to simulate the tunnel’s stability under all conditions. Also, three-dimensional solid elements (C3D8) are selected. For boundary conditions, gravity is applied to the whole model; the bottom of the model is fully constrained, the top of the model is unconstrained, and on the four side planes, only vertical deformations are allowed. For simplification, the effects of the excavation process, rock support, and reinforcement are not considered.


## 5.2. Modeling Results and Analysis

Figure18 presents the plastic zone under different dip angles; the failure region is centralized on the side wall. When it is 45°, the plastic area is the largest. Comparing the depth of the plastic zone on the side wall, the medium dip angle (45°) is the largest, followed by 60° and 30°. The sketch in the middle of Figure 18 shows the plastic zone ahead of the tunnel face, also the damage region of medium dip angle (45°) is the largest; the depth of the plastic area ahead of the tunnel is about 1.7 times the height, followed by 30° and 60° being 1.2 and 1.1 times the height, respectively.Figure 18
Tunnel plastic zone and plastic zone sketch ahead the tunnel under different dip angles.In order to investigate the stress and deformation of the tunnel, four monitoring points a, b, c, and d (Figure16) on each cross section along the tunnel axis are arranged.The displacements are plotted in Figure19, and the distance means the spacing apart from the tunnel face (Figure 16). Comparing the value, point d is the smallest, and it is the arch feet of the tunnel, so further comparison is needed for points a, b, and c. Generally, for tunnel excavation stability, when it is medium dip angle, the displacement is larger. Deformation is mainly concentrated on the vault and upside of the tunnel, especially for points a and b, displacement increases at an inconstant speed along the distance. 0–8 m is the rapid growth stage, so in the spacing of 0–8 m apart from tunnel face, upper side of tunnel excavation boundary need to control stability.Figure 19
Tunnel displacements variation of point a b, c, and d along the tunnel.
(a)
(b)
(c)
(d)Considering the characteristics of shear failure in slate, the maximum deviator stress is employed to analyse the tunnel mechanical stability. Figure20 plots the maximum deviator stress distribution along the tunnel, and the maximum points are marked with spots. As is shown, the red zone is concentrated on the side wall; it is the same as the plastic zone distribution (Figure 18), and the medium-dip angle (45°) is the largest. In this region (Figure 20), the maximum point position of 45° was selected, and the stress-strain curve at all dip angle (0°, 15°, 30°, 45°, 60°, 75°, and 90°) is plotted. In Figure 20, for medium dip angles (30°, 45°, and 60°), the curves yield at a lower stress. With strain increasing, the curve of the medium dip angles grows higher than the others.Figure 20
Stress (σ1−σ3) distribution of the tunnel and stress-strain curve of the selected point.In Figure21, the upper vertical axis is the stress (σ1−σ3) maximum of all values at all dip angles (0°, 15°, 30°, 45°, 60°, 75°, and 90°); the lower vertical axis is the dip angle that corresponds to the maximum value, and the horizontal axis is the distance from tunnel face. Indicated by the graph, the maximum is mainly centralized on the line of points c and d along the tunnel axis, and the stress (σ1−σ3) of the tunnel vault is minimum. The marked part (4–24 m) in Figure 21 illustrates that stress in the middle of tunnel length is higher, and the dip angle corresponds to marker part is medium dip angle similarly. So when medium dip angle is necessary, support is needed on the middle of sidewall.Figure 21
Maximum stress (σ1−σ3) and the corresponding dip angle variation along the tunnel.


## 6. Conclusions

In order to investigate the tunnel stability in anisotropic rock, slate in eastern Guizhou was selected to achieve anisotropic mechanical behavior and parameters. A proposed numerical model that can represent the slate properties was verified by comparing simulation and test, and it was utilized to analys e the relationship between tunnel stability and dip angle (θ = 0°, 15°, 30°, 45°, 60°, 75°, and 90°). The conclusions are as follows:(1)
Bedding of slate in eastern Guizhou cemented closely, and slate can be regarded as a special continuous material. The shear strength parameters and failure pattern vary when the bedding plane direction is changed; by comparing the failure evolution process, weak bedding planes are the fundamental causes of anisotropic behaviors in slate.(2)
As bedding properties of slate fit the prerequisite of the numerical model by Zienkiewicz, it can consider three sets of anisotropic shear strength simultaneously; so it was used to represent the anisotropic behavior in the test. The numerical results show that it is tensile strain in bedding that causes the anisotropy in the initial shear failure.(3)
When dip angle (θ) is 45°, the tunnel plastic zone is the largest, especially for the face and side wall, and the depth ahead of face is 1.7 times the tunnel height. Displacement is mainly concentrated on the tunnel vault; it increases with distance from tunnel face, 0–8 m is the rapid growth region. The stress σ1−σ3 is centralized on the side wall; similarly, stress of 45° is the greatest. When medium dip angle is necessary, in the middle of the tunnel, support is needed on the sidewall.



---

*Source: 4653401-2019-10-13.xml*