边坡稳定分析的极限分析下限解有限元法

sin 2 θ s
cos2 θ s
− sin 2 θ s
{ } { } σ s T = σ x1 ,σ y1 ,τ xy1 ,σ x2 ,σ y2 ,τ xy2
L 为边界 s 的长度，θ s 为边界 s 与 x 轴的夹角，
（ σ n1 ， σ n2 ）为作用于两端点的未知应力。 对于边坡、隧洞或地基稳定等问题，更为简便
收稿日期：2004-04-30 作者简介：王雪涛，女，1977 年生，硕士，2003 年毕业于中国水利水电科学研究院，主要从事岩土数值计算方面的工作。
增刊
王雪涛等：边坡稳定分析的极限分析下限解有限元ig. 1 3-Noded linear stress triangle
2 下限解有限元法
2.1 基本原理 下限解有限元法以塑性力学下限定理为基础。
下限定理认为，在不违反屈服准则的条件下，对应 于任何一个静力许可的应力场都可以求得一个相应 的外荷载，这个外荷载总是小于或等于极限荷载， 是极限荷载的下界（ 或‘安全’）的估计。一个静 力许可的应力场要满足应力边界条件、应力平衡及 屈服条件（应力空间的应力必须位于屈服面内）。在 下限解有限单元法的求解过程中未知量为结点应 力。
摘 要： 研制开发了极限分析下限解有限元程序。在程序的开发过程中，着重解决了下限解有限元法转化为标准线性规
划问题和线性优化方法的选择问题，同时，通过典型算例的分析对比，对程序的合理性和可行性进行了验证。
关 键 词： 极限分析；下限解；有限元法；线性规划
中图分类号：Tu 443
文献标识码：A
Lower bound limit analysis of slope stability using finite element method
Fig 3 Stress boundary conditions （a）boundary load (b) boundary condition
2.5 应力边界条件 设在边界单元上作用有相应的荷载（ q1 ， t1 ）
和（ q2 ， t2 ），如图 3 所示。
为了使应力边界条件得到满足，必须对边界上
σ n1 = σ n2 ； σ n3 = σ n4 ；τ1 = τ 2 ；τ 3 = τ 4 （5）
将式（4）代入式（5）中，平衡条件变为
[ ] { } { } Ad equil
σd
=
bd equil
（6）
式中
[ ] Ad equil
=
⎡T
⎢ ⎣
0
−T 0
0 T
0⎤
−
T
⎥ ⎦
其中
⎡ sin 2 θ
T
第 25 卷增刊 2004 年 11 月
文章编号：1000－7598－(2004)增－0134－05
岩土力学 Rock and Soil Mechanics
Vol.25 Supp.2 Nov. 2004
边坡稳定分析的极限分析下限解有限元法
王雪涛 1, 汪小刚 1
(1. 中国水利水电科学研究院，北京 100044）
一系列的约束条件，以满足应力平衡条件、应力边
界条件以及屈服准则。
2.3 应力平衡条件
对图 2 所示的单元，将式（1）代入平面应变条
件下的应力平衡条件中，可得到满足应力平衡条件
的结点应力约束条件为
[ ] { } { } Ae equil
σe
=
be equil
（3）
式中
[ ] A
e equil
=
1 2Ae
由 Sloan 等学者提出的有限单元法上、下限解 [2..3]弥补了极限平衡法和极限分析上限法的缺陷，它 借鉴了传统有限元考虑解的应力-应变关系的作法， 并从上、下限两个方向逼近问题的真实解，由于极 限分析上限解已取得一定的进展，本文通过总结极
限分析法这一研究领域的相关成果，研制开发了极 限分析下限解有限元程序。在程序的开发过程中， 着重解决了下限解有限元法转化为标准线性规划问 题和线性优化方法的选择问题，同时，通过典型算 例的分析对比，对程序的合理性和可靠性进行了验 证。
1引言
自从上世纪 20 年代瑞典圆弧法出现以来，基于 极限平衡理论的条分法始终是边坡稳定分析的主要 手段。但极限平衡法为使问题变得静定可解，引入 了大量关于条间作用力的位置和方向的假定。基于 塑性力学上限原理的极限分析方法相对于极限平衡 法而言，更具有严密的理论基础，并逐步得到了广 泛的应用。极限分析上限法[1]的要点是不关心结构 加载破坏的全过程，只是通过简单考虑土体的应力– 应变关系，直接求得结构所能承受的极限荷载，但 这一方法得到问题的上限解答与现今基于安全或偏 保守的工程设计理念有所冲突。
x32 = x3 − x2 ，
y23 = y2 − y3 ；
x13 = x1 − x3 ，
y31 = y3 − y1 ；
x21 = x2 − x1 ，
y12 = y1 − y2 。
2 A = x13 y23 − x32 y31
下限解有限元法中每一个单元的任一结点仅属
于该单元自己。
为了求得严格的下限解，必须对结点应力附加
在上述简化条件下满足 F≤0 的条件，即为
Fki ≤0，i=1，2，3 表示单元结点， k = 1,2,L, p 代 表多边形各边。此时满足屈服准则的约束条件，可
进一步表示为
[ ] { } { } Ai yield
σi
≤
bi yield
（10）
这里
{ } bi yield
= 2ci cosϕ cos(π
=
⎢ ⎢−
1
sin 2
θ
⎣2
cos 2 θ 1 sin 2 θ 2
− sin 2 θ ⎤ ⎥
cos2 θ ⎥ ⎦
{ } { } σ d T = σ x1 ,σ y1 ,τ xy1 ,L,σ x4 ,σ y4 ,τ xy4
{ } { } bd equil
T=
0
0
0
0
136
岩 土 力学
2004 年
图 3 应力边界条件（a）边界荷载；（b）边界条件
式中
Ak = cos(2π k p) + sin ϕ cos(π p) Bk = sin ϕ cos(π p) − cos(2π k p)
Ck = 2 sin(2π k p) k = 1,2,L, p
一般来讲，采用 12 个边的多边形来简化可达到
足够的精度。当然在土体摩擦角较大时，宜采用更
多边的多边形。
2.2 单元形式和离散模式 下限解有限元法采用如图 1 所示的三角形单元
形式，单元内部应力呈线性变化，与结点应力的关 系可表示为
3
3
3
∑ ∑ ∑ σ x = Niσ xi ； σ y = N iσ yi τ xy = N iτ xyi
i =1
i =1
i =1
（1） 其中 σ xi ，σ yi ，τ xyi 分别为结点应力：Ni 为形函数， 可进一步用单元结点坐标表示为
WANG Xue-tao1, WANG Xiao-gang1
(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
服准则，都位于屈服面内，即屈服函数 F≤0。 下限解的屈服准则采用 Mohr-Coulomb 准则，
为避免出现非线性的约束条件，在下限解中用 p 个 边的内接等边多边形来简化 Mohr 圆。在图 4 中， 采用 6 边内接多边形简化，对于第 k 个面的 Mohr-Coulom 准则可表示为
Fk = Akσ x + Bkσ y + Ckτ xy + 2c cosϕ cos(π p) （9）
Q
=
h∫
σ
s
n
ds
（11）
Q 为极限荷载，假定未知应力呈线性分布，并
利用式（7），上式可用矩阵表示为
Q
=
CT{σ s} =
Lh 2
(σ
n1
+ σ n2 )
（12）
增刊
王雪涛等：边坡稳定分析的极限分析下限解有限元法
137
式中
{ } CT
=
L 2
sin 2 θ s
cos2 θ s
− sin 2 θ s
[T
]
=
⎡ ⎢⎣
sin 2
1 2
sin
θ 2θ
cos 2 θ
1 2
cos
2θ
− sin 2θ ⎤ cos 2θ ⎥⎦
{ } { } bl bound
T=
q1
t1
q2
t2
{ } { } σ l T = σ x1 ,σ y1 ,τ xy1 ,σ x2 ,σ y2 ,τ xy2
2.6 屈服条件 对于严格的下限解，每一点的应力应不违反屈
p){1,1,L,1,L,1}；
{ } { } σi T = σ xi ,σ yi ,τ xyi ；
⎡ A1 B1 C1 ⎤
[ ] Ai yield
=
⎢ ⎢ ⎢ ⎢ ⎢
A2 M Ak M
B2 M Bk M
C2 M Ck M
⎥ ⎥ ⎥ ⎥ ⎥
；
⎢⎣ Ap Bp C p ⎥⎦
ci 代表不同结点处的凝聚力。
（4）
图 2 应力连续要求 Fig. 2 Statically admissible stress continuity between
adjacent triangles
图 2 为相邻单元应力不连续的典型情况。①、 ②；③、④两对结点具有相同的坐标，由于假定了 应力呈线性变化，为了满足平衡条件必须使相邻单 元共轭结点的正应力和剪应力相等。即
的方法是把单位重量作为未知量，通过寻找与最大 的单位重量相对应的静力许可的应力场，而直接获 得问题的解。在这种情况下，目标函数中的向量 C 具有更为简单的形式。
是其它条件的等式约束矩阵。
对大多数常用的线性规划问题的求解方法，都
需要将上述的一般形式变换为以下的标准形式
min s = c T x, x ∈ E n
图 4 Mohr-Coulomb 屈服面从内部线性化 Fig .4 Linearized Mohr-Coulomb yield function
2.7 目标函数 对于下限问题，我们的目标是要寻找一个静力
许可的应力场，它使得作用于某一段边界 s 上的正 应力所形成的总荷载达到最大，如图 5 所示。如果 用 h 表示厚度，则
的单元结点应力增加相应的约束条件
σ n1 = q1； τ 1 = t1；
σ n2 = q2 τ 2 = t2
（7）
将式（4）代入式（7）中可得到约束条件的矩阵表 达形式
[ ] { } { } A σ = b l
lT
l
T
bound
bound
（8）
式中
[ ] Al bound
=
⎡T ⎢⎣0
0⎤ T ⎥⎦
⎡y23 ⎢⎣ 0
0 x32
x32 y23
y31 0
0 x13
x13 y31
y12 0
0 x21
x21⎤ y12⎥⎦
{ } { } σe T = σx1,σy1,τxy1,L,σx3,σy3,τxy3
{ } { } be equil
=
0,γ
2.4 应力不连续处的平衡条件 为了保证在应力不连续处满足静力许可的条
[ ] N1 = (x2 y3 − x3 y2 ) + y23 x + x32 y / 2A； [ ] N 2 = (x3 y1 − x1 y3 ) + y31 x + x13 y / 2 A； [ ] N 3 = (x1 y2 − x2 y1 ) + y12 x + x21 y / 2A。
（2）
式中
件，必须要求作用在不连续面上的正应力和剪应力
的连续性。对于图 2 所示的情况，作用在与 x 轴夹 角为θ的平面上的正应力和剪应力可以用单元应力
分量表示为
σ n = sin 2 θσ x + cos 2 θσ y − sin 2 θτ xy
τ
=
−
1 2
sin
2
θσ
x
+
1 2
sin
2
θσ
y
+
cos 2 θτ xy