# Failure Process and Stability Analysis of Rock Blocks in a Large Underground Excavation Based on a Numerical Method

Authors: Shijie Chen; Ming Xiao; Juntao Chen
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A numerical analysis method for block failure is proposed that is based on continuum mechanics. First, a mesh model that includes marked blocks was established based on the grid-based block identification method. Then, expressions of the contact force under various contact states were derived based on the explicit contact force algorithm, and a contact simulation method between blocks and the surrounding rock was proposed. The safety factors of the blocks were calculated based on the strength reduction method. This numerical analysis method can simulate both the continuous deformation of the surrounding rock and the discontinuous failure processes of the blocks. A simple example of a sliding block was used to evaluate the accuracy and rationality of the numerical method. Finally, combined with a deep underground excavation project under complex geological conditions, the stability of the blocks and rock were analyzed. The results indicate that the key blocks are damaged after excavation, the potentially dangerous blocks loosen and undergo large deformations, and the cracks between the blocks and the rock gradually increase as the excavation proceeds. The safety factors of the blocks change during the excavation. The numerical results demonstrate the influence of the surrounding rock on the failure process and on the stability of the blocks, and an effective analysis method is provided for the stability analysis of blocks under complex geological conditions.


---

## Body

## 1. Introduction

During the excavation of large underground caverns, complex geological discontinuities (faults, cracks, joints, and other kinds of discontinuities) are inevitably encountered. These discontinuities cut the rock mass into numerous regions. The opening, closing, and shearing of the discontinuities led to deformation and destruction of the rock mass. They facilitate the formation of rock blocks and local collapse, especially for an underground powerhouse that is of large size and has a high sidewall [1]. Therefore, rock block stability in a large underground excavation must be urgently investigated.For the stability analysis of rock blocks, the key block theory proposed by Warburton [2] and Goodman and Shi [3] is an effective approach. It established a strict mathematical algorithm for block stability analysis via geometrical and topological methods. The derived software, such as Unwedge [4, 5] and GeneralBlock [6], were widely used in the block analysis of underground caverns. Block theory makes the following assumptions for its applications in engineering: The blocks are considered as rigid bodies and the movability and failure modes of blocks are determined based only on the gravity stress and the force on the sliding face. In routine analysis, there are three types of block failure modes, namely, lifting, single-plane sliding, and double-plane sliding, while the safety factor of a block is typically determined via the rigid body limit equilibrium method (RBLEM) [2, 7]. These conventional stability analysis methods yield more accurate safety factors for shallow buried blocks. However, the accuracy is insufficient for deep buried blocks. The high and complex compressive stress from the surrounding rock, including the stress from nonsliding faces, has a substantial influence on its stability [8, 9]. According to engineering theory [10], neglecting the interactions between blocks and rocks may lead to deviations in the stability assessment and might increase the support cost for blocks.For blocks in complex geological conditions, analytical methods (such as RBLEM) cannot calculate the complex compressive stress accurately. In contrast, numerical methods are more suitable for complex block problems. Based on the DDA method, Fu et al. and Zhao et al. [11, 12] studied the failure modes of blocks during the excavation process of an underground powerhouse. Zhu et al. and Cui et al. [13, 14] used the discrete element method to analyze the interactions between blocks and calculated the block safety factor using limit equilibrium theory. Scholars frequently use the discontinuum method to analyze the block problems [15, 16]. Although it can reflect the deformation and failure process of the block, it has low calculation accuracy for the intact and homogenous surrounding rock [17, 18]. The continuum method is highly mature in underground engineering and it can effectively consider the geostress, excavation load, and deformation characteristics of the surrounding rock. Few studies have been conducted based on the continuum method. Zhang et al. [19, 20] proposed an element reconstruction modeling method for geological discontinuities that could identify blocks based on a finite element grid. It enables the introduction of the continuum method into block stability analysis. Zhang et al. [21] further proposed a general method for the stability evaluation of blocks that is based on the Flac3D software. The current numerical method for block stability analysis cannot consider both the continuous deformation of the rock and the discontinuous deformation of the block. The discontinuum method and the continuum method differ in terms of their advantages and the existing analysis methods cannot combine the advantages of the two methods effectively.This paper proposes a numerical analysis method for block failure that is based on continuum mechanics. It can simulate both the deformation characteristics of the surrounding rock and the failure process of blocks and solve for the safety factors. First, a block identification method that is based on a grid is introduced. Then, based on the explicit solution of the contact force algorithm [22, 23], the contact simulation method between the block and rock is proposed. The numerical method considers the point-to-point and point-to-surface contact types in the grid model, and it can simulate three contact states between blocks and rock: bonded, sliding, and separation. Finally, the safety factor of each block is calculated via the strength reduction method (SRM). For an underground powerhouse project that is crossed by a fault, the stability of the block and the surrounding rock were analyzed. The research results can provide effective references for the block stability under complex geological conditions.


## 2. Block Identification Based on a Grid

### 2.1. Basic Steps

The rock block can be divided into a finite block and an infinite block. The infinite block is on the boundary of the calculation model. Due to gravity and the interaction of the surrounding rock, the finite block tends to loosen or slide down after excavation. A finite block that is exposed to the excavation face can be referred to as a key block. After losing stability, it may cause a chain reaction and lead to damage to entire caverns. A block that may cause a chain reaction can be referred to as a potentially dangerous block. As the volume of a potentially dangerous block increases, the construction risk of cavern excavation also increases [24]. Dangerous block regions need to be investigated. Therefore, identifying these blocks prior to numerical calculation is highly important. It is necessary to build a grid model that contains complex blocks. This paper proposes a grid-based block identification method. Combined with the example that is presented in Figure 1, the basic steps are as follows:(1)
Obtain the necessary information on the geological discontinuities via geological exploration. This information includes the occurrence, spacing, and spatial distribution of the discontinuities. Using this information, the spatial relationship between the discontinuities and the excavation face is determined, as shown in Figure1(a).(2)
Establish the finite element model, including geological discontinuities via existing mature gridding techniques. The mesh tool of ANSYS is used to generate the mesh automatically. The element shapes are mainly quadrilateral and partly triangular.(3)
Run the algorithm for generating blocks to generate the finite blocks near the excavation face, which is described in detail in Section2.2. Then, mark the blocks. According to Figure 1(b), three key blocks and four potentially dangerous blocks have been identified and marked.(4)
Separate the common nodes on the interfaces between the marked blocks and the surrounding rock and identify the initial contact point pairs. A contact model on the interface will be constructed in a later calculation.Figure 1
Simple example for block identification. (a) Multiple discontinuities; (b) finite element model and marked blocks.
(a)
(b)Using the above basic steps, a finite element model that contains marked blocks can be obtained, which will be used for numerical calculation. The above basic steps can be programmed to realize automatic processing in 2D and 3D models. In contrast to conventional block theory, it can effectively mark all key blocks and potentially dangerous blocks around the excavation face.


### 2.2. Algorithm for Generating Blocks

Generating blocks is one of the most important steps in identifying a block in a grid model. The main objectives of the block generation algorithm are to aggregate the elements in the same spatial domain and determine the finite block.Suppose the grid model in step (2) containsm elements and n discontinuities. From the geometric point of view, the discontinuity j can be described as follows:(1)Ajx+Bjy+Cjz+Dj=0,where Aj, Bj, Cj, and Dj are parameters that define a spatial plane of discontinuity j.Point (xi, yi, zi) is the centroid of element i. Then, each element’s spatial location relationships with all discontinuities can be expressed by the following formulas:(2)Pi,j=Ajxi+Bjyi+Cjzi+DjAj2+Bj2+Cj2,(3)Hi,j=1,Pi,j>0−1,Pi,j<0,i=1,…,m;j=1,…,n.IfPi,j>0, element i is judged to be at the upper wall of discontinuity j and Hi,j is marked as 1. If Pi,j<0, element i is judged to be at the footwall of discontinuity j and Hi,j is marked as −1. The elements that correspond to the same value of Hi,j for all discontinuities are aggregated to generate blocks. If elements E1,E2,E3,…,Ek−1, and Ek form block B, then these elements should satisfy(4)H1,j=H2,j=H3,j=…=Hk−1,j=Hk,jj=1,2,3,…,n,and block B can be generated:(5)B=E1∪E2∪E3∪…∪Ek−1∪Ek.After searching and aggregating the elements, multiple block systems are identified. However, some blocks may be infinite blocks or may contain excavation elements. These blocks require further treatment.The infinite blocks are on the boundary and they are stable. Combine the infinite blocks into a single group. For the blocks that are cut by the excavation faces, special treatment is required. Mark all excavation elements as one group. The excavation faces cut the blocks into one or more new blocks. In the example that is shown in Figure2, block A is identified by only considering the geological discontinuities. After considering the excavation faces, it can be cut into block B and block C. The blocks near the underground cavern are always dangerous and require additional consideration.Figure 2
Block identification in consideration of excavation.


## 2.1. Basic Steps

The rock block can be divided into a finite block and an infinite block. The infinite block is on the boundary of the calculation model. Due to gravity and the interaction of the surrounding rock, the finite block tends to loosen or slide down after excavation. A finite block that is exposed to the excavation face can be referred to as a key block. After losing stability, it may cause a chain reaction and lead to damage to entire caverns. A block that may cause a chain reaction can be referred to as a potentially dangerous block. As the volume of a potentially dangerous block increases, the construction risk of cavern excavation also increases [24]. Dangerous block regions need to be investigated. Therefore, identifying these blocks prior to numerical calculation is highly important. It is necessary to build a grid model that contains complex blocks. This paper proposes a grid-based block identification method. Combined with the example that is presented in Figure 1, the basic steps are as follows:(1)
Obtain the necessary information on the geological discontinuities via geological exploration. This information includes the occurrence, spacing, and spatial distribution of the discontinuities. Using this information, the spatial relationship between the discontinuities and the excavation face is determined, as shown in Figure1(a).(2)
Establish the finite element model, including geological discontinuities via existing mature gridding techniques. The mesh tool of ANSYS is used to generate the mesh automatically. The element shapes are mainly quadrilateral and partly triangular.(3)
Run the algorithm for generating blocks to generate the finite blocks near the excavation face, which is described in detail in Section2.2. Then, mark the blocks. According to Figure 1(b), three key blocks and four potentially dangerous blocks have been identified and marked.(4)
Separate the common nodes on the interfaces between the marked blocks and the surrounding rock and identify the initial contact point pairs. A contact model on the interface will be constructed in a later calculation.Figure 1
Simple example for block identification. (a) Multiple discontinuities; (b) finite element model and marked blocks.
(a)
(b)Using the above basic steps, a finite element model that contains marked blocks can be obtained, which will be used for numerical calculation. The above basic steps can be programmed to realize automatic processing in 2D and 3D models. In contrast to conventional block theory, it can effectively mark all key blocks and potentially dangerous blocks around the excavation face.


## 2.2. Algorithm for Generating Blocks

Generating blocks is one of the most important steps in identifying a block in a grid model. The main objectives of the block generation algorithm are to aggregate the elements in the same spatial domain and determine the finite block.Suppose the grid model in step (2) containsm elements and n discontinuities. From the geometric point of view, the discontinuity j can be described as follows:(1)Ajx+Bjy+Cjz+Dj=0,where Aj, Bj, Cj, and Dj are parameters that define a spatial plane of discontinuity j.Point (xi, yi, zi) is the centroid of element i. Then, each element’s spatial location relationships with all discontinuities can be expressed by the following formulas:(2)Pi,j=Ajxi+Bjyi+Cjzi+DjAj2+Bj2+Cj2,(3)Hi,j=1,Pi,j>0−1,Pi,j<0,i=1,…,m;j=1,…,n.IfPi,j>0, element i is judged to be at the upper wall of discontinuity j and Hi,j is marked as 1. If Pi,j<0, element i is judged to be at the footwall of discontinuity j and Hi,j is marked as −1. The elements that correspond to the same value of Hi,j for all discontinuities are aggregated to generate blocks. If elements E1,E2,E3,…,Ek−1, and Ek form block B, then these elements should satisfy(4)H1,j=H2,j=H3,j=…=Hk−1,j=Hk,jj=1,2,3,…,n,and block B can be generated:(5)B=E1∪E2∪E3∪…∪Ek−1∪Ek.After searching and aggregating the elements, multiple block systems are identified. However, some blocks may be infinite blocks or may contain excavation elements. These blocks require further treatment.The infinite blocks are on the boundary and they are stable. Combine the infinite blocks into a single group. For the blocks that are cut by the excavation faces, special treatment is required. Mark all excavation elements as one group. The excavation faces cut the blocks into one or more new blocks. In the example that is shown in Figure2, block A is identified by only considering the geological discontinuities. After considering the excavation faces, it can be cut into block B and block C. The blocks near the underground cavern are always dangerous and require additional consideration.Figure 2
Block identification in consideration of excavation.


## 3. Contact Simulation Method between Blocks and the Surrounding Rock

Due to excavation unloading, a block may slide or separate from the surrounding rock. During the failure process of a block, complex contact actions occur at the interface between the block and the rock. Therefore, based on the explicit solution of the contact force algorithm and considering the various contact types, a contact simulation method between a block and the surrounding rock is proposed.


### 3.1. Explicit Integration Format that Considers the Contact Force

On the basis of the basic equation of the finite element method, the contact force of the interface and the inertial force are considered. After the finite element discretization and considering the contact force, the equilibrium equation of the nodes of blocks and the surrounding rock can be defined as follows:(6)Mu¨+Cu˙+Ku=F+R,where M, C, and K are the mass, damping, and stiffness matrices, respectively, of the interface node; u¨, u˙, and u are the acceleration, velocity, and displacement vectors, respectively; F is the external load vector; R is the contact force vector; and R=N+T, in which N and T are the normal and tangential components of R.The velocity and acceleration of the node are calculated via the central difference method and can be expressed as follows:(7)u˙t=12Δtut+Δt−ut−Δt,(8)u¨t=1Δt2ut+Δt−2ut+ut−Δt.Substituting (7) and (8) into (6) yields the node displacement at time step t+Δt:(9)ut+Δt=u¯t+Δt+Δut+Δt,(10)u¯t+Δt=ut+Δtu˙t+Δt2M−1Ft−Kut2,(11)Δut+Δt=Δt2M−1Rt,where t is the time, Δt is the time step size, u¯t+Δt is the nodal displacement vector without considering the contact force, and Δut+Δt is the vector of the additional displacement that is caused by the contact force.According to equation (9) to equation (11), the block displacement at t+Δt is determined by the motion state and the contact force at t, which is the main feature of the explicit algorithm. The motion state at t is known, while the contact force is unknown. The contact type and the contact force Rt can be judged and calculated according to the contact types and the contact states at t and t+Δt.


### 3.2. Calculation of the Contact Force for Several Contact Types

During the block identification, the initial contact point pairs were established between the block and the rock. Prior to excavation, the interface between the block and the rock is well cemented and the contacts are of point-to-point type. After the block slides, the interface node of the block slides relative to the interface node of the rock. At this time, the interface node is in contact with a surface of an element; this contact type is called point-to-surface type. Figure3 shows sketches of point-to-point and point-to-surface contact types. In each type, the contact state between the block and the surrounding rock can be bonded, sliding, or separated. The following presents the calculation method for the contact force under different contact types.Figure 3
Sketch of contact types between a block and a rock. (a) Point-to-point contact type; (b) point-to-surface contact type.
(a)
(b)


#### 3.2.1. Contact Force Calculation for the Point-to-point Types

Assume that the block and rock are in a bonded state at time stept+Δt. Therefore, for the contact node pair l and l′, they should satisfy the deformation coordination conditions, namely, nonmutual embedding in the normal direction and nonrelative slip in the tangent direction:(12)nlTul′t+Δt−ult+Δt=0,(13)τlTul′t+Δt−ult+Δt=τlTul′t−ult,where nl is the unit normal vector of the contact node pair, which points to l from l′, and τl is the corresponding unit tangent vector.By substituting equation (9) into equations (12) and (13) and according to Rlt=−Rlt, the normal contact force and tangent contact force are obtained:(14)Nlt=2MlMl′Ml+Ml′Δt2Δ1lnl,(15)Tlt=2MlMl′Ml+Ml′Δt2Δ2lτl,(16)Δ1l=nlTu¯l′t+Δt−u¯lt+Δt,(17)Δ2l=τlTu¯l′t+Δt−u¯lt+Δt−ul′t−ult,where Δ1l and Δ2l are the relative normal and tangential displacements, respectively, of the contact point pair without considering the contact force; Ml and Ml′ are the lumped masses of l and l′, respectively; and Nlt and Tlt are the normal and tangential components of RltRlt=Nlt+Tlt.The above formulas are obtained under the assumption of a bonded contact state. However, the contact point pairs may slide or separate during the failure process of the block. Therefore, it is necessary to judge the contact state and correct the contact force after each step.IfΔ1l > 0, the contact point pairs have a tendency to be embedded with each other. Hence, the contact state may be a bonded state or a sliding state. In addition, the tangential contact force is compared with the shear strength of the interface. If Tlt>μs⋅Nlt+cA, the contact point is in a sliding contact state. Therefore, the contact force should be corrected as follows:(18)Tlt=μd⋅Nlt⋅TltTlt.IfΔ1l < 0, the contact point pairs tend to be separated from each other. Therefore, the contact state may be a bonded state or a separation state. In addition, the contact force and the cohesive force are compared. If Tlt2+Nlt2>cA, the contact point is in a separation state. Hence, the contact force should be(19)Nlt=0,Tlt=0,where μs and μd are the static friction coefficient and the kinetic friction coefficient, respectively; A is the control area of node l; and c is the cohesion force of the interface. If the contact point has not been in a sliding or separation state prior to time step t+Δt, then c > 0; otherwise, c = 0.


#### 3.2.2. Contact Force Calculation for the Point-to-Surface Type

When the contact type is the point-to-surface type, it is again assumed that the block and the rock are in a bonded state at time stept+Δt. Hence, the contact point on the contact surface that corresponds to the contact node l is l′, and l′ is on the surface that belongs to element E. The displacement field ul′t and the equivalent lumped mass Ml′ of node l′ can be calculated via the shape function interpolation of the finite element:(20)ul′t=∑jϕjujt,mj=ϕjMj∑iϕi2,Ml′=∑jmj,where ϕj is the shape function value of node j in E at l′, Mj is the lumped mass of j, and mj is the mass contribution of j at l′.According to the deformation coordination conditions andRlt=−Rlt, the contact force expressions of l can also be obtained as equations (14) and (15).In addition, it is necessary to correct the contact state and the contact force after each step. The cohesion under the point-to-surface type is zero. Therefore, the tensile strength of the interface is zero and the shear strength isμs⋅Nlt.IfΔ1l > 0 and Tlt>μs⋅Nlt, the contact point is in the sliding contact state; hence, the contact force can be calculated via equation (18). If Δ1l < 0, the contact point is in the separation state; therefore, the contact force can be calculated via equation (19).


### 3.3. Basic Steps of the Contact Simulation Method

Prior to excavation, the nodes on the interface are bonded in the point-to-point contact type. When the contact point pairs break the cohesive force and enter into a sliding or separating state, the elements and the element surfaces that are in contact with these points should be identified via contact search. The contact type and the contact force must be determined at each time step. The basic steps of the calculation method at time stept+Δt are illustrated in Figure 4.Figure 4
Flowchart of the contact simulation method.The contact simulation method need not determine the spring stiffness. It can realize the discontinuity characteristic in the finite element model.


## 3.1. Explicit Integration Format that Considers the Contact Force

On the basis of the basic equation of the finite element method, the contact force of the interface and the inertial force are considered. After the finite element discretization and considering the contact force, the equilibrium equation of the nodes of blocks and the surrounding rock can be defined as follows:(6)Mu¨+Cu˙+Ku=F+R,where M, C, and K are the mass, damping, and stiffness matrices, respectively, of the interface node; u¨, u˙, and u are the acceleration, velocity, and displacement vectors, respectively; F is the external load vector; R is the contact force vector; and R=N+T, in which N and T are the normal and tangential components of R.The velocity and acceleration of the node are calculated via the central difference method and can be expressed as follows:(7)u˙t=12Δtut+Δt−ut−Δt,(8)u¨t=1Δt2ut+Δt−2ut+ut−Δt.Substituting (7) and (8) into (6) yields the node displacement at time step t+Δt:(9)ut+Δt=u¯t+Δt+Δut+Δt,(10)u¯t+Δt=ut+Δtu˙t+Δt2M−1Ft−Kut2,(11)Δut+Δt=Δt2M−1Rt,where t is the time, Δt is the time step size, u¯t+Δt is the nodal displacement vector without considering the contact force, and Δut+Δt is the vector of the additional displacement that is caused by the contact force.According to equation (9) to equation (11), the block displacement at t+Δt is determined by the motion state and the contact force at t, which is the main feature of the explicit algorithm. The motion state at t is known, while the contact force is unknown. The contact type and the contact force Rt can be judged and calculated according to the contact types and the contact states at t and t+Δt.


## 3.2. Calculation of the Contact Force for Several Contact Types

During the block identification, the initial contact point pairs were established between the block and the rock. Prior to excavation, the interface between the block and the rock is well cemented and the contacts are of point-to-point type. After the block slides, the interface node of the block slides relative to the interface node of the rock. At this time, the interface node is in contact with a surface of an element; this contact type is called point-to-surface type. Figure3 shows sketches of point-to-point and point-to-surface contact types. In each type, the contact state between the block and the surrounding rock can be bonded, sliding, or separated. The following presents the calculation method for the contact force under different contact types.Figure 3
Sketch of contact types between a block and a rock. (a) Point-to-point contact type; (b) point-to-surface contact type.
(a)
(b)


### 3.2.1. Contact Force Calculation for the Point-to-point Types

Assume that the block and rock are in a bonded state at time stept+Δt. Therefore, for the contact node pair l and l′, they should satisfy the deformation coordination conditions, namely, nonmutual embedding in the normal direction and nonrelative slip in the tangent direction:(12)nlTul′t+Δt−ult+Δt=0,(13)τlTul′t+Δt−ult+Δt=τlTul′t−ult,where nl is the unit normal vector of the contact node pair, which points to l from l′, and τl is the corresponding unit tangent vector.By substituting equation (9) into equations (12) and (13) and according to Rlt=−Rlt, the normal contact force and tangent contact force are obtained:(14)Nlt=2MlMl′Ml+Ml′Δt2Δ1lnl,(15)Tlt=2MlMl′Ml+Ml′Δt2Δ2lτl,(16)Δ1l=nlTu¯l′t+Δt−u¯lt+Δt,(17)Δ2l=τlTu¯l′t+Δt−u¯lt+Δt−ul′t−ult,where Δ1l and Δ2l are the relative normal and tangential displacements, respectively, of the contact point pair without considering the contact force; Ml and Ml′ are the lumped masses of l and l′, respectively; and Nlt and Tlt are the normal and tangential components of RltRlt=Nlt+Tlt.The above formulas are obtained under the assumption of a bonded contact state. However, the contact point pairs may slide or separate during the failure process of the block. Therefore, it is necessary to judge the contact state and correct the contact force after each step.IfΔ1l > 0, the contact point pairs have a tendency to be embedded with each other. Hence, the contact state may be a bonded state or a sliding state. In addition, the tangential contact force is compared with the shear strength of the interface. If Tlt>μs⋅Nlt+cA, the contact point is in a sliding contact state. Therefore, the contact force should be corrected as follows:(18)Tlt=μd⋅Nlt⋅TltTlt.IfΔ1l < 0, the contact point pairs tend to be separated from each other. Therefore, the contact state may be a bonded state or a separation state. In addition, the contact force and the cohesive force are compared. If Tlt2+Nlt2>cA, the contact point is in a separation state. Hence, the contact force should be(19)Nlt=0,Tlt=0,where μs and μd are the static friction coefficient and the kinetic friction coefficient, respectively; A is the control area of node l; and c is the cohesion force of the interface. If the contact point has not been in a sliding or separation state prior to time step t+Δt, then c > 0; otherwise, c = 0.


### 3.2.2. Contact Force Calculation for the Point-to-Surface Type

When the contact type is the point-to-surface type, it is again assumed that the block and the rock are in a bonded state at time stept+Δt. Hence, the contact point on the contact surface that corresponds to the contact node l is l′, and l′ is on the surface that belongs to element E. The displacement field ul′t and the equivalent lumped mass Ml′ of node l′ can be calculated via the shape function interpolation of the finite element:(20)ul′t=∑jϕjujt,mj=ϕjMj∑iϕi2,Ml′=∑jmj,where ϕj is the shape function value of node j in E at l′, Mj is the lumped mass of j, and mj is the mass contribution of j at l′.According to the deformation coordination conditions andRlt=−Rlt, the contact force expressions of l can also be obtained as equations (14) and (15).In addition, it is necessary to correct the contact state and the contact force after each step. The cohesion under the point-to-surface type is zero. Therefore, the tensile strength of the interface is zero and the shear strength isμs⋅Nlt.IfΔ1l > 0 and Tlt>μs⋅Nlt, the contact point is in the sliding contact state; hence, the contact force can be calculated via equation (18). If Δ1l < 0, the contact point is in the separation state; therefore, the contact force can be calculated via equation (19).


## 3.2.1. Contact Force Calculation for the Point-to-point Types

Assume that the block and rock are in a bonded state at time stept+Δt. Therefore, for the contact node pair l and l′, they should satisfy the deformation coordination conditions, namely, nonmutual embedding in the normal direction and nonrelative slip in the tangent direction:(12)nlTul′t+Δt−ult+Δt=0,(13)τlTul′t+Δt−ult+Δt=τlTul′t−ult,where nl is the unit normal vector of the contact node pair, which points to l from l′, and τl is the corresponding unit tangent vector.By substituting equation (9) into equations (12) and (13) and according to Rlt=−Rlt, the normal contact force and tangent contact force are obtained:(14)Nlt=2MlMl′Ml+Ml′Δt2Δ1lnl,(15)Tlt=2MlMl′Ml+Ml′Δt2Δ2lτl,(16)Δ1l=nlTu¯l′t+Δt−u¯lt+Δt,(17)Δ2l=τlTu¯l′t+Δt−u¯lt+Δt−ul′t−ult,where Δ1l and Δ2l are the relative normal and tangential displacements, respectively, of the contact point pair without considering the contact force; Ml and Ml′ are the lumped masses of l and l′, respectively; and Nlt and Tlt are the normal and tangential components of RltRlt=Nlt+Tlt.The above formulas are obtained under the assumption of a bonded contact state. However, the contact point pairs may slide or separate during the failure process of the block. Therefore, it is necessary to judge the contact state and correct the contact force after each step.IfΔ1l > 0, the contact point pairs have a tendency to be embedded with each other. Hence, the contact state may be a bonded state or a sliding state. In addition, the tangential contact force is compared with the shear strength of the interface. If Tlt>μs⋅Nlt+cA, the contact point is in a sliding contact state. Therefore, the contact force should be corrected as follows:(18)Tlt=μd⋅Nlt⋅TltTlt.IfΔ1l < 0, the contact point pairs tend to be separated from each other. Therefore, the contact state may be a bonded state or a separation state. In addition, the contact force and the cohesive force are compared. If Tlt2+Nlt2>cA, the contact point is in a separation state. Hence, the contact force should be(19)Nlt=0,Tlt=0,where μs and μd are the static friction coefficient and the kinetic friction coefficient, respectively; A is the control area of node l; and c is the cohesion force of the interface. If the contact point has not been in a sliding or separation state prior to time step t+Δt, then c > 0; otherwise, c = 0.


## 3.2.2. Contact Force Calculation for the Point-to-Surface Type

When the contact type is the point-to-surface type, it is again assumed that the block and the rock are in a bonded state at time stept+Δt. Hence, the contact point on the contact surface that corresponds to the contact node l is l′, and l′ is on the surface that belongs to element E. The displacement field ul′t and the equivalent lumped mass Ml′ of node l′ can be calculated via the shape function interpolation of the finite element:(20)ul′t=∑jϕjujt,mj=ϕjMj∑iϕi2,Ml′=∑jmj,where ϕj is the shape function value of node j in E at l′, Mj is the lumped mass of j, and mj is the mass contribution of j at l′.According to the deformation coordination conditions andRlt=−Rlt, the contact force expressions of l can also be obtained as equations (14) and (15).In addition, it is necessary to correct the contact state and the contact force after each step. The cohesion under the point-to-surface type is zero. Therefore, the tensile strength of the interface is zero and the shear strength isμs⋅Nlt.IfΔ1l > 0 and Tlt>μs⋅Nlt, the contact point is in the sliding contact state; hence, the contact force can be calculated via equation (18). If Δ1l < 0, the contact point is in the separation state; therefore, the contact force can be calculated via equation (19).


## 3.3. Basic Steps of the Contact Simulation Method

Prior to excavation, the nodes on the interface are bonded in the point-to-point contact type. When the contact point pairs break the cohesive force and enter into a sliding or separating state, the elements and the element surfaces that are in contact with these points should be identified via contact search. The contact type and the contact force must be determined at each time step. The basic steps of the calculation method at time stept+Δt are illustrated in Figure 4.Figure 4
Flowchart of the contact simulation method.The contact simulation method need not determine the spring stiffness. It can realize the discontinuity characteristic in the finite element model.


## 4. Calculation Method for the Block Safety Factor

The strength reduction method (SRM) is applied to calculate the block safety factor and the block stability is quantitatively evaluated. Since the mechanical properties of the interface between the block and the rock play an important role in the stability of the block, the critical instability state can be calculated by reducing the strength parameters of the interface:(21)c′=cFs,tanφ′=tanφFs,where Fs is the strength reduction factor, c and c′ are the cohesion before and after the strength reduction, and φ and φ′ are friction angles before and after the strength reduction.First, the calculation is conducted underFs = 1. A stable block can reach a force-balanced state and a stable displacement value in finite steps, of which the strength reduction factor Fs is considered to be within [1, +∞]. Otherwise, the block is in an unstable state and Fs is within [0, 1). In the corresponding interval, Fs is gradually increased with a gradient of 0.01 until the block is in the critical instability state. The critical state can be judged by the displacement mutation criterion [25]. The reduction factor at the critical state is considered to be the safety factor of the block.For a falling block in the top arch, the safety factor is close to zero. For an unmovable block, the safety factor approaches positive infinity. IfFs < 0.1, the block falls freely, and if Fs > 100, the block is unmovable.


## 5. Example Evaluation

The contact simulation method is embedded into the explicit finite element calculation framework in the Fortran language to form a numerical analysis program for block failure.To evaluate the accuracy and the rationality of the numerical analysis program and the solution method for the safety factor, a simple yet representative example of a sliding block is used to compare the numerical results and the analytical results. As shown in Figure5, the model consists of a sliding block and a base. The dip angle of the base is 40°. The size of the sliding block is 0.25 m × 0.25 m. The materials are both elastic. The elastic moduli are both 20 MPa, the Poisson ratio is 0.3, and the density is 2.0 g·cm−3. Two discontinuities are set on the left and bottom of the sliding block. The cohesion of the discontinuities is 0 MPa, and the friction angle is represented as φ. The bottom and lateral sides of the base are fixed. The block slides along the slope from rest under gravity.Figure 5
Sliding block model.During the numerical calculation process, the sliding block separates from the left discontinuity and the contact state on the bottom discontinuity changes from the bonded contact to sliding contact. The point-to-point contact type and the point-to-surface contact type occur alternately.The displacement of the block can be calculated analytically as follows:(22)st=12gsin40°−cos40°tanφt2,where st is the displacement of the sliding block, g is the gravitational acceleration, t is time, and φ is the friction angle of the discontinuity. If φ < 40°, the resting block will slide.In this example, the numerical calculations are conducted under the friction angles of 25°, 30°, and 35°. The displacement time-history of the sliding block that is calculated via this numerical method is compared with the analytical result, as shown in Figure6. The numerical solution accords with the analytical solution, thereby demonstrating the accuracy of the contact simulation method between a block and the surrounding rock. Overall, this numerical analysis program can simulate various contact types with satisfactory accuracy.Figure 6
Displacement results that were obtained numerically and analytically.The safety factor of the sliding block is calculated by reducing the strength parameters of the sliding surface, as shown in Figure7. In the case of a friction angle of 35°, when the reduction factor exceeds 0.83, the block is in the critical state. Therefore, the safety factor Fs−35 of the block is 0.83; similarly, the safety factors under the friction angles of 30° and 25° are 0.68 and 0.55, respectively. The expected safety factors according to RBLEM are 0.834, 0.688, and 0.556, respectively. The safety factors that are obtained via the numerical method are very close to those that are obtained via RBLEM. Therefore, the numerical method for calculating the safety factor of a block is accurate and reliable.Figure 7
Safety factor of the sliding block.


## 6. Engineering Example

### 6.1. Engineering Profile and Calculation Model

Hydropower stations in the southeast of China are always large in scale and the underground powerhouses are complex in structure. Moreover, the geological structure is complex and various geological discontinuities have developed. An underground powerhouse in the southeast of China has a size of 154.0 m × 26 m × 54.8 m and a depth of approximately 200 m. In the middle of the powerhouse, a fault F1 with a thickness of 5.0 m is obliquely crossed and it intersects with a thin fault F2. There are also two distinct fractures that intersect with fault F1 at the top of the powerhouse. Under the cut of multiple geological discontinuities, block failure tends to occur at the top and the left sidewall of the powerhouse. Figure8 illustrates the relative spatial positions between discontinuities and the excavation faces, the calculation range, and the boundary conditions. The underground cavern is excavated in 8 stages from top to bottom. In this study, a 2D model is established for analyzing the stability of the blocks and the surrounding rock.Figure 8
Sketch of the model for an underground powerhouse.The finite element model, including geological discontinuities, is shown in Figure9. After using the block identification method, four blocks are identified near the excavation faces: blocks B1 and B2 at the top and blocks B3 and B4 in the sidewall.Figure 9
Calculation model and marked blocks.


### 6.2. Calculation Conditions

The mechanical parameters of the rock, fault, and interface are provided in Table1. The materials are based on the Mohr–Coulomb criterion. The initial geostress field is obtained via stress inversion of measured points. The horizontal lateral pressure coefficient is approximately 1.10. According to the maximum principle, the stress that surrounds the caverns is between −6.0 MPa and −7.6 MPa. The bottom of the model is fixed and the left and right sidewalls are fixed with a normality constraint. The top is a free boundary and is subject to the gravitational stress of the overlying rock mass.Table 1
Mechanical parameters of the materials.
Material
Deformation modulus (GPa)
Poisson ratio
Cohesion force (MPa)
Frictional angle (°)
Density (g·cm-3)
Tensile strength (MPa)
Rock
15
0.3
0.6
35.0
2.7
1.5
Fault
2
0.4
0.5
25.0
2.0
0.5
Interface
0.0
30.0The numerical analysis program for block failure is adopted to analyze the failure process and the stability of the blocks under 8 stages of excavation. The initial geostress balance was carried out in 930 time steps. The simulation of 5000 time steps for each stage excavation could reach a stable deformation of the surrounding rock.


### 6.3. Results and Analysis

#### 6.3.1. Failure Process of Blocks during Excavation

A monitoring scheme, which is illustrated in Figure10, is used to monitor the deformation process of the blocks. Figure 11 shows the deformation of the surrounding rock after the 2nd, 4th, 6th, and 8th stages of excavation of the powerhouse (the displacement is magnified by a factor of 15). Figure 12 shows the monitored displacement. The monitored displacement gradually increases as the excavation progresses. The point P1 in block B1 undergoes a large displacement and continues to increase after the first stage of excavation. Hence, block B1 falls freely during excavation. The displacement of point P2 in block B2 is 14 mm after the first stage due to the excavation unloading. It slowly increases as the excavation proceeds and stabilizes at approximately 30 mm. Blocks B3 and B4 are completely exposed after the 5th stage of excavation and their displacements are increased by 22 mm and 18 mm, respectively. According to Figures 9(b)–9(d), blocks B3 and B4 undergo large deformations and cracks are generated between the blocks and the rock. The depth and width of the cracks gradually increase as the excavation proceeds. Points P4 and P5 are a contact point pair prior to excavation. After the 5th excavation stage, they have a relative displacement of approximately 10 mm and a crack occurs in their interface. The relative displacement eventually reaches approximately 25 mm. Hence, the interface deforms discontinuously. After the excavation has been completed, except for block B1 falling off, the blocks does not suffer from instability or damage. However, the blocks have loosened from the surrounding rock and there are high security risks.Figure 10
Monitoring schemes.Figure 11
Deformation process during excavation (the displacement is amplified by a factor of 15): (a) 2nd stage; (b) 4th stage; (c) 6th stage; (d) 8th stage.
(a)
(b)
(c)
(d)Figure 12
Monitored displacement.If using key block theory, only blocks B1 and B3, which are exposed on the excavation face, can be identified and it cannot simulate the failure process of the blocks. The proposed numerical method can effectively identify the key blocks, timely discover potentially dangerous blocks, and simulate the failure process of the blocks.


#### 6.3.2. Deformation Characteristics of the Blocks and the Surrounding Rock

To reflect the deformation law of blocks and surrounding rock, an additional calculation in which contact is not considered is conducted for comparison with the result in which contact is considered, as shown in Figure13. The case of without contact means the blocks and the surrounding rock are assumed to be tied together during modeling, and the common nodes on the interfaces are not separated. Without consideration of contact, the failure of the blocks cannot be simulated. Only a large displacement of the rock through which the fault passes is observed. With the consideration of contact, a block falls or separates from the surrounding rock and a discontinuous displacement field is observed near the interface of the fault and the rock. The maximum displacement when contact is considered approximately 35 mm larger than without contact. In the upper contact interface of block B4, the displacement changes substantially, which further indicates the loosening of block B4.Figure 13
Deformation of the surrounding rock. (a) Without contact and (b) considering contact.
(a)
(b)When contact between the blocks and the surrounding rock is considered, the simulation results are more consistent with the learned empirical rules. Therefore, this numerical analysis method better accords with the actual situation.


#### 6.3.3. Stress Evolution of the Interfaces between the Blocks and Rock

To analyze the stress evolution process of the interfaces between the blocks and rock, block B2 is selected as the research object. Three monitored points, namely, F1, F2, and F3, are arranged in the middles of the interfaces, as illustrated in Figure10. The normal stress and shear stress evolutions at the three points are plotted in Figure 14. During excavation, the interfaces are relieved of the initial stress. The stresses of interfaces F2 and F3 gradually increase as the excavation progresses, while the stress of interface F1 decreases. Nonsliding face F3 has the largest stress with block B2. Hence, the block is subjected to multiple stresses from the surrounding rock during its failure process, in addition to the stress from the sliding face.Figure 14
Stress evolutions of the interfaces. (a) Norm stress and (b) shear stress.
(a)
(b)


#### 6.3.4. Safety Factors of the Blocks

During the calculation of each excavation step, the safety factors of the blocks are obtained via reduction of the shear strength parameters of the interfaces. In the 1st stage of the excavation, the displacement of block B1 is large and continuously growing, even if the strength reduction factor is less than 0.1. Block B1 falls freely from the top of the cavern. Block B2 is fully exposed after the 1st stage of the excavation and blocks B3 and B4 are fully exposed after the 5th stage of the excavation. Figure15 plots the displacements of the blocks as functions of the strength reduction factor during the 5th stage of the excavation. For block B2, when the reduction factor exceeds 4.51, the monitored displacement suddenly increases; hence, the safety factor of block B2 is 4.51. Similarly, the safety factors of blocks B3 and B4 are 1.45 and 1.49, respectively.Figure 15
Safety factors of three blocks after the 5th stage excavation.The safety factors of the blocks during each stage of the excavation are plotted in Figure16. The changes in the safety factors vary among the stages. As the excavation progresses, the safety factor of block B2, which is in the top arch, gradually increases, while the safety factors of the blocks in the sidewall gradually decrease. However, the magnitudes of the changes among the stages are small. The block stability is greatly affected by geostress. Although the excavation has adjusted the stress of the surrounding rock, its influence degree is not as great as that of the geostress. So, block depth and excavation process are two important factors affecting the safety degree of blocks, and the influences of the excavation process to the safety factors of deep buried blocks are relatively small.Figure 16
Variations of the safety factors of three blocks.According to the safety factors of the blocks, in addition to strengthening the support of key blocks B1 and B3, it is also necessary to focus on potentially dangerous blocks B2 and B4. The potentially dangerous block in the left sidewall is more dangerous than that in the top arch.


## 6.1. Engineering Profile and Calculation Model

Hydropower stations in the southeast of China are always large in scale and the underground powerhouses are complex in structure. Moreover, the geological structure is complex and various geological discontinuities have developed. An underground powerhouse in the southeast of China has a size of 154.0 m × 26 m × 54.8 m and a depth of approximately 200 m. In the middle of the powerhouse, a fault F1 with a thickness of 5.0 m is obliquely crossed and it intersects with a thin fault F2. There are also two distinct fractures that intersect with fault F1 at the top of the powerhouse. Under the cut of multiple geological discontinuities, block failure tends to occur at the top and the left sidewall of the powerhouse. Figure8 illustrates the relative spatial positions between discontinuities and the excavation faces, the calculation range, and the boundary conditions. The underground cavern is excavated in 8 stages from top to bottom. In this study, a 2D model is established for analyzing the stability of the blocks and the surrounding rock.Figure 8
Sketch of the model for an underground powerhouse.The finite element model, including geological discontinuities, is shown in Figure9. After using the block identification method, four blocks are identified near the excavation faces: blocks B1 and B2 at the top and blocks B3 and B4 in the sidewall.Figure 9
Calculation model and marked blocks.


## 6.2. Calculation Conditions

The mechanical parameters of the rock, fault, and interface are provided in Table1. The materials are based on the Mohr–Coulomb criterion. The initial geostress field is obtained via stress inversion of measured points. The horizontal lateral pressure coefficient is approximately 1.10. According to the maximum principle, the stress that surrounds the caverns is between −6.0 MPa and −7.6 MPa. The bottom of the model is fixed and the left and right sidewalls are fixed with a normality constraint. The top is a free boundary and is subject to the gravitational stress of the overlying rock mass.Table 1
Mechanical parameters of the materials.
Material
Deformation modulus (GPa)
Poisson ratio
Cohesion force (MPa)
Frictional angle (°)
Density (g·cm-3)
Tensile strength (MPa)
Rock
15
0.3
0.6
35.0
2.7
1.5
Fault
2
0.4
0.5
25.0
2.0
0.5
Interface
0.0
30.0The numerical analysis program for block failure is adopted to analyze the failure process and the stability of the blocks under 8 stages of excavation. The initial geostress balance was carried out in 930 time steps. The simulation of 5000 time steps for each stage excavation could reach a stable deformation of the surrounding rock.


## 6.3. Results and Analysis

### 6.3.1. Failure Process of Blocks during Excavation

A monitoring scheme, which is illustrated in Figure10, is used to monitor the deformation process of the blocks. Figure 11 shows the deformation of the surrounding rock after the 2nd, 4th, 6th, and 8th stages of excavation of the powerhouse (the displacement is magnified by a factor of 15). Figure 12 shows the monitored displacement. The monitored displacement gradually increases as the excavation progresses. The point P1 in block B1 undergoes a large displacement and continues to increase after the first stage of excavation. Hence, block B1 falls freely during excavation. The displacement of point P2 in block B2 is 14 mm after the first stage due to the excavation unloading. It slowly increases as the excavation proceeds and stabilizes at approximately 30 mm. Blocks B3 and B4 are completely exposed after the 5th stage of excavation and their displacements are increased by 22 mm and 18 mm, respectively. According to Figures 9(b)–9(d), blocks B3 and B4 undergo large deformations and cracks are generated between the blocks and the rock. The depth and width of the cracks gradually increase as the excavation proceeds. Points P4 and P5 are a contact point pair prior to excavation. After the 5th excavation stage, they have a relative displacement of approximately 10 mm and a crack occurs in their interface. The relative displacement eventually reaches approximately 25 mm. Hence, the interface deforms discontinuously. After the excavation has been completed, except for block B1 falling off, the blocks does not suffer from instability or damage. However, the blocks have loosened from the surrounding rock and there are high security risks.Figure 10
Monitoring schemes.Figure 11
Deformation process during excavation (the displacement is amplified by a factor of 15): (a) 2nd stage; (b) 4th stage; (c) 6th stage; (d) 8th stage.
(a)
(b)
(c)
(d)Figure 12
Monitored displacement.If using key block theory, only blocks B1 and B3, which are exposed on the excavation face, can be identified and it cannot simulate the failure process of the blocks. The proposed numerical method can effectively identify the key blocks, timely discover potentially dangerous blocks, and simulate the failure process of the blocks.


### 6.3.2. Deformation Characteristics of the Blocks and the Surrounding Rock

To reflect the deformation law of blocks and surrounding rock, an additional calculation in which contact is not considered is conducted for comparison with the result in which contact is considered, as shown in Figure13. The case of without contact means the blocks and the surrounding rock are assumed to be tied together during modeling, and the common nodes on the interfaces are not separated. Without consideration of contact, the failure of the blocks cannot be simulated. Only a large displacement of the rock through which the fault passes is observed. With the consideration of contact, a block falls or separates from the surrounding rock and a discontinuous displacement field is observed near the interface of the fault and the rock. The maximum displacement when contact is considered approximately 35 mm larger than without contact. In the upper contact interface of block B4, the displacement changes substantially, which further indicates the loosening of block B4.Figure 13
Deformation of the surrounding rock. (a) Without contact and (b) considering contact.
(a)
(b)When contact between the blocks and the surrounding rock is considered, the simulation results are more consistent with the learned empirical rules. Therefore, this numerical analysis method better accords with the actual situation.


### 6.3.3. Stress Evolution of the Interfaces between the Blocks and Rock

To analyze the stress evolution process of the interfaces between the blocks and rock, block B2 is selected as the research object. Three monitored points, namely, F1, F2, and F3, are arranged in the middles of the interfaces, as illustrated in Figure10. The normal stress and shear stress evolutions at the three points are plotted in Figure 14. During excavation, the interfaces are relieved of the initial stress. The stresses of interfaces F2 and F3 gradually increase as the excavation progresses, while the stress of interface F1 decreases. Nonsliding face F3 has the largest stress with block B2. Hence, the block is subjected to multiple stresses from the surrounding rock during its failure process, in addition to the stress from the sliding face.Figure 14
Stress evolutions of the interfaces. (a) Norm stress and (b) shear stress.
(a)
(b)


### 6.3.4. Safety Factors of the Blocks

During the calculation of each excavation step, the safety factors of the blocks are obtained via reduction of the shear strength parameters of the interfaces. In the 1st stage of the excavation, the displacement of block B1 is large and continuously growing, even if the strength reduction factor is less than 0.1. Block B1 falls freely from the top of the cavern. Block B2 is fully exposed after the 1st stage of the excavation and blocks B3 and B4 are fully exposed after the 5th stage of the excavation. Figure15 plots the displacements of the blocks as functions of the strength reduction factor during the 5th stage of the excavation. For block B2, when the reduction factor exceeds 4.51, the monitored displacement suddenly increases; hence, the safety factor of block B2 is 4.51. Similarly, the safety factors of blocks B3 and B4 are 1.45 and 1.49, respectively.Figure 15
Safety factors of three blocks after the 5th stage excavation.The safety factors of the blocks during each stage of the excavation are plotted in Figure16. The changes in the safety factors vary among the stages. As the excavation progresses, the safety factor of block B2, which is in the top arch, gradually increases, while the safety factors of the blocks in the sidewall gradually decrease. However, the magnitudes of the changes among the stages are small. The block stability is greatly affected by geostress. Although the excavation has adjusted the stress of the surrounding rock, its influence degree is not as great as that of the geostress. So, block depth and excavation process are two important factors affecting the safety degree of blocks, and the influences of the excavation process to the safety factors of deep buried blocks are relatively small.Figure 16
Variations of the safety factors of three blocks.According to the safety factors of the blocks, in addition to strengthening the support of key blocks B1 and B3, it is also necessary to focus on potentially dangerous blocks B2 and B4. The potentially dangerous block in the left sidewall is more dangerous than that in the top arch.


## 6.3.1. Failure Process of Blocks during Excavation

A monitoring scheme, which is illustrated in Figure10, is used to monitor the deformation process of the blocks. Figure 11 shows the deformation of the surrounding rock after the 2nd, 4th, 6th, and 8th stages of excavation of the powerhouse (the displacement is magnified by a factor of 15). Figure 12 shows the monitored displacement. The monitored displacement gradually increases as the excavation progresses. The point P1 in block B1 undergoes a large displacement and continues to increase after the first stage of excavation. Hence, block B1 falls freely during excavation. The displacement of point P2 in block B2 is 14 mm after the first stage due to the excavation unloading. It slowly increases as the excavation proceeds and stabilizes at approximately 30 mm. Blocks B3 and B4 are completely exposed after the 5th stage of excavation and their displacements are increased by 22 mm and 18 mm, respectively. According to Figures 9(b)–9(d), blocks B3 and B4 undergo large deformations and cracks are generated between the blocks and the rock. The depth and width of the cracks gradually increase as the excavation proceeds. Points P4 and P5 are a contact point pair prior to excavation. After the 5th excavation stage, they have a relative displacement of approximately 10 mm and a crack occurs in their interface. The relative displacement eventually reaches approximately 25 mm. Hence, the interface deforms discontinuously. After the excavation has been completed, except for block B1 falling off, the blocks does not suffer from instability or damage. However, the blocks have loosened from the surrounding rock and there are high security risks.Figure 10
Monitoring schemes.Figure 11
Deformation process during excavation (the displacement is amplified by a factor of 15): (a) 2nd stage; (b) 4th stage; (c) 6th stage; (d) 8th stage.
(a)
(b)
(c)
(d)Figure 12
Monitored displacement.If using key block theory, only blocks B1 and B3, which are exposed on the excavation face, can be identified and it cannot simulate the failure process of the blocks. The proposed numerical method can effectively identify the key blocks, timely discover potentially dangerous blocks, and simulate the failure process of the blocks.


## 6.3.2. Deformation Characteristics of the Blocks and the Surrounding Rock

To reflect the deformation law of blocks and surrounding rock, an additional calculation in which contact is not considered is conducted for comparison with the result in which contact is considered, as shown in Figure13. The case of without contact means the blocks and the surrounding rock are assumed to be tied together during modeling, and the common nodes on the interfaces are not separated. Without consideration of contact, the failure of the blocks cannot be simulated. Only a large displacement of the rock through which the fault passes is observed. With the consideration of contact, a block falls or separates from the surrounding rock and a discontinuous displacement field is observed near the interface of the fault and the rock. The maximum displacement when contact is considered approximately 35 mm larger than without contact. In the upper contact interface of block B4, the displacement changes substantially, which further indicates the loosening of block B4.Figure 13
Deformation of the surrounding rock. (a) Without contact and (b) considering contact.
(a)
(b)When contact between the blocks and the surrounding rock is considered, the simulation results are more consistent with the learned empirical rules. Therefore, this numerical analysis method better accords with the actual situation.


## 6.3.3. Stress Evolution of the Interfaces between the Blocks and Rock

To analyze the stress evolution process of the interfaces between the blocks and rock, block B2 is selected as the research object. Three monitored points, namely, F1, F2, and F3, are arranged in the middles of the interfaces, as illustrated in Figure10. The normal stress and shear stress evolutions at the three points are plotted in Figure 14. During excavation, the interfaces are relieved of the initial stress. The stresses of interfaces F2 and F3 gradually increase as the excavation progresses, while the stress of interface F1 decreases. Nonsliding face F3 has the largest stress with block B2. Hence, the block is subjected to multiple stresses from the surrounding rock during its failure process, in addition to the stress from the sliding face.Figure 14
Stress evolutions of the interfaces. (a) Norm stress and (b) shear stress.
(a)
(b)


## 6.3.4. Safety Factors of the Blocks

During the calculation of each excavation step, the safety factors of the blocks are obtained via reduction of the shear strength parameters of the interfaces. In the 1st stage of the excavation, the displacement of block B1 is large and continuously growing, even if the strength reduction factor is less than 0.1. Block B1 falls freely from the top of the cavern. Block B2 is fully exposed after the 1st stage of the excavation and blocks B3 and B4 are fully exposed after the 5th stage of the excavation. Figure15 plots the displacements of the blocks as functions of the strength reduction factor during the 5th stage of the excavation. For block B2, when the reduction factor exceeds 4.51, the monitored displacement suddenly increases; hence, the safety factor of block B2 is 4.51. Similarly, the safety factors of blocks B3 and B4 are 1.45 and 1.49, respectively.Figure 15
Safety factors of three blocks after the 5th stage excavation.The safety factors of the blocks during each stage of the excavation are plotted in Figure16. The changes in the safety factors vary among the stages. As the excavation progresses, the safety factor of block B2, which is in the top arch, gradually increases, while the safety factors of the blocks in the sidewall gradually decrease. However, the magnitudes of the changes among the stages are small. The block stability is greatly affected by geostress. Although the excavation has adjusted the stress of the surrounding rock, its influence degree is not as great as that of the geostress. So, block depth and excavation process are two important factors affecting the safety degree of blocks, and the influences of the excavation process to the safety factors of deep buried blocks are relatively small.Figure 16
Variations of the safety factors of three blocks.According to the safety factors of the blocks, in addition to strengthening the support of key blocks B1 and B3, it is also necessary to focus on potentially dangerous blocks B2 and B4. The potentially dangerous block in the left sidewall is more dangerous than that in the top arch.


## 7. Conclusions

This paper proposed a numerical analysis method for block failure that is based on continuum mechanics. The following conclusions can be drawn from this study:(1)
It considered the contact action between the blocks and the surrounding rock and combined the advantages of the discontinuum method and the continuum method. It can simulate the continuous deformation characteristics of the surrounding rock and the discontinuous failure processes of the blocks. A simple siding block example demonstrates that the simulation results are consistent with the results that were obtained via RBLEM under specified conditions. The main advantages of this numerical method are that it can effectively consider the influence of excavation unloading on the stability of the blocks and simulate the bonded, sliding, and separation behaviors of the interfaces between the blocks and the rock.(2)
The numerical example of the underground powerhouse demonstrates that the key block is damaged after excavation, the potentially dangerous blocks undergo large deformations, and the safety factor is low in the local area, which increases the safety risk during excavation. The numerical results also demonstrate the influence of the surrounding rock to the stability of the blocks, and the discontinuous deformation characteristics and various failure modes of the complex blocks. The basic law is consistent with the real scenario. This numerical method provides an effective analysis method for the stability analysis of blocks under complex geological conditions.



---

*Source: 4280428-2020-01-13.xml*