A General Block Stability Analysis Algorithm for Arbitrary Block Shapes

In rock engineering, block theory is a fundamental theory that aims to analyze the finiteness, removability, and mechanical stability of convex blocks under different engineering conditions. In practice, the possible combinations of the fractures and joint sets that may generate key blocks can be identified by stereographic projection graphs of block theory. However, classic key block theory does not provide solutions for nonconvex blocks, which are very common in civil projects, such as those with underground edges, corners, and portals. To enhance the availability of block theory, a general algorithm that can analyze the removability and stability of blocks of arbitrary shapes is proposed in this paper. In the proposed algorithm, the joint pyramid for blocks of arbitrary shapes can be computed, and the faces of the blocks are grouped according to their normal vectors such that parallel or nonadjacent sliding faces with the same normal vector can be immediately identified when the sliding mode is determined. With this algorithm, blocks of arbitrary shapes can be analyzed, and users do not need to have experience interpreting graphs of block theory to take advantage of its accuracy and effectiveness. The proposed algorithm was verified by several benchmarking examples, and it was further applied to investigate the stability of the left bank rock slope of a dam. The results showed that the proposed algorithm is correct, effective, and feasible for use in the design and support of excavation in complex rock masses.


INTRODUCTION
Block theory Shi (1977), Shi (1982) is an important analysis method in computational rock mechanics and has been applied in many rock engineering projects in the past 30 years (Zhang et al., 2014;Zheng et al., 2014;Sun et al., 2015). The core purpose of classic block theory is to analyze the finiteness, removability, and mechanical stability of convex blocks Goodman and Shi (1985),  under different engineering conditions. There are three assumptions in block theory. First, it is assumed that all fractures are infinite planes. Second, all the blocks are assumed to be rigid bodies. Third, it is assumed that blocks move in only translation and do not rotate. Under these assumptions, block theory judges finiteness, removability and stability of blocks by considering joint planes and space planes that the faces of the convex blocks are located. The mechanical stability of blocks is computed based on the removability, frictional angles of joint planes and resultant force applied on the block. Then, suggestions for the design and construction of rock engineering are proposed based on the analysis results. The joint pyramid (JP) is the core concept of block theory; it is a convex pyramid intersected by inward half spaces that faces of blocks are located. The evacuation pyramid (EP) is similar to the JP and corresponds to the intersection of the evacuation half spaces. A block pyramid (BP) corresponds to the intersection of joint pyramid (JP) and evacuation pyramid (EP). As shown in Figure 1A, four key blocks are generated considering the intersection of orange and green joint sets and the tunnel boundary. In Figure 1B, joint pyramid codes of the blocks at the surface of a tunnel are shown; the label "01" means above an orange joint plane and below a green joint plane. The removability of blocks can be analyzed by the joint pyramid of convex blocks.
With block theory, key blocks of different joint pyramids can be estimated before excavation, and the shape and volume of the maximum key block can be analyzed. The concept of a maximum block means that the key blocks belonging to the same joint pyramid cannot be removed if they have larger volumes. The maximum key block must correspond to the largest supporting force; thus, block theory provides conservative and safe support suggestions, and it is a useful and efficient tool in engineering. As shown in Figure 2, block theory is used to analyze an underground cavern. There are three main sets of joints in the cavern region whose geometrical parameters (dip direction ∠dip) are 0°∠50°, 120°∠60°and 240°∠70°. And their friction angles are 20°, 30°and 40°, respectively. Different combinations of joint sets form different joint pyramids and key blocks with different excavation surfaces. Block theory requires only a few joint planes; thus, the orientation and location of planes can be directly measured by geological compasses or other equipment on site (Gonzalez-Palacio et al., 2005). Subsequently, the volumes, weight and shapes of the maximum key blocks can be obtained. In Figure 2, the solid circles denote the projections of joints, and the label "101" is the joint pyramid code, which means below joint sets 1 and 3 and above joint set 2. The label "+0.58" (below the joint pyramid code) corresponds to the residual sliding force on a unit volume of the corresponding joint pyramid "101". The corresponding maximum key blocks for joint pyramids "010", "111", and "101" are shown in Figure 2.
In practice, the orientations and locations of fractures and joint sets can be obtained by field survey, and blocks can be generated in two ways: 1) Obtain intersection lines by intersecting the fractures and joint disks. Then, intersection lines are intersected to obtain all the intersection points and final intersection lines. Then, search all the closed loops and save them as polygons. Finally, search the polyhedrons by determining the correlation between the intersection lines and searched polygons (Zhang and Wu, 2005;Zhang and Wu, 2007). 2) Extend all fractures and joint sets into infinite planes. Then, perform the same procedure as in the first approach to find all the polyhedrons. Then, return the infinite planes back to their original sizes. Finally, polyhedrons that share the same imaginary faces are combined (Yu et al., 2005). These two methods are commonly used to obtain approximate blocks, but the second method requires more computation. They are both based on computational geometry and require a discrete fracture network (DFN) model to be established first. Because the orientations of joints are complex and difficult to obtain, stochastic approaches are usually adopted to generate joints (Mauldon, 1993;Ahn and Lee, 2004;Dowd et al., 2009). After block generation, the convex blocks can be analyzed by classic block theory to obtain more accurate results. To make block theory more applicable, researchers have focused on nonconvex blocks Shi and Goodman (1989), Shi (2006), Menéndez-Díaz et al., (2009), Elmouttie et al. (2010 and more complex application scenarios (Fu et al., 2016;Zheng et al., 2016;Wu et al., 2018). Compared to the much more complex stability analysis He et al. (2018), Zhang and Shi (2021) of discontinuous deformation analysis Shi (1992), block theory is a static method that has three features: 1) the removability of blocks is analyzed by geometrical and topological methods; 2) stability is determined by simple mechanical analysis; and 3) only translation movement at the initial moment is considered. Thus, block theory requires much less computational cost and is suitable to find key blocks and provide supporting advice during the survey and construction period of civil engineering projects. To facilitate the utilization of block theory Cheng et al. (2017), threedimensional visualization for stereographic projection  and equilibrium region graphs Xue and Song (2017) have been implemented. However, a general algorithm with which to compute the stability of blocks of arbitrary shapes is still needed.
We focus on finite blocks in this paper, and a general algorithm to analyze the stability of blocks of arbitrary shapes is proposed, as described in Proposed Algorithm. Joint pyramids for blocks of arbitrary shapes can computed by proposed algorithm that based on classical block theory, by which the stability of blocks can be analyzed. Validations and engineering applications are given in Validations and Engineering Applications. Conclusions are given in Conclusion. By means of the method proposed in this paper, dangerous key blocks of arbitrary shapes, such as those at the edges, corners, and portals of underground chambers, can be identified by engineers immediately.

PROPOSED ALGORITHM
All blocks formed by fractures of a finite size in a given rock mass domain can be identified using 3D block-cutting analysis. Among these blocks, there may exist blocks that have parallel faces or concave angles formed by fracture surfaces. However, a general key block judgment algorithm for blocks of arbitrary shapes is seldom studied. To make block theory suitable for blocks of any shape, this paper proposes a general algorithm for analyzing block stability. It simplifies faces of arbitrary shapes into joint planes and computes a joint pyramid for the joint planes, even those with concave angles. The algorithm can utilize block theory to create industrial products that everyone can use.
The algorithm has been designed and implemented as a library, whose architecture is shown as Figure 3. The input module reads blocks generated by the 3D cutting process. Then, the analysis module invokes the classes "BlockFacesTagger", "JointPyramidAnalyzer" and "BlockStabilityAnalyzer" to analyze the stability of each block. Finally, the output module saves and displays the analysis result. The analysis takes the following steps: 1) "BlockFacesTagger" classifies and tag faces of a block and give each face an id according to the face normal vectors (see Grouping and Tagging Block Faces). 2) "JointPyramidAnalyzer" computes the joint pyramid based on the normal vectors from the first step to build the relation Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 723320 between the faces of the block and the joint pyramid (see Computing the Joint Pyramid). 3) "BlockStabilityAnalyzer" computes block stability according to the joint pyramid, e.g., sliding mode (see Analyze Stability of Joint Pyramid).

Grouping and Tagging Block Faces
To analyze the stability of blocks, engineers must simplify the faces of blocks into joint planes first, which is tedious and complex work. To make block theory easier to utilize, the proposed algorithm groups and tags the faces of the block by the face normal vectors and generates joint planes. It has the following steps: 1) Prepare a normal array N and normal id array NI.
2) Iterate through the faces of the block, and check whether each normal already exists in N. 3) If the normal does not exist, identify the current face normal with a size of N and add a new normal into N; otherwise, just use the id of the existing normal.
The key steps are described in Algorithm 1.

Computing the Joint Pyramid
According to block theory, "a convex block is removable if and only if its block pyramid is empty and joint pyramid is not empty" (Goodman and Shi, 1985). Because blocks in the real world are all finite, the block pyramid is always empty; see Eq. 1. Therefore, we only need to consider the joint pyramid, and then the theorem can be simplified as "finite convex blocks are removable if and only if its joint pyramid is not empty".
In classic block theory, joint pyramids are only used to analyze convex blocks. However, blocks have various shapes in practice, and we must find a general way to compute the joint pyramid. "A finite nonconvex block is a composite block that consists of some finite convex blocks. If all its convex subblocks are removable, the nonconvex block is removable" (Li et al., 2012). If a block is removable, the movement direction must be in the joint pyramid of this block, and the joint pyramid is not empty. Thus, we need to compute the intersection of joint pyramids of all the convex subblocks only. Because the number of joint planes of these joint pyramids is finite, we compute the intersection of all half spaces of the joint planes of this block.
Therefore, a nonconvex block is removable if and only if its joint pyramid is not empty. The joint pyramid of a nonconvex block is the minimum pyramid that belongs to any other combination of half spaces of joint planes of the block. Thus, we need to obtain the orientations of the joint planes first, i.e., dip angles and dip direction angles.
There are four steps to compute the joint pyramid for any number of joint planes, even those with concave edges: 1) Determine all the inward normal vectors of the joint planes. 2) Compute the cross-product for every two normal vectors and save the cross-product vectors if the dot-product between itself and any other normal is positive, then do the same procedure for the reverse vector of the cross-product vector. 3) Simultaneously reserve the face pairs of all the cross-product vectors saved in step 2. 4) Iterate all face pairs, sort them into end-to-end order, and reverse any face pairs if needed.
Algorithm 2 provides a detailed explanation of this computation.

Analyze Stability of Joint Pyramid
According to block theory, a block is in a stable mode if its joint pyramid is empty, which means that there is no space available to move into. Otherwise, there are three types of translational failure modes: lifting, single-face sliding and double-face sliding. If a block is in lifting mode, the motion direction d → is parallel to the direction of resultant force r → : where d i → is the orthographic projection of r → on face i.
If a block slides along one face of its joint pyramid, i.e., sliding along one face or parallel faces of the block. If the index of the sliding face is i, then its translational direction d → is: where n i is normal vector of face i. For the single-face sliding mode, Eqs. 4, 5must be satisfied: If a block slides along two nonparallel faces i and j of a joint pyramid simultaneously, its sliding direction is: Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 723320 d → n i × n j n i × n j sign n i × n j · r → where sign[x] is equal to -1, 0 or 1 when x is < 0, 0, or >0, respectively. Thus, Eqs. 7-9 must be satisfied at the same time: To compute the stability of the joint pyramid, all the sliding vectors need to be computed to find the final sliding mode. This requires the following steps, as shown in Figure 4: 1) If the joint pyramid is empty, the block is in a stable mode; go to step 7; otherwise, go to the next step. 2) Check if the resultant force is in the joint pyramid. If it is in, the block is in lifting mode, and go to step 7; otherwise, go to the next step. 3) Compute all the single-face sliding vectors and double-face sliding vectors. 4) Filter all the sliding vectors that are out of the joint pyramid. 5) Find the sliding vector that has the smallest angle with the resultant force, which corresponds to the final sliding direction, from the filtered vectors. 6) Judge the sliding mode and compute the sliding force for the final sliding direction. 7) Process ends.

VALIDATIONS
In this section, the correctness and feasibility of the proposed algorithm are verified by block with parallel faces and concave angles. This block is on the edge of a tunnel, and its faces are tagged with f 1 , f 2 , f 3 , f 4 , f 5 , f 6 , and f 7 , and the normal ids of the faces are tagged with 1, 2, 3, 4, 4, 5, and 6. Among these faces, f 4 and f 5 are parallel, while f 1 and f 3 form a concave edge . Faces f 6 and f 7 of the block are located on the boundary of tunnel, so they act as space planes, and all the other faces are used as joint planes, as shown in Table 1. Under gravity force, i.e. (0, 0, −1), the concave block is in the double-face sliding mode along faces f 3 , f 4 and f 5 , and the sliding force of unit weight is 0.806, as shown in Figure 5.

ENGINEERING APPLICATIONS
The proposed algorithm is used to study the stability of the left bank of a dam in China. The left bank slope is excavated on a large scale, and the strike of the excavated slope changes significantly, showing a different angle with the strike of the rock stratum, as shown in Figure 6. According to the direction of the excavation slope and the different engineering positions, the left bank slope is divided into four slope sections: 1) the upstream slope of the dam abutment, 2) the front and the downstream slope of the dam abutment, 3) the tailrace slope, and 4) the downstream slope of the tailrace canal.  The average occurrence of the strata is 268∠39°(dip direction ∠dip), with some interlayer shear zones developed. There are three main sets of joints developed in the rock mass. The joints are straight and closed, with a calcite film attached. The geometric and mechanical features of the discontinuities considered in the analysis are shown in Table 2.
When adopting block theory, it is necessary to analyze only the blocks formed by the combination of the strata and the joints. The blocks formed by a combination of joints are omitted due to their small volume, and the stability of this kind of block can be ensured by a general slope hanging net and anchor rod support system. Generally, gravity is considered the only force applied to the blocks.   Using full-space stereographic projection analysis, the removable joint pyramids and their sliding modes of different slope sections are obtained. See Table 3. Among them, the stereographic projection of the slope near section 1 is shown in Figure 7. 1) For the same slope section, if the result of the stereographic projection analysis shows that there are two or more types of   Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 723320 7 removable blocks, then only blocks with a poor stability, large scale, or high support requirement need to be analyzed in the block stability and engineering support analysis. When the stability of this type of block is ensured, other types of removable blocks are stabilized. 2) When only the effect of gravity is considered, the double-face sliding mode has a smaller safety factor if both faces are upward. Therefore, if there are joint pyramids "00" and "01" at the same slope section, only the blocks formed by joint pyramid "00" need to be analyzed.
The strata on the upstream slope of the dam abutment incline outward, and the strata direction (strike 178°, dip direction 268°) is close to the slope direction (strike 162°, dip direction 252°). According to the qualitative analysis of engineering geology, the slope is a near-forward slope, and the blocks exhibit a single-face sliding mode along the strata. According to the analysis of the full-space stereographic projection, when combining interlayer shear zone II1-7005 with the fractures of the NWW group, block 1 in Figure 8 can be formed, and the block is in a single-face sliding mode along the interlayer shear zone. The block traces the NWW group and forms a shear zone. Since zone II1-7005 is exposed at the ∇173 m platform, the block with the largest body is formed when the shear surface passes through the intersection of zone II1-7005 and ∇173 m. However, the block may be unremovable due to the shape restriction. In fact, it is more reasonable to cut the rock at a certain position on the inner side of the slope to form a shear surface, i.e., Block 2 in Figure 8 and Block 3 in Figure 9.
The shear surface on the inner side of the slope follows other cracks or shear surfaces that cut through the rock mass. With the different orientations of the inner shear surfaces of the slope, the shape and stability of the blocks are different. Thus, to determine the final support design, the orientation of the inner shear surface of the slope with the smallest safety factor must be identified. Within the possible variation range of the dip direction, the shape characteristics of blocks are shown in Figure 10, and the stability is shown in Table 4. The analysis results show that when the shear surface orientation is 350°∠85°, the safety factors of the blocks in Figure 10E, F and have the same value of 0.222, which is the smallest value calculated in this case, and the sliding mode of the blocks is single-face sliding. However, the block in Figure 10F is more dangerous because of its larger volume.
The conclusion obtained by the analysis of block theory is similar to the qualitative analysis of engineering geology. The block undergoes single-face sliding with the interlayer shear zone  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 723320 as the sliding face, but the block theory clearly identifies the shear surface through the analysis of block morphology. When the orientation of the shear surface is different, the sliding mode of the block can transition from double-face sliding to single-face sliding. Accordingly, the safety factor of the block changes. The spatial position of the block can be analyzed after determining the orientation of the shear surface. In different positions, the size and shape of the blocks are different, and the stability is different as well. Therefore, it is necessary to compute the safety factors of blocks of different sizes and shapes within the possible range of block sizes to fully describe the stability characteristics of this type of block. The analysis results show that when a block is as large as possible, its safety factor is smaller; the largest possible block is shown in Figure 11. The volume of this block is 45,808 m 3 .

CONCLUSION
In this paper, we propose a general block stability analysis algorithm based on the traditional block theory analysis framework. It has 2 advantages: first, parallel or nonadjacent sliding faces with the same normal vector can be identified immediately, which may provide more detailed support advice for engineers; second, the joint pyramids and removability of blocks of arbitrary shapes can be computed, which greatly enhances the applicability of block theory. The feasibility of the proposed algorithm in treating blocks with parallel faces and concave angles is verified first. Then, the proposed algorithm is used to locate the key block on the left bank slope of a dam in China, and the final key block identified using the proposed algorithm is similar to that identified via a qualitative engineering geology analysis. Using this proposed algorithm, engineers can analyze block stability and find key blocks of arbitrary shapes at any position without the additional cost of interpreting graphs of block theory, which could greatly promote the development and application of block theory.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
XC developed the algorithm and wrote the manuscript; LL performed the experiments and participated in writing the manuscript; JX contributed significantly to analysis and manuscript preparation; QZ performed the engineering application; JX helped perform the analysis with constructive discussions. YW contributed to the conception of the study.