第9章 有限元法在边坡稳定分析中的应用
9,1 概述
随着计算机软件硬件的飞速发展采用理论体系更为严格的方法进行边坡稳定分析已经成为可能有限单元法全面满足了静力许可应变相容和应力应变之间的本构关系同时因为是采用数值分析方法可以不受边坡几何形状的不规则和材料的不均匀性的限制因此应该是比较理想的分析边坡应力变形和稳定性态的手段
与传统的极限平衡法相比边坡稳定分析的有限元法的优点可总结如下
(1) 破坏面的形状或位置不需要事先假定破坏自然地发生在土的抗剪强度不能抵抗剪应力的地带
(2) 由于有限元法引入变形协调的本构关系因此也不必引入假定条件保持了严密的理论体系
(3) 有限元解提供了应力变形的全部信息
有关有限元的基本原理及其在边坡稳定分析中的应用已有很多的文献沈珠江1990
陈胜宏2001介绍作者也另有专著陈祖煜等2003介绍了采用有限元法求解第6章中所列渗流固结和应力应变分析控制方程的详细方法本章9.2节简要介绍这一计算方法的解题思路以及三角形四边形单元的离散方法并介绍有限元法在小浪底和务平工程中的应用
如何将有限元计算成果与传统的安全系数挂钩成为直接用于边坡设计的判别依据也是广泛受到重视的课题近期在这方面获得了有意义的成果本章将简要介绍这方面的成果
9,2 求解渗流和应力应变控制方程的有限元方法
9,2,1 基本原理
第6章所述求解渗流和应力应变的数学提法是在体力{fb}和边界上面力{T }的作用下求满足微分方程式(6.27)式(6.28)式(6.31)式(6.40)和边界条件式(6.30)式(6.43)
式(6.44)式(6.46)和式(6.48)的位移场{W}和水头场h用有限元求解这些偏微分方程组边值问题其基本途径是找出一个泛函π({W} h)当{W}和h为该偏微分方程边值问题解时
π获得极值
可以证明所寻找的泛函由下式确定
240 土质边坡稳定分析?原理
方法
程序
(9.1) {}()
7654321
,ππππππππ ++++++=hW
其中
∫∫
′′=
VV
T
VV
0
1
d}{}{d}{}{
2
1
εσεσπ (9.2)
(9.3) {}
∫
=
V
T
w
Vayh
2
d}{)( εγπ
{} {}{}
∫
=
V
T
w
Vhkh
3
d**
2
1
γπ (9.4)
∫
=
V
w
Vhhh
Q
0
2
4
d)2(
2
1 γ
π (9.5)
(9.6) {} {}
∫
=
V
T
b
VWf
5
dπ
{} {}
∫
=
1
d
6
s
T
sWTπ 9.7)
(9.8) shq
s
w
d
4
7
∫
= γπ
这一命题的理论证明可参见文献陈祖煜等2003
9,2,2 有限元解法
采用数值分析方法求解使获得极值的水头和位移场包括以下步骤
1,离散化
将所研究的域V分成n个单元共m个节点这m个节点的{W}用{v}来代表即
(9.9) ),,,,,,(}{
2211
T
ymxmyxyx
WWWWWW LL=ν
同样用{?}代表这m个节点的{h}值
(9.10) ),,,(}{
21
T
m
hhh LL=?
任一点的{W}和h可用该点所属单元的节点{W}
e
和{h}
e
近似表示本质上也就是可用
{ν}和{?}代表于是π可以近似地用{ν}和{?}来代表即
(9.11) }){},({)},({?νππ =hW
根据里兹法的原理使π取得极值的{ν}和{?}满足
Ο
Ο
对于某一行向量{a}
T
我们定义C/{a}
T
为一列向量{b}它满足
第9章 有限元法在边坡稳定分析中的应用 241
{}
{}0
2
1
1
=
=
ym
x
y
x
T
W
W
W
W
π
π
π
π
ν
π
M
M
(9.12)
{}
}0{
2
1
=
=
m
T
h
h
h
h
π
π
π
π
M
M
(9.13)
其中{0}为元素均为零的向量由式(9.12)和式(9.13)可得3m个线性方程可用来求解由
{ν}和{h}所包含的3m个未知数这就是有限单元法的基本原理
2,应用形状函数表达单元内的物理量
单元内任一点的{W}和h可用该单元节点的{W}
e
和{h}
e
来近似表达
(9.14)
e
W
WNW }]{[}{ =
9.15)
T
h
eTe
h
NhhNh ][}{}]{[ ==
因此
(9.16)
e
W
e
W
TT
WBWNW }]{[}]{[][}{][}{ =?=?=ε
9.17)
e
h
e
h
hBhNh }]{[}]{}[{}{ =?=?
式中
(9.18) ][][][
W
T
W
NB?=
9.19) [] ]}[{
hh
NB?=
[NW] [Nh]称形状函数或插值函数对三角形和四边形单元具有不同表达形式将在第9.2.3节中详述同时我们记
C={a}
T
{b}或 C={b}
T
{a}
242 土质边坡稳定分析?原理
方法
程序
(9.20) ][}{][
W
TT
BaB =
将式(9.14)式(9.15)式(9.16)式(9.17)和式(9.20)代入式(9.1)得
{}[]{}
∑
∫
∑
∫
∑
∫
==
=
++
+
+
+
′=
n
i
e
h
T
h
eTe
W
T
W
S
i
s
eT
e
h
T
h
eTweT
h
T
h
eTw
e
W
T
b
e
h
T
h
e
eT
w
eTT
h
eT
e
W
T
n
i
V
e
W
T
W
eT
shNNqsWNNT
VhNNh
Q
hNNh
Q
WfhNKNh
WByWBNh
WBWBCBW
w
1
s
1
0
22
T
0
1
d}{*][][}{*d}{*][][}{
d)}{*][][}{ }{*][][}{
2
1
}{*][}{-}{*]][[][*}{
2
1
}{*][}{*][][}{
*}]{][[][*}{
2
1
(
41
γ
γγ
γ
γγ
σπ
ω
ω
N (9.21)
注意边界上{T }和q也被离散化为
e
W
TNT }]{[}{ = (9.22)
9.23)
e
h
qNq }]{[=
如果将式(9.21)中的单元节点位移{W}
e
和水头{h}
e
改写成系统整体的位移{ν}和水头
{?}可表达为
[]{} { } [] {}
[] [] []
{ } {} { } {} { } }{****
}{*}{}{*}{
2
1
}{**}{
2
1
}{* }{*}{}{*}{
2
1
212
4
2
4
2
2
4311
γνν
γ
γ
γ
νγν?γνννπ
T
w
TT
T
o
TwT
w
T
w
TT
w
TT
PPM
k
Q
k
Q
k
MkMk
w
++?
++
+=
(9.24)
其中
(9.25) []
∑
∫
=
=
n
i
V
W
T
W
VBCBk
1
1
d]][[][
(9.26) []
∑
∫
=
=
n
i
v
hh
VNkNk
1
2
d]][][[
(9.27) []
∑
∫
=
=
n
i
V
h
VNBk
1
3
d]][[
(9.28) []
∑
∫
=
=
n
i
V
h
T
h
VNNk
1
4
d][][
(9.29) {}
∫
=
V
T
w
VBM
01
d}{][ σ
(9.30) {}
∑
∫
=
=
n
i
V
T
bw
VfNM
1
2
d}]{[
第9章 有限元法在边坡稳定分析中的应用 243
(9.31) {}
∑
∫
=
=
n
i
V
VByM
1
4
d][
{}
∑
∫
=
=
S
i
S
e
u
T
u
sTNNP
1
1
1
d}]{[][ (9.32)
(9.33) {}
∑
∫
=
=
S
i
S
e
h
T
h
sqNNP
1
2
4
d}]{[][
式(9.25)~式(9.33)系指各矩阵按其在整体结构矩阵中的位置叠加
3,应用里兹法求解泛函的极小值
对式(9.24)所表达的 π 进行式(9.12)式(9.13)的运算可以得到最终的线性方程组
(9.34) }{}{}{}{}]{[}]{[
141231
MMPMkk
w
++?=+γν
}]{[}{}]){[][(}{][
04
2
2423
γ
γ?
γ
γνγ k
Q
Pk
Q
kgk
w
w
w
w
T
w
+?=+?+ (9.35)
求解这个方程组即获得了用有限元法得到的固结问题的解实际计算时常采用增量法具体的数值分析步骤可参见文献陈祖煜等2003
9,2,3 结构矩阵的形成和计算
1,三角形单元
(1) 形状函数三角形单元的自变量用以下矩阵表达
(9.36) ),,,,,(}{
332211 yxyxyx
eT
WWWWWWW =
9.37) ),,(}{
321
hhhh
eT
=
单元内任一点(x,y)的{W}和h可用该单元的{W}
e
和{h}
e
表示即式(9.14)和(9.15)
其中
(9.38) ),,(][
321
NNNN
h
=
(9.39)
=
321
321
000
000
][
NNN
NNN
N
W
而
(9.40)?2/)(
1111
ycxbaN ++=
9.41)
23321
yxyxa?=
(9.42)
321
yyb?=
(9.43)
231
xxc?=
244 土质边坡稳定分析?原理
方法
程序
33
22
11
1
1
1
det2
yx
yx
yx
=? (9.44)
参考图9.1在脚标为2 3时N a b c的相应表达式可依次类推式(9.18)和式(9.19)
中
图 9,1 三角形单元
=
=
=
321
321
321
321
321
2
1
),,(][
aaa
bbb
y
N
y
N
y
N
x
N
x
N
x
N
NNN
y
x
B
h
(9.45)
=
=
=
332211
321
321
332211
321
321
321
321
000
000
2
1
000
000
000
000
0
0
][
bcbcbc
ccc
bbb
x
N
y
N
x
N
y
N
x
N
y
N
y
N
y
N
y
N
x
N
x
N
x
N
NNN
NNN
xy
y
x
B
W
(9.46)
(2) 单元矩阵的积分将式(9.38)式(9.39)式(9.45)式(9.46)代入式(9.25)~式(9.33)
可以得到相应各单元矩阵一般来说它们都具有这样的表达形式
∫
v
cba
dVNNNk
321
其中k
为不包含x,y的常数
可以证明对一个三角形单元
∫
+++
=
V
cba
cba
cba
dvNNN
321
2
)!2(
!!!
(9.47)
式中N1 N2 N3为形函数!表示阶乘运算?为三角形单元面积a b c为指数
这样就算得式(9.25)~ 式(9.33)各单元矩阵系数的数值下一步具体求解线性方程式(9.34)和式(9.35)即可得到问题的最后解
第9章 有限元法在边坡稳定分析中的应用 245
2,四节点四边形单元
(1) 形状函数对任一单元在该单元形心处建立局部座标(s,t)如图9.2所示将单元内任一点的x,y,{W},u均表达成s和t的函数
4321
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
xtsxstxtsxtsx +?++++?++= (9.48)
4321
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
ytsystytsytsy +?++++?++= (9.49)
采用式(9.48)式(9.49)这样的表达式当四边形节点(x,y)值分别为(x1,y1),(x2,y2),
(x3,y3) 和 (x4,y4)时使得相应节点的(s,t坐标分别为(-1,-1),(+1,-1),(+1,+1),
(-1,+1)
用矩阵来表达式(9.48)和式(9.49)可写成
(9.50)
=
4
3
2
1
4321
),,,(
x
x
x
x
NNNNx
(9.51)
=
4
3
2
1
4321
),,,(
y
y
y
y
NNNNy
图 9,2 四节点四边形单元
其中
+?=
++=
+=
=
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
4
3
2
1
tsN
tsN
tsN
tsN
(9.52)
和式(9.38)式(9.39)类似有
(9.53) ),,,(][
4321
NNNNN
h
=
246 土质边坡稳定分析?原理
方法
程序
(9.54)
=
4321
4321
0000
0000
][
NNNN
NNNN
N
w
而相应的{W}
e
和{h}
e
表达式应为
(9.55)
T
yxyxyxyx
e
WWWWWWWWW ),,,,,,,(}{
44332211
=
(9.56)
Te
hhhhh ),,,(}{
4321
=
计算应变和水头梯度时首先需要建立对整体坐标微分和对局部坐标微分之间的关系
=
y
x
t
s
J (9.57)
其中
[][](YX
t
s
t
y
t
x
s
y
s
x
,
][
][
=
=J ) (9.58)
=
=
=
=
),,,(
,
,
,
][
),,,(
,
,
,
][
4321
4321
4321
4321
tttt
ssss
NNNN
t
N
t
N
t
N
t
N
t
NNNN
s
N
s
N
s
N
s
N
s
(9.59)
(9.60)
=
=
),,,(][
),,,(][
4321
4321
yyyyY
xxxxX
T
T
而
=+?=
+=+=
+?=?=
==
)1(
4
1
,)1(
4
1
)1(
4
1
,)1(
4
1
)1(
4
1
,)1(
4
1
)1(
4
1
,)1(
4
1
44
33
22
11
sNtN
sNtN
sNtN
sNtN
ts
ts
ts
ts
(9.61)
第9章 有限元法在边坡稳定分析中的应用 247
式(9.57)可变换成
=
t
s
y
x
1
J (9.62)
和式(9.45)式(9.46)相似
=
=
=
tttt
ssss
h
NNNN
NNNN
NNNN
t
s
NNNN
y
x
B
4321
43211
4321
1
4321
),,,(
),,,(
][
J
J (9.63)
=
4321
4321
0000
0000
0
0
][
NNNN
NNNN
xy
y
x
B
W
(9.64)
假如将式(9.62)中的J
-1
写成
(9.65)
=
2221
12111
QQ
QQ
J
代入后有
t
Q
s
Q
x
1211
+
=
(9.66)
t
Q
s
Q
y
2221
+
=
9.67)
将式(9.66)式(9.67)代入式(9.64)得
248 土质边坡稳定分析?原理
方法
程序
+
+
+
+
=
1424132312221121
24232221
14131211
4321
4321
12112221
2221
1211
0000
0000
=
0000
0000
0
0
][
mmmmmmmm
mmmm
mmmm
NNNN
NNNN
t
Q
s
Q
t
Q
s
Q
t
Q
s
Q
t
Q
s
Q
B
W
(9.68)
其中
(i
itisi
NQNQm
12111
+=
=1,2,3,4) (9.69)
(i
itisi
NQNQm
22212
+=
=
1,2,3,4) (9.70)
(2) 单元矩阵的积分将式(9.55)式(9.56)式(9.57)式(9.58)分别代入式(9.25)~
式(9.33)可以得到相应各单元矩阵注意到dxdy = det[J]dtds可得到一般的表达式
∫∫
1
1
1
1
dd),( tstsF用高斯积分法来计算式(9.25)~式(9.33)各单元矩阵系数的值一般用四点法此时
(9.71)
∫∫
∑∑
==
≈
1
1
1
1
2
1
2
1
),(),(
ij
jiji
tsFdsdttsF αα
式中 57735.0
11
== ts ; ;57735.0
22
== ts 0.1
21
==αα
9,2,4 本构关系
进行固结和结构分析的一个重要工作是确定土的应力应变关系即矩阵[C]这里涉及到土的本构关系已有许多这方面的研究成果现结合土体的应力应变分析作一简要的回顾
1,弹性模型
建立在广义定律基础上的弹性理论对式(6.32)中的[C]= [C
e
]的表达式为
+
=
ν
νν
νν
νν
2
1
00
01
01
)21)(1(
][
E
C
e
(9.72)
式中E为弹性模量ν为泊桑比
2,非线性弹性模型
土的变形特征至少与周围应力σ3有直接的关系因此相应不同的σ3采用不同的E ν
值是比较符合实际的
第9章 有限元法在边坡稳定分析中的应用 249
图 9,3 通过固结仪的单向压缩曲线整理压缩系数
在筑填土形成的边坡中Penman(1973)
等人建议利用固结仪单向压缩所得轴向应变
εv和应力σ
v
关系曲线 整理相应坝体某处上覆荷重的平均压缩系数m
v
见图9.3 E和ν
分别由下式确定
]
1
)1)(21(
[
1
0
00
K
KK
m
E
v
+
+
= (9.73)
(9.74) )1/(
00
KK +=ν
式中K0为侧压力系数即大小主应力比
σ3/σ1需要对常规固结仪进行改造以测定
σ3从而确定K0
土石坝中大部分土单元在施工期σ1和σ3同时增长K0大致保持不变这和固结仪中单向压缩条件下的应力路径很相似查理斯(Charles,1976)等曾专门著文论述此问题认为此简单的模型仍能得出合理的成果Eisensfein 和 Simmons (1975) 曾使用此法分析加拿大高250m的Mica坝Cathie 和 Dungar (1978)曾系统报导过一些成功的实例在确定材料的力学特性时他们采用了一种可以测定σ3的单向压缩仪
3,双曲线弹性模型
邓肯和张提出双曲线?指数模型用常规三轴试验确定土的非线性参数用双曲线函数拟合轴向应力σa和轴向应变εa的关系用指数函数拟合体积模量K和周围应力σ3的关系据此可按下式确定E K
(9.75)
if
ESRE
2
)1(?=
9.76)
m
aab
PPKK )/(
3
σ=
其中S为应力水平按下式定义
(9.77)
f
S )/()(
3131
σσσσ=
(9.78) )sin1/()cos2sin2()(
331
φφφσσσ?+=? c
f
式中c,φ分别为土的凝聚力的内摩擦角P
a
为大气压E
i
为相应某一固结应力σ
3
的初始模量即
(9.79)
n
aai
PPKE )/(
3
σ′=
采用非线性弹性模型进行有限元计算时存在的另一个问题是如何处于确定土体处于卸荷状态以及如何使用卸荷模量邓肯建议卸荷模量按下式确定
(9.80)
n
aaurur
PPKE )/(
3
σ=
参数c,φ,R
f
,K′,n,K
ur
的整理方法见土工试验规程
4,弹塑性理论模型
土的弹塑性理论是把土的总变形分成弹性变形和塑性变形两部分用虎克定律计算弹性
250 土质边坡稳定分析?原理
方法
程序
变形部分用塑性理论来计算塑性变形部分对于塑性变形部分要作三方面的假定即破坏准则和屈服准则硬化规律和流动法则
土的有效应力弹塑性本构关系用有效应力增量d与应变增量可表达为 σ′ εd
(9.81) }d]{[}d{ εσ
ep
C=′
其中
A
g
C
f
C
fg
C
CC
eT
eTe
eep
+
′?
′?
′?
′?
=
}]{[}{
][}}{]{[
][][
σσ
σσ
(9.82)
式中[C
ep
]为弹塑性矩阵[C
e
]为弹性矩阵f是以H为硬化参数的屈服函数即)(HFf =
g为势函数A是反映硬化特性的一个变量与硬化参数的选择有关
当f = g时称为相关联流动法则当f ≠ g 时称为不相关联流动法则
(1) 加载过程中的弹塑性模型剑桥理论英国剑桥大学K,H,Roscoe (1958)提出了状态边界面临界状态线的概念其屈服函数为
0
2
2
)1(),( p
M
pqpff ′=+′=′=
η
(9.83)
式中pq ′=η相应的屈服面见图9.4
剑桥模型是对正常固结重塑粘土进行三轴试验得到的为了将此模型用于土石坝中填筑土等情况时邓肯(Duncan,1981)在下列几个方面做了修正
(a) 对土的强度包线的处理剑桥模型假定正常固结土的强度包线在 p′–q 平面上是通过原点的但填筑土具有粘聚力如图9.5所示邓肯建议采用强度包线不通过原点而在p′ 轴上相交于离原点距离为 p′r 的点上因此上面推导的所有公式中有关 p′ 的量都应加上 p′r 的值
(b) 对等向压缩曲线的处理邓肯认为等向压缩曲线的斜率 λ 不是常数可随 p′ 改变而改变因此输入的 λ 是与p′ 有关的一系列数值在计算时通过内插来决定处于某一应力状态的λ值
(c) 对土的泊松比的处理邓肯认为在弹性部分泊松比v也随 p′ 改变而改变因此也需输入一组v和p′ 的数值供计算时内插取值这样可更好地反映土的应力应变关系通过三轴和单向压缩试验确定和拟合参数作者在另著陈祖煜等2003中结合小浪底工程介绍了拟合过程
(2) 处于极限平衡状态时的弹塑性模型参见图9.4和图9.5土体在加载过程中可能处于A点或B点在A点土体处于正常固结加载状态其屈服面通常可用剑桥模型来描述在B点土体处于破坏面上塑性变形呈剪胀模型此时需要采用摩尔库仑或
Drucker-Prager准则来描述分述如下
在有限元计算中当某一点进入破坏面即图9.4和图9.5中B点时土体不再按加载第9章 有限元法在边坡稳定分析中的应用 251
的模型变形通常会出现明显的剪胀特性此时应使用OB线所代表的屈服面即摩尔库仑准则其屈服函数为
φφθθφ cossinsin
3
1
cossin
22
′+′+′= JJpf (9.84)
其中
)(
3
1
zyx
p σσσ ′+′+′=′ (9.85)
2
2Jq ′= (9.86)
}666)()(){(
6
1
222222
2 zxyzxyxzzyyx
J τττσσσσσσ +++′?′+′?′+′?′=′ (9.87)
2/3
2
3
)(2
33
J
J
′
′?
=θ (9.88)
9.89)
3
233
ppJIJ ′?′+=′
在σ1,σ2,σ3坐标中摩尔库仑准则的屈服面如图9.6示这是一个带角点的六边形
为了克服对有限元计算带来的困难通常采用Drnckel-Prager准则此时屈服面为
eJpf ′?′+′′=
2
sin3 φ (9.90)
相应的屈服面如图9.5所示在以下的讨论9.4.2节中还将介绍进一步的修正
图 9,4 剑桥模型 图 9,5 修正的剑桥模型
252 土质边坡稳定分析?原理
方法
程序
图 9,6 摩尔?库仑和Drucker?Prager准则
(3) 平面应变条件下弹塑性结构矩阵的构建从本节推导的平面应变条件下的应力和应变之间的关系式中可以看到屈服准则通常是用应力不变量p′
和q来代表的但是建立在平面应变条件下的刚度矩阵却是以,,τ
x
σ′
y
σ′ xy 为自变量构筑的
在使用式(9.82)计算[C
ep
]时需要知道}
σ′?
f
{此时使用了以下的关系式
T
xy
z
y
x
xy
z
y
x
xy
z
y
x
q
f
p
f
q
q
q
q
p
p
p
p
f
f
f
f
′?
′?
′?
′?
′?
′?
′?
′?
′?
′?
′?
′?
=
′?
′?
′?
′?
τ
σ
σ
σ
τ
σ
σ
σ
τ
σ
σ
σ
(9.91)
故对式(9.83)代表的剑桥模型可以得到
]
3
1
2
[
3
1
22
2
2
pM
q
pM
f
kji
i
′
′
′?′?′
+=
′?
σσσ
σ
(9.92)
式中i = x,y,z在i代表x时j,k分别代表y,z依此类推
xy
xy
Mp
f
τ
τ
2
6
′
=
(9.93)
作者在其它文献陈祖煜等,2003中介绍了按公式(9.91)构筑的子程序CAMCLY
第9章 有限元法在边坡稳定分析中的应用 253
9,2,5 实例分析
1,小浪底心墙固结分析
在第6章中介绍了该工程的概况和孔隙水压计算成果图9.7 是竣工蓄水时坝体填筑完毕上游水位升至250m高程大坝的累积位移矢量图从整体看由于心墙刚度相对较低坝体的变形呈现心墙部分大坝体两侧和基础部分小的态势其变形在坝体中间部位偏向下游侧而在坝体两侧基本呈竖直向
由图中可以看出竣工蓄水时坝体的最大沉降约为2.14m竣工后由于固结的作用使坝体沉降继续增大至两年半形成稳定渗流以后坝体的变形不再明显蓄水5年后坝体的累积最大沉降值达到2.62m
大坝竣工蓄水时的大主应力小主应力等值线如图9.8图9.9所示从总体上看高程较低部分的应力较大下游部位比上游部位的应力要大在同一高程上无论大主应力还是小主应力其最大值都出现在下游坝壳与心墙毗邻的部位且心墙部位的大主应力明显小于两侧坝壳的大主应力上下游坝壳处的剪应力较大而斜心墙部位的剪应力非常小斜心墙基本以竖直方向的压缩为主
在整个坝体中除了高塑性土区个别单元出现小主应力小于零外其它单元的小主应力全部大于零未出现拉应力
对比分析结果表明心墙部位的应力随着时间推移有一定的增长这是心墙固结过程中孔隙压力逐渐消散的结果在同一高程上大小主应力在心墙部位有较明显的增长心墙底部单元的应力在蓄水以后仍有相当的增长分析心墙内部的应力状态可以发现各个单元的有效小主应力全部大于零不会出现拉伸破坏发生水力劈裂的可能性较小
图 9,7 <250m水位>竣工蓄水时的总位移矢量图
图 9,8 <250m水位>竣工蓄水时大主应力等值线单位10kPa
254 土质边坡稳定分析?原理
方法
程序
图 9,9 <250m水位>竣工蓄水时小主应力等值线单位10kPa
2,务平大坝坝基固结分析
有关务平工程的概况和孔隙水压分析成果同样也已在第6章作过介绍坝体竣工时的总位移等值线图9.10很好地反映基础和坝体的位移与沉降从整体上看软基的位移量比较小竣工时最大值为0.36m固结完成时最大值为0.4m位移最大值出现在软基最深厚处的软基表面由于软基厚度的减小整个位移有向河谷岸边减小的趋势坝体的位移与沉降呈现中间大两侧小的态势竣工时坝体的最大位移为0.82m出现在心墙1/3高度处固结完成时坝体的最大位移为0.9m
大坝固结完成时的大小主应力的等值线如图9.11和图9.12所示总体上看高程较低部分应力较大由于没有蓄水心墙部分的大小应力呈对称状态由于弹性模量大小的差异同一高度处心墙部分的大主应力明显小于两侧坝壳的大主应力小主应力则相差不大坝体的最大大主应力约为500kPa而坝基的最大大主应力达900kPa心墙和坝壳的剪应力都不大说明坝体主要以竖直方向的压缩为主坝体的应力水平较低一般都在0.7以下
软基部分的应力呈现显明的波浪形说明在碎石桩上存在应力集中正是这种应力集中使得处理后的软基能够满足工程的需求从图9.11和图9.12中可以看出软土的大主应力一般比同一高度的碎石桩上的值要小约100kPa随着位置的不同软基部分的大主应力的范围为300~900kPa
图 9,10 竣工时的总位移等值线图单位m
第9章 有限元法在边坡稳定分析中的应用 255
图 9,11 竣工时的大主应力等值线图单位kPa
图 9,12 竣工时的小主应力等值线图单位kPa
9,3 建立在滑裂面应力分析基础上的边坡稳定有限元分析
9,3,1 概述
几乎在有限元法开发的同时研究者就开始了其与边坡稳定分析中传统的条分法关系的研究早期的工作Wright,1973; Yamagam,1988等是在边坡中定义一个潜在的滑裂面根据有限元得出的应力分布计算滑裂面上各点的应力水平然后根据加权平均的原则定义安全系数他们发现虽然有限元计算仍处于弹性分析的水平但按照这一方法计算的安全系数和极限平衡方法十分接近
在这一工作的基础上研究者进一步研究在诸多滑裂面中寻找一个临界滑裂面的可能性这工作使用了与极限平衡法相同的计算步骤同时也获得了与极限平衡法接近的最小安全系数和临界滑裂面近期研究沿用传统的缩小强度指标的作法将非线性有限元数值计算分析收敛条件崩溃作为边坡失稳的判据此时强度指标被缩小的倍数即为安全系数
Griffith(1999)赵尚毅和郑颖人等(2002)在这方面作出了有益的探索按照这一方法求解首先需要具备一个成熟的非线性有限元计算程序否则无法判断数值计算发散的原因是边坡系统的物理原因还是有限元程序本身的缺陷
本节简要介绍上述有关的研究成果
9,3,2 安全系数的定义
在极限平衡法中整个滑面的安全系数定义为
256 土质边坡稳定分析?原理
方法
程序
τ
τ
f
F = (9.94)
式中τf为土的抗剪强度τ为土的剪应力安全系数即为抗滑力与滑动力的比值
根据沿滑裂面剪应力的计算方法安全系数大致有如下三种定义法
1,基于应力水平的定义法
保持不变
3
σ′作一个与摩尔库仑强度包线相切的应力圆称为破坏应力圆相应的应力圆直径为
f
)(
31
σσ ′?′
f
)/()(
3131
σσσσ ′?′′?′为应力水平如式(9.77)定义它表示当前应力圆直径与破坏应力圆直径之比反映了强度的发挥程度
整个滑面的安全系数定义为
∫
∫
′?′
′?′
=
l
l
F
f
FE
d
)(
d
31
31
1
σσ
σσ
(9.95)
2,基于剪应力的定义法
滑裂面上某点的应力状态为σx σy τxy法向应力σ′n和切向应力τ可由下式确定
ατασστ 2cos2sin)(
2
1
xyxy
+?= (9.96)
9.97) ατασασσ 2sincossin
22
xyyxn
+=′
式中α为此处滑面切线与x轴夹角
该点的抗剪强度由下式确定
(9.98) φστ ′′+′= tan
nf
c
式中c′
为材料的内聚力φ′ 为材料的内摩擦角
整个滑面的安全系数定义为
∫
∫
′′+′
=
l
lc
F
n
FE
d
d)tan(
2
τ
φσ
(9.99)
3,基于应力水平加权强度的定义法
Donald和Tan(1985)提出以如下公式计算一个滑裂面的安全系数方法式中各符号的含义同上
∫
∫
′′+′
′?′
′?′
′′+′
=
lc
lc
F
n
f
n
FE
d)tan(
)(
d)tan(
31
31
3
φσ
σσ
σσ
φσ
(9.100)
第9章 有限元法在边坡稳定分析中的应用 257
9,3,3 有关成果介绍
1,Donald 等的模式搜索法
Giam和Donald于1988年提出了一种由已知的应力场确定临界滑裂面及最小安全系数的方法称为CRISS法CRISS法的基本原理示于图9.13和图9.14由有限元法计算所得应力水平较高点P出发形成一个由顶到底的破坏面其过程如下将边坡分成若干土条
(stage)见图9.13由P点出发分别向上向下在相邻两土条内构筑局部滑裂面图9.14
给出了由P点向上在第N+1个土条内构筑局部滑裂面的过程即以P点为基点呈放射状按1°倾角变化形成m个条块每一条块分为长度为?L的五小段每一小段以其中点为控制点这些控制点的应力
xx
σ′
xy
σ′
xy
τ′由有限元计算结果而定这些控制点的正应力和剪应力计算如下
ατασσσσσ 2sin2cos)(
2
1
)(
2
1
xyxxyyxxyyni
′?′+′+′=′ (9.101)
ατασστ 2cos2sin)(
2
1
xyxxyyni
+′?′= (9.102)
由P点出发第i个条块对应的安全系数为
L
Lc
F
ni
ni
i
∑
′+′∑
=
)(
)tan(
τ
φσ
(9.103)
找出F
i
中的最小值F
imin
(
i
=
1,…,m
) Fimin对应的条块PQ即为第N+1个土条内的局部滑裂面位置由Q点为基点可以找到第N
+2个土条内的局部滑裂面位置依此类推到边坡顶部由P点向下找寻第N,N?1,…,1个土条内的滑裂面过程类同每一条块内的F
imin
对应的局部滑裂面即构成这一边坡的滑裂面其安全系数为
L
Lc
F
n
n
FE
d
d)tan(
2
τ
φσ
∫
∫
′′+′
= (9.104)
258 土质边坡稳定分析?原理
方法
程序
图 9,13 CRISS法求解过程简图
图 9,14 CRISS法的一个求解步
[例9.1] 简单边坡的稳定分析
图9.15为用CRISS法的F
FE2
计算结果此例为第11章11.4.2节介绍的ACAD边坡稳定分析程序考题1(a)详细数据可参见该节曾经使用各种不同的极限平衡分析程序进行计算得到一致的解答最小安全系数为1.0表9.1汇总了对应于不同安全系数定义法的有限元法计算结果可以看出一次填筑和分层填筑边坡的计算值非常接近开挖边坡的计算值略高所有计算结果和极限平衡法的解答一致
表 9,1 安全系数比较
安全系数
类型
FFE1 FFE2 FFE3
分层填筑1.001 1.001 1.001
一次填筑 1.000 1.003 1.000
开挖1.082 1.044 1.044
第9章 有限元法在边坡稳定分析中的应用 259
图 9,15 [例9.1] 简单边坡的CRISS法分析结果
2,Zou等(1995)的动态规划法
Zou等通过有限单元法获得的应力分布规律确定滑裂面的范围和初始滑裂面然后利用动态规划的数值方法(IDPM)搜索最危险滑裂面及相应的安全系数
在这一方法中滑裂面安全系数定义为
Li
Li
L
L
F
i
Ns
i
fi
Ns
i
B
A
f
B
A
s
≈=
∑
∑
∫
∫
=
=
τ
τ
τ
τ
1
1
d
d
(9.105)
式中τ为滑裂面AB上任一点的剪应力τ
f
为在这些点的材料抗剪强度L为滑裂面的长
度N
s
为沿滑裂面的分段数?L
i
为每一段的长度对于第i段
(9.106)φστ tan
sfi
c +=
(9.107) ατασστ 2cos2sin)(5.0
xyxyi
+?=
(9.108) ατασασσ 2sincossin
22
xyyxs
+=
式中c为材料的内聚力φ是材料的内摩擦角σ
s
为滑裂面上第i段对应法向应力α为
滑裂面第i段与水平面夹角σ
x
σ
y
τ
xy
为第i段对应的有效应力
引入附加函数G
F
则
(9.109) LFG
isfi
Ns
i
F
=
∑
=
)(
1
ττ
利用动态规划法通过求G
F
的极小值从而确定F
s
的最小值在选定的搜索区域内计算每一可能的滑裂面直至找到F
s
[例9.2] 一均质粘土边坡如图9.16其材料参数见表9.2表中最后一列最小安全系数为简化毕肖普法的计算结果
260 土质边坡稳定分析?原理
方法
程序
图 9,16 [例9.2] 均质粘土边坡示意图
表 9,2 材料参数
方案 凝聚力
(kPa)
摩擦角(°) 容重
(kN/m
3
)
杨氏模量(MPa) 泊松比 最小安全系数
1 3 20 20 200 0.25 1.00
2 10 20 20 200 0.25 1.39
利用平面有限单元法计算边坡的应力分布采用8结点等参单元边界条件为边坡两侧垂直边界只约束水平向位移边坡底部水平边界为固定边界材料假定为各向同性遵守
Moho?Coulomb准则的理想弹塑性材料图9.17给出了计算所得应力水平分布图
图 9,17 [例9.2] 边坡应力水平分布
根据图9.17的应力水平分布图Zou等采用图9.18图9.19所示两种状态点格栅图图
9.18分条方向为竖直向图9.19分条方向大致垂直于滑裂面分条厚为1m状态点间距
0.2m通过IDMP法计算发现,分条方向对计算结果影响不大两个方案对应安全系数分别为1.01和1.39与表9.2中所示简化毕肖普法的计算结果几乎一样图9.20图9.21给出第9章 有限元法在边坡稳定分析中的应用 261
了对应于表9.2中两个方案的IDMP分析结果及极限平衡条分法GWEDGEM的分析结果两种方法的结果也吻合
图 9,18 假定的滑裂面及竖直向分条的搜索范围
图 9,19 [例9.2] 假定的滑裂面及垂直于滑裂面分条的搜索范围
图 9,20 [例9.2] 方案1临界滑裂面
IDPM动态规划法GWEDGEM通用的极限平衡条分法FS安全系数
262 土质边坡稳定分析?原理
方法
程序
图 9,21 [例9.2] 方案2临界滑裂面
9,4 建立在强度缩小有限元分析基础上的边坡稳定分析
9,4,1 Griffiths和Lane的工作
建立在强度缩小有限元分析基础上的边坡稳定分析的基本原理是将边坡强度参数粘聚力c′
和内摩擦角φ′
同时除以一个折减系数F得到一组新的c′
和φ′
值然后作为一组新的材料参数输入再进行试算当计算不收敛时对应的F被称为坝坡的最小安全系数此时坝坡达到极限状态发生剪切破坏同时可得到临界滑动面
1,有限单元模型的简单介绍
本小节所提及的程序是建立于Smith和Griffiths 1998年开发的程序基础上的与其它程序相比主要的不同在于能够模拟更普遍的几何形状和土性的变化包括水位和孔压的变化并增加了图形化功能程序适用于摩尔库仑破坏准则下理想弹塑性土体的二维平面应变分析采用的是简化积分的八节点四边形单元每个单元有四个高斯点在重力荷载作用下刚度矩阵生成和应力再分配的算法中都采用这种单元假定土体开始为弹性的模型在网格内所有高斯点生成正应力和剪应力然后将这些应力与摩尔库仑破坏准则相比如果在特定高斯点上的应力位于摩尔库仑破坏圆内则该点仍然是弹性的如果位于圆上或圆外则该点处于屈服状态利用粘弹性算法屈服应力在网格中被重分配当足够数目的高斯点发生屈服使机制发生变化时即认为发生了整体剪切破坏本小节中提到的分析不能模拟拉裂
2,土体模型
本次研究采用的土体模型包括六个参数如表9.3所示
表 9,3 六个参数的土体模型
φ′ c′ ψ′ E′ v′ γ
摩擦角 凝聚力 剪胀角 杨氏模量 泊松比 容重
第9章 有限元法在边坡稳定分析中的应用 263
剪胀角在屈服过程中影响土体的体积变化众所周知在屈服过程中真实的土体体积是变化的例如一种中等密度的材料在剪切过程中也许初始体积是减小的(ψ
<0)随后是剪胀阶段(ψ
>0)最终是在体积为常量条件下的屈服(ψ 0)
随之而来的问题是ψ应取何值如果取ψ φ则塑性流动法则是相关联的并可与经典塑性理论进行对比当流动法则是相关联时应力和速度是一致的因此用有限元预测的破坏机制与用极限平衡法得到的临界破坏面基本一致
尽管使用相关联的流动法则有这些优点但我们知道非粘性土模型的相关联流动法则预测的剪胀要比现实中观察的大很多导致预测的破坏荷载偏大特别是在一些有限制的问题中如承载力问题(Griffiths 1982)
在稳定分析中研究对象常常是没有约束的天然边坡因此剪胀角的选择就变得不太重要这里研究的主要目的是准确地预测边坡的安全系数Griffiths等采用的是ψ 0的折中值对应屈服时体积改变量为零属非关联流动研究者发现ψ的取值可以使模型给出可靠的安全系数并合理地预测可能破坏面的位置和形状
参数c′
和φ′
指土体的有效凝聚力和摩擦角虽然有许多模拟土的强度的破坏准则但摩尔?库仑准则仍然是岩土工程实践中应用最广泛的一个也是此文所采用的根据主应力和以压为负的规定该准则可写成如下的形式
'cos
2
sin
2
3131
φ
σσ
φ
σσ
cF ′?
′?′
′
′+′
= (9.110)
式中
1
σ′和分别是最大和最小主应力
3
σ′
破坏函数解释如下F
0<F 应力在屈服圆内弹性
0=F 应力在屈服圆上屈服
0>F 应力在屈服圆外屈服并且必须重分配
弹性参数E′
和v′
是土体的杨氏模量和泊松比如果已经假定了泊松比(0.2< v′
< 0.3)
则杨氏模量与在一维固结中得到的土的压缩率有关由第6章式(6.17)可知
)1(
)21)(1(
ν
νν
ν
′?
′?′+
=′
m
E (9.111)
式中
ν
m是体积压缩率系数
虽然目前弹性参数的值对计算破坏前的变形有很大影响但对边坡稳定分析中的安全系数的影响却很小该文作者认为没有E′
和v′
的值时可以给一个象征性的值如E′
10
5
kN/m
2
v′
0.3
在计算过程中土体的总单位容重γ与由重力产生的节点自重荷载是成比例的
总而言之在有限元边坡稳定分析中采用的最重要的参数与传统的方法是一样的即容重γ强度参数c′
和φ′以及问题的几何形状
3,重力荷载
264 土质边坡稳定分析?原理
方法
程序
每一单元上由土体自重产生的力由以下积分式计算
(9.112)
∫
=
e
V
Te
vNp
)(
dγ
式中N是单元的形函数上标e代表单元号这个积分结果是将每个单元的面积与土的容重的乘积分配到各个节点上
4,安全系数的确定
一个土质边坡的安全系数是这样定义的为了使边坡达到破坏给原始的剪切强度参数除以一个数这个数就是安全系数安全系数这种定义和极限平衡法中的定义是相同的则经过折减的剪切强度参数和变为
f
c′
f
φ′
(9.113) Fcc
f
/′=′
′
=′
F
f
φ
φ
tan
arctan (9.114)
5,破坏的定义
在指定的最大迭代次数内如果算法不能收敛就意味着没有发现同时既能满足摩尔?
库仑破坏准则又能满足整体平衡的应力分布如果算法不能满足这些准则则说明破坏已经发生了边坡破坏和数值上的不收敛同时发生并且伴随着网格中节点位移的显著增加本文中大部分的计算使用的是最高限度为1000次的迭代结果以F与Eδmax/γH
2
无量纲的位移的图形形式显示其中δmax是收敛点的最大节点位移H是边坡的高度依靠位移的网格图和向量图来显示安全系数和破坏机制的特性
6,边坡稳定实例及验证
Griffiths提供了几个有限元边坡稳定分析的实例并在可能的条件下与传统的边坡稳定分析进行对比验证
[例9.3] 无地基层的均质边坡(D
=
1)
图9.22显示的均质边坡的土性如下20φ′=
o
/0.cHγ′ = 05边坡与水平面成
26.57°(2:1)左边界是水平向约束的而底部边界是完全固定的D的定义可参见图9.25在迭代限度内不断改变安全系数进行计算直到计算不收敛结果如表9.4所示
第9章 有限元法在边坡稳定分析中的应用 265
图 9,22 [例9.3] 坡角为26.57°的均质边坡的网格
表9.4显示的是6次试算得到的成果安全系数从0.8增大到1.4每一次的强度参数都除以相应安全系数值由于在每次分析中重力荷载和整体刚度矩阵都是相同的所以只要生成一次就可以了因而效率较高
如表9.4示当安全系数为0.8和1.0时只需进行2次和10次迭代计算就收敛越接近真实的安全系数算法就越难达到收敛当F 1.4时无量纲位移有一个突然增加
2
max
/EHδγ′
此时算法在允许迭代次数内不能收敛图9.23是根据表9.4的数据得到的安全系数与无量纲位移关系当F 1.4时位移突然增加计算不收敛意味着边坡破坏从该图可以看出对同样的问题有限元结果与由Bishop & Morgenstern得到的安全系数很接近
图9.24(a)和9.24(b)给出节点位移向量和对应于不收敛情况F 1.4时的变形的网格
图 9,23 [例9.3]和[例9.4] 安全系数与无量纲位移关系
266 土质边坡稳定分析?原理
方法
程序
图 9,24 [例9.3] 对应于不收敛解F 1.4的网格变形图
(a) 节点位移向量(b) 变形网格
表 9,4 [例9.3] 结果
F
2
max
/ HE γδ′ 迭代次数
0.80 0.379 2
1.00 0.381 10
1.20 0.422 20
1.30 0.453 41
1.35 0.544 792
1.40 1.476 1000
[例9.4] 有地基层的均质边坡(D
=
1.5)
图9.25显示的是在例9.3的边坡底面增加了一个厚度为H/2地基层的例子而其它的土性及几何形状都不变该例坡角为26.57° c′/γH
= 0.05 φ′
20° D = 1.5
第9章 有限元法在边坡稳定分析中的应用 267
图 9,25 [例9.4] 有地基层的均质边坡
(a) 没变形的网格(b) 与不收敛解F = 1.4相对应的网格
图9.25(a)和9.25(b)分别为初始网格和破坏时的变形网格从图9.25(b)可以清楚地获得
坡脚破坏类型的机制虽然由于可压缩性土体积增大位移也增加了但临界安全系数本质上并没有改变仍为F= 1.4有限元结果表明增加的地基层并没有明显改变边坡安全系数
Bishop 和 Morgenstern (1960)给出安全系数F
=
1.752作为该例D
=
1.5 c′/γH
=
0.05 φ′
20°坡度=2:1的一个可能解假定破坏时的圆弧与地基的底部相切用一个比较通用的滑弧程序也可以得到F 1.7这个可能的结果但要获得正确的安全系数1.4必须使滑弧通过坡脚这个例子说明了有限元边坡稳定分析法与传统边坡稳定分析方法相比的一个主要优点即有限元方法不要求假定临界面位置或形状破坏自然地发生在土的剪切强度不足以抵抗剪应力的地带因而不需要第4章所论述的搜索过程
[例9.5] 带有软弱夹层的不排水粘土边坡
下面的例子是关于不排水粘土边坡的稳定性分析在本例中采用总应力分析使用
Tresca破坏准则(φ
u
=
0)
图9.26显示的是在不排水粘土地基层上的边坡这个边坡包括一个软弱夹层这个软弱夹层一部分与边坡平行在地基的部分呈水平最后地面出露的部分与坡踵夹45°的角虽然这个例子看起来有点人为的痕迹但它有点像一个有软弱衬垫的废渣填埋系统的情况在计算该边坡的安全系数时保持不变25.0/
1
=Hc
u
γ只改变软弱薄层不排水抗剪强度
2u
c
268 土质边坡稳定分析?原理
方法
程序
图 9,26 [例9.5] 在地基层中有软弱夹层的不排水粘土边坡
有限元结果见图9.27对均质边坡(c
u2
/c
u1
=1)计算的安全系数与Taylor(1937)的F
=
1.47
的解很接近并且给出了期望的圆弧基本破坏机制随着薄层强度的逐渐减小在c
u2
/c
u1
≈0.6
时结果有很明显的改变
图 9,27 [例9.5] 不同c
u2
/c
u1
比值得到的安全系数F
该图还显示了用Janbu法得到的极限平衡结果该法同时采用了沿软弱层的圆弧和三线楔形破坏机制当c
u2
/c
u1
≈0.6时在由软弱层支配的圆弧和非圆弧之间的过渡阶段很清楚的显示了不连续性当c
u2
/c
u1
>0.6时为圆弧破坏机制支配安全系数本质上不受软弱薄层强度的影响当c
u2
/c
u1
<0.6非圆弧的薄层机制支配了破坏模式
图9.28可以很清楚地解释破坏模式图中显示了三种不同的c
u2
/c
u1
比值下破坏时的变形网格图9.28(a)示均质情况(c
u2
/c
u1
1)显示了一个本质上圆弧状的破坏机制图9.28(c)
中薄层的强度只有其周围土强度的20 (c
u2
/c
u1
0.2)从图中可以看出非圆弧机制与软弱薄层的路径几乎重合图9.28(b)中薄层的强度为其周围土强度(c
u2
/c
u1
0.6)的60图形复杂第9章 有限元法在边坡稳定分析中的应用 269
且模糊但是很明显至少有两种互相冲突的机制第一种是在坡脚以上与软弱层重合的破坏机制第二种是沿着与边坡和坡脚出露处平行的软弱层破坏的机制
图 9,28 [例9.5] 对应三种不同的c
u2
/c
u1
比值破坏时网格变形图
a c
u2
/c
u1
1.0; b c
u2
/c
u1
0.6; c c
u2
/c
u1
0.2
如果没有前面的关于两种不同机制的知识传统的极限平衡搜索出的安全系数可能会过高图9.27已经证明了这一点例如当c
u2
/c
u1
0.2时圆弧机制得到的F为1.3而正确的
F约为0.6
9,4,2 赵尚毅和郑颖等人的工作
上述Griffith等的研究成果受到国内外学者的重视赵尚毅和郑颖人等(2002)近期也获得有意义的成果他们采用了广义米赛斯准则
κα =+=
21
JIF (9.115)
式中I1 J2分别为应力张量的第一不变量和应力偏张量的第二不变量α κ是与岩土材料内摩擦角φ和内聚力c有关的常数不同的α κ 在 π
平面上代表不同的圆图9.29
这是一个通用的表达式通过变换α κ的表达式就可以在有限元中实现不同的屈服准则各准则的参数换算关系见表9.5
270 土质边坡稳定分析?原理
方法
程序
图 9,29 各屈服准则在π平面上的曲线
传统的极限平衡法采用摩尔?库仑准则但由于摩尔?库仑准则的屈服面为不规则的六角形截面的角锥体表面存在尖顶和菱角给数值计算带来困难为了与传统方法进行比较这里采用徐干成郑颖人等(1990)提出的摩尔库仑等面积圆屈服准则(DP4)代替传统的摩尔库仑准则
表 9,5 各准则参数换算表
编号 准则种类 α κ
DP1 外角点外接D-P圆
()φ
φ
sin33
sin2
()φ
φ
sin33
cos6
c
DP2 内角点外接D-P圆
()φ
φ
sin33
sin2
+
()φ
φ
sin33
cos6
+
c
DP3 内切D-P圆
φ
φ
2
sin33
sin
+
φ
φ
2
sin33
cos3
+
c
DP4 等面积D-P圆
()φπ
φ
2
sin932
sin32
()φπ
φ
2
sin932
cos36
c
[例9.6] 均质土坡稳定分析算例
此例为一均质土坡坡高H
=
20m土容重
=
25kN/m
3
粘聚力c
=
42kPa内摩擦角
=
17°
表 9,6 用不同方法求得的安全系数
坡角(°) 30 35 40 45 50
有限元法外接圆屈服准则 1.78 1.62 1.48 1.36 1.29
有限元法摩尔-库仑等面积圆屈服准则
1.47 1.34 1.22 1.12 1.06
Bishop 1.39 1.26 1.15 1.06 0.99
Spencer 1.46 1.32 1.21 1.12 1.04
第9章 有限元法在边坡稳定分析中的应用 271
计算结果见表9.7可见DP4准则与简化Bishop法所得的稳定安全系数最为接近其误差的平均值为5.7%最大误差小于8%且离散度很小而DP1准则的平均误差为29.5
采用DP2和DP3准则所得计算结果的离散度也很大因此在数值分析中可以用DP4准则代替摩尔库仑准则
表 9,7 不同屈服准则所得最小安全系数
H
=
20m β=45
o
c
=
42kPa
φ (
o
) 0.1 10 25 35 45
DP1 0.525 1.044 1.769 2.254 3.051
DP2 0.525 0.930 1.332 1.530 1.887
DP3 0.454 0.848 1.279 1.499 1.870
DP4 0.477 0.879 1.396 1.689 2.182
简化Bishop法 0.494 0.846 1.316 1.623 2.073
(DP1-Bishop)/Bishop 0.063 0.234 0.344 0.355 0.472
(DP2-Bishop)/Bishop 0.063 0.099 0.012?0.080?0.090
(DP3-Bishop)/Bishop?0.081 0.002?0.028?0.099?0.098
(DP4-Bishop)/Bishop?0.034 0.059 0.061 0.041 0.053
研究者还发现边界范围的大小在有限元法中对计算结果的影响比在传统极限平衡法中表现的更为敏感当坡角到左端边界的距离为坡高的1.5倍坡顶到右端边界的距离为坡高的2.5倍且上下边界总高不低于2倍坡高时计算精度最为理想另外如果网格划分太粗将会造成很大的误差计算时必须考虑适当的网格密度
从表9.7中计算结果可以看出采用外接圆屈服准则计算的安全系数比传统的方法大许多而采用摩尔?库仑等面积圆屈服准则计算的结果与传统的极限平衡方法Spencer法计算的结果十分接近说明采用摩尔?库仑等面积圆屈服准则代替摩尔?库仑不等角六边形屈服准则是可行的
另外通过4组方案的算例的比较分析表明用摩尔?库仑等面积圆屈服准则求得的安全系数与Bishop法的误差为3%~7%与Spencer法的误差为1%~3%说明有限元强度折减法完全在边坡稳定分析中应用的可行性
9,4,3 NDM结点位移法
Tan和Donald于1985年提出了一种利用有限元求解得到的结点位移来确定安全系数的简单图解法称为结点位移法这一方法是将结构的组成材料的粘聚力c和内摩擦系数乘以一个Nφtan=f
(N<1.0)
的折减系数进而通过有限元计算跟踪某一结点的位移δ
T
的变化随N的减小δ
T
不断增大如图9.30所示尖角点给出了N的临界值即为安全系数
T
F
δ
当各参数设计值乘以N之后结构物就可达到极限平衡 1/N即为设计安全系数为识别可能的滑裂面必须清楚在潜在破坏区域内结点的位移状况这可以很容易地从位移矢量图中观察到
272 土质边坡稳定分析?原理
方法
程序
这种方法包含了如下假定
1) 材料强度参数的降低不影响它的其它性质
2) 尽管有限元分析是建立在小变形理论基础上的但在出现尖角点前后大变形δ
T
的计算是有效的
图 9,30 [例9.7] 节点位移随1/N变化过程图
[例 9.7] Talbingo 坝
图9.31所示高为162米的Talbingo心墙坝强度参数从100%的设计值 (1/N
=
1.0)降到
40%(1/N
=
2.50)图9.32示一典型结点的位移变化过程F
δT
的范围大约在2.05~2.14之间对应于1/N
=
2.06和2.11的位移矢量参见图9.33在1/N
=
2.06时无明显大变形δ
T
出现坝中无明显滑裂面在1/N
=
2.11时在坝的上游坡堆石料区过渡料区和心墙料区可以很明显看到大变形出现由图 9.33(b)可以看到一明显的滑裂面用条分法计算所得通过心墙的滑裂面的安全系数为2.06
图 9,31 Talbingo坝网格及观测点
第9章 有限元法在边坡稳定分析中的应用 273
图 9,32 [例9.7] 151号节点位移变化过程
根据位移矢量图9.33(b)假定一圆弧滑裂面表9.8列出了各种方法计算的稳定分析结果可以看出与F
T
F
δ
FE1 FFE3很接近但略小于FFLE和FFE2各安全系数之间的一致性可以证明用δT与F来确定土石坝的安全系数是可行的
T
δ
图 9,33 [例9.7] Talbingo坝位移矢量图
(a) 1/N=2.06 (b) 1/N=2.11
表 9,8 安全系数比较
临界圆弧信息 总应力安全系数
x y r F
FE1
F
FE2
F
FE3
F
FLE
T
F
δ
100.3 290.0 278.8 2.10 2.23 2.11 2.21 2.11
274 土质边坡稳定分析?原理
方法
程序
注 坐标原点为上游坡与坝底的交点
参考文献
1 Biot,M,A,Theory of elasticity and consolidation for a porous anisotropic solid,Journal
of Applied Physics 1955,26,182-185
2 Bishop,A,W,and Morgenstern,N,R,Stability coefficients for earth slopes,Geotechnique,
1960,10:129-150
3 Cathie,D,N.,and Dungar,R,Evaluation of finite element predictions for constructional
behavior of a rockfill dam,Proc,Instn,Civ,Engrs.,Part 2,1978,65:551-568
4 Charles,J,A,The use of one dimensional compression tests and elastic in prediction
deformation of embankments,Canadian Geotechnioal Journal.1976,13,189-200
5 Chen L.H,etc,Consolidation and stability behavior of a high rockfill dam built on soft
clay foundation Proceedings of 3rd International Conference on Soft Soil Engineering ICSSE3,
2001,DEC,Hong Kong
6 陈胜宏,高坝复杂岩石地基及岩石高边坡稳定分析,北京中国水利水电出版社,2001
7 陈祖煜,用有限元进行渗流结构和固结计算中的基本原理,水利水电科学研究院,1985
8 陈祖煜陈立宏李新强,第7章 土工数值分析(二)渗流固结和应力应变的有限元方法高等土力学2003清华大学出版社,
9 Donald,I,B,Tan,C,P,And Goh,T,C,A,Stability of geomechanical structures assessed
by finite element method,Proc,2nd Int,Conf,In Civil Engr,Hangzhou,Science Press,Beijing,
1985,845-856
10 Duncan,J,M,et,al,CON2D:A finite element computer program for analysin of consolidation,
Report No,VCB1GT/81-01,University of California,Berkeley,1981
11 Duncan,J,M.,and Chang,C,Y,Nonlinear analysis of stress and strain in soils,J,Soil
Mech,and Found,Div,ASCE,1970,96(5),1629-1653
12 Donald,I,B,Finite element assessment of slope stability,Austrlian-China Landslide Seminar,
1993,3-18 July,1-42.D1-D18
13 Eisensein,Z,and Simmons,J,V,Three –dimensional analysis of Mica Dam,Proc,Int,Symp,
Criteria and Assumptions for Numer,Anal,of Dams,Quadrant,Swansea,U,K.1975,1051-1070
14 Griffith,D,V,and Lane,P,A,Slope stability analysis by finite elements,Geotechnique,
1999,49(3),387-403
15 Giam,S,K,and Donald,I,B,Determination of critical slip surfaces for slopes via
stress-strain calculations,5th A,N,Z,Conf,Geomechanics,Sydney,1988,461-464
16 Penman,A,D,M.,Burland,J,B.,and Charles,J,A,Observed and Predicted deformations in
a large embankment cam during con-struction,Proc,Ins,Civ,Engrs.,London,1971,49:1-21
17 Penman,A.,and Charles,A,Construction deformations in rock-fill dam,J,Soil Mech,And
Found,Div.,ASCE,1973.99(2),139-163
18 沈珠江,计算土力学,上海科技出版社,1990
19 Tan,C,P,and Donald,I,B,Finite element calculation of dam stability,Proc,11th Int,
Conf,Soil Mech,and Fnd,Engr,San Francisco,1980
20 Taylor,D,W,Stability of earth slopes,J,Boston Soc,Civ,Eng,1937,24:197-246
21 Wright,S,G,Kulhawy,F,H,and Duncan,J,M,Accuracy of equilibrium slope stability analysis,
J,Soil Mech,And FDN,Div,1973,99(SM10),783-791
第9章 有限元法在边坡稳定分析中的应用 275
22 Wang,C,and Duncan,J,M,Analysin of consolidation of earth and rockfill dans,Report No,
TE77-3,University Of California,Berkeley,1977
23 徐干城郑颖人,岩土工程中屈服准则应用的研究,岩土工程学报,1990,12(2):93-99
24 Yamagami,T,and Ueta,Y,Search for critical slip lines in finite element stress fields
by dynamic programming,Proc,6th International Conference on Numerical Methods in
Geomechanics,Innsbruck,Australia,1988,1335-1339
25 Zou,J,Z.,Williams,D,J.,and Xiong,W,L,Search for critical slip surfaces based on
finite element method,Canadian Geotechnical Journal,1995,32,233-246
26 赵尚毅郑颖人时卫民王敬林,有限元强度折减法求边坡稳定安全系数,岩土工程学报,2002,
129(3),343-346
9,1 概述
随着计算机软件硬件的飞速发展采用理论体系更为严格的方法进行边坡稳定分析已经成为可能有限单元法全面满足了静力许可应变相容和应力应变之间的本构关系同时因为是采用数值分析方法可以不受边坡几何形状的不规则和材料的不均匀性的限制因此应该是比较理想的分析边坡应力变形和稳定性态的手段
与传统的极限平衡法相比边坡稳定分析的有限元法的优点可总结如下
(1) 破坏面的形状或位置不需要事先假定破坏自然地发生在土的抗剪强度不能抵抗剪应力的地带
(2) 由于有限元法引入变形协调的本构关系因此也不必引入假定条件保持了严密的理论体系
(3) 有限元解提供了应力变形的全部信息
有关有限元的基本原理及其在边坡稳定分析中的应用已有很多的文献沈珠江1990
陈胜宏2001介绍作者也另有专著陈祖煜等2003介绍了采用有限元法求解第6章中所列渗流固结和应力应变分析控制方程的详细方法本章9.2节简要介绍这一计算方法的解题思路以及三角形四边形单元的离散方法并介绍有限元法在小浪底和务平工程中的应用
如何将有限元计算成果与传统的安全系数挂钩成为直接用于边坡设计的判别依据也是广泛受到重视的课题近期在这方面获得了有意义的成果本章将简要介绍这方面的成果
9,2 求解渗流和应力应变控制方程的有限元方法
9,2,1 基本原理
第6章所述求解渗流和应力应变的数学提法是在体力{fb}和边界上面力{T }的作用下求满足微分方程式(6.27)式(6.28)式(6.31)式(6.40)和边界条件式(6.30)式(6.43)
式(6.44)式(6.46)和式(6.48)的位移场{W}和水头场h用有限元求解这些偏微分方程组边值问题其基本途径是找出一个泛函π({W} h)当{W}和h为该偏微分方程边值问题解时
π获得极值
可以证明所寻找的泛函由下式确定
240 土质边坡稳定分析?原理
方法
程序
(9.1) {}()
7654321
,ππππππππ ++++++=hW
其中
∫∫
′′=
VV
T
VV
0
1
d}{}{d}{}{
2
1
εσεσπ (9.2)
(9.3) {}
∫
=
V
T
w
Vayh
2
d}{)( εγπ
{} {}{}
∫
=
V
T
w
Vhkh
3
d**
2
1
γπ (9.4)
∫
=
V
w
Vhhh
Q
0
2
4
d)2(
2
1 γ
π (9.5)
(9.6) {} {}
∫
=
V
T
b
VWf
5
dπ
{} {}
∫
=
1
d
6
s
T
sWTπ 9.7)
(9.8) shq
s
w
d
4
7
∫
= γπ
这一命题的理论证明可参见文献陈祖煜等2003
9,2,2 有限元解法
采用数值分析方法求解使获得极值的水头和位移场包括以下步骤
1,离散化
将所研究的域V分成n个单元共m个节点这m个节点的{W}用{v}来代表即
(9.9) ),,,,,,(}{
2211
T
ymxmyxyx
WWWWWW LL=ν
同样用{?}代表这m个节点的{h}值
(9.10) ),,,(}{
21
T
m
hhh LL=?
任一点的{W}和h可用该点所属单元的节点{W}
e
和{h}
e
近似表示本质上也就是可用
{ν}和{?}代表于是π可以近似地用{ν}和{?}来代表即
(9.11) }){},({)},({?νππ =hW
根据里兹法的原理使π取得极值的{ν}和{?}满足
Ο
Ο
对于某一行向量{a}
T
我们定义C/{a}
T
为一列向量{b}它满足
第9章 有限元法在边坡稳定分析中的应用 241
{}
{}0
2
1
1
=
=
ym
x
y
x
T
W
W
W
W
π
π
π
π
ν
π
M
M
(9.12)
{}
}0{
2
1
=
=
m
T
h
h
h
h
π
π
π
π
M
M
(9.13)
其中{0}为元素均为零的向量由式(9.12)和式(9.13)可得3m个线性方程可用来求解由
{ν}和{h}所包含的3m个未知数这就是有限单元法的基本原理
2,应用形状函数表达单元内的物理量
单元内任一点的{W}和h可用该单元节点的{W}
e
和{h}
e
来近似表达
(9.14)
e
W
WNW }]{[}{ =
9.15)
T
h
eTe
h
NhhNh ][}{}]{[ ==
因此
(9.16)
e
W
e
W
TT
WBWNW }]{[}]{[][}{][}{ =?=?=ε
9.17)
e
h
e
h
hBhNh }]{[}]{}[{}{ =?=?
式中
(9.18) ][][][
W
T
W
NB?=
9.19) [] ]}[{
hh
NB?=
[NW] [Nh]称形状函数或插值函数对三角形和四边形单元具有不同表达形式将在第9.2.3节中详述同时我们记
C={a}
T
{b}或 C={b}
T
{a}
242 土质边坡稳定分析?原理
方法
程序
(9.20) ][}{][
W
TT
BaB =
将式(9.14)式(9.15)式(9.16)式(9.17)和式(9.20)代入式(9.1)得
{}[]{}
∑
∫
∑
∫
∑
∫
==
=
++
+
+
+
′=
n
i
e
h
T
h
eTe
W
T
W
S
i
s
eT
e
h
T
h
eTweT
h
T
h
eTw
e
W
T
b
e
h
T
h
e
eT
w
eTT
h
eT
e
W
T
n
i
V
e
W
T
W
eT
shNNqsWNNT
VhNNh
Q
hNNh
Q
WfhNKNh
WByWBNh
WBWBCBW
w
1
s
1
0
22
T
0
1
d}{*][][}{*d}{*][][}{
d)}{*][][}{ }{*][][}{
2
1
}{*][}{-}{*]][[][*}{
2
1
}{*][}{*][][}{
*}]{][[][*}{
2
1
(
41
γ
γγ
γ
γγ
σπ
ω
ω
N (9.21)
注意边界上{T }和q也被离散化为
e
W
TNT }]{[}{ = (9.22)
9.23)
e
h
qNq }]{[=
如果将式(9.21)中的单元节点位移{W}
e
和水头{h}
e
改写成系统整体的位移{ν}和水头
{?}可表达为
[]{} { } [] {}
[] [] []
{ } {} { } {} { } }{****
}{*}{}{*}{
2
1
}{**}{
2
1
}{* }{*}{}{*}{
2
1
212
4
2
4
2
2
4311
γνν
γ
γ
γ
νγν?γνννπ
T
w
TT
T
o
TwT
w
T
w
TT
w
TT
PPM
k
Q
k
Q
k
MkMk
w
++?
++
+=
(9.24)
其中
(9.25) []
∑
∫
=
=
n
i
V
W
T
W
VBCBk
1
1
d]][[][
(9.26) []
∑
∫
=
=
n
i
v
hh
VNkNk
1
2
d]][][[
(9.27) []
∑
∫
=
=
n
i
V
h
VNBk
1
3
d]][[
(9.28) []
∑
∫
=
=
n
i
V
h
T
h
VNNk
1
4
d][][
(9.29) {}
∫
=
V
T
w
VBM
01
d}{][ σ
(9.30) {}
∑
∫
=
=
n
i
V
T
bw
VfNM
1
2
d}]{[
第9章 有限元法在边坡稳定分析中的应用 243
(9.31) {}
∑
∫
=
=
n
i
V
VByM
1
4
d][
{}
∑
∫
=
=
S
i
S
e
u
T
u
sTNNP
1
1
1
d}]{[][ (9.32)
(9.33) {}
∑
∫
=
=
S
i
S
e
h
T
h
sqNNP
1
2
4
d}]{[][
式(9.25)~式(9.33)系指各矩阵按其在整体结构矩阵中的位置叠加
3,应用里兹法求解泛函的极小值
对式(9.24)所表达的 π 进行式(9.12)式(9.13)的运算可以得到最终的线性方程组
(9.34) }{}{}{}{}]{[}]{[
141231
MMPMkk
w
++?=+γν
}]{[}{}]){[][(}{][
04
2
2423
γ
γ?
γ
γνγ k
Q
Pk
Q
kgk
w
w
w
w
T
w
+?=+?+ (9.35)
求解这个方程组即获得了用有限元法得到的固结问题的解实际计算时常采用增量法具体的数值分析步骤可参见文献陈祖煜等2003
9,2,3 结构矩阵的形成和计算
1,三角形单元
(1) 形状函数三角形单元的自变量用以下矩阵表达
(9.36) ),,,,,(}{
332211 yxyxyx
eT
WWWWWWW =
9.37) ),,(}{
321
hhhh
eT
=
单元内任一点(x,y)的{W}和h可用该单元的{W}
e
和{h}
e
表示即式(9.14)和(9.15)
其中
(9.38) ),,(][
321
NNNN
h
=
(9.39)
=
321
321
000
000
][
NNN
NNN
N
W
而
(9.40)?2/)(
1111
ycxbaN ++=
9.41)
23321
yxyxa?=
(9.42)
321
yyb?=
(9.43)
231
xxc?=
244 土质边坡稳定分析?原理
方法
程序
33
22
11
1
1
1
det2
yx
yx
yx
=? (9.44)
参考图9.1在脚标为2 3时N a b c的相应表达式可依次类推式(9.18)和式(9.19)
中
图 9,1 三角形单元
=
=
=
321
321
321
321
321
2
1
),,(][
aaa
bbb
y
N
y
N
y
N
x
N
x
N
x
N
NNN
y
x
B
h
(9.45)
=
=
=
332211
321
321
332211
321
321
321
321
000
000
2
1
000
000
000
000
0
0
][
bcbcbc
ccc
bbb
x
N
y
N
x
N
y
N
x
N
y
N
y
N
y
N
y
N
x
N
x
N
x
N
NNN
NNN
xy
y
x
B
W
(9.46)
(2) 单元矩阵的积分将式(9.38)式(9.39)式(9.45)式(9.46)代入式(9.25)~式(9.33)
可以得到相应各单元矩阵一般来说它们都具有这样的表达形式
∫
v
cba
dVNNNk
321
其中k
为不包含x,y的常数
可以证明对一个三角形单元
∫
+++
=
V
cba
cba
cba
dvNNN
321
2
)!2(
!!!
(9.47)
式中N1 N2 N3为形函数!表示阶乘运算?为三角形单元面积a b c为指数
这样就算得式(9.25)~ 式(9.33)各单元矩阵系数的数值下一步具体求解线性方程式(9.34)和式(9.35)即可得到问题的最后解
第9章 有限元法在边坡稳定分析中的应用 245
2,四节点四边形单元
(1) 形状函数对任一单元在该单元形心处建立局部座标(s,t)如图9.2所示将单元内任一点的x,y,{W},u均表达成s和t的函数
4321
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
xtsxstxtsxtsx +?++++?++= (9.48)
4321
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
ytsystytsytsy +?++++?++= (9.49)
采用式(9.48)式(9.49)这样的表达式当四边形节点(x,y)值分别为(x1,y1),(x2,y2),
(x3,y3) 和 (x4,y4)时使得相应节点的(s,t坐标分别为(-1,-1),(+1,-1),(+1,+1),
(-1,+1)
用矩阵来表达式(9.48)和式(9.49)可写成
(9.50)
=
4
3
2
1
4321
),,,(
x
x
x
x
NNNNx
(9.51)
=
4
3
2
1
4321
),,,(
y
y
y
y
NNNNy
图 9,2 四节点四边形单元
其中
+?=
++=
+=
=
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
)1)(1(
4
1
4
3
2
1
tsN
tsN
tsN
tsN
(9.52)
和式(9.38)式(9.39)类似有
(9.53) ),,,(][
4321
NNNNN
h
=
246 土质边坡稳定分析?原理
方法
程序
(9.54)
=
4321
4321
0000
0000
][
NNNN
NNNN
N
w
而相应的{W}
e
和{h}
e
表达式应为
(9.55)
T
yxyxyxyx
e
WWWWWWWWW ),,,,,,,(}{
44332211
=
(9.56)
Te
hhhhh ),,,(}{
4321
=
计算应变和水头梯度时首先需要建立对整体坐标微分和对局部坐标微分之间的关系
=
y
x
t
s
J (9.57)
其中
[][](YX
t
s
t
y
t
x
s
y
s
x
,
][
][
=
=J ) (9.58)
=
=
=
=
),,,(
,
,
,
][
),,,(
,
,
,
][
4321
4321
4321
4321
tttt
ssss
NNNN
t
N
t
N
t
N
t
N
t
NNNN
s
N
s
N
s
N
s
N
s
(9.59)
(9.60)
=
=
),,,(][
),,,(][
4321
4321
yyyyY
xxxxX
T
T
而
=+?=
+=+=
+?=?=
==
)1(
4
1
,)1(
4
1
)1(
4
1
,)1(
4
1
)1(
4
1
,)1(
4
1
)1(
4
1
,)1(
4
1
44
33
22
11
sNtN
sNtN
sNtN
sNtN
ts
ts
ts
ts
(9.61)
第9章 有限元法在边坡稳定分析中的应用 247
式(9.57)可变换成
=
t
s
y
x
1
J (9.62)
和式(9.45)式(9.46)相似
=
=
=
tttt
ssss
h
NNNN
NNNN
NNNN
t
s
NNNN
y
x
B
4321
43211
4321
1
4321
),,,(
),,,(
][
J
J (9.63)
=
4321
4321
0000
0000
0
0
][
NNNN
NNNN
xy
y
x
B
W
(9.64)
假如将式(9.62)中的J
-1
写成
(9.65)
=
2221
12111
J
代入后有
t
Q
s
Q
x
1211
+
=
(9.66)
t
Q
s
Q
y
2221
+
=
9.67)
将式(9.66)式(9.67)代入式(9.64)得
248 土质边坡稳定分析?原理
方法
程序
+
+
+
+
=
1424132312221121
24232221
14131211
4321
4321
12112221
2221
1211
0000
0000
=
0000
0000
0
0
][
mmmmmmmm
mmmm
mmmm
NNNN
NNNN
t
Q
s
Q
t
Q
s
Q
t
Q
s
Q
t
Q
s
Q
B
W
(9.68)
其中
(i
itisi
NQNQm
12111
+=
=1,2,3,4) (9.69)
(i
itisi
NQNQm
22212
+=
=
1,2,3,4) (9.70)
(2) 单元矩阵的积分将式(9.55)式(9.56)式(9.57)式(9.58)分别代入式(9.25)~
式(9.33)可以得到相应各单元矩阵注意到dxdy = det[J]dtds可得到一般的表达式
∫∫
1
1
1
1
dd),( tstsF用高斯积分法来计算式(9.25)~式(9.33)各单元矩阵系数的值一般用四点法此时
(9.71)
∫∫
∑∑
==
≈
1
1
1
1
2
1
2
1
),(),(
ij
jiji
tsFdsdttsF αα
式中 57735.0
11
== ts ; ;57735.0
22
== ts 0.1
21
==αα
9,2,4 本构关系
进行固结和结构分析的一个重要工作是确定土的应力应变关系即矩阵[C]这里涉及到土的本构关系已有许多这方面的研究成果现结合土体的应力应变分析作一简要的回顾
1,弹性模型
建立在广义定律基础上的弹性理论对式(6.32)中的[C]= [C
e
]的表达式为
+
=
ν
νν
νν
νν
2
1
00
01
01
)21)(1(
][
E
C
e
(9.72)
式中E为弹性模量ν为泊桑比
2,非线性弹性模型
土的变形特征至少与周围应力σ3有直接的关系因此相应不同的σ3采用不同的E ν
值是比较符合实际的
第9章 有限元法在边坡稳定分析中的应用 249
图 9,3 通过固结仪的单向压缩曲线整理压缩系数
在筑填土形成的边坡中Penman(1973)
等人建议利用固结仪单向压缩所得轴向应变
εv和应力σ
v
关系曲线 整理相应坝体某处上覆荷重的平均压缩系数m
v
见图9.3 E和ν
分别由下式确定
]
1
)1)(21(
[
1
0
00
K
KK
m
E
v
+
+
= (9.73)
(9.74) )1/(
00
KK +=ν
式中K0为侧压力系数即大小主应力比
σ3/σ1需要对常规固结仪进行改造以测定
σ3从而确定K0
土石坝中大部分土单元在施工期σ1和σ3同时增长K0大致保持不变这和固结仪中单向压缩条件下的应力路径很相似查理斯(Charles,1976)等曾专门著文论述此问题认为此简单的模型仍能得出合理的成果Eisensfein 和 Simmons (1975) 曾使用此法分析加拿大高250m的Mica坝Cathie 和 Dungar (1978)曾系统报导过一些成功的实例在确定材料的力学特性时他们采用了一种可以测定σ3的单向压缩仪
3,双曲线弹性模型
邓肯和张提出双曲线?指数模型用常规三轴试验确定土的非线性参数用双曲线函数拟合轴向应力σa和轴向应变εa的关系用指数函数拟合体积模量K和周围应力σ3的关系据此可按下式确定E K
(9.75)
if
ESRE
2
)1(?=
9.76)
m
aab
PPKK )/(
3
σ=
其中S为应力水平按下式定义
(9.77)
f
S )/()(
3131
σσσσ=
(9.78) )sin1/()cos2sin2()(
331
φφφσσσ?+=? c
f
式中c,φ分别为土的凝聚力的内摩擦角P
a
为大气压E
i
为相应某一固结应力σ
3
的初始模量即
(9.79)
n
aai
PPKE )/(
3
σ′=
采用非线性弹性模型进行有限元计算时存在的另一个问题是如何处于确定土体处于卸荷状态以及如何使用卸荷模量邓肯建议卸荷模量按下式确定
(9.80)
n
aaurur
PPKE )/(
3
σ=
参数c,φ,R
f
,K′,n,K
ur
的整理方法见土工试验规程
4,弹塑性理论模型
土的弹塑性理论是把土的总变形分成弹性变形和塑性变形两部分用虎克定律计算弹性
250 土质边坡稳定分析?原理
方法
程序
变形部分用塑性理论来计算塑性变形部分对于塑性变形部分要作三方面的假定即破坏准则和屈服准则硬化规律和流动法则
土的有效应力弹塑性本构关系用有效应力增量d与应变增量可表达为 σ′ εd
(9.81) }d]{[}d{ εσ
ep
C=′
其中
A
g
C
f
C
fg
C
CC
eT
eTe
eep
+
′?
′?
′?
′?
=
}]{[}{
][}}{]{[
][][
σσ
σσ
(9.82)
式中[C
ep
]为弹塑性矩阵[C
e
]为弹性矩阵f是以H为硬化参数的屈服函数即)(HFf =
g为势函数A是反映硬化特性的一个变量与硬化参数的选择有关
当f = g时称为相关联流动法则当f ≠ g 时称为不相关联流动法则
(1) 加载过程中的弹塑性模型剑桥理论英国剑桥大学K,H,Roscoe (1958)提出了状态边界面临界状态线的概念其屈服函数为
0
2
2
)1(),( p
M
pqpff ′=+′=′=
η
(9.83)
式中pq ′=η相应的屈服面见图9.4
剑桥模型是对正常固结重塑粘土进行三轴试验得到的为了将此模型用于土石坝中填筑土等情况时邓肯(Duncan,1981)在下列几个方面做了修正
(a) 对土的强度包线的处理剑桥模型假定正常固结土的强度包线在 p′–q 平面上是通过原点的但填筑土具有粘聚力如图9.5所示邓肯建议采用强度包线不通过原点而在p′ 轴上相交于离原点距离为 p′r 的点上因此上面推导的所有公式中有关 p′ 的量都应加上 p′r 的值
(b) 对等向压缩曲线的处理邓肯认为等向压缩曲线的斜率 λ 不是常数可随 p′ 改变而改变因此输入的 λ 是与p′ 有关的一系列数值在计算时通过内插来决定处于某一应力状态的λ值
(c) 对土的泊松比的处理邓肯认为在弹性部分泊松比v也随 p′ 改变而改变因此也需输入一组v和p′ 的数值供计算时内插取值这样可更好地反映土的应力应变关系通过三轴和单向压缩试验确定和拟合参数作者在另著陈祖煜等2003中结合小浪底工程介绍了拟合过程
(2) 处于极限平衡状态时的弹塑性模型参见图9.4和图9.5土体在加载过程中可能处于A点或B点在A点土体处于正常固结加载状态其屈服面通常可用剑桥模型来描述在B点土体处于破坏面上塑性变形呈剪胀模型此时需要采用摩尔库仑或
Drucker-Prager准则来描述分述如下
在有限元计算中当某一点进入破坏面即图9.4和图9.5中B点时土体不再按加载第9章 有限元法在边坡稳定分析中的应用 251
的模型变形通常会出现明显的剪胀特性此时应使用OB线所代表的屈服面即摩尔库仑准则其屈服函数为
φφθθφ cossinsin
3
1
cossin
22
′+′+′= JJpf (9.84)
其中
)(
3
1
zyx
p σσσ ′+′+′=′ (9.85)
2
2Jq ′= (9.86)
}666)()(){(
6
1
222222
2 zxyzxyxzzyyx
J τττσσσσσσ +++′?′+′?′+′?′=′ (9.87)
2/3
2
3
)(2
33
J
J
′
′?
=θ (9.88)
9.89)
3
233
ppJIJ ′?′+=′
在σ1,σ2,σ3坐标中摩尔库仑准则的屈服面如图9.6示这是一个带角点的六边形
为了克服对有限元计算带来的困难通常采用Drnckel-Prager准则此时屈服面为
eJpf ′?′+′′=
2
sin3 φ (9.90)
相应的屈服面如图9.5所示在以下的讨论9.4.2节中还将介绍进一步的修正
图 9,4 剑桥模型 图 9,5 修正的剑桥模型
252 土质边坡稳定分析?原理
方法
程序
图 9,6 摩尔?库仑和Drucker?Prager准则
(3) 平面应变条件下弹塑性结构矩阵的构建从本节推导的平面应变条件下的应力和应变之间的关系式中可以看到屈服准则通常是用应力不变量p′
和q来代表的但是建立在平面应变条件下的刚度矩阵却是以,,τ
x
σ′
y
σ′ xy 为自变量构筑的
在使用式(9.82)计算[C
ep
]时需要知道}
σ′?
f
{此时使用了以下的关系式
T
xy
z
y
x
xy
z
y
x
xy
z
y
x
q
f
p
f
q
q
q
q
p
p
p
p
f
f
f
f
′?
′?
′?
′?
′?
′?
′?
′?
′?
′?
′?
′?
=
′?
′?
′?
′?
τ
σ
σ
σ
τ
σ
σ
σ
τ
σ
σ
σ
(9.91)
故对式(9.83)代表的剑桥模型可以得到
]
3
1
2
[
3
1
22
2
2
pM
q
pM
f
kji
i
′
′
′?′?′
+=
′?
σσσ
σ
(9.92)
式中i = x,y,z在i代表x时j,k分别代表y,z依此类推
xy
xy
Mp
f
τ
τ
2
6
′
=
(9.93)
作者在其它文献陈祖煜等,2003中介绍了按公式(9.91)构筑的子程序CAMCLY
第9章 有限元法在边坡稳定分析中的应用 253
9,2,5 实例分析
1,小浪底心墙固结分析
在第6章中介绍了该工程的概况和孔隙水压计算成果图9.7 是竣工蓄水时坝体填筑完毕上游水位升至250m高程大坝的累积位移矢量图从整体看由于心墙刚度相对较低坝体的变形呈现心墙部分大坝体两侧和基础部分小的态势其变形在坝体中间部位偏向下游侧而在坝体两侧基本呈竖直向
由图中可以看出竣工蓄水时坝体的最大沉降约为2.14m竣工后由于固结的作用使坝体沉降继续增大至两年半形成稳定渗流以后坝体的变形不再明显蓄水5年后坝体的累积最大沉降值达到2.62m
大坝竣工蓄水时的大主应力小主应力等值线如图9.8图9.9所示从总体上看高程较低部分的应力较大下游部位比上游部位的应力要大在同一高程上无论大主应力还是小主应力其最大值都出现在下游坝壳与心墙毗邻的部位且心墙部位的大主应力明显小于两侧坝壳的大主应力上下游坝壳处的剪应力较大而斜心墙部位的剪应力非常小斜心墙基本以竖直方向的压缩为主
在整个坝体中除了高塑性土区个别单元出现小主应力小于零外其它单元的小主应力全部大于零未出现拉应力
对比分析结果表明心墙部位的应力随着时间推移有一定的增长这是心墙固结过程中孔隙压力逐渐消散的结果在同一高程上大小主应力在心墙部位有较明显的增长心墙底部单元的应力在蓄水以后仍有相当的增长分析心墙内部的应力状态可以发现各个单元的有效小主应力全部大于零不会出现拉伸破坏发生水力劈裂的可能性较小
图 9,7 <250m水位>竣工蓄水时的总位移矢量图
图 9,8 <250m水位>竣工蓄水时大主应力等值线单位10kPa
254 土质边坡稳定分析?原理
方法
程序
图 9,9 <250m水位>竣工蓄水时小主应力等值线单位10kPa
2,务平大坝坝基固结分析
有关务平工程的概况和孔隙水压分析成果同样也已在第6章作过介绍坝体竣工时的总位移等值线图9.10很好地反映基础和坝体的位移与沉降从整体上看软基的位移量比较小竣工时最大值为0.36m固结完成时最大值为0.4m位移最大值出现在软基最深厚处的软基表面由于软基厚度的减小整个位移有向河谷岸边减小的趋势坝体的位移与沉降呈现中间大两侧小的态势竣工时坝体的最大位移为0.82m出现在心墙1/3高度处固结完成时坝体的最大位移为0.9m
大坝固结完成时的大小主应力的等值线如图9.11和图9.12所示总体上看高程较低部分应力较大由于没有蓄水心墙部分的大小应力呈对称状态由于弹性模量大小的差异同一高度处心墙部分的大主应力明显小于两侧坝壳的大主应力小主应力则相差不大坝体的最大大主应力约为500kPa而坝基的最大大主应力达900kPa心墙和坝壳的剪应力都不大说明坝体主要以竖直方向的压缩为主坝体的应力水平较低一般都在0.7以下
软基部分的应力呈现显明的波浪形说明在碎石桩上存在应力集中正是这种应力集中使得处理后的软基能够满足工程的需求从图9.11和图9.12中可以看出软土的大主应力一般比同一高度的碎石桩上的值要小约100kPa随着位置的不同软基部分的大主应力的范围为300~900kPa
图 9,10 竣工时的总位移等值线图单位m
第9章 有限元法在边坡稳定分析中的应用 255
图 9,11 竣工时的大主应力等值线图单位kPa
图 9,12 竣工时的小主应力等值线图单位kPa
9,3 建立在滑裂面应力分析基础上的边坡稳定有限元分析
9,3,1 概述
几乎在有限元法开发的同时研究者就开始了其与边坡稳定分析中传统的条分法关系的研究早期的工作Wright,1973; Yamagam,1988等是在边坡中定义一个潜在的滑裂面根据有限元得出的应力分布计算滑裂面上各点的应力水平然后根据加权平均的原则定义安全系数他们发现虽然有限元计算仍处于弹性分析的水平但按照这一方法计算的安全系数和极限平衡方法十分接近
在这一工作的基础上研究者进一步研究在诸多滑裂面中寻找一个临界滑裂面的可能性这工作使用了与极限平衡法相同的计算步骤同时也获得了与极限平衡法接近的最小安全系数和临界滑裂面近期研究沿用传统的缩小强度指标的作法将非线性有限元数值计算分析收敛条件崩溃作为边坡失稳的判据此时强度指标被缩小的倍数即为安全系数
Griffith(1999)赵尚毅和郑颖人等(2002)在这方面作出了有益的探索按照这一方法求解首先需要具备一个成熟的非线性有限元计算程序否则无法判断数值计算发散的原因是边坡系统的物理原因还是有限元程序本身的缺陷
本节简要介绍上述有关的研究成果
9,3,2 安全系数的定义
在极限平衡法中整个滑面的安全系数定义为
256 土质边坡稳定分析?原理
方法
程序
τ
τ
f
F = (9.94)
式中τf为土的抗剪强度τ为土的剪应力安全系数即为抗滑力与滑动力的比值
根据沿滑裂面剪应力的计算方法安全系数大致有如下三种定义法
1,基于应力水平的定义法
保持不变
3
σ′作一个与摩尔库仑强度包线相切的应力圆称为破坏应力圆相应的应力圆直径为
f
)(
31
σσ ′?′
f
)/()(
3131
σσσσ ′?′′?′为应力水平如式(9.77)定义它表示当前应力圆直径与破坏应力圆直径之比反映了强度的发挥程度
整个滑面的安全系数定义为
∫
∫
′?′
′?′
=
l
l
F
f
FE
d
)(
d
31
31
1
σσ
σσ
(9.95)
2,基于剪应力的定义法
滑裂面上某点的应力状态为σx σy τxy法向应力σ′n和切向应力τ可由下式确定
ατασστ 2cos2sin)(
2
1
xyxy
+?= (9.96)
9.97) ατασασσ 2sincossin
22
xyyxn
+=′
式中α为此处滑面切线与x轴夹角
该点的抗剪强度由下式确定
(9.98) φστ ′′+′= tan
nf
c
式中c′
为材料的内聚力φ′ 为材料的内摩擦角
整个滑面的安全系数定义为
∫
∫
′′+′
=
l
lc
F
n
FE
d
d)tan(
2
τ
φσ
(9.99)
3,基于应力水平加权强度的定义法
Donald和Tan(1985)提出以如下公式计算一个滑裂面的安全系数方法式中各符号的含义同上
∫
∫
′′+′
′?′
′?′
′′+′
=
lc
lc
F
n
f
n
FE
d)tan(
)(
d)tan(
31
31
3
φσ
σσ
σσ
φσ
(9.100)
第9章 有限元法在边坡稳定分析中的应用 257
9,3,3 有关成果介绍
1,Donald 等的模式搜索法
Giam和Donald于1988年提出了一种由已知的应力场确定临界滑裂面及最小安全系数的方法称为CRISS法CRISS法的基本原理示于图9.13和图9.14由有限元法计算所得应力水平较高点P出发形成一个由顶到底的破坏面其过程如下将边坡分成若干土条
(stage)见图9.13由P点出发分别向上向下在相邻两土条内构筑局部滑裂面图9.14
给出了由P点向上在第N+1个土条内构筑局部滑裂面的过程即以P点为基点呈放射状按1°倾角变化形成m个条块每一条块分为长度为?L的五小段每一小段以其中点为控制点这些控制点的应力
xx
σ′
xy
σ′
xy
τ′由有限元计算结果而定这些控制点的正应力和剪应力计算如下
ατασσσσσ 2sin2cos)(
2
1
)(
2
1
xyxxyyxxyyni
′?′+′+′=′ (9.101)
ατασστ 2cos2sin)(
2
1
xyxxyyni
+′?′= (9.102)
由P点出发第i个条块对应的安全系数为
L
Lc
F
ni
ni
i
∑
′+′∑
=
)(
)tan(
τ
φσ
(9.103)
找出F
i
中的最小值F
imin
(
i
=
1,…,m
) Fimin对应的条块PQ即为第N+1个土条内的局部滑裂面位置由Q点为基点可以找到第N
+2个土条内的局部滑裂面位置依此类推到边坡顶部由P点向下找寻第N,N?1,…,1个土条内的滑裂面过程类同每一条块内的F
imin
对应的局部滑裂面即构成这一边坡的滑裂面其安全系数为
L
Lc
F
n
n
FE
d
d)tan(
2
τ
φσ
∫
∫
′′+′
= (9.104)
258 土质边坡稳定分析?原理
方法
程序
图 9,13 CRISS法求解过程简图
图 9,14 CRISS法的一个求解步
[例9.1] 简单边坡的稳定分析
图9.15为用CRISS法的F
FE2
计算结果此例为第11章11.4.2节介绍的ACAD边坡稳定分析程序考题1(a)详细数据可参见该节曾经使用各种不同的极限平衡分析程序进行计算得到一致的解答最小安全系数为1.0表9.1汇总了对应于不同安全系数定义法的有限元法计算结果可以看出一次填筑和分层填筑边坡的计算值非常接近开挖边坡的计算值略高所有计算结果和极限平衡法的解答一致
表 9,1 安全系数比较
安全系数
类型
FFE1 FFE2 FFE3
分层填筑1.001 1.001 1.001
一次填筑 1.000 1.003 1.000
开挖1.082 1.044 1.044
第9章 有限元法在边坡稳定分析中的应用 259
图 9,15 [例9.1] 简单边坡的CRISS法分析结果
2,Zou等(1995)的动态规划法
Zou等通过有限单元法获得的应力分布规律确定滑裂面的范围和初始滑裂面然后利用动态规划的数值方法(IDPM)搜索最危险滑裂面及相应的安全系数
在这一方法中滑裂面安全系数定义为
Li
Li
L
L
F
i
Ns
i
fi
Ns
i
B
A
f
B
A
s
≈=
∑
∑
∫
∫
=
=
τ
τ
τ
τ
1
1
d
d
(9.105)
式中τ为滑裂面AB上任一点的剪应力τ
f
为在这些点的材料抗剪强度L为滑裂面的长
度N
s
为沿滑裂面的分段数?L
i
为每一段的长度对于第i段
(9.106)φστ tan
sfi
c +=
(9.107) ατασστ 2cos2sin)(5.0
xyxyi
+?=
(9.108) ατασασσ 2sincossin
22
xyyxs
+=
式中c为材料的内聚力φ是材料的内摩擦角σ
s
为滑裂面上第i段对应法向应力α为
滑裂面第i段与水平面夹角σ
x
σ
y
τ
xy
为第i段对应的有效应力
引入附加函数G
F
则
(9.109) LFG
isfi
Ns
i
F
=
∑
=
)(
1
ττ
利用动态规划法通过求G
F
的极小值从而确定F
s
的最小值在选定的搜索区域内计算每一可能的滑裂面直至找到F
s
[例9.2] 一均质粘土边坡如图9.16其材料参数见表9.2表中最后一列最小安全系数为简化毕肖普法的计算结果
260 土质边坡稳定分析?原理
方法
程序
图 9,16 [例9.2] 均质粘土边坡示意图
表 9,2 材料参数
方案 凝聚力
(kPa)
摩擦角(°) 容重
(kN/m
3
)
杨氏模量(MPa) 泊松比 最小安全系数
1 3 20 20 200 0.25 1.00
2 10 20 20 200 0.25 1.39
利用平面有限单元法计算边坡的应力分布采用8结点等参单元边界条件为边坡两侧垂直边界只约束水平向位移边坡底部水平边界为固定边界材料假定为各向同性遵守
Moho?Coulomb准则的理想弹塑性材料图9.17给出了计算所得应力水平分布图
图 9,17 [例9.2] 边坡应力水平分布
根据图9.17的应力水平分布图Zou等采用图9.18图9.19所示两种状态点格栅图图
9.18分条方向为竖直向图9.19分条方向大致垂直于滑裂面分条厚为1m状态点间距
0.2m通过IDMP法计算发现,分条方向对计算结果影响不大两个方案对应安全系数分别为1.01和1.39与表9.2中所示简化毕肖普法的计算结果几乎一样图9.20图9.21给出第9章 有限元法在边坡稳定分析中的应用 261
了对应于表9.2中两个方案的IDMP分析结果及极限平衡条分法GWEDGEM的分析结果两种方法的结果也吻合
图 9,18 假定的滑裂面及竖直向分条的搜索范围
图 9,19 [例9.2] 假定的滑裂面及垂直于滑裂面分条的搜索范围
图 9,20 [例9.2] 方案1临界滑裂面
IDPM动态规划法GWEDGEM通用的极限平衡条分法FS安全系数
262 土质边坡稳定分析?原理
方法
程序
图 9,21 [例9.2] 方案2临界滑裂面
9,4 建立在强度缩小有限元分析基础上的边坡稳定分析
9,4,1 Griffiths和Lane的工作
建立在强度缩小有限元分析基础上的边坡稳定分析的基本原理是将边坡强度参数粘聚力c′
和内摩擦角φ′
同时除以一个折减系数F得到一组新的c′
和φ′
值然后作为一组新的材料参数输入再进行试算当计算不收敛时对应的F被称为坝坡的最小安全系数此时坝坡达到极限状态发生剪切破坏同时可得到临界滑动面
1,有限单元模型的简单介绍
本小节所提及的程序是建立于Smith和Griffiths 1998年开发的程序基础上的与其它程序相比主要的不同在于能够模拟更普遍的几何形状和土性的变化包括水位和孔压的变化并增加了图形化功能程序适用于摩尔库仑破坏准则下理想弹塑性土体的二维平面应变分析采用的是简化积分的八节点四边形单元每个单元有四个高斯点在重力荷载作用下刚度矩阵生成和应力再分配的算法中都采用这种单元假定土体开始为弹性的模型在网格内所有高斯点生成正应力和剪应力然后将这些应力与摩尔库仑破坏准则相比如果在特定高斯点上的应力位于摩尔库仑破坏圆内则该点仍然是弹性的如果位于圆上或圆外则该点处于屈服状态利用粘弹性算法屈服应力在网格中被重分配当足够数目的高斯点发生屈服使机制发生变化时即认为发生了整体剪切破坏本小节中提到的分析不能模拟拉裂
2,土体模型
本次研究采用的土体模型包括六个参数如表9.3所示
表 9,3 六个参数的土体模型
φ′ c′ ψ′ E′ v′ γ
摩擦角 凝聚力 剪胀角 杨氏模量 泊松比 容重
第9章 有限元法在边坡稳定分析中的应用 263
剪胀角在屈服过程中影响土体的体积变化众所周知在屈服过程中真实的土体体积是变化的例如一种中等密度的材料在剪切过程中也许初始体积是减小的(ψ
<0)随后是剪胀阶段(ψ
>0)最终是在体积为常量条件下的屈服(ψ 0)
随之而来的问题是ψ应取何值如果取ψ φ则塑性流动法则是相关联的并可与经典塑性理论进行对比当流动法则是相关联时应力和速度是一致的因此用有限元预测的破坏机制与用极限平衡法得到的临界破坏面基本一致
尽管使用相关联的流动法则有这些优点但我们知道非粘性土模型的相关联流动法则预测的剪胀要比现实中观察的大很多导致预测的破坏荷载偏大特别是在一些有限制的问题中如承载力问题(Griffiths 1982)
在稳定分析中研究对象常常是没有约束的天然边坡因此剪胀角的选择就变得不太重要这里研究的主要目的是准确地预测边坡的安全系数Griffiths等采用的是ψ 0的折中值对应屈服时体积改变量为零属非关联流动研究者发现ψ的取值可以使模型给出可靠的安全系数并合理地预测可能破坏面的位置和形状
参数c′
和φ′
指土体的有效凝聚力和摩擦角虽然有许多模拟土的强度的破坏准则但摩尔?库仑准则仍然是岩土工程实践中应用最广泛的一个也是此文所采用的根据主应力和以压为负的规定该准则可写成如下的形式
'cos
2
sin
2
3131
φ
σσ
φ
σσ
cF ′?
′?′
′
′+′
= (9.110)
式中
1
σ′和分别是最大和最小主应力
3
σ′
破坏函数解释如下F
0<F 应力在屈服圆内弹性
0=F 应力在屈服圆上屈服
0>F 应力在屈服圆外屈服并且必须重分配
弹性参数E′
和v′
是土体的杨氏模量和泊松比如果已经假定了泊松比(0.2< v′
< 0.3)
则杨氏模量与在一维固结中得到的土的压缩率有关由第6章式(6.17)可知
)1(
)21)(1(
ν
νν
ν
′?
′?′+
=′
m
E (9.111)
式中
ν
m是体积压缩率系数
虽然目前弹性参数的值对计算破坏前的变形有很大影响但对边坡稳定分析中的安全系数的影响却很小该文作者认为没有E′
和v′
的值时可以给一个象征性的值如E′
10
5
kN/m
2
v′
0.3
在计算过程中土体的总单位容重γ与由重力产生的节点自重荷载是成比例的
总而言之在有限元边坡稳定分析中采用的最重要的参数与传统的方法是一样的即容重γ强度参数c′
和φ′以及问题的几何形状
3,重力荷载
264 土质边坡稳定分析?原理
方法
程序
每一单元上由土体自重产生的力由以下积分式计算
(9.112)
∫
=
e
V
Te
vNp
)(
dγ
式中N是单元的形函数上标e代表单元号这个积分结果是将每个单元的面积与土的容重的乘积分配到各个节点上
4,安全系数的确定
一个土质边坡的安全系数是这样定义的为了使边坡达到破坏给原始的剪切强度参数除以一个数这个数就是安全系数安全系数这种定义和极限平衡法中的定义是相同的则经过折减的剪切强度参数和变为
f
c′
f
φ′
(9.113) Fcc
f
/′=′
′
=′
F
f
φ
φ
tan
arctan (9.114)
5,破坏的定义
在指定的最大迭代次数内如果算法不能收敛就意味着没有发现同时既能满足摩尔?
库仑破坏准则又能满足整体平衡的应力分布如果算法不能满足这些准则则说明破坏已经发生了边坡破坏和数值上的不收敛同时发生并且伴随着网格中节点位移的显著增加本文中大部分的计算使用的是最高限度为1000次的迭代结果以F与Eδmax/γH
2
无量纲的位移的图形形式显示其中δmax是收敛点的最大节点位移H是边坡的高度依靠位移的网格图和向量图来显示安全系数和破坏机制的特性
6,边坡稳定实例及验证
Griffiths提供了几个有限元边坡稳定分析的实例并在可能的条件下与传统的边坡稳定分析进行对比验证
[例9.3] 无地基层的均质边坡(D
=
1)
图9.22显示的均质边坡的土性如下20φ′=
o
/0.cHγ′ = 05边坡与水平面成
26.57°(2:1)左边界是水平向约束的而底部边界是完全固定的D的定义可参见图9.25在迭代限度内不断改变安全系数进行计算直到计算不收敛结果如表9.4所示
第9章 有限元法在边坡稳定分析中的应用 265
图 9,22 [例9.3] 坡角为26.57°的均质边坡的网格
表9.4显示的是6次试算得到的成果安全系数从0.8增大到1.4每一次的强度参数都除以相应安全系数值由于在每次分析中重力荷载和整体刚度矩阵都是相同的所以只要生成一次就可以了因而效率较高
如表9.4示当安全系数为0.8和1.0时只需进行2次和10次迭代计算就收敛越接近真实的安全系数算法就越难达到收敛当F 1.4时无量纲位移有一个突然增加
2
max
/EHδγ′
此时算法在允许迭代次数内不能收敛图9.23是根据表9.4的数据得到的安全系数与无量纲位移关系当F 1.4时位移突然增加计算不收敛意味着边坡破坏从该图可以看出对同样的问题有限元结果与由Bishop & Morgenstern得到的安全系数很接近
图9.24(a)和9.24(b)给出节点位移向量和对应于不收敛情况F 1.4时的变形的网格
图 9,23 [例9.3]和[例9.4] 安全系数与无量纲位移关系
266 土质边坡稳定分析?原理
方法
程序
图 9,24 [例9.3] 对应于不收敛解F 1.4的网格变形图
(a) 节点位移向量(b) 变形网格
表 9,4 [例9.3] 结果
F
2
max
/ HE γδ′ 迭代次数
0.80 0.379 2
1.00 0.381 10
1.20 0.422 20
1.30 0.453 41
1.35 0.544 792
1.40 1.476 1000
[例9.4] 有地基层的均质边坡(D
=
1.5)
图9.25显示的是在例9.3的边坡底面增加了一个厚度为H/2地基层的例子而其它的土性及几何形状都不变该例坡角为26.57° c′/γH
= 0.05 φ′
20° D = 1.5
第9章 有限元法在边坡稳定分析中的应用 267
图 9,25 [例9.4] 有地基层的均质边坡
(a) 没变形的网格(b) 与不收敛解F = 1.4相对应的网格
图9.25(a)和9.25(b)分别为初始网格和破坏时的变形网格从图9.25(b)可以清楚地获得
坡脚破坏类型的机制虽然由于可压缩性土体积增大位移也增加了但临界安全系数本质上并没有改变仍为F= 1.4有限元结果表明增加的地基层并没有明显改变边坡安全系数
Bishop 和 Morgenstern (1960)给出安全系数F
=
1.752作为该例D
=
1.5 c′/γH
=
0.05 φ′
20°坡度=2:1的一个可能解假定破坏时的圆弧与地基的底部相切用一个比较通用的滑弧程序也可以得到F 1.7这个可能的结果但要获得正确的安全系数1.4必须使滑弧通过坡脚这个例子说明了有限元边坡稳定分析法与传统边坡稳定分析方法相比的一个主要优点即有限元方法不要求假定临界面位置或形状破坏自然地发生在土的剪切强度不足以抵抗剪应力的地带因而不需要第4章所论述的搜索过程
[例9.5] 带有软弱夹层的不排水粘土边坡
下面的例子是关于不排水粘土边坡的稳定性分析在本例中采用总应力分析使用
Tresca破坏准则(φ
u
=
0)
图9.26显示的是在不排水粘土地基层上的边坡这个边坡包括一个软弱夹层这个软弱夹层一部分与边坡平行在地基的部分呈水平最后地面出露的部分与坡踵夹45°的角虽然这个例子看起来有点人为的痕迹但它有点像一个有软弱衬垫的废渣填埋系统的情况在计算该边坡的安全系数时保持不变25.0/
1
=Hc
u
γ只改变软弱薄层不排水抗剪强度
2u
c
268 土质边坡稳定分析?原理
方法
程序
图 9,26 [例9.5] 在地基层中有软弱夹层的不排水粘土边坡
有限元结果见图9.27对均质边坡(c
u2
/c
u1
=1)计算的安全系数与Taylor(1937)的F
=
1.47
的解很接近并且给出了期望的圆弧基本破坏机制随着薄层强度的逐渐减小在c
u2
/c
u1
≈0.6
时结果有很明显的改变
图 9,27 [例9.5] 不同c
u2
/c
u1
比值得到的安全系数F
该图还显示了用Janbu法得到的极限平衡结果该法同时采用了沿软弱层的圆弧和三线楔形破坏机制当c
u2
/c
u1
≈0.6时在由软弱层支配的圆弧和非圆弧之间的过渡阶段很清楚的显示了不连续性当c
u2
/c
u1
>0.6时为圆弧破坏机制支配安全系数本质上不受软弱薄层强度的影响当c
u2
/c
u1
<0.6非圆弧的薄层机制支配了破坏模式
图9.28可以很清楚地解释破坏模式图中显示了三种不同的c
u2
/c
u1
比值下破坏时的变形网格图9.28(a)示均质情况(c
u2
/c
u1
1)显示了一个本质上圆弧状的破坏机制图9.28(c)
中薄层的强度只有其周围土强度的20 (c
u2
/c
u1
0.2)从图中可以看出非圆弧机制与软弱薄层的路径几乎重合图9.28(b)中薄层的强度为其周围土强度(c
u2
/c
u1
0.6)的60图形复杂第9章 有限元法在边坡稳定分析中的应用 269
且模糊但是很明显至少有两种互相冲突的机制第一种是在坡脚以上与软弱层重合的破坏机制第二种是沿着与边坡和坡脚出露处平行的软弱层破坏的机制
图 9,28 [例9.5] 对应三种不同的c
u2
/c
u1
比值破坏时网格变形图
a c
u2
/c
u1
1.0; b c
u2
/c
u1
0.6; c c
u2
/c
u1
0.2
如果没有前面的关于两种不同机制的知识传统的极限平衡搜索出的安全系数可能会过高图9.27已经证明了这一点例如当c
u2
/c
u1
0.2时圆弧机制得到的F为1.3而正确的
F约为0.6
9,4,2 赵尚毅和郑颖等人的工作
上述Griffith等的研究成果受到国内外学者的重视赵尚毅和郑颖人等(2002)近期也获得有意义的成果他们采用了广义米赛斯准则
κα =+=
21
JIF (9.115)
式中I1 J2分别为应力张量的第一不变量和应力偏张量的第二不变量α κ是与岩土材料内摩擦角φ和内聚力c有关的常数不同的α κ 在 π
平面上代表不同的圆图9.29
这是一个通用的表达式通过变换α κ的表达式就可以在有限元中实现不同的屈服准则各准则的参数换算关系见表9.5
270 土质边坡稳定分析?原理
方法
程序
图 9,29 各屈服准则在π平面上的曲线
传统的极限平衡法采用摩尔?库仑准则但由于摩尔?库仑准则的屈服面为不规则的六角形截面的角锥体表面存在尖顶和菱角给数值计算带来困难为了与传统方法进行比较这里采用徐干成郑颖人等(1990)提出的摩尔库仑等面积圆屈服准则(DP4)代替传统的摩尔库仑准则
表 9,5 各准则参数换算表
编号 准则种类 α κ
DP1 外角点外接D-P圆
()φ
φ
sin33
sin2
()φ
φ
sin33
cos6
c
DP2 内角点外接D-P圆
()φ
φ
sin33
sin2
+
()φ
φ
sin33
cos6
+
c
DP3 内切D-P圆
φ
φ
2
sin33
sin
+
φ
φ
2
sin33
cos3
+
c
DP4 等面积D-P圆
()φπ
φ
2
sin932
sin32
()φπ
φ
2
sin932
cos36
c
[例9.6] 均质土坡稳定分析算例
此例为一均质土坡坡高H
=
20m土容重
=
25kN/m
3
粘聚力c
=
42kPa内摩擦角
=
17°
表 9,6 用不同方法求得的安全系数
坡角(°) 30 35 40 45 50
有限元法外接圆屈服准则 1.78 1.62 1.48 1.36 1.29
有限元法摩尔-库仑等面积圆屈服准则
1.47 1.34 1.22 1.12 1.06
Bishop 1.39 1.26 1.15 1.06 0.99
Spencer 1.46 1.32 1.21 1.12 1.04
第9章 有限元法在边坡稳定分析中的应用 271
计算结果见表9.7可见DP4准则与简化Bishop法所得的稳定安全系数最为接近其误差的平均值为5.7%最大误差小于8%且离散度很小而DP1准则的平均误差为29.5
采用DP2和DP3准则所得计算结果的离散度也很大因此在数值分析中可以用DP4准则代替摩尔库仑准则
表 9,7 不同屈服准则所得最小安全系数
H
=
20m β=45
o
c
=
42kPa
φ (
o
) 0.1 10 25 35 45
DP1 0.525 1.044 1.769 2.254 3.051
DP2 0.525 0.930 1.332 1.530 1.887
DP3 0.454 0.848 1.279 1.499 1.870
DP4 0.477 0.879 1.396 1.689 2.182
简化Bishop法 0.494 0.846 1.316 1.623 2.073
(DP1-Bishop)/Bishop 0.063 0.234 0.344 0.355 0.472
(DP2-Bishop)/Bishop 0.063 0.099 0.012?0.080?0.090
(DP3-Bishop)/Bishop?0.081 0.002?0.028?0.099?0.098
(DP4-Bishop)/Bishop?0.034 0.059 0.061 0.041 0.053
研究者还发现边界范围的大小在有限元法中对计算结果的影响比在传统极限平衡法中表现的更为敏感当坡角到左端边界的距离为坡高的1.5倍坡顶到右端边界的距离为坡高的2.5倍且上下边界总高不低于2倍坡高时计算精度最为理想另外如果网格划分太粗将会造成很大的误差计算时必须考虑适当的网格密度
从表9.7中计算结果可以看出采用外接圆屈服准则计算的安全系数比传统的方法大许多而采用摩尔?库仑等面积圆屈服准则计算的结果与传统的极限平衡方法Spencer法计算的结果十分接近说明采用摩尔?库仑等面积圆屈服准则代替摩尔?库仑不等角六边形屈服准则是可行的
另外通过4组方案的算例的比较分析表明用摩尔?库仑等面积圆屈服准则求得的安全系数与Bishop法的误差为3%~7%与Spencer法的误差为1%~3%说明有限元强度折减法完全在边坡稳定分析中应用的可行性
9,4,3 NDM结点位移法
Tan和Donald于1985年提出了一种利用有限元求解得到的结点位移来确定安全系数的简单图解法称为结点位移法这一方法是将结构的组成材料的粘聚力c和内摩擦系数乘以一个Nφtan=f
(N<1.0)
的折减系数进而通过有限元计算跟踪某一结点的位移δ
T
的变化随N的减小δ
T
不断增大如图9.30所示尖角点给出了N的临界值即为安全系数
T
F
δ
当各参数设计值乘以N之后结构物就可达到极限平衡 1/N即为设计安全系数为识别可能的滑裂面必须清楚在潜在破坏区域内结点的位移状况这可以很容易地从位移矢量图中观察到
272 土质边坡稳定分析?原理
方法
程序
这种方法包含了如下假定
1) 材料强度参数的降低不影响它的其它性质
2) 尽管有限元分析是建立在小变形理论基础上的但在出现尖角点前后大变形δ
T
的计算是有效的
图 9,30 [例9.7] 节点位移随1/N变化过程图
[例 9.7] Talbingo 坝
图9.31所示高为162米的Talbingo心墙坝强度参数从100%的设计值 (1/N
=
1.0)降到
40%(1/N
=
2.50)图9.32示一典型结点的位移变化过程F
δT
的范围大约在2.05~2.14之间对应于1/N
=
2.06和2.11的位移矢量参见图9.33在1/N
=
2.06时无明显大变形δ
T
出现坝中无明显滑裂面在1/N
=
2.11时在坝的上游坡堆石料区过渡料区和心墙料区可以很明显看到大变形出现由图 9.33(b)可以看到一明显的滑裂面用条分法计算所得通过心墙的滑裂面的安全系数为2.06
图 9,31 Talbingo坝网格及观测点
第9章 有限元法在边坡稳定分析中的应用 273
图 9,32 [例9.7] 151号节点位移变化过程
根据位移矢量图9.33(b)假定一圆弧滑裂面表9.8列出了各种方法计算的稳定分析结果可以看出与F
T
F
δ
FE1 FFE3很接近但略小于FFLE和FFE2各安全系数之间的一致性可以证明用δT与F来确定土石坝的安全系数是可行的
T
δ
图 9,33 [例9.7] Talbingo坝位移矢量图
(a) 1/N=2.06 (b) 1/N=2.11
表 9,8 安全系数比较
临界圆弧信息 总应力安全系数
x y r F
FE1
F
FE2
F
FE3
F
FLE
T
F
δ
100.3 290.0 278.8 2.10 2.23 2.11 2.21 2.11
274 土质边坡稳定分析?原理
方法
程序
注 坐标原点为上游坡与坝底的交点
参考文献
1 Biot,M,A,Theory of elasticity and consolidation for a porous anisotropic solid,Journal
of Applied Physics 1955,26,182-185
2 Bishop,A,W,and Morgenstern,N,R,Stability coefficients for earth slopes,Geotechnique,
1960,10:129-150
3 Cathie,D,N.,and Dungar,R,Evaluation of finite element predictions for constructional
behavior of a rockfill dam,Proc,Instn,Civ,Engrs.,Part 2,1978,65:551-568
4 Charles,J,A,The use of one dimensional compression tests and elastic in prediction
deformation of embankments,Canadian Geotechnioal Journal.1976,13,189-200
5 Chen L.H,etc,Consolidation and stability behavior of a high rockfill dam built on soft
clay foundation Proceedings of 3rd International Conference on Soft Soil Engineering ICSSE3,
2001,DEC,Hong Kong
6 陈胜宏,高坝复杂岩石地基及岩石高边坡稳定分析,北京中国水利水电出版社,2001
7 陈祖煜,用有限元进行渗流结构和固结计算中的基本原理,水利水电科学研究院,1985
8 陈祖煜陈立宏李新强,第7章 土工数值分析(二)渗流固结和应力应变的有限元方法高等土力学2003清华大学出版社,
9 Donald,I,B,Tan,C,P,And Goh,T,C,A,Stability of geomechanical structures assessed
by finite element method,Proc,2nd Int,Conf,In Civil Engr,Hangzhou,Science Press,Beijing,
1985,845-856
10 Duncan,J,M,et,al,CON2D:A finite element computer program for analysin of consolidation,
Report No,VCB1GT/81-01,University of California,Berkeley,1981
11 Duncan,J,M.,and Chang,C,Y,Nonlinear analysis of stress and strain in soils,J,Soil
Mech,and Found,Div,ASCE,1970,96(5),1629-1653
12 Donald,I,B,Finite element assessment of slope stability,Austrlian-China Landslide Seminar,
1993,3-18 July,1-42.D1-D18
13 Eisensein,Z,and Simmons,J,V,Three –dimensional analysis of Mica Dam,Proc,Int,Symp,
Criteria and Assumptions for Numer,Anal,of Dams,Quadrant,Swansea,U,K.1975,1051-1070
14 Griffith,D,V,and Lane,P,A,Slope stability analysis by finite elements,Geotechnique,
1999,49(3),387-403
15 Giam,S,K,and Donald,I,B,Determination of critical slip surfaces for slopes via
stress-strain calculations,5th A,N,Z,Conf,Geomechanics,Sydney,1988,461-464
16 Penman,A,D,M.,Burland,J,B.,and Charles,J,A,Observed and Predicted deformations in
a large embankment cam during con-struction,Proc,Ins,Civ,Engrs.,London,1971,49:1-21
17 Penman,A.,and Charles,A,Construction deformations in rock-fill dam,J,Soil Mech,And
Found,Div.,ASCE,1973.99(2),139-163
18 沈珠江,计算土力学,上海科技出版社,1990
19 Tan,C,P,and Donald,I,B,Finite element calculation of dam stability,Proc,11th Int,
Conf,Soil Mech,and Fnd,Engr,San Francisco,1980
20 Taylor,D,W,Stability of earth slopes,J,Boston Soc,Civ,Eng,1937,24:197-246
21 Wright,S,G,Kulhawy,F,H,and Duncan,J,M,Accuracy of equilibrium slope stability analysis,
J,Soil Mech,And FDN,Div,1973,99(SM10),783-791
第9章 有限元法在边坡稳定分析中的应用 275
22 Wang,C,and Duncan,J,M,Analysin of consolidation of earth and rockfill dans,Report No,
TE77-3,University Of California,Berkeley,1977
23 徐干城郑颖人,岩土工程中屈服准则应用的研究,岩土工程学报,1990,12(2):93-99
24 Yamagami,T,and Ueta,Y,Search for critical slip lines in finite element stress fields
by dynamic programming,Proc,6th International Conference on Numerical Methods in
Geomechanics,Innsbruck,Australia,1988,1335-1339
25 Zou,J,Z.,Williams,D,J.,and Xiong,W,L,Search for critical slip surfaces based on
finite element method,Canadian Geotechnical Journal,1995,32,233-246
26 赵尚毅郑颖人时卫民王敬林,有限元强度折减法求边坡稳定安全系数,岩土工程学报,2002,
129(3),343-346
