# Strength Reduction Method for Stability Analysis of Local Discontinuous Rock Mass with Iterative Method of Partitioned Finite Element and Interface Boundary Element

Authors: Tongchun Li; Jinwen He; Lanhao Zhao; Xiaona Li; Zhiwei Niu
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 Figure1, a rock slope is composed of slip mass Ω
1, bed rock Ω
2, and the discontinuous joint with the gap of δ. ξ
η
ζ and x
y
z 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 F and contact internal force f. f
1 and f
2 are the contact forces on two nodes of each nodal pair and f
ξ
,
f
η
,
a
n
d
f
ζ are the three stress components in the local coordinate system, assuming tension is positive in all forces.Figure 1
Mechanical model.If the contact force incrementΔ
f is given, the stress and deformation of Ω
1 and Ω
2 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Δ
u
t
i at point i on interface boundary will be changed by the effect of contact force increment Δ
f. Δ
u
t
i includes the incremental displacement of deformable body Δ
u
i and the rigid incremental displacement of rigid body Δ
u
g
i; that is,(1)
Δ
u
t
i
=
Δ
u
i
+
Δ
u
g
i
.
(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 Table1, 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.Table 1
Definition of contact state on interface boundary.
Contact state
Displacement criterion
Stress criterion
open
δ
=
0
f
1
=
f
2
=
0
close
δ
>
0
f
1
=
-
f
2
,
f
ς
<
σ
t
A
,
f
ξ
2
+
f
η
2
<
-
μ
f
ζ
+
A
c
slip
δ
=
0
f
ξ
2
+
f
η
2
=
-
μ
f
ζ
+
A
c
,
f
ς
>
σ
t
A
A is the area controlled by nodal pair; σ
t, c, and μ are tensile strength, cohesion, and friction coefficient, respectively.Table 2
Iterative operation.
Serial number
State transition
Before iteration of state
After iteration of state
Update
1
open—close
δ
>
0
δ
≤
0
δ
=
0
2
close or slip—open
f
ζ
<
0
f
ς
>
σ
t
A
f
1
=
f
2
=
0
3
close or slip—slip
f
ζ
<
0
f
ξ
+
f
η
>
-
μ
f
ζ
+
A
c
α
1
=
|
f
ξ
|
/
f
ξ
+
f
η
;
α
2
=
|
f
η
|
/
f
ξ
+
f
η
f
ξ
=
α
1
(
-
μ
f
ζ
+
A
c
)
;
f
η
=
α
2
(
-
μ
f
ζ
+
A
c
)


### 2.2. Partitioned Finite Element Equations

For contact bodiesΩ
1 and Ω
2, if the analysis of nth step has been finished and the contact force increment Δ
f
n been solved, the static incremental equilibrium equation at step (
n
+
1
) without considering rigid motion is stated as follows:(2)
K
Δ
u
n
=
F
n
+
1
-
∫
B
T
σ
n
d
Ω
-
f
n
+
Δ
f
n
,where K is the global stiffness matrix; Δ
u
n represents the vector of incremental displacement at step n; F
n
+
1 is the vector of total external load at step n
+
1; ∫
B
T
σ
n
d
Ω
-
f
n stands for the external load at step n; f
n is the vector of contact internal force at step n; σ
n is the stress matrix.C
=
K
-
1 is introduced to (2); the following equation can be gotten: (3)
Δ
u
n
=
Δ
u
-
n
+
C
Δ
f
n
,where Δ
u
-
n
=
C
[
F
n
+
1
-
(
∫
B
T
σ
n
d
Ω
-
f
n
)
]; an arbitrary component c
i
j in the matrix C represents the flexibility coefficient corresponding to the displacement at the freedom i due to a unit force at the freedom j. 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 theith block is composed of translation components and rotation components. Assuming the centroid displacement increment Δ
γ
i of ith block is {
Δ
u
-
,
Δ
v
-
,
Δ
w
-
,
Δ
θ
x
,
Δ
θ
y
,
Δ
θ
z
}
i
T, the rigid incremental displacement of the kth node with a relative coordinate value of (
Δ
x
,
Δ
y
,
Δ
z
) to centroidin ith block is given as follows:(4)
Δ
u
g
=
1
0
0
0
Δ
z
-
Δ
y
0
1
0
-
Δ
z
0
Δ
x
0
0
1
Δ
y
-
Δ
x
0
Δ
u
-
Δ
v
-
Δ
w
-
Δ
θ
x
Δ
θ
y
Δ
θ
z
=
ω
i
k
Δ
γ
n
i
,where ω
i
k is the transform matrix of kth node to the centroid of ith block.Equations (3) and (4) plugged into (1), the total incremental displacement of node k on any nodal pair of contact blocks is (5)
Δ
u
t
n
k
=
Δ
u
-
n
k
+
C
k
Δ
f
n
k
+
ω
i
k
Δ
γ
n
i
.Equation (5) can be written in another form; that is, (6)
C
k
Δ
f
n
k
+
ω
i
k
Δ
γ
n
i
=
Δ
u
t
n
k
-
Δ
u
-
n
k
,where C
k is the matrix relevant to node k in the total flexibility matrix; Δ
u
-
n
k 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. Δ
u
t
n
k 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 Δ
f
n
k and the rigid displacement increment Δ
γ
n
i of ith block.For each contact body, the equilibrium equation related to centroid point at stepn for all forces acted on it could be written as follows: (7)
W
i
γ
¨
n
i
+
α
W
i
γ
˙
n
i
=
∑
k
=
1
p
i
ω
i
k
T
f
n
k
+
∑
l
=
1
q
i
ω
i
l
T
F
n
l
,where W
i
=
m
i
0
0
I
i; m
i, I
i are mass and rotational inertia related to the centroid of ith block, respectively; α represents the mass-damping coefficient; γ
¨
n
i
,
γ
˙
n
i stand for the rigid acceleration and rigid velocity of ith block, respectively; p
i is the number of contact nodes on ith block; q
i is the number of total nodes on ith block; f
n
k is the contact force of node k; F
n
l is the external load acted on the lth node of ith block.For static problem, left side of (7) is zero, and the incremental formula is(8)
∑
K
=
1
p
i
ω
i
k
T
Δ
f
n
k
+
∑
L
=
1
q
i
ω
i
l
T
Δ
F
n
l
=
0
.Combing (6) and (8), the static iterative equations on discontinuous interface boundaries of ith block can be gotten as follows:(9)
C
i
ω
Γ
i
ω
Γ
i
T
0
Δ
f
n
i
Δ
γ
n
i
=
I
0
0
-
ω
Ω
i
T
Δ
u
t
n
-
Δ
u
-
n
i
Δ
F
n
i
,where subscript Γ is the interface boundary and the dimension of [
ω
Γ
] is (
p
i
∗
d
i
m
)
×
[
(
d
i
m
-
1
)
∗
3
]; subscript Ω stands for the whole region of block and the dimension of [
ω
Ω
]
T is [
(
d
i
m
-
1
)
∗
3
]
×
(
q
i
∗
d
i
m
); 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 inZ-direction is far less than the size in X and Y directions, three stress components, σ
z, τ
z
x, and τ
z
y, can be used to describe nodal pair. The incremental constitutive law can be described as follows: (10)
d
σ
j
=
K
0
d
u
j
=
k
z
z
k
z
x
k
z
y
d
u
d
v
d
ω
,where {
d
σ
j
}
=
{
d
σ
z
,
d
τ
z
x
,
d
τ
z
y
}
T; [
K
0
] is the stiffness matrix of weak structural plane.AssumingE 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 k
z
z, such as (500~1000)E, will be given. If the elastic module E
j and thickness B
j of weak structural plane are given, real normal stiffness can be calculated by k
z
z
=
E
j
/
B
j according to the principle of equivalent deformation under uniaxial compression [18]. If a pair is at the open state, a small k
z
z
_
m
i
n, 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]: (11)
k
z
x
k
z
y
=
k
z
z
2
1
+
μ
.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:(12)
k
z
x
k
z
y
=
k
z
z
2
1
+
μ
1
-
R
2
,where R is damage ratio; a value close to 1.00 is suggested.Equation (10) multiplies A at both sides. A
{
d
σ
j
} is the contact force increment at point j; {
d
u
j
} represents the relative incremental displacement at nodal pair j. Therefore, (13) can be gotten; that is,(13)
K
0
-
1
A
Δ
f
j
=
Δ
u
j
.Assuming[
c
¯
j
]
=
[
K
0
]
-
1
/
A represent the influence flexibility matrix of nodal pair j on weak structural plane, [
C
¯
]
i can be assembled by the nodal pairs around block i and add it in (9); the static iterative equations on discontinuous interface boundaries considering the influence of weak structural planes are(14)
C
i
+
C
¯
i
ω
Γ
i
ω
Γ
i
T
0
Δ
f
n
i
Δ
γ
n
i
=
I
0
0
-
ω
Ω
i
T
Δ
u
t
n
-
Δ
u
-
n
i
Δ
F
n
i
,where Δ
u
t
n is the total incremental displacement on interface boundaries of block i, 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:(15)
A
Δ
f
n
1
⋯
Δ
f
n
m
Δ
γ
n
1
⋯
Δ
γ
n
m
T
=
B
,where (16)
A
=
ω
Γ
1
C
+
C
¯
⋯
ω
Γ
m
ω
Γ
1
T
⋯
ω
Γ
m
T
,
B
=
I
-
ω
Ω
1
T
⋯
-
ω
Ω
m
T
Δ
u
t
n
-
Δ
u
-
n
1
⋯
Δ
u
t
n
-
Δ
u
-
n
m
Δ
F
n
1
⋯
Δ
F
n
m
.


### 2.5. Numerical Implementation

Combining the principles and equations above, the numerical implementation of the iterative method can be summarized as follows (the superscripti represents the contact iteration step while the subscript n stands for the loading increment step).Step 1.
Assemble and decompose the total stiffness matrixK, solve the flexibility matrix C on possible contact interface boundaries by the unit force method, form the initial influence flexibility matrix c
¯
j on weak structural planes, and get the transform matrix ω on each node.Step 2.
If the load increment at stepn has been finished, get into step (
n
+
1
) 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Δ
u
-
n by back substitution of the equation Δ
u
-
n
=
K
-
1
(
F
n
+
1
+
f
n
-
∫
B
T
σ
n
d
Ω
), add Δ
u
-
n, rigid incremental displacement ω
T
Δ
γ
n, and the incremental deformation of weak structural planes in total incremental displacement Δ
u
t
n, form ω
T
Δ
F
n and assemble vector B, and assemble and decompose matrix A in (15).Step 4 (the iteration of contact force and rigid displacement).
If the iteration of contact force and rigid displacement at stepi has been finished, get into step (
i
+
1
), solve vector {
Δ
f
n
}
1
i
+
1
⋯
{
Δ
f
n
}
m
i
+
1
{
Δ
γ
n
}
1
i
+
1
⋯
{
Δ
γ
n
}
m
i
+
1
T by using (15), and judge the convergence of iteration by using Δ
f
n
i
+
1
=
f
n
i
+
1
-
f
n
i
<
ε.
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Δ
u
-
n
i
+
1 and achieve the rigid incremental displacement ω
T
Δ
γ
n
i
+
1 on interface boundaries, add them to Δ
u
t
n, get the external load increment ω
T
Δ
F
n
i
+
1 on each body, and assemble vector B in (15).
If iteration is convergent, the mutual transition of the three states will be carried out as follows according to Table2.(1)
If the current state is at the open state, new gap can be obtained by usingδ
n
+
1
i
+
1
=
δ
n
+
1
i
+
(
Δ
u
t
n
1
,
i
-
Δ
u
t
n
2
,
i
)
·
ξ
⃑. If δ
n
+
1
i
+
1 is less than zero, the current state is updated to the close state and δ
n
+
1
i
+
1 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(
σ
n
i
+
1
)
ξ
=
σ
n
i
+
1
·
ξ
⃑. If (
σ
n
i
+
1
)
ξ is less than σ
t
i, calculate f
ξ
η
=
f
ξ
2
+
f
η
2. If (
σ
n
i
+
1
)
ξ is greater than σ
t
i, the current state is changed to the open state and contact internal force f
n
+
1
i
+
1 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 usingf
ξ
η
=
f
ξ
2
+
f
η
2. If f
ξ
η
>
A
c
+
N
μ is satisfied, the current state is updated to the slip state, and a new condition of f
ξ
2
+
f
η
2
=
-
μ
f
ξ
+
A
c 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 c
¯
j 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 matrixC, and go to Step 3.Step 6.
If the contact state is convergent, go to Step2 and get into the next load increment until all load increment steps are completed.


## 2.1. Mechanical Model

In Figure1, a rock slope is composed of slip mass Ω
1, bed rock Ω
2, and the discontinuous joint with the gap of δ. ξ
η
ζ and x
y
z 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 F and contact internal force f. f
1 and f
2 are the contact forces on two nodes of each nodal pair and f
ξ
,
f
η
,
a
n
d
f
ζ are the three stress components in the local coordinate system, assuming tension is positive in all forces.Figure 1
Mechanical model.If the contact force incrementΔ
f is given, the stress and deformation of Ω
1 and Ω
2 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Δ
u
t
i at point i on interface boundary will be changed by the effect of contact force increment Δ
f. Δ
u
t
i includes the incremental displacement of deformable body Δ
u
i and the rigid incremental displacement of rigid body Δ
u
g
i; that is,(1)
Δ
u
t
i
=
Δ
u
i
+
Δ
u
g
i
.
(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 Table1, 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.Table 1
Definition of contact state on interface boundary.
Contact state
Displacement criterion
Stress criterion
open
δ
=
0
f
1
=
f
2
=
0
close
δ
>
0
f
1
=
-
f
2
,
f
ς
<
σ
t
A
,
f
ξ
2
+
f
η
2
<
-
μ
f
ζ
+
A
c
slip
δ
=
0
f
ξ
2
+
f
η
2
=
-
μ
f
ζ
+
A
c
,
f
ς
>
σ
t
A
A is the area controlled by nodal pair; σ
t, c, and μ are tensile strength, cohesion, and friction coefficient, respectively.Table 2
Iterative operation.
Serial number
State transition
Before iteration of state
After iteration of state
Update
1
open—close
δ
>
0
δ
≤
0
δ
=
0
2
close or slip—open
f
ζ
<
0
f
ς
>
σ
t
A
f
1
=
f
2
=
0
3
close or slip—slip
f
ζ
<
0
f
ξ
+
f
η
>
-
μ
f
ζ
+
A
c
α
1
=
|
f
ξ
|
/
f
ξ
+
f
η
;
α
2
=
|
f
η
|
/
f
ξ
+
f
η
f
ξ
=
α
1
(
-
μ
f
ζ
+
A
c
)
;
f
η
=
α
2
(
-
μ
f
ζ
+
A
c
)


## 2.2. Partitioned Finite Element Equations

For contact bodiesΩ
1 and Ω
2, if the analysis of nth step has been finished and the contact force increment Δ
f
n been solved, the static incremental equilibrium equation at step (
n
+
1
) without considering rigid motion is stated as follows:(2)
K
Δ
u
n
=
F
n
+
1
-
∫
B
T
σ
n
d
Ω
-
f
n
+
Δ
f
n
,where K is the global stiffness matrix; Δ
u
n represents the vector of incremental displacement at step n; F
n
+
1 is the vector of total external load at step n
+
1; ∫
B
T
σ
n
d
Ω
-
f
n stands for the external load at step n; f
n is the vector of contact internal force at step n; σ
n is the stress matrix.C
=
K
-
1 is introduced to (2); the following equation can be gotten: (3)
Δ
u
n
=
Δ
u
-
n
+
C
Δ
f
n
,where Δ
u
-
n
=
C
[
F
n
+
1
-
(
∫
B
T
σ
n
d
Ω
-
f
n
)
]; an arbitrary component c
i
j in the matrix C represents the flexibility coefficient corresponding to the displacement at the freedom i due to a unit force at the freedom j. 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 theith block is composed of translation components and rotation components. Assuming the centroid displacement increment Δ
γ
i of ith block is {
Δ
u
-
,
Δ
v
-
,
Δ
w
-
,
Δ
θ
x
,
Δ
θ
y
,
Δ
θ
z
}
i
T, the rigid incremental displacement of the kth node with a relative coordinate value of (
Δ
x
,
Δ
y
,
Δ
z
) to centroidin ith block is given as follows:(4)
Δ
u
g
=
1
0
0
0
Δ
z
-
Δ
y
0
1
0
-
Δ
z
0
Δ
x
0
0
1
Δ
y
-
Δ
x
0
Δ
u
-
Δ
v
-
Δ
w
-
Δ
θ
x
Δ
θ
y
Δ
θ
z
=
ω
i
k
Δ
γ
n
i
,where ω
i
k is the transform matrix of kth node to the centroid of ith block.Equations (3) and (4) plugged into (1), the total incremental displacement of node k on any nodal pair of contact blocks is (5)
Δ
u
t
n
k
=
Δ
u
-
n
k
+
C
k
Δ
f
n
k
+
ω
i
k
Δ
γ
n
i
.Equation (5) can be written in another form; that is, (6)
C
k
Δ
f
n
k
+
ω
i
k
Δ
γ
n
i
=
Δ
u
t
n
k
-
Δ
u
-
n
k
,where C
k is the matrix relevant to node k in the total flexibility matrix; Δ
u
-
n
k 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. Δ
u
t
n
k 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 Δ
f
n
k and the rigid displacement increment Δ
γ
n
i of ith block.For each contact body, the equilibrium equation related to centroid point at stepn for all forces acted on it could be written as follows: (7)
W
i
γ
¨
n
i
+
α
W
i
γ
˙
n
i
=
∑
k
=
1
p
i
ω
i
k
T
f
n
k
+
∑
l
=
1
q
i
ω
i
l
T
F
n
l
,where W
i
=
m
i
0
0
I
i; m
i, I
i are mass and rotational inertia related to the centroid of ith block, respectively; α represents the mass-damping coefficient; γ
¨
n
i
,
γ
˙
n
i stand for the rigid acceleration and rigid velocity of ith block, respectively; p
i is the number of contact nodes on ith block; q
i is the number of total nodes on ith block; f
n
k is the contact force of node k; F
n
l is the external load acted on the lth node of ith block.For static problem, left side of (7) is zero, and the incremental formula is(8)
∑
K
=
1
p
i
ω
i
k
T
Δ
f
n
k
+
∑
L
=
1
q
i
ω
i
l
T
Δ
F
n
l
=
0
.Combing (6) and (8), the static iterative equations on discontinuous interface boundaries of ith block can be gotten as follows:(9)
C
i
ω
Γ
i
ω
Γ
i
T
0
Δ
f
n
i
Δ
γ
n
i
=
I
0
0
-
ω
Ω
i
T
Δ
u
t
n
-
Δ
u
-
n
i
Δ
F
n
i
,where subscript Γ is the interface boundary and the dimension of [
ω
Γ
] is (
p
i
∗
d
i
m
)
×
[
(
d
i
m
-
1
)
∗
3
]; subscript Ω stands for the whole region of block and the dimension of [
ω
Ω
]
T is [
(
d
i
m
-
1
)
∗
3
]
×
(
q
i
∗
d
i
m
); 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 inZ-direction is far less than the size in X and Y directions, three stress components, σ
z, τ
z
x, and τ
z
y, can be used to describe nodal pair. The incremental constitutive law can be described as follows: (10)
d
σ
j
=
K
0
d
u
j
=
k
z
z
k
z
x
k
z
y
d
u
d
v
d
ω
,where {
d
σ
j
}
=
{
d
σ
z
,
d
τ
z
x
,
d
τ
z
y
}
T; [
K
0
] is the stiffness matrix of weak structural plane.AssumingE 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 k
z
z, such as (500~1000)E, will be given. If the elastic module E
j and thickness B
j of weak structural plane are given, real normal stiffness can be calculated by k
z
z
=
E
j
/
B
j according to the principle of equivalent deformation under uniaxial compression [18]. If a pair is at the open state, a small k
z
z
_
m
i
n, 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]: (11)
k
z
x
k
z
y
=
k
z
z
2
1
+
μ
.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:(12)
k
z
x
k
z
y
=
k
z
z
2
1
+
μ
1
-
R
2
,where R is damage ratio; a value close to 1.00 is suggested.Equation (10) multiplies A at both sides. A
{
d
σ
j
} is the contact force increment at point j; {
d
u
j
} represents the relative incremental displacement at nodal pair j. Therefore, (13) can be gotten; that is,(13)
K
0
-
1
A
Δ
f
j
=
Δ
u
j
.Assuming[
c
¯
j
]
=
[
K
0
]
-
1
/
A represent the influence flexibility matrix of nodal pair j on weak structural plane, [
C
¯
]
i can be assembled by the nodal pairs around block i and add it in (9); the static iterative equations on discontinuous interface boundaries considering the influence of weak structural planes are(14)
C
i
+
C
¯
i
ω
Γ
i
ω
Γ
i
T
0
Δ
f
n
i
Δ
γ
n
i
=
I
0
0
-
ω
Ω
i
T
Δ
u
t
n
-
Δ
u
-
n
i
Δ
F
n
i
,where Δ
u
t
n is the total incremental displacement on interface boundaries of block i, 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:(15)
A
Δ
f
n
1
⋯
Δ
f
n
m
Δ
γ
n
1
⋯
Δ
γ
n
m
T
=
B
,where (16)
A
=
ω
Γ
1
C
+
C
¯
⋯
ω
Γ
m
ω
Γ
1
T
⋯
ω
Γ
m
T
,
B
=
I
-
ω
Ω
1
T
⋯
-
ω
Ω
m
T
Δ
u
t
n
-
Δ
u
-
n
1
⋯
Δ
u
t
n
-
Δ
u
-
n
m
Δ
F
n
1
⋯
Δ
F
n
m
.


## 2.5. Numerical Implementation

Combining the principles and equations above, the numerical implementation of the iterative method can be summarized as follows (the superscripti represents the contact iteration step while the subscript n stands for the loading increment step).Step 1.
Assemble and decompose the total stiffness matrixK, solve the flexibility matrix C on possible contact interface boundaries by the unit force method, form the initial influence flexibility matrix c
¯
j on weak structural planes, and get the transform matrix ω on each node.Step 2.
If the load increment at stepn has been finished, get into step (
n
+
1
) 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Δ
u
-
n by back substitution of the equation Δ
u
-
n
=
K
-
1
(
F
n
+
1
+
f
n
-
∫
B
T
σ
n
d
Ω
), add Δ
u
-
n, rigid incremental displacement ω
T
Δ
γ
n, and the incremental deformation of weak structural planes in total incremental displacement Δ
u
t
n, form ω
T
Δ
F
n and assemble vector B, and assemble and decompose matrix A in (15).Step 4 (the iteration of contact force and rigid displacement).
If the iteration of contact force and rigid displacement at stepi has been finished, get into step (
i
+
1
), solve vector {
Δ
f
n
}
1
i
+
1
⋯
{
Δ
f
n
}
m
i
+
1
{
Δ
γ
n
}
1
i
+
1
⋯
{
Δ
γ
n
}
m
i
+
1
T by using (15), and judge the convergence of iteration by using Δ
f
n
i
+
1
=
f
n
i
+
1
-
f
n
i
<
ε.
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Δ
u
-
n
i
+
1 and achieve the rigid incremental displacement ω
T
Δ
γ
n
i
+
1 on interface boundaries, add them to Δ
u
t
n, get the external load increment ω
T
Δ
F
n
i
+
1 on each body, and assemble vector B in (15).
If iteration is convergent, the mutual transition of the three states will be carried out as follows according to Table2.(1)
If the current state is at the open state, new gap can be obtained by usingδ
n
+
1
i
+
1
=
δ
n
+
1
i
+
(
Δ
u
t
n
1
,
i
-
Δ
u
t
n
2
,
i
)
·
ξ
⃑. If δ
n
+
1
i
+
1 is less than zero, the current state is updated to the close state and δ
n
+
1
i
+
1 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(
σ
n
i
+
1
)
ξ
=
σ
n
i
+
1
·
ξ
⃑. If (
σ
n
i
+
1
)
ξ is less than σ
t
i, calculate f
ξ
η
=
f
ξ
2
+
f
η
2. If (
σ
n
i
+
1
)
ξ is greater than σ
t
i, the current state is changed to the open state and contact internal force f
n
+
1
i
+
1 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 usingf
ξ
η
=
f
ξ
2
+
f
η
2. If f
ξ
η
>
A
c
+
N
μ is satisfied, the current state is updated to the slip state, and a new condition of f
ξ
2
+
f
η
2
=
-
μ
f
ξ
+
A
c 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 c
¯
j 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 matrixC, and go to Step 3.Step 6.
If the contact state is convergent, go to Step2 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 c
,
φ are divided by reduction coefficient F
s, new parameters of shear strength c
′
,
φ
′ can be gotten. Calculating by using (17) repeatedly until the critical failure state, the reduction coefficient is safety factor:(17)
c
′
=
c
F
s
,
φ
′
=
a
r
c
t
a
n
t
a
n
φ
F
s
.As shown in Figure2, 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 ρ
=
1.25 g/cm3, elastic module E
=
2.4
×
1
0
8 Pa, Poisson’s ratio μ
=
0.163, vertical uniform load N
=
25 kN/m, horizontal uniform load Q
=
125 kN/m, cohesion C
=
0.1 MPa, and friction coefficient f
=
0.5. 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 R are given 1000
E and 0.992, respectively. If any pair is at the open or slip state, the nodal pair is failure.Figure 2
Grid and blocks of slope system.Cohesion and friction coefficient reduced gradually at the same level;X-direction displacement of point A and the number of failure pairs are shown in Table 3. Curve of reduction coefficient and X-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: K
=
(
f
∑
N
+
c
A
)
/
∑
Q. Thus, the states of all nodal pairs at the open or slip state for the first time can be regarded as the failure criterion.Table 3
Results in difference reduction coefficient.
Reductioncoefficients
X-displacementof point A/mm
The number offailure pairs
0.500
2.309
0
0.600
2.309
0
0.650
3.109
1
0.700
6.485
3
0.750
6.577
3
0.800
6.657
3
0.850
6.728
3
0.900
6.791
3
0.950
6.847
3
0.999
6.859
3
1.000
16.070
5Figure 3
Curve of reduction coefficient andX-direction displacement of point A.


## 4. Numerical Examples

### 4.1. Stability Analysis of Wedge

Figure4 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.Table 4
Attitude of surfaces around wedge.
Surface
Dip direction angle/°
Dip angle/°
Top surface
180
10
Front surface
180
60
L-WSP
115
45
R-WSP
245
45Figure 4
Model of slope system.
(a)
Model of slope
(b)
Model of wedgeThe method presented in this paper is applied to solve the safety factors of a material symmetric wedge and a material asymmetric wedge. Table5 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 R are given 1000
E and 0.995, respectively.Table 5
Mechanical parameters of wedge.
Blocks
Natural density
Shear strength
Elastic module
Poisson’s ratio
kg/m3
c/kPa
φ/°
MPa
L-WSP
2000
50
20
50
0.3
R-WSP
2000
30
10
50
0.3
Bedrock/wedge
2600
1000
45
1000
0.15Cohesion 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.Figures5(a) and 5(b) show the displacement vector of slip surfaces on X-Y 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.Figure 5
Displacement vector of slip surface onX-Y plane.
(a)
Material symmetric wedge
(b)
Material asymmetric wedge


### 4.2. Stability Analysis of Fractured Rock Slope

In Figure6, 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.Figure 6
Grid and blocks of slope system.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 Table6 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 1000
E and damage ratio R0.99, respectively.Table 6
Mechanical parameters of material.
Material
Natural density
Elastic module
Poisson’s ratio
Shear strength
Tensile strength
kg/m3
MPa
c/kPa
φ/°
kPa
Rock
2000
1000
0.25
1000
50
1000
Weak structural plane
2000
1000
0.25
50
30
0Reduction coefficient Fs is set to 2.0; computation is divergent. DecreaseF
s 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.Figure 7
Failure pairs on interface joints (△ is slip; □ is open).
(a)
Fs = 1.90
(b)
Fs = 1.95Figure 8
Deformation of critical failure state.
(a)
Fs = 1.90
(b)
Fs = 1.95In 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%.


## 4.1. Stability Analysis of Wedge

Figure4 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.Table 4
Attitude of surfaces around wedge.
Surface
Dip direction angle/°
Dip angle/°
Top surface
180
10
Front surface
180
60
L-WSP
115
45
R-WSP
245
45Figure 4
Model of slope system.
(a)
Model of slope
(b)
Model of wedgeThe method presented in this paper is applied to solve the safety factors of a material symmetric wedge and a material asymmetric wedge. Table5 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 R are given 1000
E and 0.995, respectively.Table 5
Mechanical parameters of wedge.
Blocks
Natural density
Shear strength
Elastic module
Poisson’s ratio
kg/m3
c/kPa
φ/°
MPa
L-WSP
2000
50
20
50
0.3
R-WSP
2000
30
10
50
0.3
Bedrock/wedge
2600
1000
45
1000
0.15Cohesion 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.Figures5(a) and 5(b) show the displacement vector of slip surfaces on X-Y 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.Figure 5
Displacement vector of slip surface onX-Y plane.
(a)
Material symmetric wedge
(b)
Material asymmetric wedge


## 4.2. Stability Analysis of Fractured Rock Slope

In Figure6, 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.Figure 6
Grid and blocks of slope system.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 Table6 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 1000
E and damage ratio R0.99, respectively.Table 6
Mechanical parameters of material.
Material
Natural density
Elastic module
Poisson’s ratio
Shear strength
Tensile strength
kg/m3
MPa
c/kPa
φ/°
kPa
Rock
2000
1000
0.25
1000
50
1000
Weak structural plane
2000
1000
0.25
50
30
0Reduction coefficient Fs is set to 2.0; computation is divergent. DecreaseF
s 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.Figure 7
Failure pairs on interface joints (△ is slip; □ is open).
(a)
Fs = 1.90
(b)
Fs = 1.95Figure 8
Deformation of critical failure state.
(a)
Fs = 1.90
(b)
Fs = 1.95In 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

Figure9 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.Figure 9
Cross section of slope.Figure 10
Grid and blocks of slope system.Figure 11
Weak structural planes.In Figure10, the simulated region of bedrock covers the whole area with a dimension of 290.0
m
×
274.5
m
×
345.0
m. 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.Table 7
Mechanical material parameters.
Material
Natural density
Elastic module
Poisson’s ratio
Shear strength
Tensile strength
Thickness
g/cm3
GPa
c/kPa
φ/°
kPa
cm
DR2, A1
2.68
7.0
0.30
0.90
31.5
FGC
2.40
3.0
0.30
0.45
30.0
B1, B2
2.67
3.0
0.35
0.30
24.0
bedrock
2.70
10.0
0.27
0.80
50.0
J1
2.67
3.0
0.35
0.30
24.0
0.3
0.1
J2~J4
2.67
3.0
0.35
0.01
19.5
0.01
0.1
L1
1.90
0.9
0.35
0.09
20.0
0.3
1~5
JC1~JC3
2.00
3.0
0.30
0.45
30.0
0.4
0.2
JC4
1.95
3.0
0.30
0.015
25
0.02
0.2Cohesion and friction coefficients are reduced gradually at the same level. All nodal pairs on L1 are failure at the reduction coefficient 1.0. Figures12(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.Figure 12
Failure pairs on Jc1~Jc4 weak structural plane (red cube is at the close state).
(a)
Fs = 1.03
(b)
Fs = 1.05Figure 13
Failure pairs on J3, J4 weak structural plane (red cube is at the close state).
(a)
Fs = 1.03
(b)
Fs = 1.05Increase reduction coefficient Fs to 1.05; the results of the key weak structural planes are shown in Figures12(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.



---

*Source: 872834-2015-05-20.xml*