Slope Stability Analysis by Strength Reduction FEM

Abstract

An analysis method for slope safety factor through c - cp reduction algorithm by finite elements is presented. When the system reaches instability, the numerical non-convergence occurs simultaneously. The safety factor is then obtained by c - ψ reduction algorithm. The same time, the critical failure surface is found automatically. The traditional limit equilibrium method can't get the safety factor and failure surface of jointed rock slope. Strength Reduction FEM (SRFEM) presents a powerful alternative approach for slope stability analysis, especially to jointed rock slope. This paper analyzes the precision and error caused by different soil yield criterions、FEM itself、slope height、slope angle´cohesion and friction angle thoroughly. Through a series of case studies, the results show average error of safety factor between Strength Reduction FEM and traditional limit equilibrium method (Bishop simplified method) is 3 %～8 %，the error between SRFEM and Spencer's method is 1 %～4 % . The applicability of the proposed method was clearly exhibited.

Content

《1 引言》

1 引言

西部开发是我国实现地区平衡发展和可持续发展的重大战略举措。然而, 我国西部地区山高坡陡、沟壑纵横, 城市建筑依山而立, 公路、铁路翻山越岭, 复杂多变的地形地貌决定了我国西部开发将面临大量滑 (边) 坡工程, 滑坡与边坡事故日益增多。例如2001年, 重庆市云阳县就发生两次大型滑坡, 重庆市武隆边坡失稳坍塌造成79人死亡。重庆市近20年来累计发生地质灾害3.33万处, 仅2000年就发生6371处, 受灾19.33万人, 倒塌房屋8.68万间, 直接经济损失7.67亿元。重庆市已经成为地质灾害的重灾区, 尤其是随着三峡库区蓄水和新兴城市的建设, 有可能诱发更大的地质灾害, 隐患无穷。现在, 国务院已经决定拨款40亿元, 用于三峡库区地质灾害治理, 仅重庆三峡库区计划的地质灾害治理工程就有143个。频发的地质灾害以及大量灾害治理资金的投入, 使得滑 (边) 坡稳定性问题成为西部开发中的热点与难点问题。

边坡稳定分析是经典土力学最早试图解决而至今仍未圆满解决的课题, 各种稳定分析方法在国内外水平大致相当。对于均质土坡, 传统方法主要有:极限平衡法、极限分析法和滑移线场法等。就目前工程应用而言, 主要还是极限平衡法, 但需要事先知道滑动面位置和形状。对于均质土坡, 可以通过各种优化方法来搜索危险滑动面。但是, 对于岩质边坡, 由于实际岩体中含有大量不同构造、产状和特性等不连续结构面 (比如层面、节理、裂隙、软弱夹层、岩脉和断层破碎带等) , 给岩质边坡的稳定分析带来了巨大的困难。传统极限平衡方法尚不能搜索出危险滑动面以及相应的稳定安全系数。而目前的各种数值分析方法, 一般只是得出边坡应力、位移、塑性区, 也无法得到边坡危险滑动面以及相应的安全系数。随着计算机技术的发展, 尤其是岩土材料的非线性弹塑性有限元计算技术的发展, 有限元强度折减法近来在国内外受到关注, 对于均质土坡已经得到了较好的结论, 但尚未在工程中实用。笔者采用有限元强度折减法 ^[1], 对均质土坡进行了系统分析, 证实了用于工程的可行性, 得到了节理岩质边坡坡体的危险滑动面和相应的稳定安全系数。该方法可以对贯通和非贯通的节理岩质边坡进行稳定分析, 同时可以考虑地下水、施工过程对边坡稳定性的影响, 可以考虑各种支挡结构与岩土材料的共同作用, 为岩质边坡稳定分析开辟了新的途径。

《2 有限元强度折减系数法基本原理》

2 有限元强度折减系数法基本原理

有限元强度折减系数法的基本原理是将坡体强度参数:粘聚力c和内摩擦角φ值同时除以一个折减系数F_trial, 得到一组新的c′、φ′值, 然后作为新的资料参数输入, 再进行试算, 当计算不收敛时, 对应的F_trial被称为坡体的最小稳定安全系数, 此时坡体达到极限状态, 发生剪切破坏, 同时可得到坡体的破坏滑动面。

$c^{'} = c / F_{t r i a l}, φ^{'} = \arctan (\tan φ / F_{t r i a l})$ $c^{'} = c / F_{t r i a l}, φ^{'} = \arctan (\tan φ / F_{t r i a l})$

《3 有限元强度折减系数法精度分析》

3 有限元强度折减系数法精度分析

《3.1屈服准则的选用》

3.1屈服准则的选用

安全系数大小与程序采用的屈服准则密切相关, 不同的准则得出不同的安全系数。目前流行的大型有限元软件ANSYS, 以及美国MSC公司的MARC、PATRAN、NASTRAN均采用了广义米赛斯准则:

$F = α Ι_{1} + \sqrt{J_{2}} = k ‚$ $F = α Ι_{1} + \sqrt{J_{2}} = k ‚$

式中I₁, J₂分别为应力张量的第一不变量和应力偏张量的第二不变量。α、k是与岩土材料内摩擦角φ和内聚力c有关的常数, 这是一个通用表达式, 通过变换α、k的表达式就可以在有限元中实现不同的屈服准则, 各准则的参数换算关系见表1。传统的极限平衡法采用摩尔-库仑准则, 但是由于摩尔-库仑准则的屈服面为不规则六角形截面的角锥体表面, 存在尖顶和菱角, 给数值计算带来困难。为了与传统方法进行比较, 本文采用了徐干成、郑颖人 (1990) 提出的摩尔-库仑等面积圆屈服准则 (DP4) 代替传统摩尔-库仑准则 ^[2]。

表1各准则参数换算表

Table 1 The relationship of different yield criterions

《表1》

编号	准则种类	α	k
DP1	外角点外接D-P圆	$\frac{2 \sin φ}{\sqrt{3} (3 - \sin φ)}$ $\frac{2 \sin φ}{\sqrt{3} (3 - \sin φ)}$	$\frac{6 c c o s φ}{\sqrt{3} (3 - \sin φ)}$ $\frac{6 c c o s φ}{\sqrt{3} (3 - \sin φ)}$
DP2	内角点外接D-P圆	$\frac{2 \sin φ}{\sqrt{3} (3 + \sin φ)}$ $\frac{2 \sin φ}{\sqrt{3} (3 + \sin φ)}$	$\frac{6 c c o s φ}{\sqrt{3} (3 + \sin φ)}$ $\frac{6 c c o s φ}{\sqrt{3} (3 + \sin φ)}$
DP3	内切D-P圆	$\frac{\sin φ}{\sqrt{3} \sqrt{3 + \sin^{2} φ}}$ $\frac{\sin φ}{\sqrt{3} \sqrt{3 + \sin^{2} φ}}$	$\frac{3 c c o s φ}{\sqrt{3} \sqrt{3 + \sin^{2} φ}}$ $\frac{3 c c o s φ}{\sqrt{3} \sqrt{3 + \sin^{2} φ}}$
DP4	等面积D-P圆	$\frac{2 \sqrt{3} \sin φ}{\sqrt{2} \sqrt{3} π (9 - \sin^{2} φ)}$ $\frac{2 \sqrt{3} \sin φ}{\sqrt{2} \sqrt{3} π (9 - \sin^{2} φ)}$	$\frac{6 \sqrt{3} c c o s φ}{\sqrt{2} \sqrt{3} π (9 - \sin^{2} φ)}$ $\frac{6 \sqrt{3} c c o s φ}{\sqrt{2} \sqrt{3} π (9 - \sin^{2} φ)}$

算例分析表明 (表2) :DP4准则与简化Bishop法所得稳定安全系数最为接近。通过对误差进行统计分析可知, 当选用DP4准则时, 误差的平均值为5.7%, 最大误差小于8%, 且离散度很小, 而DP1的平均误差为29.5%, 同时采用DP2、DP3准则所得计算结果的离散度也非常大。因此, 在数值分析中可用DP4准则代替摩尔-库仑准则。

表2不同屈服准则所得最小安全系数

Table 2 Safety factor by different method

《图1》

《3.2有限元计算精度分析》

3.2有限元计算精度分析

边界范围的大小在有限元法中对计算结果的影响比传统极限平衡法表现得更为敏感, 当坡角到左端边界的距离为坡高的1.5倍, 坡顶到右端边界的距离为坡高的2.5倍, 且上下边界总高不低于2倍坡高时, 计算精度最为理想。另外如果网格划分太粗, 将会造成很大的误差, 计算时必须考虑适当的网格密度。

《4 均质土坡稳定分析》

4 均质土坡稳定分析

均质土坡, 坡高h=20 m, 土容重γ=25 kN/m³, 粘聚力c=42 kPa, 内摩擦角φ=17°, 求坡角β=30°, 35°, 40°, 45°, 50°时边坡的稳定安全系数。

如图1, 按照平面应变建立有限元模型, 边界条件为左右两侧水平约束, 下部固定, 上部为自由边界, 计算结果见表3。

从表3看出, 采用外接圆屈服准则计算的安全系数比传统的方法大许多, 而采用摩尔-库仑等面积圆屈服准则计算的结果与传统极限平衡方法 (Spencer法) 计算的结果十分接近, 说明采用摩尔-库仑等面积圆屈服准则来代替摩尔-库仑不等角六边形屈服准则是可行的。

《图2》

图1有限元单元网格划分

Fig.1 FEM model

表3用不同方法求得的稳定安全系数

Table 3 Safety factor by different method

《表2》

方法	坡角 / (°)
方法	30	35	40	45	50
DP1	1.78	1.62	1.48	1.36	1.29
DP4	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

另外, 通过4组计算方案 (改变内摩擦角φ、粘聚力c、坡角β、坡高h的值) 共计106个算例的比较分析表明, 用摩尔-库仑等面积圆屈服准则求得的安全系数与Bishop法的误差为3%～8%, 与Spencer法的误差为1%～4%, 说明了有限元强度折减法完全可以实用于土坡工程。

《5 岩质边坡稳定分析》

5 岩质边坡稳定分析

根据岩体中结构面的贯通情况, 可以将结构面分为贯通性、半贯通性、非贯通性三种类型。根据结构面的胶结和充填情况, 可以将结构面分为硬性结构面 (无充填结构面) 和软弱结构面。由于岩体结构的复杂性, 要十分准确地反映岩体结构的特征并使之模型化是不可能的, 也没有必要使问题复杂化。基于这种考虑, 对于一个实际工程来说, 往往根据现场地质资料, 根据结构面的长度、密度、贯通率, 展布方向等, 着重考虑2～3组对边坡稳定起主要控制作用的节理组或其他主要结构面。

岩体是弱面体, 起控制作用的是结构面强度, 对于软弱结构面, 可采用低强度实体单元模拟, 按照连续介质处理;对于无充填的硬性结构面可以采用无厚度的接触单元来模拟, 安全系数的求解与均质土坡相同。结构面单元的设置会影响计算精度, 一般可先以较大间距设置结构面, 求出滑动面大概位置后再在滑动面附近将结构面加密, 增加结构面单元以提高计算精度。

《5.1具有一组平行节理面的岩质边坡算例》

5.1具有一组平行节理面的岩质边坡算例

如图2和图3所示, 一组软弱结构面倾角40°, 间距10 m, 岩体和结构面采用平面6节点三角形单元模拟。岩体以及结构面材料物理力学参数取值见表4。采用不同方法的计算结果见表5, 其中极限平衡方法计算结果是在滑动面确定的情况下算出的。

《图3》

图2几何模型

Fig.2 Geometry model

《图4》

图3有限元模型

Fig.3 FEM model

表4计算采用的物理力学参数

Table 4 Material properties

《表3》

材料名称	容重 /kN·m^-3	弹性模量/Pa	泊松比	内聚力 /MPa	内摩擦角/ (°)
岩体	25	1E+10	0.2	1.0	38
结构面	17	1E+7	0.3	0.12	24

通过有限元强度折减计算, 当有限元计算不收敛时, 程序自动找出了滑动面, 如图4。在一组平行的结构面中, 只出现了一条滑动面, 其余结构面没有出现塑性区和滑动。

表5计算结果

Table 5 Safety factor by different method

《表4》

计算方法	安全系数
有限元法 (外接圆屈服准则)	1.26
有限元法 (等面积圆屈服准则)	1.03
极限平衡方法 (解析解)	1.06
极限平衡方法 (Spencer)	1.06

《图5》

图4坡体达到极限状态时形成的滑动面

Fig.4 The failure surface at limited state

《5.2具有两组节理面的岩质边坡算例》

5.2具有两组节理面的岩质边坡算例

如图5所示, 两组方向不同的节理, 贯通率100%, 第一组软弱结构面倾角30°, 平均间距10 m;第二组软弱结构面倾角75°, 平均间距10 m, 岩体以及结构面计算物理力学参数见表6。

《图6》

图5几何模型

Fig.5 Geometry model

按照二维平面应变问题建立有限元模型, 在单元划分的过程中, 在两个滑动面的交汇处形成了尖角, 在建模时应在曲率突然变化的区域使用更细的网格以保证计算的准确性, 或者将尖角抹圆。计算步骤同上, 通过有限元强度折减, 求得的滑动面如图6a所示, 它是最先贯通的塑性区。塑性区贯通并不等于破坏, 当塑性区贯通后塑性发展到一定程度, 岩体发生整体破坏, 同时出现第二、三条贯通的塑性区, 如图6b。程序还可以动画模拟边坡失去稳定的过程, 从动画演示过程可以看出边坡的破坏过程也就是塑性区逐渐发展、最后贯通形成整体破坏的过程。求得的稳定安全系数见表7。其中, 极限平衡方法计算结果是根据最先贯通的那一条滑动面求得的。

表6 物理力学参数计算取值Table 6 Material properties

《图7》

《图8》

图6 极限状态后产生的滑动面和塑性区Fig.6 Failure surface and plastic zone

表7计算结果

Table 7 Safety factor by different method

《表5》

计算方法	安全系数
有限元法 (外接圆屈服准则)	1.62
有限元法 (等面积圆屈服准则)	1.33
极限平衡方法 (Spencer)	1.36

《6 结论》

6 结论

1) 有限元强度折减法不需要对滑动面形状和位置做假定, 通过强度折减使边坡达到不稳定状态时, 非线性有限元静力计算将不收敛, 此时的折减系数就是稳定安全系数, 同时可得到边坡破坏时的滑动面。本文对其计算精度进行了分析, 算例表明采用徐干成、郑颖人 (1990) 提出的摩尔-库仑等面积圆屈服准则求得的稳定安全系数与简化Bishop法的误差为3 %～8%, 与Spencer法的误差为1 %～4 %, 证实了其实用于工程的可行性。

2) 目前对复杂节理岩质边坡的稳定分析尚没有好的办法, 传统的极限平衡方法无法得到岩质边坡的滑动面及其稳定安全系数, 而各种数值分析方法只能算出应力、位移、塑性区等, 无法判断边坡的稳定安全系数以及相应的滑移面。本文利用非线性有限元强度折减系数法由程序自动求得边坡的危险滑动面以及相应的稳定安全系数。通过算例分析表明了此法的可行性, 为岩质边坡稳定分析开辟了新的途径。

3) 屈服条件的选择、计算模型的建立, 包括计算范围、边界条件、网格划分密度和收敛标准等, 应满足有限元计算的精度要求。如果使有限元计算保持足够的计算精度, 那么有限元法较传统的方法具有如下优点:

a.能够对具有复杂地貌、地质的边坡进行计算;

b.考虑了土体的非线性弹塑性本构关系, 以及变形对应力的影响;

c.能够模拟土坡的失稳过程及其滑移面形状。如图7、图8, 可见滑移面大致在水平位移突变的地方, 也是在塑性区塑性发展最充分的地方, 呈条带状;

d.能够模拟土体与支护的共同作用, 图9为无锚杆 (锚杆单元被杀死) 时边坡稳定安全系数为1.1, 图10为有锚杆支护时安全系数为1.5, 且塑性区后移。

e.求解安全系数时, 可以不需要假定滑移面的形状, 也无需进行条分。

《图9》

图7水平方向位移等值云图

Fig.7 Continuous contours of X displacement

《图10》

图8均质土坡达到极限状态时形成的滑动面

Fig.8 The failure surface at limited state of soil slope

《图11》

图9不加锚杆ω=1.1时的塑性区

Fig.9 The plastic strain zone without support

《图12》

图10加锚杆ω=1.5时的塑性区

Fig.10 The plastic strain zone with support

Show More