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.
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 (Figure 1). 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.
(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 Figure 2(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.
(a)
(b)
The cubic specimens (100 mm × 100 mm × 100 mm) were prepared (Figure 3). 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.
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.
3. Anisotropic Properties Analysis
3.1. Experimental Results
The digital ultrasonic instrument (RSM-SY5) used for acoustic testing was in Figure 5, and results of the acoustic wave test (Table 1) indicate that (average velocity of P-wave along the bedding plane) is 6058.4 m/s, and (average velocity of P-wave across the bedding plane) is 5274.7 m/s, is 12.9% less than . 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.
(a)
(b)
(c)
3.2. Failure Pattern Analysis
The failure patterns (Figure 7) 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.
(a)
(b)
(c)
In condition c, failure appeared along the bedding direction on the loading plane (Figure 7(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 (Figure 8), 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.
(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: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].
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. Figure 10 shows the failure evolution of condition b and c, and the observed zone is plotted.
(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 . According to the research of Zienkiewicz and Pande in 1977 [26], a particular bedding plane is considered by the normal vector , where , 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, and , and similarly, . So and 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
(a)
(b)
The shear stress magnitude is defined as
The local strain components are the normal strain across the bedding plane:and the shear strain in direction () of the bedding plane iswhere is the strain tensor.
The linear strain rate in the bedding plane system can be divided into elastic and plastic parts:where is the total strain rate, is the elastic strain rate, and is the plastic strain rate. According to Hooke’s law,where D is the stiffness matrix of the elastic constant, and the plastic strain rate is linear:
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 aswhere is the component of direct elastic strain across the bedding plane and is the component of hypothetical direct elastic strain across the bedding plane that is calculated aswhere E is the Young’s modulus of the material, is the Poisson’s ratio, and 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 bywhere is the friction angle of bedding system a, is the cohesion for system a, is the shear stress on the bedding plane, and is the pressure stress normal to the bedding plane; so it is the transformation of Mohr−Coulomb law.
As long as , the bedding plane in system a does not slip. When , the bedding plane in system a slips. The plastic strain on the system is then given bywhere is the plastic shear strain rate in direction ( are the orthogonal directions) on the bedding plane, is the magnitude of the plastic strain rate, is the component of shear stress on the bedding plane, is the dilation angle for the bedding system, and 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 (Figure 12), and three-dimensional solid finite elements C3D8 were chosen as the model element algorithm.
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 () were selected from the approximate result of condition a; even the friction angle may be small (), 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.
4.2.2. Simulation Results
Figure 13 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.
(a)
(b)
(c)
Figures 13(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 Figure 14, 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.
5. Stability Analysis of Slate Tunnel
5.1. Modeling and Parameters
A tunnel construction is needed in Tongren city (Figure 15). 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:where i is the number of joint fractures per meter, is the rock cohesion (MPa), and 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 16 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
Figure 18 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.
In order to investigate the stress and deformation of the tunnel, four monitoring points a, b, c, and d (Figure 16) on each cross section along the tunnel axis are arranged.
The displacements are plotted in Figure 19, 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.
(a)
(b)
(c)
(d)
Considering the characteristics of shear failure in slate, the maximum deviator stress is employed to analyse the tunnel mechanical stability. Figure 20 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.
In Figure 21, the upper vertical axis is the stress () 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 () 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.
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 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge financial support from the China Guizhou Science and Technology Planning Project (grant no. 20172803).