第25卷增刊 岩 土 力 学 Vol.25 Supp.2 2004年11月 Rock and Soil Mechanics Nov. 2004
文章编号：1000－7598－(2004)增－0134－05
边坡稳定分析的极限分析下限解有限元法
王雪涛, 汪小刚1
1
(1. 中国水利水电科学研究院，北京 100044）
摘 要： 研制开发了极限分析下限解有限元程序。在程序的开发过程中，着重解决了下限解有限元法转化为标准线性规划问题和线性优化方法的选择问题，同时，通过典型算例的分析对比，对程序的合理性和可行性进行了验证。 关 键 词： 极限分析；下限解；有限元法；线性规划 中图分类号：Tu 443 文献标识码：A
Lower bound limit analysis of slope stability using finite element method
WANG Xue-tao, WANG Xiao-gang1
1
(China Institute of Water Resources and Hydropower Research, Beijing, 100044, China)
Abstract ：The procedure of the finite-element method of lower bound. is deve loped. The two main problems of the realization of the finite-element method. one is the transformation of standard linear programming and the other is the choice of optimization, are solved well in this procedure The rationality and reliability of this procedure are examined and certified by the typical examples.
Key words：limit analysis; lower bound; finite element method; linear proqramming
1 引 言
自从上世纪20年代瑞典圆弧法出现以来，基于极限平衡理论的条分法始终是边坡稳定分析的主要手段。但极限平衡法为使问题变得静定可解，引入了大量关于条间作用力的位置和方向的假定。基于塑性力学上限原理的极限分析方法相对于极限平衡法而言，更具有严密的理论基础，并逐步得到了广泛的应用。极限分析上限法[1]的要点是不关心结构加载破坏的全过程，只是通过简单考虑土体的应力–应变关系，直接求得结构所能承受的极限荷载，但这一方法得到问题的上限解答与现今基于安全或偏保守的工程设计理念有所冲突。
由Sloan 等学者提出的有限单元法上、下限解
[2..3]
限分析法这一研究领域的相关成果，研制开发了极限分析下限解有限元程序。在程序的开发过程中，着重解决了下限解有限元法转化为标准线性规划问题和线性优化方法的选择问题，同时，通过典型算例的分析对比，对程序的合理性和可靠性进行了验证。
2 下限解有限元法
2.1 基本原理
下限解有限元法以塑性力学下限定理为基础。下限定理认为，在不违反屈服准则的条件下，对应于任何一个静力许可的应力场都可以求得一个相应的外荷载，这个外荷载总是小于或等于极限荷载，是极限荷载的下界（ 或‘安全’）的估计。一个静力许可的应力场要满足应力边界条件、应力平衡及屈服条件（应力空间的应力必须位于屈服面内）。在下限解有限单元法的求解过程中未知量为结点应力。
弥补了极限平衡法和极限分析上限法的缺陷，它
借鉴了传统有限元考虑解的应力-应变关系的作法，并从上、下限两个方向逼近问题的真实解，由于极限分析上限解已取得一定的进展，本文通过总结极收稿日期：2004-04-30
作者简介：王雪涛，女，1977年生，硕士，2003年毕业于中国水利水电科学研究院，主要从事岩土数值计算方面的工作。
增刊 王雪涛等：边坡稳定分析的极限分析下限解有限元法
135
式中
⎡y 0
[A ]=21
A ⎢⎣0x
e
equil e T
23e
32
x 32y 310x 13y 120x 21⎤y 230x 13y 310x 21y 12⎥⎦
{σ}={σ, σ, τ
{b }={0, γ}
x1
e equil
y1xy1
, L , σx3, σy3, τxy3}
2.4 应力不连续处的平衡条件
为了保证在应力不连续处满足静力许可的条件，必须要求作用在不连续面上的正应力和剪应力的连续性。对于图2所示的情况，作用在与x 轴夹角为θ的平面上的正应力和剪应力可以用单元应力
图1 下限解采用的单元形式 Fig. 1 3-Noded linear stress triangle
分量表示为
σn =sin 2θσx +cos 2θσy −sin 2θτxy
（4）
11222
τ=−sin θσx +sin θσy +cos θτxy
22
2.2 单元形式和离散模式
下限解有限元法采用如图1所示的三角形单元形式，单元内部应力呈线性变化，与结点应力的关系可表示为
3
3
3
σx =∑N i σx i ；σy =∑N i σy i τxy =∑N i τxy i
i =1
i =1
i =1
（1）
其中σx i ，σy i ，τxy i 分别为结点应力：N i 为形函数，可进一步用单元结点坐标表示为
N 1=[(x 2y 3−x 3y 2) +y 23x +x 32y ]/2A ；N 3=[(x 1y 2−x 2y 1) +y 12x +x 21y ]/2A 。
图2 应力连续要求
Fig. 2 Statically admissible stress continuity between
adjacent triangles
（2） N 2=[(x 3y 1−x 1y 3) +y 31x +x 13y ]/2A ；
式中 图2为相邻单元应力不连续的典型情况。①、
x 32=x 3−x 2， y 23=y 2−y 3； x 13=x 1−x 3 ， y 31=y 3−y 1； x 21=x 2−x 1， y 12=y 1−y 2。
2A =x 13y 23−x 32y 31
下限解有限元法中每一个单元的任一结点仅属于该单元自己。
为了求得严格的下限解，必须对结点应力附加一系列的约束条件，以满足应力平衡条件、应力边界条件以及屈服准则。 2.3 应力平衡条件
对图2所示的单元，将式（1）代入平面应变条件下的应力平衡条件中，可得到满足应力平衡条件的结点应力约束条件为
②；③、④两对结点具有相同的坐标，由于假定了应力呈线性变化，为了满足平衡条件必须使相邻单元共轭结点的正应力和剪应力相等。即
σn1=σn2 ；σn3=σn4 ；τ1=τ2 ；τ3=τ4 （5）
将式（4）代入式（5）中，平衡条件变为
[A ]{σ}={b } （6）
d
equil
d
d equil
式中
T
[A ]=⎡⎢0
d
equil
−T 0
0T
⎣
0⎤−T ⎥⎦cos 2θ1
sin 2θ2
−sin 2θ⎤
⎥2
cos θ⎥
⎦
其中
⎡sin 2θT =⎢1
⎢−sin 2θ⎣2
[A ]{σ}={b } （3）
e equil
e
e equil
{σ}={σ{b }={0
d T T d
equil
x1
, σy1, τxy1, L , σx4, σy4, τxy4}000}
136
岩 土 力 学 2004年
C k =2sin(2πk p ) k =1, 2, L , p
一般来讲，采用12个边的多边形来简化可达到足够的精度。当然在土体摩擦角较大时，宜采用更多边的多边形。
在上述简化条件下满足F ≤0的条件，即为
F k i ≤0，i =1，2，3表示单元结点，k =1, 2, L , p 代
图3 应力边界条件（a ）边界荷载；（b ）边界条件
Fig 3 Stress boundary conditions
（a ）boundary load (b) boundary condition
表多边形各边。此时满足屈服准则的约束条件，可进一步表示为
[A ]{σ}≤{b } （10）
i yield
i
i yield
2.5 应力边界条件
设在边界单元上作用有相应的荷载（q 1，t 1）和（q 2，t 2），如图3所示。
为了使应力边界条件得到满足，必须对边界上的单元结点应力增加相应的约束条件
这里
{b }=2c cos ϕcos(πp ) {1, 1, L , 1, L , 1}；
i
yield
i
{σi }T ={σx i , σy i , τxy i }；
⎡A 1⎢A 2⎢M =⎢A ⎢k ⎢M ⎢⎣A p
B 1
B 2M B k M B p
C 1⎤C 2⎥M ⎥； C k ⎥⎥M ⎥C p ⎥⎦
σn1=q 1；σ=q 2
n2 （7）
τ1=t 1；τ2=t 2
将式（4）代入式（7）中可得到约束条件的矩阵表达形式
[A ]
i
yield
[A
式中
l
bound
]{σ}={b
l T
T l
bound
} （8）
c i 代表不同结点处的凝聚力。
T l
[A bound ]=⎡
⎢0
⎣⎣0⎤T ⎥⎦
1
sin 2θ[T ]=⎡1
⎢sin 2θ
l
{b bound }T ={q 1
cos 2θ−sin 2θ⎤
cos 2θcos 2θ⎥⎦
t 1q 2t 2}
{σ}={σ
l T
x1
, σy1, τxy1, σx2, σy2, τxy2}
图4 Mohr-Coulomb屈服面从内部线性化
Fig .4 Linearized Mohr-Coulomb yield function
2.6 屈服条件
对于严格的下限解，每一点的应力应不违反屈服准则，都位于屈服面内，即屈服函数F ≤0。
下限解的屈服准则采用Mohr-Coulomb 准则，为避免出现非线性的约束条件，在下限解中用p 个边的内接等边多边形来简化Mohr 圆。在图4中，采用6边内接多边形简化，对于第k 个面的Mohr-Coulom 准则可表示为
F k =A k σx +B k σy +C k τxy +2c cos ϕcos(πp ) （9）
2.7 目标函数
对于下限问题，我们的目标是要寻找一个静力许可的应力场，它使得作用于某一段边界s 上的正应力所形成的总荷载达到最大，如图5所示。如果用h 表示厚度，则
Q =h ∫s σn d s （11）
Q 为极限荷载，假定未知应力呈线性分布，并
式中
A k =cos(2πk p ) +sin ϕcos(πp ) B k =sin ϕcos(πp ) −cos(2πk p )
利用式（7），上式可用矩阵表示为
Q =C T {σs }=
Lh
(σn1+σn2) （12） 2
增刊 王雪涛等：边坡稳定分析的极限分析下限解有限元法
137
式中
C T =
L
sin 2θs 2
是其它条件的等式约束矩阵。
{
cos 2θs −sin 2θs sin 2θs cos 2θs −sin 2θs
}
对大多数常用的线性规划问题的求解方法，都需要将上述的一般形式变换为以下的标准形式
min s =c T x ,
{σ}={σ
s T
x1
, σy1, τxy1, σx2, σy2, τxy2}
x ∈E n
L 为边界s 的长度，θs 为边界s 与x 轴的夹角，
满足于
A x =b , χ≥0，b ≥0， (15)
式（14）与式（15）相比，不同之处在于：①结点应力分量{σ}不是非负变量，可正可负；②约束条件中存在不等式。
因此，将式（14）转换为式（16）所示的标准
{σ−}和松弛变量形式时，需要引入非负变量{σ+}、
{s }，变换后的形式如下式所示：
（σn1，σn2）为作用于两端点的未知应力。
对于边坡、隧洞或地基稳定等问题，更为简便的方法是把单位重量作为未知量，通过寻找与最大的单位重量相对应的静力许可的应力场，而直接获得问题的解。在这种情况下，目标函数中的向量C 具有更为简单的形式。
最小值 {c *}
∗∗
T
{x *}
∗
∗
[A ]{x }={b }
约束条件 （17）
{x }≥0
式中各向量定义如下：
{c *}T ={−{c }T {c }T {0}T }；{x *}T ={{σ+}T {σ−}T {s }T }； {b *}T ={{b 1}T {b 2}T }。
图5 作用在边界上的荷载
Fig. 5 Load in a direction normal to a boundary edge
[A ]−[A ][I ]⎤
[A ]=⎡⎢[A ]−[A ][0]⎥
∗
1
1
⎣
22
⎦
其中
2.8 下限解有限元法对应的线性规划问题
(1) 线性规划问题的标准化
{σ}={σ+}−{σ−}
{σ+}{, σ−}, {s }≥{0}
一旦建立了单元结点应力的约束条件和下限解问题的目标函数，下限解的求解就可归结为如下的线性规划问题：
约束条件 最大值 C T σs
{s }是施加给不等式约束的松弛变量；[I ]是单位矩
阵。
(2) 优化方法的选择
{}
经过标准化之后，就可以进行优化了，本次研究分别尝试了两种优化方法，一种是目前在优化计
e
equil
[A ]{σ}={b }；[A ]{σ}={b }；
e equil d equil
e d
d equil
算中常用的单纯形法；一种是直接采用数学软件MATLAB 中的优化方法。单纯形法是一种很成熟并
（13）
l l
[A bound ]{σl }T ={b bound }；i i
}。[A yield ]{σi }≤{b yield
且应用很广的优化方法，其具体算法可参见文献[4]。下面将数学软件MATLAB 中的优化方法作一
简单介绍。
针MATLAB 是一个大型的计算分析工具软件，对上述标准的线性规划问题，可直接调用函数linprog 来求解。本次研究采用的linprog 函数形式
上述线性规划问题可转换为如下的一般形式： 最小值 约束条件
−{c }
T
{σ}
[A 1]{σ}≤{b 1}；
(14)
[A 2]{σ}={b 2}；
式中，[A 1]为屈服条件的不等式约束矩阵，[A 2]
如下：
函数 linprog
格式 [x，fval] = linprog(f，Aeq ，beq ，lb ，ub ，options) % lb和ub 为指定x 的范围lb ≤x ≤ub ，Aeq
138
岩 土 力 学 2004年
和beq 为等式约束矩阵，fval 为返回目标函数最优值，即fval=f x ，options 为指定的优化参数。
'
经由软件MATLAB 的优化方法得到的极限荷
载为109.4 365 kPa，则其加载系数为
3 算 例
3.1 算例1
η=
T −T 0T 0
=
109. 4365−111. 437
=−0. 0179
111. 437
本算例是一个承受均布垂直荷载的半无限、无重量的条形地基。地基为均一材料，材料参数为：凝聚力c ＝30 kPa，土的重度γ＝0.0 kN/m ，摩擦角分别取ϕ＝0°，5°，10°，15°，20° 5种情况进行讨论。
具体有限元网格划分如图6，单元总数为12个，总结点编号为36个，独立的坐标结点总数为11个，
简化摩尔圆的等边多边形的边数p ＝6。
3
其中T 0为索科洛夫斯基（Sokolovski ）给出的理论解，其值为111.437 kPa。
而单纯形法则认为，该问题为无解问题。
4 结 论
本文在认真总结极限分析方法这一研究领域相关成果的基础上，以下限解有限元法为研究对象，研制开发了相应的计算机程序。并选取土力学中两个典型的算例对程序的可靠性以及优化方法的适用性进行了验证。计算结果表明，单纯形法由于要预先设定初值，很难适用于下限解有限元法，而采用数学软件MATLAB 中的优化方法得到了令人满意的结果。
当然，由于目前优化方法发展很快，是否能够找到更加适用于下限解有限元法的优化方法是今后值得进一步深入研究的问题。
此外，本次开发的程序仅仅通过了两个典型算例的验证，如何将该程序推广应用于更为复杂的实
图6 算例1单元划分模式
Fig .6 Element meshes for example 1
计算结果见表1。
表1 算例1的计算结果
ϕ /(°) 5 10 15 20 14.3
际工程问题，也是今后应进一步研究的重点。
参 考 文 献
[1] Donald,I and Chen,Z.Y ..1997. Slope stability analysis by
the upper bound approach: fundamentals and methods.[J] Canadian Geotechnical Journal, 1997,34:853-862 [2] Sloan S.W.. Lower bound limit analysis using finite
elements and linear programming[J]. International journarl for Numerical and Analytical Methods in Geomechanics, 1989,12：61-67.
[3] Sloan S. W.. Upper bound limit analysis using finite
elements and linear programming[J]. International journal for Numerical and Analytical Methods in Geomechanics, 1989, 13: 263-282.
[4] 陈祖煜，边坡稳定分析的塑性力学上限解和下限解
[A]，中国水利电力技术发展和成就[C]. 北京：中国电力出版社，1997，198-208.
[5] 龚晓南主编，土工计算机分析[M]，中国建筑工业出版
社，2000.
[6] M. J. Best and K. Ritter ,Linear Programming: Active Set
Analysis and Computer Programs[M]. Prentice-Hall, New Jersey，1985.
N cN 5.14 6.49 8.35 11.0
N cL 5.0 6.25 7.94 10.20 13.32 N ris 5.02 7.12 10.77 20.75 114.387 2 N MAT
4.9627 6.213 7.915 10.156 13.238
注：N cN 和N cL 分别为极限分析上下限解的理论值，N ris 为下限解有限元采用单纯形法优化方法得到的下限解，N MAT 为下限解有限元采用数学软件MATLAB 中的优化方法得到的下限解
图7 算例2单元划分模式
Fig. 7 Element meshes for example 2
3.2 算例2
本算例为作用有垂直荷载的无重量的均质边坡，材料参数为：凝聚力c ＝98 kPa，土的重度
γ＝0.0 kN/m ，摩擦角ϕ＝30°。将假定的滑坡体
分成12个单元，如图7所示。
3