Abstract
SRM (strength reduction method) with iterative method of PFE (partitioned fnite element) and IBE (interface boundary element) is proposed to solve the safety factor of local discontinuous rock mass. Slope system is divided into several continuous bodies and local discontinuous interface boundaries. Each block is treated as a partition of the system and contacted by discontinuous joints. The displacements of blocks are chosen as basic variables and the rigid displacements in the centroid of blocks are chosen as motion variables. The contact forces on interface boundaries and the rigid displacements to the centroid of each body are chosen as mixed variables and solved iteratively using the interface boundary equations. Flexibility matrix is formed through PFE according to the contact states of nodal pairs and spring flexibility is used to reflect the influence of weak structural plane so that nonlinear iteration is only limited to the possible contact region. With cohesion and friction coefficient reduced gradually, the states of all nodal pairs at the open or slip state for the first time are regarded as failure criterion, which can decrease the effect of subjectivity in determining safety factor. Examples are used to verify the validity of the proposed method.
1. Introduction
SRM [1] is one of frequently used methods to solve safety factor. There are two important factors in stability analysis of rock slope with SRM, that is, the validity of numerical simulation for local discontinuous joints and the correctness of failure criteria in determining safety factor. Several methods are used to simulate local discontinuous problem, such as FEM (finite element method) [2, 3], DEM (discrete element method) [4, 5], and DDA (discontinuous deformation analysis) [6–9], but there are more or less disadvantages in computational efficiency. Although there are some reports about stability analysis by using SRM based on these methods, the failure criteria need to be improved [10–12].
TLEM (thin layer element method of FEM) [13] belongs to continuum method. In this method, layer elements with elastic-plastic constitutive are used to simulate weak structural planes and plastic yield is just along the layer elements. For compatibility requirement needed in FEM, TLEM is unable to simulate the open state of discontinuous joint effectively. Moreover, the yield criterion chosen is empirical and different safety factors are always gotten using different yield criteria.
DEM and DDA belong to noncontinuum method. In DEM, stability analysis of slope is considered by a number of explicit contact equilibrium equations, and contact forces are adjusted continuously to keep them balanced, but the choice of time step is strictly limited. In DDA, the contact states can be taken into the total stiffness matrix and the equilibrium equations built up are implicit. The uniqueness and convergence of results can be guaranteed [14]. Although the time step of DDA can be assigned larger than DEM and will not cause the problem of numerical convergence [15], this method belongs to the whole iteration scheme and will cause computational inefficiency for rock mass with few weak structural planes. Due to dynamic relaxation method conducted in DEM and DDA, calculation without convergence at the last several steps is always taken as a failure criterion in determining safety factor with SRM based on DEM or DDA [16]. Besides, the remarkable growth of displacement is also used as a failure criterion [17]. But the effect of subjectivity in a set of convergence criteria or the choice of inflection point of displacement is unavoidable.
In this paper, SRM based on iterative method of PFE and IBE is proposed for stability analysis of local discontinuous rock slope. The displacements of blocks are chosen as basic variables. The contact forces and the rigid displacements to the centroid of each body are chosen as mixed variables and solved iteratively using the interface boundary equations. The displacements of blocks can be solved by FEM after the contact forces obtained and added in the total load vector. Flexibility matrix is formed through PFE according to contact state and spring flexibility is used to reflect the influence of weak structural plane so that nonlinear iteration is only limited to the possible contact region. With cohesion and friction coefficient reduced gradually, the states of all nodal pairs at the open or slip state for the first time are regarded as a failure criterion in determining safety factor of rock slope.
2. Iterative Method of Partitioned Finite Element and Interface Boundary Element
2.1. Mechanical Model
In Figure 1, a rock slope is composed of slip mass , bed rock , and the discontinuous joint with the gap of . and are the local coordinate system and the global coordinate system, respectively. The system of forces acting on blocks is divided into two parts: external force and contact internal force . and are the contact forces on two nodes of each nodal pair and are the three stress components in the local coordinate system, assuming tension is positive in all forces.
If the contact force increment is given, the stress and deformation of and can be solved by FEM. In fact, the contact force increment and the incremental displacement on interface boundary are unknown, but there is a nonlinear relationship between them. The nonlinear relationship can be described as follows.(1)The total incremental displacement at point on interface boundary will be changed by the effect of contact force increment . includes the incremental displacement of deformable body and the rigid incremental displacement of rigid body ; that is,(2)Contact state will be changed by the difference of total incremental displacement on interface boundary, and the updated state will make contact force change further.(3)Total incremental displacement and contact force increment can be solved iteratively by using the interface boundary element equations.
In numerical calculation, three states are defined in Table 1, that is, open, close, and slip. Mutual transition among states is listed in Table 2 and iterative operation is executed sequentially by a serial number. Two principles must be fulfilled, that is, the principle of normal impenetrability and the principle of normal pressure.
2.2. Partitioned Finite Element Equations
For contact bodies and , if the analysis of th step has been finished and the contact force increment been solved, the static incremental equilibrium equation at step without considering rigid motion is stated as follows:where is the global stiffness matrix; represents the vector of incremental displacement at step ; is the vector of total external load at step ; stands for the external load at step ; is the vector of contact internal force at step ; is the stress matrix.
is introduced to (2); the following equation can be gotten: where ; an arbitrary component in the matrix represents the flexibility coefficient corresponding to the displacement at the freedom due to a unit force at the freedom . If the nodes are located on interface boundaries, (3) is the deformable displacement increment of boundary nodes of blocks.
2.3. Interface Boundary Element Equations
The centroid displacement of the th block is composed of translation components and rotation components. Assuming the centroid displacement increment of th block is , the rigid incremental displacement of the th node with a relative coordinate value of to centroidin th block is given as follows:where is the transform matrix of th node to the centroid of th block.
Equations (3) and (4) plugged into (1), the total incremental displacement of node on any nodal pair of contact blocks is
Equation (5) can be written in another form; that is, where is the matrix relevant to node in the total flexibility matrix; stands for the deformable incremental displacement only due to the external load increment; it has nothing to do with the current contact state and can be obtained by back substitution directly. is the total displacement increment at current iterative step induced by both the external load increment and the contact force increment. It can be solved iteratively by the contact forces increment and the rigid displacement increment of th block.
For each contact body, the equilibrium equation related to centroid point at step for all forces acted on it could be written as follows: where ; , are mass and rotational inertia related to the centroid of th block, respectively; represents the mass-damping coefficient; stand for the rigid acceleration and rigid velocity of th block, respectively; is the number of contact nodes on th block; is the number of total nodes on th block; is the contact force of node ; is the external load acted on the th node of th block.
For static problem, left side of (7) is zero, and the incremental formula is
Combing (6) and (8), the static iterative equations on discontinuous interface boundaries of th block can be gotten as follows:where subscript is the interface boundary and the dimension of is ; subscript stands for the whole region of block and the dimension of is ; dim is the dimension of space.
2.4. Influence of Weak Structural Plane
Fractured rock mass is cut by weak structural planes and the influence of weak structural planes can be reflected by adding spring stiffness on nodal pairs without establishing solid elements. If the thickness in -direction is far less than the size in and directions, three stress components, , , and , can be used to describe nodal pair. The incremental constitutive law can be described as follows: where ; is the stiffness matrix of weak structural plane.
Assuming is the elastic module of continuous body, if a pair is at the close or slip state, the properties of weak structural plane are controlled by shear strength and the deformation of it can be ignored on the condition that the thickness is very thin (less than 1~3 mm). So a big normal stiffness , such as (500~1000), will be given. If the elastic module and thickness of weak structural plane are given, real normal stiffness can be calculated by according to the principle of equivalent deformation under uniaxial compression [18]. If a pair is at the open state, a small , such as 5 kPa, will be given to simulate the infinite open state.
If a pair is at the close state, the tangent stiffness of the pair is given as follows [19]:
If a pair is at the slip or open state, considering residual strength of weak structural plane and taking tangential stress-strain hyperbolic model for reference, tangent stiffness can be calculated in the following form:where is damage ratio; a value close to 1.00 is suggested.
Equation (10) multiplies at both sides. is the contact force increment at point ; represents the relative incremental displacement at nodal pair . Therefore, (13) can be gotten; that is,
Assuming represent the influence flexibility matrix of nodal pair on weak structural plane, can be assembled by the nodal pairs around block and add it in (9); the static iterative equations on discontinuous interface boundaries considering the influence of weak structural planes arewhere is the total incremental displacement on interface boundaries of block , including three parts: the incremental displacement of deformable body, the rigid incremental displacement of rigid body, and the incremental deformation of weak structural plane. Taking all blocks and all discontinuous joints into account, the total static interface boundary equations can be assembled as follows:where
2.5. Numerical Implementation
Combining the principles and equations above, the numerical implementation of the iterative method can be summarized as follows (the superscript represents the contact iteration step while the subscript stands for the loading increment step).
Step 1. Assemble and decompose the total stiffness matrix , solve the flexibility matrix on possible contact interface boundaries by the unit force method, form the initial influence flexibility matrix on weak structural planes, and get the transform matrix on each node.
Step 2. If the load increment at step has been finished, get into step and update the contact states and the contact forces of nodal pairs, which will be set as the initial condition for the iteration of contact state.
Step 3 (the iteration of contact state). Obtain deformable incremental displacement by back substitution of the equation , add , rigid incremental displacement , and the incremental deformation of weak structural planes in total incremental displacement , form and assemble vector , and assemble and decompose matrix in (15).
Step 4 (the iteration of contact force and rigid displacement). If the iteration of contact force and rigid displacement at step has been finished, get into step , solve vector by using (15), and judge the convergence of iteration by using .
If iteration is not convergent, add contact forces in the whole load vector and solve the equations of FEM, back substitution of deformable incremental displacement and achieve the rigid incremental displacement on interface boundaries, add them to , get the external load increment on each body, and assemble vector in (15).
If iteration is convergent, the mutual transition of the three states will be carried out as follows according to Table 2.(1)If the current state is at the open state, new gap can be obtained by using . If is less than zero, the current state is updated to the close state and is set to zero according to the principle of normal impenetrability.(2)If the current state is at the close or slip state, solve the new normal stress by using . If is less than , calculate . If is greater than , the current state is changed to the open state and contact internal force is set to zero according to the principle of normal stress pressure.(3)After Step (1) and Step (2), if the current state is still at the close or slip state, calculate tangential force by using . If is satisfied, the current state is updated to the slip state, and a new condition of is set. Mark the number of freedoms on condition that three freedoms exist in the close state and one in the slip state. Update influence flexibility matrix on weak structural planes.
Step 5 (convergence judgment for contact state). If the contact state of any nodal pair at current step is not the same as last, the contact state is not convergent. Update new contact state, assemble and decompose new flexibility matrix , and go to Step 3.
Step 6. If the contact state is convergent, go to Step 2 and get into the next load increment until all load increment steps are completed.
3. Strength Reduction Method
It is not the safety factor of slope, but stress, deformation, the rigid displacements of blocks, and the contact states of nodal pairs can be obtained directly by iterative method of PFE and IBE. Duncan [20] pointed out that the safety factor of slope can be defined as the reduction degree on shear strength of soil when the slope is at the critical failure state. If are divided by reduction coefficient , new parameters of shear strength can be gotten. Calculating by using (17) repeatedly until the critical failure state, the reduction coefficient is safety factor:
As shown in Figure 2, a slope system is divided into slip mass, slide bed, and the interface boundary, which is simulated by 5 nodal pairs. The slip mass and the slide bed are considered as continuous elastic blocks and nonlinear iteration is only conducted on joint. Density g/cm3, elastic module Pa, Poisson’s ratio , vertical uniform load kN/m, horizontal uniform load kN/m, cohesion MPa, and friction coefficient . The bottom of the mesh is fixed and the horizontal displacements on both sides of the slide bed are constrained. Only considering the strength properties of the discontinuous joint, normal contact stiffness and damage ratio are given and 0.992, respectively. If any pair is at the open or slip state, the nodal pair is failure.
Cohesion and friction coefficient reduced gradually at the same level; -direction displacement of point A and the number of failure pairs are shown in Table 3. Curve of reduction coefficient and -direction displacement of point A is shown in Figure 3. It can be seen that the number of the failure nodal pairs is increased at the reduction coefficients 0.65 and 1.00, and inflection points appear two times. The conventional failure critical to the inflection point of displacement is not effective to determine the safety factor. But all nodal pairs are failure for the first time at the reduction coefficient 1.00, which is the same as the theoretical solution with the rigid body limit equilibrium method: . Thus, the states of all nodal pairs at the open or slip state for the first time can be regarded as the failure criterion.
4. Numerical Examples
4.1. Stability Analysis of Wedge
Figure 4 shows the analytical model of a typical slope with a wedge. The slope system consists of bedrock, wedge, L-WSP (left weak structural plane), and R-WSP (right weak structural plane). The maximum height of the wedge is equal to 64.89 m. There are 2212 elements in the bedrock, 504 elements in the wedge, and 239 nodal pairs used to simulate the interfaces between the wedge and the bedrock. The number of nodes is 3592. The mesh size of the solid elements is about 10.00 m~30.00 m. The wedge has four surfaces, attitudes of which are listed in Table 4. L-WSP and R-WSP are the sliding planes.
(a) Model of slope
(b) Model of wedge
The method presented in this paper is applied to solve the safety factors of a material symmetric wedge and a material asymmetric wedge. Table 5 shows the mechanical parameters of slope. The wedge and the bedrock are considered as elastic-plastic bodies. In analysis of the material symmetric wedge, the mechanical parameters of L-WSP and R-WSP used are the same as L-WSP. The bottom of the slope is fixed and the normal displacements on the four vertical sides of the bedrock are constrained. Only considering the strength properties of the sliding planes, normal contact stiffness and damage ratio are given and 0.995, respectively.
Cohesion and friction coefficients are reduced gradually at the same level until all nodal pairs are at the open or slip state for the first time. The safety factor of the material symmetric wedge is 1.554 with this paper’s method. It is in good agreement with the theoretical solution 1.556 from Hoek and Bray [21]. The safety factor of the material asymmetric wedge with this paper’s method is 1.16. It is 2.0% less than 1.18 from Hoek and Bray. But the paper [22] gives a range from 1.143 to 1.167 according to the different assumption of the shear stress distribution on WSP. The safety factor 1.16 in this paper is also in that range.
Figures 5(a) and 5(b) show the displacement vector of slip surfaces on - plane in the material symmetric wedge and in the material asymmetric wedge, respectively. It can be clearly seen that the wedge slides along the intersection between L-WSP and R-WSP but the displacement in R-WSP is larger than L-WSP in the material asymmetric wedge. Due to the shear strength parameters of R-WSP being lower than L-WSP, R-WSP is easier to slide. Obviously, the actual characteristics of deformation can be reflected with this paper’s method.
(a) Material symmetric wedge
(b) Material asymmetric wedge
4.2. Stability Analysis of Fractured Rock Slope
In Figure 6, the width of the slope is equal to 27.00 m and the height 15.46 m. The slope system consists of 40 blocks and 75 interface joints, which are cut by 8 vertical fractures and 6 transverse fractures. The vertical spacing of fractures is 2.00 m and the transverse spacing of fractures is 3.00 m. 39 sliding blocks are on the top of the slide bed with dip angle 20°. The number of total nodes is 1515. 1100 solid elements are used to simulate the partitioned blocks and 374 nodal pairs are used to simulate the interface joints. The mesh size of the solid elements is from 0.50 m to 0.70 m.
The method presented in this paper is applied to solve the safety factor and to explore the failure mode. The mechanical parameters are listed in Table 6 and all blocks are regarded as elastic-plastic bodies. The bottom of the mesh is fixed and the horizontal displacements on both sides of the slide bed are fixed. Only considering the strength properties of the interface joints, normal contact stiffness is given and damage ratio 0.99, respectively.
Reduction coefficient Fs is set to 2.0; computation is divergent. Decrease to 1.90; the results of deformation and failure pairs are shown in Figures 7(a) and 8(a), from which the nodal pairs located in the bottom of the first and the second lists are at the slip state, but these located in the bottom side are still at the close state and the slope can keep stable. Gradually increase Fs to 1.95; Figure 7(b) shows the distribution of the failure nodal pairs and Figure 8(b) shows the deformation of the slope. It can be seen that the nodal pairs located in the bottom of the first to the third lists are at the slip state and these located in the bottom side of the first list are all at the open state; but the nodal pairs located in the bottom side of the second and the third lists are still at the close state. Due to the failure pairs of the first list blocks run through for the first time, the first list blocks cannot keep stable, and the toppling failure will occur.
(a) Fs = 1.90
(b) Fs = 1.95
(a) Fs = 1.90
(b) Fs = 1.95
In paper [23], the results are achieved by UDEC software, which show that the safety factor is 1.96 and the failure mode is the same as the toppling failure occurring in the first list blocks. The difference of safety factors is 0.5%.
5. Engineering Example
Figure 9 shows a cross section of an engineering slope, the height of which is equal to 240 m. There is a nearly vertical cliff from the 930 m elevation to the top. Slope system consists of 6 blocks (Figure 10) and 9 weak structural planes (Figure 11). Six blocks are the upper dangerous rock masses DR2 and A1, the middle mudstone blocks B1 and B2, and bedrock and the left accumulation body FGC with the thickness of 20~50 m. Nine weak structural planes are L1, J1~J4, and Jc1~Jc4, respectively. L1 cut through J1~J3 until J4, which is the interface between B2 and bedrock. J3 cuts the mudstone into B1 and B2 blocks. J2 is the interface between A1 and B1. J1 cuts the dangerous rock mass into DR2 and A1 blocks. Jc1~Jc4 are located in the bottom of FGC. Sliding channel may be formed by L1, J1~J4, and Jc1~Jc4.
In Figure 10, the simulated region of bedrock covers the whole area with a dimension of . The width of the dangerous rock mass DR2 is 175.0 m. There are 17062 nodes, 12142 solid elements, and 2584 nodal pairs used to simulate the interface boundaries. The bottom of the mesh is fixed and normal displacements on the four sides of the bedrock are constrained. The elastic-plastic constitutive is assign to the blocks B1 and B2, for geological exploration shows that there is a certain compression deformation in the mudstone area. Mechanical material parameters are listed in Table 7. Not only the strength but also the deformation properties of the interface joints are considered. The method presented in this paper is applied to solve the safety factor and to search the sliding channel.
Cohesion and friction coefficients are reduced gradually at the same level. All nodal pairs on L1 are failure at the reduction coefficient 1.0. Figures 12(a) and 13(a) show the failure pairs on the key weak structural planes at the reduction coefficient 1.03. These figures show that the nodal pairs on Jc4 are all failure and the failure nodal pairs on J3 occupy 2/5 of the area of J3. Although most of the nodal pairs on J4 are at the failure state, the failure pairs are not run through and the rock mass can keep stable.
(a) Fs = 1.03
(b) Fs = 1.05
(a) Fs = 1.03
(b) Fs = 1.05
Increase reduction coefficient Fs to 1.05; the results of the key weak structural planes are shown in Figures 12(b) and 13(b), from which the number of the failure nodal pairs on J3 continues increasing and the nodal pairs on Jc3, Jc4, and J4 are all in failure. The failure pairs run through from L1, Jc4 to J4 for the first time, and the rock mass cannot keep stable. Thus, the safety factor of slope is 1.05 and the possible sliding channel is along Jc4, J4, and L1.
6. Conclusion
RSM based on iterative method of PFE and IBE is studied and used to the stability analysis of local discontinuous rock mass in this paper. The contact forces and the total incremental displacements on joints can be solved iteratively by IBE, and the displacements of continuous bodies can be calculated by FEM after the contact forces are solved and added in the total load vector. Not only the deformation of blocks but also the large deformation of slope system can be considered. Flexibility matrix is formed through PFE according to the contact states and spring flexibility is used to reflect the influence of weak structural planes. Nonlinear iteration is only limited to the possible contact region so that the computational efficiency can be improved. There is a great potential use for local discontinuous rock mass with iterative method of PFE and IBE.
With cohesion and friction coefficient reduced gradually, the states of all nodal pairs at the open or slip state for the first time were regarded as a failure criterion, which can decrease the effect of subjectivity in determining safety factor. Stability analysis of the numerical and engineer slopes shows that SRM based on iterative method of PFE and IBE is an efficient tool for the stability analysis of local discontinuous rock mass.
The process from stability to collapsing can be simulated, but the movement of slope after all nodal pairs are at the open state cannot be achieved with this method at present.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by Chinese National Natural Science Foundation (nos. 51209076 and 51279050), Research and Innovation Program for Graduate Students in Jiangsu Province (no. CXZZ13_0254), and Special and Fund for Scientific Research in the Public Interest of Ministry of Water Resources (no. 201001035).