基于变形分析的边坡潜在滑面的确定

岩石力学与工程学报 23(5)：709～716

2004年3月 Chinese Journal of Rock Mechanics and Engineering March ，2004

2003年5月7日收到初稿，2003年7月31日收到修改稿。

作者 郑 宏 简介：男，40岁，博士，现任三峡大学特聘教授、中国科学院武汉岩土力学研究所研究员，主要从事岩土力学数值方法的教学和科研工作。E-mail ：hzheng@ctgu.edu.cn 。

基于变形分析的边坡潜在滑面的确定

郑 宏1，

2 刘德富1 罗先启1

(1三峡大学 宜昌 443002) (2中国科学院武汉岩土力学研究所 武汉 430071)

摘要 基于弹塑性有限元分析的计算结果，提出了二维情况下边坡潜在滑移线应满足的一个常微分方程组初值问题，给出了该初值问题的预测-校正算法以及确保其收敛的充分必要条件，讨论了潜在滑移线的自动搜寻技术。最后，分别与极限分析法和极限平衡法的计算结果进行了对比，验证了方法的有效性。 关键词 边坡稳定性，有限元法，潜在滑移线

分类号 P 2.22 文献标识码 A 文章编号 1000-6915(2004)05-0709-08

DETERMINATION OF POTENTIAL SLIDE LINE OF SLOPES BASED ON

DEFORMATION ANALYSIS

Zheng Hong

1，2

，Liu Defu 1，Luo Xianqi 1

(1China Three Gorges University ， Yichang 443002 China )

(2Institute of Rock and Soil Mechanics ，The Chinese Academy of Sciences ， Wuhan 430071 China )

Abstract Based on the computational results from elasto-plasticity finite element analysis ，an initial value problem related to a system of ordinary differential equations(ODEs) is formulated to define the potential slide lines(PSL) for two dimensional cases. A prediction-correction algorithm for the ODEs and a necessary and sufficient condition that assures the convergence of the algorithm are presented. The skills for searching PSL with presented method are discussed ，and the procedures of this method have been validated against traditional limit analysis methods and limit equilibrium method.

Key words slope stability ，finite element method ，potential slide lines (PSL)

1 引 言

与经典的极限平衡法相比，基于变形分析的有限元法具有许多独特的优势。首先，可以得到极限 状态下的失效形式，这在很多场合是非常重要的，例如，如果采用削坡方式来改善滑坡的稳定性，就需要充分了解滑带的滑动部分和阻滑部分，滑动部分是对应于滑带中有较大剪切变形的那部分，开挖体只有位于滑带的滑动部分之上，才能起到积极作用，否则，会适得其反[1]。其次，可以了解边坡随

强度的恶化而呈现出的渐进失稳过程[2]，这样，就

可将有限的加固措施置于最紧要的部位。第三， 可考虑不同的施工工序对边坡最终安全度的影响

[3]

。第四，可以考虑影响边坡稳定的某些更复杂的

因素，如，模拟降雨过程，采用更贴合实际的本构模型，等等。

尽管在有限元分析中可以考虑更复杂的本构模型[4]，但目前在工程分析中最普遍的还是理想弹塑性模型，因为理想弹塑性分析的结果与工程师们熟悉的极限平衡法分析的结果最具可比性。文[5]也证明了，如果对理想弹塑性结构施以简单加载，则极

限载荷与相同强度参数的刚塑性体的极限载荷相等。

工程设计通常要求有限元法分析也能象极限平衡法那样为边坡提供一个安全系数和潜在滑面[6]。但不幸的是，利用强度折减方法仅能提供近似的安全系数，而且迄今为止，尚没有一个确定的数学模型去定义潜在滑面的空间位置。

目前在工程计算中，一般是根据临界平衡态的塑性区或变形图来大致估计潜在滑面的，或者是根据分析者的经验，通过手工指定一系列线段和圆弧的组合作为可能的滑移路径，并从中搜寻出安全系数最小的滑移路径作为潜在滑面。

文[7]利用贯通的广义塑性剪应变的等色图来定义滑动面。对于人工边坡，文[8]建议采用位移增量等值线来确定潜在滑面。文[2]通过使用非关联流动法则，将剪胀角ψ取为0，发现在变形后的网格中会出现一条明显的畸变带，他们将这条畸变带就定义为潜在滑面。取ψ为0意味着完全忽略了岩土材料的剪胀特性，而仅突出其剪切变形。果然，在使用他们的程序时发现，取非零的ψ后这条畸变带并不明显；而且，即使将ψ取为0，所使用的网格也必须相当规则，否则，也难以出现畸变带。根据应力计算结果，文[9]利用动态规划法，文[10]利用人工智能型优化方法的蚂蚁算法，讨论了基于瑞典法安全系数(抗滑力比滑动力)概念的潜在滑面的搜寻问题。

最近，文[4]提出的基于能量准则的岩体稳定性分析方法是一个在概念上有别于经典极限平衡法的分析方法。该方法是通过干扰能量等值线图来确定潜在滑面的。为了得到单元的干扰能量，需考虑系统的几何非线性，计算中会涉及到几何刚度矩阵的求逆和特征值问题的求解，其运算量和存贮量甚至高于应力分析本身。虽然一般情况下，所求得的安全系数及其潜在滑面都不同于经典的极限平衡法，但这仍然是一个有前景的研究方法。

文[3，5，11～17]从不同角度探讨了将有限元法用于边坡稳定性分析的一些技巧和问题，如文[12]发现如果采用在π平面上与Mohr-Coulomb屈服面有相同面积的Drucker-Prager准则，则可求得与经典极限平衡法相接近的安全系数。

本文基于弹塑性有限元分析，给出了二维情况下边坡的潜在滑面所满足的常微分方程初值问题和数值求解该初值问题的预测-校正算法，从而使这类问题有了牢固的数学力学基础。

2 滑移线的初值问题

为了不失一般性，可假设所研究的边坡处于极限平衡状态，岩土材料满足Mohr-Coulomb强度准则。因为若该边坡未处于极限平衡状态，总可以通过强度折减，进行弹塑性计算，使其达到极限平衡状态。顺便指出，文[18]证明了对于任何满足Mohr-

Coulomb强度准则的弹塑性材料，必有称之为φ-ν不等式的关系成立：sinφ≥ν2

1−，式中，φ为内摩擦角，ν为泊松比。在进行强度折减的过程中，通常要根据这一不等式调整ν，否则，极限状态下的塑性区可能偏大，而安全系数可能偏小。另外，为了使边坡处于真正的极限平衡状态，可以采用文[19]所推荐的位移控制法。

为简便起见，仅研究压剪型破坏。因材料满足Mohr-Coulomb准则，所以，边坡内任一点M在2

个方向T1和T2上的抗滑能力最弱，其中，T

1

与第1主应力1

σ(拉正压负)的方向所成的角为

2

4

πs

1

φ

μ

α−

=

=(1) 式中：sφ为极限状态下的内摩擦角。T2与第1主应力1

σ的方向所成的角为

μ

α−

=

2

(2) 因此，若点M在以弧长s为参数的某一滑移线S：

⎩

⎨

⎧

=

=

)

(

)

(

s

y

y

s

x

x(3) 上，则S在点M的切线应沿着T1或者T2的方向，即S应为下列2个常微分方程组之一的1条积分曲线

⎪

⎪

⎭

⎪⎪

⎬

⎫

=

=

=

)2

1

()

(

d

d

)

(

d

d

，

，

，

i

y

x

q

s

y

y

x

p

s

x

i

i

(4)

⎪⎭

⎪

⎬

⎫

=

+

=

+

=

)2

1

(

)

sin(

)

(

)

cos(

)

(

，

，

，

i

y

x

q

y

x

p

i

i

i

i

α

β

α

β

(5)

式中：β为第1主应力1σ的方向角。方程(4)定义了2簇积分曲线：对应于1

=i的曲线属于C1簇，对应于2

=

i的曲线属于C2簇。因为任一滑移线S在通过某一特定介质时与第1主应力1σ的方向所成的角度为常数，所以，也将这种方法称为等倾线法。

第23卷 第5期 郑 宏等. 基于变形分析的边坡潜在滑面的确定 • 711·

在边坡内任取一点M 0，均有2条分别属于C 1

簇的S 1线和C 2簇的S 2线经过点M 0，S 1线与S 2线在点M 0的切线角分别为μβ+和μβ−。因此，这2簇曲线中的任何一簇都充满着整个边坡，而位于塑性区内部的这2簇曲线恰恰构成极限状态下的特征线[20]。

在常微分方程的几何理论中，将满足，

0(x p i =)0y 0)(00=y x q i ，的点)(00y x ，称为方程(4)的奇点，且有如下定理[21]成立：若点)(00y x ，不是奇点，则对于点)(00y x ，的任一充分小的邻域U ，都存在拓扑变换T ：将U 变为矩形V ，将U 内的 1C 簇(或 2C 簇)都变为V 内的平行线。利用这一定理可以得出如下重要结论：若在所研究的区域内的任一点)(y x ，都非奇点，则由方程(4)定义的属于同一簇的积分曲线必不相交。由式(5)可知在边坡内处处都有1)()(22

=+y x q y x p i

i

，即边坡内不存在奇点，所以任意两条 1S 线(或2S )必不相交。利用这一特性，则在沿着某一曲线求其抗滑力和滑动力时，就不必担心路径的分叉问题。

由上述讨论可知，求过一指定点)(000y x M ，的滑移线S ，需指定S 在M 0的切线方向，也就是指定式(5)中的i ，然后再求解下列初值问题：

⎪⎪

⎪⎭

⎪⎪

⎪⎬⎫

====0000)()()(d d )(d d y s y x s x y x q s y

y x p s

x i i ， (6) 由于积分曲线上的每一点都要求有唯一的切线，即要求)(y x p i ，和)(y x q i ，是连续的，而当边坡中含有结构面时，)(y x p i ，和)(y x q i ，在跨越结构面时通常是不连续的，此时滑移线应是由分段光滑曲线组成的。但另一方面，考虑到任何分段光滑曲线都可由光滑曲线来逼近，如果不作任何修改地用式(6)来求解滑移线，得到的曲线将是真实滑移线的平均位置，这条滑移线在穿越变形较大的结构面时会发生明显的扭曲。

3 确定滑移线的数值方法

通过上一节分析即可知道滑移线属于初值问题(6)所定义的2条积分曲线中的1条，但是一般情况下，滑移线所经过的初始点)(00y x ，以及它属于2

簇曲线中的哪1簇却是未知的，好在这2点都可以凭借工程经验，通过设置初始搜索线来加以解决。以图1所示的边坡为例，显然，滑移线应与初始搜索线AB 相交。

图1 设置初始搜索线搜索潜在滑移线 Fig.1 Set initial searching line to search PSL

将AB 剖分成若干个小区间，小区间的端点设为)(00j j y x ，m j ，

L 2 1=，以每1个端点)(00j j y x ，作为方程(4)的初值，都可以得到2条积分

曲线j i S ，2 1，

=i 。而对于每1条j i S ，都可以求得其稳定系数)(w i f j

(抗滑力比滑动力)：

∫

∫=

j

i j

i

S i

S i j

s

s

i f d d )(w ττσ (7)

式中：i τ为沿j i S 切向的剪应力分量；i

στ为对应于

法向应力i σ的抗剪强度，即

s s tan φστσi i

c += (8)

在这2m 条积分曲线中，找出具有最小稳定系

数)(w

i f k 的曲线k i S ，若)(w i f k

接近于1，则k i S 就是要寻找的潜在滑移线，对应的安全系数就是强度折减系数i z 。否则，应重新设置初始搜索线。

无论采用何种数值方法，如Euler 法、预测-校正法(PC 法)或Runge-Kutta 法，来求解式(6)，都要涉及到求域内任一点P 的应力分量，而通过有限元能够得到的仅是单元的Gauss 点应力，因此，需通过某一插值技术来得到点P 的应力分量。对这种插值技术的要求之一是插得的应力点不得位于屈服面之外，一些应力磨平技术是不满足这一要求的。作者推荐下列插值方法来求点P 的应力ij σ：

∑∑===

g

g

1

1

n k k

n k k ij

i

ij w

w σ

σ (9)

式中：g n 为包含点P 的单元内的Gauss 点数目；k

ij

σ1 2

3

A

B

• 712 • 岩石力学与工程学报 2004年

为第k 个Gauss 点的应力；k w 为权系数，权系数可取为

⎩

⎨

⎧=∞≠=−)0( )

0( 2k k k k r r r w ， (10) 式中：k r 为点P 到第k 个Gauss 点的距离。显然，

式(9)是k

ij

σ的一个凸组合，由屈服面的凸性可知由式(9)得到的应力点必位于屈服面内。

若距点P 最近的那个Gauss 点，如第2个Gauss 点，处于塑性状态，则可认为点P 也应处于塑性状态，那么，就应该将上述插值所求得的应力ij σ调整到屈服面上。这种调整的目的在于，当潜在滑移线穿过与Gauss 点2处相对应的一小块塑性区时（如图2所示），使得在相应的这一小段滑移线上的抗 滑力等于滑动力，从而确保当整个潜在滑移线都位

于塑性区内时，其稳定系数k

f w

等于1。调整的方法不是唯一的，本项研究所采用的调整策略是调整后的应力ij

σ′与插值求得的应力ij σ有相同的静水压力。为此令

ij ij ij

s λσδσ+=′ (11) 式中：σ为ij σ的平均应力，3

ii

σσ=

；ij s 为ij σ的

偏应力，ij ij ij s σδσ−=；λ为待定系数。

将式(11)代入Mohr-Coulomb 屈服面方程：

0)()(=+=′ij ij ij

s F F λσδσ (12) 从而可求得参数λ。

除非所选择的步长s 非常小，否则，采用Euler

Fig.2 If the Gauss point 2 is in the plastic state ，adjust the

interpolated point P onto the yield surface

Runge-Kutta 法，因需要做太多的应力插值和调整

运算，影响了搜索速度。所以，可以采用折中方案 ——PC 法来求解式(6)。为此，将式(6)写成向量形式：

⎪⎭

⎪⎬⎫==00)()(d d r s s r r f r (13) 式中：T

T )()(i i q p y x ，==f r 。用梯形积分公式来求(13)式在区间][1+n n s s ，上的积分，得

)]()([211++++=n n n n s

r f r f r r (14)

原则上，可以利用求解非线性方程组的各种方

法来从式(14)中获得1+n r ，但因难以获得Jacobian

矩阵r

f d d ，所以，采用简单迭代法：

⎪⎭

⎪⎬⎫

=++=+=++++) 2 1 0()]()([2)

(11

1

01L ，

m s s m

n n m n n n r f r f r r r f r r n n (15)

式(15)的第1式为预测(prediction)步，第2式则为 校正(correction)步。

为了给出式(15)收敛的充要条件，引入应力旋转度(stress rotation degree)概念。根据对二维应力场的认识，在复杂场合(包括有裂纹存在)的情况，主应力矢量场本身尽管可能存在奇异性，但主应力在方向上的改变却比较均匀。基于上述认识，可以假设主应力1σ与+x 轴之间的夹角)(r μ是一个

函数，即对于域内任何2个位矢r ′和r ′′，L ——应力旋转度，使得下列不等式成 )()(r r ′−′′μμ≤1r r ′−′′L (16)

1为范数，定义为y x +=1

r ，L 的量纲为

。

L ：

e e

L L max = (17)

e L 为单元e 的应力旋转度，可由其内Gauss

1

1

21

22

1max

ig ig ig ig ig ig e

L r r −−=≠μμ (18)

第23卷 第5期 郑 宏等. 基于变形分析的边坡潜在滑面的确定 • 713·

下面讨论使式(15)收敛的充要条件，为此令

)]()([2)(r f r f r r R ++=n n s

(19)

容易证明，当步长s 满足下列不等式：

sL <1 (20)

时，由式(19)定义的)(r R 是一个压缩映象。

事实上，因为 =−=

−112112)()(2

)()(r f r f r R r R s

++−+)cos()cos([2

12i i s

αβαβ ])sin()sin(12i i αβαβ+−+≤ =−−′′+′1212 )cos sin (2

ββββββs s

11212 r r −−sL s μμ

在上述推导过程中，β′和β′′是位于1β和2β之间的2个角度。因此，当步长s 满足式(20)时，)(r R 是

一个压缩映象。

显然，式(15)是)(r R 的一个不动点迭代格式，由压缩映象原理可知，当步长s 满足式(20)时，式 (15)收敛。实算表明，预测后仅需做一次校正即可

满足精度要求。

4 算 例

对于本文所有的算例，假设岩土为理想弹塑性的Mohr-Coulomb 材料，满足关联的流动法则，非线性有限元方程组的求解器是文[19]所建议的位移控制法，弹塑性本构积分为自适应方法[22]。 4.1 均质边坡计算结果与极限分析解的对比

假设有一均质边坡，坡高10 m ，坡度45°，抗剪强度参数kPa 38.12=c ，=φ20°，重度=γ20 kN/m 3

。对于本例，文[23]基于上限法和对潜在滑移

线的对数螺旋线假定，给出的安全系数为1。

根据φ-ν不等式：φsin ≥ν21−，可知ν应满足ν≥0.5)sin 1(φ−≈0.33。为了使φ-ν不等式在对

c 和φ 打折扣的过程中始终成立，作者取 ν = 0.35。

为了验证方法关于网格的非敏感性，作者用了2套网格：Mesh-a 和Mesh-b ，见图3，4所示。2

套网格所取得的安全系数均为1.06，所对应的迭代次数分别为22和27。图5，6分别为2套网格所 对应的潜在滑移线和在极限状态下的广义塑性应变

等值线图。尽管等值线关于网格有一定的依赖性，但2条潜在滑移线却非常接近。从图5，6还可看 出，达到临界状态时，坡角的破损程度最大而坡顶的破损程度最小。

Fig.4 Irregular mesh for the homogeneous slope (mesh-b)

如果使用Mesh-a 并取泊松比=ν0.22( <0.33，不满足φ-ν不等式)，则经过多达54次的迭代后 所求得的安全系数为 1.05，虽然数值能够令人满 意，但其所对应的临界状态下的塑性区(如图7所示)，几乎充满了整个剖面，显然被夸大了，而且此时的潜在滑移线也不再唯一。

Fig.5 PSL and the contours of the generalized plastic strain

(mesh-a)

≤ ≤

Fig.7 Contours of the generalized plastic strain corresponding

to ν = 0.22

为临界平衡态的标准求解了本例，给出的安全系数为1.03，

小于本项研究的结果，这大概是因为FLAC 中所定义的非平衡力未被完全抵消，从而使得边坡并未达到真正的临界平衡态。

4.2 均质边坡计算结果与极限平衡法的计算结果

的对比

图8为均质边坡的有限元模型。假定抗剪强度参数=c 0.058 86 MPa ，=φ11.31°，重度=γ19.62 kN/m 3，弹模=E 80 MPa ，泊松比=ν0.43。文[24]

图8 均质边坡的有限元网格

Fig.8 Finite element mesh of a homogenerous slope

图9 边坡的潜在滑移线及其安全系数

4.3 工程算例之一——水布垭水利工程马崖高边坡

图10是正在建设中的水布垭水利工程的马崖岩质边坡剖面图，该边坡含一陡倾角的大断层F158。该边坡的岩性十分复杂，软硬相间，最硬岩

石的弹性模量是最软岩石的将近40倍，表1列出 了岩性参数。采用强度折减法，同时对E 进行调 整，所算得该边坡的安全系数为1.40。

图11给出了2滑移路径ABC 和ABD ，二者在

层F158的左侧相重合，稳定系数分别为=w

ABC

f 1.00和=w

ABD f 1.25，因此，ABC 才是潜在滑移线。

这是因为断层F158从边坡中部将其切断，如果发

表1 岩体力学参数

Table 1 Rock mass mechanics parameters

抗剪断强度 抗剪断强度 岩性

变模 E / GPa

泊松比

μ

c /MPa

f

重度 /kN ·m

－3

抗拉强度R t / MPa

岩性 变模 E / GPa

泊松比μ

c /MPa

f

重度 /kN ·m

－3

抗拉强度R t / MPa

1 15 0.25 1.

2 1.2 27 1.0 4′

0.4 0.35 0.1 0.4 20 0.1

2 15 0.25 1.0 1.2 27 1.0

2″ 8 0.3 0.6 0.8 26 0

3 5 0.30 0.6 0.8 25 0.5 3″ 2 0.35 0.3 0.5 25 0

4 1 0.3

5 0.3 0.5 22 0.2 2′ 10 0.2

6 0.8 1.0 26 0.5 5 3 0.30 0.5 0.

7 25 0.4 3′ 4 0.3 0.5 0.7 25 0.3 6 10 0.27 0.

8 1.0 25 0.7 F158 — — 0.3 0.5 — — 注：*′表示位于卸荷带内的岩性*，*″表示位于极强卸荷带内的岩性*。

SF = 1.35

第23卷第5期郑宏等. 基于变形分析的边坡潜在滑面的确定• 715·

图11 两条潜在的滑移路径ABC和ABD

Fig.11 Two PSLs ⎯⎯ABC and ABD

生滑坡，首先是F158以下的岩体沿着AB进行滑动，滑动时将无法得到来自F158以上岩体的阻滑作用。图11还显示出，滑移路径ABD基本上沿着强卸荷带与新鲜岩体的交界面，上凸下凹，而拐点正好位于剪切变形较大的剪切带，图12所示的极

Fig.12 Deformation plot of the rock slope in critical

equilibrium state

限状态的变形图可以清楚地反映这一点。有关单位利用经典的极限平衡法所求得的安全系数为1.59，远远大于本文得出的1.40，究其原因就在于，像本例这样复杂的滑移路径，是很难再根据经验而想象得到的。

• 716 • 岩石力学与工程学报 2004年

5 结论

通过本项研究并结合文[18]所建议的基于φ-ν不等式的强度折减技术，建立了完整的基于变形分析的边坡稳定性分析系统，该系统包含边坡稳定性分析时所涉及的2个过程——安全系数的求解和潜在滑面的确定。

本文再次验证在应用基于理想弹塑性模型(Mohr-Coulomb屈服准则)的有限元-强度折减法分析边坡稳定性时，应确保φ-ν不等式成立，否则可能会高估极限状态下的塑性区，而低估安全系数。

诸多算例的考核结果表明，本系统所建立的搜索潜在滑面的方法是非常稳定而可靠的。

算例还表明，对于简单的均质边坡，潜在滑移线更接近于对数螺旋线而不是圆弧，但对于复杂边坡，潜在滑移线的形状非常复杂，会出现多个拐点。致谢本文的工作得到Goodman R E教授宝贵的指导意见及高度的评价：φ-ν不等式体现了“Simplicity is beauty”这一科学意义上的审美原则，在此深表感谢！

参考文献

1郑宏，冯强，罗先启等. 石榴树包滑坡机制的有限元分析[J].

岩石力学与工程学报，(待刊)

2Griffiths D V，Lane P A. Slope stability analysis by finite elements[J].

Geotechnique，1999，49(3)：387～403

3Lechman J B，Griffiths D V. Analysis of the progression of failure of earth slopes by finite elements[A]. In：Griffith D V，Fenton Gordon A，Martin Timothy R，ed. Slope Stability 2000：Proceedings of Sessions of Geo-Denver 2000[C]. New York：Geotechnical Special Publication，ASCE，2000，250～265

4邵国建，卓家寿，章青. 岩体稳定性分析与评价准则研究[J]. 岩石力学与工程学报，2003，22(5)：691～696

5Johnson C. Existence thermos for plasticity problems[J]. J. Math.

Pures Appl.，1976，55：431～444

6Duncan J M. State of the art: limit equilibrium and finite-element analysis of slopes[J]. J. Geotech. Engrg. Div.，ASCE，1996，122(7)：577～596

7连镇营，韩国城，孔宪京. 强度折减有限元法研究开挖边坡的稳

定性[J]. 岩土工程学报，2001，23(4)：407～411

8宋二祥. 土工结构安全系数的有限元计算[J]. 岩土工程学报，1997，19(2)：1～7

9Zou J Z，Williams D J. Search for critical slip surface based on finite element method[J]. Canadian Geotechnical Journal，1995，32(1)：233～246

10王成华，夏绪勇，李广信. 基于应力场的土坡临界滑动面的蚂蚁算法搜索技术[J]. 岩石力学与工程学报，2003，22(5)：813～819

11史恒通，王成华. 土坡有限元稳定分析若干问题的探讨[J]. 岩土力学，2000，21(2)：152～155

12赵尚毅，郑颖人，时卫民等. 用有限元强度折减法求边坡稳定安全系数[J]. 岩土工程学报，2002，24(3)：343～346

13Matsui T，San K C. Finite element slope stability analysis by shear strength reduction technique[J]. Soils and Foundations，1992，32(1)：59～70

14Dawson E M，Roth W H，Drescher A. Slope stability analysis by strength reduction[J]. Geotechnique，1999，49(6)：835～840

15Dawson E，Motamed F，Nesarajah S，et al. Geotechnical stability analysis by strength reduction[A]. In：Griffiths D V，Fenton Gordon A，Martin Timothy R，ed. Slope Stability 2000：Proceedings of Sessions of Geo-Denver 2000[C]. New York：Geotechnical Special Publication，ASCE，2000，99～113

16Dawson E，You K，Park Y. Strength-reduction stability analysis of rock slopes using the Hoek-Brown failure criterion[A]. In：Labuz J F，Glaser S D，Dawson E，ed. Trends in Rock Mechanics[C]. New York：Geotechnical Special Publication，ASCE，2000，65～77

17Jeremić B. Finite element methods for 3D slope stability analysis[A].

In：Griffith D V，Fenton Gordon A，Martin Timothy R，ed. Slope Stability 2000：Proceedings of Sessions of Geo-Denver 2000[C]. New York：Geotechnical Special Publication，ASCE，2000，224～239 18郑宏，李春光，李焯芬等. 求解安全系数的有限元法[J]. 岩土工程学报，2002，24(5)：626～628

19郑宏，李焯芬，谭国焕等. 有限元分析的位移控制法及其应用[J].

岩土工程学报，2002，24(1)：81～85

20沈珠江. 理论土力学[M]. 北京：中国水利水电出版社，2000

21张锦炎. 常微分方程几何理论与分支问题[M]. 北京：北京大学出版社，1981

22郑宏. 岩土力学中的几类非线性问题[博士学位论文][D]. 武汉：中国科学院武汉岩土力学研究所，2000

23Chen W F. Limit analysis and soil plasticity[M]. Amsterdam：Elsevier，1975

24邓建辉. CAD技术在边坡稳定性分析中的应用[硕士学位论文][D].

武汉：中国科学院武汉岩土力学研究所，1992