有限元强度折减系数法计算土坡稳定安全系数的精度研究(1.后勤工程学院 军事土木工程系)
自弗伦纽期于1927年提出圆弧滑动法以来，至今已出现数十种土坡稳定分析方法，有极限平衡法、极限分析法、有限元法等。不少研究表明，各种方法所得稳定安全系数都比较接近，可以说，这些方法已经达到了相当高的精度。 近年来，由于计算机技术的长足发展，基于有限元的折减系数法在边坡稳定分析中的应用备受重视。与极限平衡法相比，它不需要任何假设，便能够自动地求得任意形状的临界滑移面以及对应的最小安全系数，同时它还可以真实的反映坡体失稳及塑性区的开展过程。到目前为止，已有很多学者对折减系数法进行了较为深入的研究[1,2,3]，并在一些算例中得到了与极限平衡法十分接近的结果。但总体说来，此法仍未在工程界得到确认和推广，究其原因在于影响该法计算精度的因素很多，除了有限元法引入的误差外，还依赖于所选用的屈服准则。
此论文的目的有两点：(1)力图全面分析屈服条件和有限元法本身对折减系数法计算精度的影响，并提出应选用何种屈服准则以及提高有限元法计算精度的具体措施；(2)结合工程实例，分析对边坡稳定安全系数影响最大的4个主要参数(H坡高、β坡角、C粘聚力、Φ摩擦角)对折减系数法计算精度的影响。从以往的计算结果来看，严格法(Spencer)所得稳定安全系数比简化Bishop法平均高出约2%～3%，而通过106个算例的比较分析，表明：折减系数法所得稳定安全系数比简化Bishop法平均高出约5.7%，且误差离散度极小，可以认为是正确的解答[4]。这有力地说明了将有限元折减系数法用于分析土坡稳定问题是可行的，但必须合理地选用屈服条件以及严格地控制有限元法的计算精度，同时也表明：有限元折减系数法所得安全系数稍微偏高，其原因有待进一步研究。
1 折减系数法的基本原理
Bishop等将土坡稳定安全系数F定义为沿整个滑移面的抗剪强度与实际抗剪强度之比，工程中广为采用的各种极限平衡条分法便是以此来定义坡体稳定安全系数。有限元强度折减系数法的基本思想与此一致，两者均可称之为强度储备安全度。因后者无法直接用公式计算安全系数，而需根据某种破坏判据来判定系统是否进入极限平衡状态，这样不可避免地会带来一定的人为误差。尽管如此，仍发展了一些切实可行的平衡判据，如：限定求解迭代次数，当超过限值仍未收敛则认为破坏发生；或限定节点不平衡力与外荷载的比值大小；或利用可视化技术，当广义剪应变等值线自坡角与坡顶贯通则定义坡体破坏[3]。文中平衡判据取：当节点不平衡力与外荷载的比值大于10-3时便认为坡体破坏。
有限元折减系数法的基本原理是将土体参数 C、Φ 值同时除以一个折减系数 Ftrial，得到一组新的C′、Φ′值，然后作为新的材料参数带入有限元进行试算，当计算正好收敛时，也即 Ftrial再稍大一些(数量级一般为10-3)，计算便不收敛，对应的Ftrial被称为坡体的最小安全系数，此时土体达到临界状态，发生剪切破坏，具体计算步骤可参考文献[2]，文中如无特别说明，计算结果均指达到临界状态时的 折减系数。
| (1) |
| (2) |
2 屈服准则的影响
用折减系数法求解实际边坡稳定问题时，通常将土体假设成理想弹塑性体，其中本构模型常选用摩尔-库仑准则(M-C)、Drucker-Prager准则以及摩尔-库仑等面积圆[5]准则。
摩尔-库仑准则可用不变量I1,J2,θσ表述成如下形式：
| (3) |
Drucker-prager准则：
| (4) |
式中：I1为应力张量第一不变量；J2为应力偏量第二不变量；θσ是应力洛德角。
|
图1 各屈服准则在π平面上的曲线
|
M-C准则较为可靠，它的缺点在于三维应力空间中的屈服面存在尖顶和棱角的不连续点，导致数值计算不收敛，所以有时也采用抹圆了的M-C修正准则[6]，它是用光滑连续曲线来逼进摩尔-库仑准则，此法虽然方便了数值计算，但不可避免地会引入一定的误差；而D-P准则在偏平面上是一个圆，更适合数值计算。通常取M-C准则的外角点外接圆、内角点外接圆或其内切圆作为屈服准则，以利数值计算。各准则的参数换算关系见表1。
由徐干成、郑颖人(1990)提出的摩尔库仑等效面积圆准则[5]实际上是将M-C准则转化成近似等效的D-P准则形式。该准则要求偏平面上的摩尔-库仑不等边六角形与D-P圆面积相等。计算表明它与摩尔-库仑准则十分接近。
见图1，r1为外角外接圆半径；r2为内角外接圆半径；r3为内切圆半径；摩尔-库仑准则构成的六角形面积为
| (5) |
对半径为r的圆面积S=πr2，令S=Smorl得
| (6) |
| (7) |
式(7)与式(4)对应项相等，可得
| (8) |
|
表1 各准则参数换算 编号 | 准则种类 | αφ | k | DP1 | 外角点外接D-P圆 | 
| 
| DP2 | 内角点外接D-P圆 | 
| 
| DP3 | 内切D-P圆 | 
| 
| DP4 | 等面积D-P圆 | 
| 
| 注：表中αφ、k是与D-P有关的材料参数。 | | 表2 不同屈服准则所得最小安全系数 | φ/° | | 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.896 | 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 | 注：H=20mm；β=45°；C=42kPa。 | |
算例分析表明(表2、图2)：DP4准则与简化Bishop法所得稳定安全系数最为接近。对有效算例(Φ≠0)的误差进行统计分析可知，当选用DP4准则时，误差的平均值为5.7%，且离散度很小(图3)。而DP1的平均误差为29.5%，同时采用DP2、DP3准则所得计算结果的离散度非常大，均不可用。因此在数值分析中可用DP4准则代替摩尔-库仑准则。
| |
图2 Φ～折减系数曲线
|
图3 DP4准则的计算误差
|
3 不同流动法则的影响
有限元计算中，采用关联还是非关联流动法则，取决于ψ值(剪胀角)：ψ=φ，为关联流动法则；ψ≠0，为非关联流动法则。总体说来，采用非关联流动法则所得破坏荷载比同一类型材料而采用关联流动法则所得破坏荷载小，如忽略剪胀角(ψ=0)，将会得到较为保守的结果。值得注意的是：当ψ=0时，正好与郑颖人等提出的广义塑性力学理论相符[7]，这时对应的塑性势面与q轴垂直。
表3 不同流动法则的影响 | φ=10° | φ=17° | φ=25° | 非关联 | 0.871 | 1.105 | 1.363 | 关联 | 0.887 | 1.137 | 1.425 | 相对误差 | 0.018 | 0.029 | 0.045 | β=45°；C=40kPa；H=20m；DP4准则。 | | 表4 网格疏密对计算结果的影响 | 节点数 | | 577 | 1111 | 2250 | DP4 | 0.661 | 0.618 | 0.593 | 简化Dishop法 | 0.583 | 0.583 | 0.583 | (DP4-Bishop)/Bishop | 0.134 | 0.060 | 0.017 | 注：H=20m；β=45°；φ=45°；c=10000Pa。 | |
笔者对采用不同流动法则的算例进行了初步分析，表3的计算结果表明：对同一边坡，不论采用关联流动法则还是非关联流动法则，计算结果相差不大。这是因为它们只与坡体的体积 变形有关，而在边坡稳定分析中，坡体常常为无约束天然坡体，体积变形对坡体稳定影响并不明显。然而，从破坏时位移大小及塑性区的分布来看，还是会有一些差异，有时并不能简单的忽略这种差异[8]。文中所有的算例均取ψ=0，即满足非关联流动法则，算例结果显示出较好的精度。
4 有限元法引入的误差
如前所述，本构模型的选择合理与否会对有限元折减系数法的计算精度造成较大影响，除此之外，有限元法本身也是误差的主要来源之一。
4.1 网格的疏密 网格疏密对单元精度的影响甚至大于单元类型的影响，对于精度较低的单元，可通过加密网格来达到较高的精度。表4列出了不同疏密的网格对计算结果的影响，由表4可知，对于折减系数法，有限元网格不能太稀，否则结果将不可用。通过大量算例证实，对于4节点矩形单元，当单元密度达到每10m2不少于3个节点时，计算精度较为理想，如果再增加节点，计算精度应还能提高，但此时耗费的机时也将成倍增长。
4.2 边界范围 边界范围的大小在有限元法中对计算结果的影响比在传统极限平衡法中表现的更为敏感，在极限平衡法中只要所求滑移面在边界之内就不会对计算结果有影响，安全系数只与划分的土条有关，而与土条外的区域无关，有限元法则不然，边界的大小直接影响到应力-应变的分布。
表5 边界条件对折减系数的影响 |
|
| 相对边距比 |
| 0 | 0.5 | 1.0 | 1.5 | 2.0 | 2.5 | 3.0 |
L/H
R/H
B/H | 1.129
1.097
1.106 | 1.124
1.078
1.117 | 1.124
1.121
1.120 | 1.120
1.122
1.131 | 1.122
1.122
1.124 | 1.121
1.120
1.132 | 1.129
1.123
1.131 |
注：表中所有模型均选用DP4屈服准则。L为坡脚到左端边界的距离（左
边距）；R为坡顶到右端边界的距离（右边距）；B为边坡到底端边
界的距离（底边 距）； H为坡高。 |
为了得到能使计算结果趋于稳定的边界范围，分别对左端、右端、底端三条边界范围的取值大小进行了分析(表5、图4)，计算时，令3个边距中的1个变化，其余2个不变。由图4可知：左边界对计算结果的影响最不敏感，不同的取值相差不到1%，底端边界次之，最大相差在1%左右，右端边界对计算精度的影响最大，达到5%。经对比分析得：当坡角到左端边界的距离为坡高的1.5倍，坡顶到右端边界的距离为坡高的2.5倍，且上下边界总高不低于2倍坡高时，计算精度最为理想。
5 可行性研究
5.1 模型的建立及研究方案 某一处于施工阶段的土质边坡，不考虑孔隙水压力的影响，土坡天然容重γ=25kN/m3。 坡体底边界为固定约束，左右边界为水平约束，其它边界为自由端。计算程序采用商业有限元软件Ansys，计算单元采用平面4节点矩形单元，当坡高保持不变时(20m)，节点1111个，对于坡高H与坡度β变化的情况，为保证足够的精度，必须重新划分单元。同时还应在临界滑移面可能区域对单元网格进行局部加密，见图5。有限元计算所需参数共6个，见表6，其中所有算例均取ψ=0。为避免个别特例在计算结果上存在的巧合，文中对影响坡体稳定安全系数的4个主要参数(坡高H、坡角β、粘聚力C、摩擦角Φ)进行了对比分析，计算方案如下：保持4个参数中的3个为常量，只取1个参数为变量，共4组计算方案，计算结果见表2、表7、表8、表9。
| 
|
图4 边界距离与坡高比～折减系数曲线
|
图5 有限元单元网格划分
|
|
表6 有限元计算参数 Drucker-Prager准则 | 弹性模量E/kPa | 1000 | 帕松比ν | 0.3 | 土体密度γ/(kN/m3) | 25 | 磨擦角ø/° | 变量 | 粘聚力C/Pa | 变量 | | 表7c为变量时的最小安全系数(节点数1111个) | C/kPa | | 0.1 | 20 | 40 | 60 | 90 | DP4 | 0.304 | 0.793 | 1.101 | 1.379 | 1.781 | 简化Bishop法 | 0.254 | 0.752 | 1.036 | 1.302 | 1.685 | (DP4-Bishop)/Bishop | 0.197 | 0.055 | 0.063 | 0.059 | 0.057 | 注：H=20m；β=45°；ø=17°。 | |
|
表8H为变量时的最小安全系数(节点数≥1190个) | H/m | | 10 | 20 | 30 | 40 | 50 | DP4 | 1.733 | 1.128 | 0.923 | 0.820 | 0.735 | 简化Bishop法 | 1.612 | 1.064 | 0.867 | 0.764 | 0.698 | (DP4-Bishop)/Bishop | 0.075 | 0.060 | 0.065 | 0.073 | 0.053 | 注：β=45°；c=42kPa；ø=17°。 | | 表9β为变量时的最小安全系数(节点数≥1210个) | 坡角β/° | | 30 | 35 | 40 | 45 | 50 | DP4 | 1.455 | 1.323 | 1.214 | 1.128 | 1.044 | 简化Bishop法 | 1.398 | 1.269 | 1.156 | 1.064 | 0.987 | (DP4-Bishop)/Bishop | 0.041 | 0.043 | 0.050 | 0.060 | 0.058 | 注：H=20m；c=42kPa；ø=17°。 | |
表2、表7、表8、表9列出了各种情况下的最小安全系数。作为对比，文中还给出了极限平衡法(简化Bishop法)的计算结果，考虑到算例坡体土质均匀，可采用圆弧滑移面，在此情况下简化Bishop法已具有足够的精度[4]。误差定义详见各表。
5.2 结果分析 屈服准则对计算精度的影响很明显。在相同网格密度下，DP4的计算精度明显好于DP1、DP2、DP3(图2)。因为各种方案所显示的规律相同，所以文中只给出了方案1的详细解答(表2)，而其它方案只给出DP4的计算结果。
Φ值的大小对计算精度的影响是明显的。如图2所示，Φ值增大，误差也呈增大的趋势。在文中所给出的4个屈服准则中，DP4精度最高，其与简化Bishop法最接近，平均误差为5.7%。
C、β、H对计算精度的影响不明显，图6分别给出了不同C、β、H值对计算结果的影响，它们对结果的影响约在1%左右。
|
图6C、H、β～折减系数曲线
|
值得注意的是：当Φ、C分别为零时，DP4的误差较大，这是因D-P类本构只适用于摩擦型材料，当ø=0时，计算是不收敛的，在这里笔者代入ø=0.1以近似ø=0；当C为零时，处理方法同Φ。因计算误差较大，所以此时折减系数法将不再适用。
需要指出的是，极限平衡法采用的圆弧滑动面也会对计算精度造成一定的影响，特别是当β很大时，如果采用任意滑动面的Spencer法，则计算精度还能提高。
6 结论
(1)通过4组计算方案共计106个算例的比较分析，表明：折减系数法所得稳定安全系数比简化Bishop法平均高出约5.7%，且误差离散度极小。这有力地说明了将有限元折减系数法用于分析土坡稳定问题是可行的，但必须合理地选用屈服条件并严格地控制有限元法的计算精度。(2)边界范围的大小在有限元法中对计算结果的影响比在传统极限平衡法中表现的更为敏感，当坡角到左端边界的距离为坡高的1.5倍，坡顶到右端边界的距离为坡高的2.5倍，且上下边界总高不低于2倍坡高时，计算精度最为理想。(3)Φ值的大小对有限元折减系数法的计算精度有较大影响，且随Φ值增大误差随之增大，增大幅度因准则类型不同而不同，其中准则DP4计算误差随Φ值增大影响相对较小。C、β、H对计算精度的影响不明显，约在1%左右。(4)有限元法的优点不仅仅在于求出折减系数，如果此法可行，那么对于具有复杂地貌、地质的边坡则可自动求出任意形状的临界滑移面，并能模拟出土坡失稳及施工开挖的自然过程，这是传统极限平衡法无法做到的。因此本文对该法的可行性分析是重要且必要的。目前的工作是个基础，如何在有限元折减系数法中考虑多种土层边坡、孔隙水的影响、非自重外荷以及岩质边坡中存在的大量节理等仍需做深入研究。
参 考 文 献：
[1] Jiang G L， Magnan J P. Stability analysis of embankments: comparison of limit analysis with methods of slices [J]. Geotechnique 1997,47(4):857-872.
[2] Griffiths D V， Lane P A. Slope stability analysis by finite elements [J]. Geotechnique, 1999,49(3):387-403.
[3] 连镇营，韩国城，孔宪京.强度折减有限元法开挖边坡的稳定性[J].岩土工程学报，2001，23(4)：407-411.
[4] 龚晓南.土工计算机分析[M].北京：中国建筑工业出版社，2000.
[5] 徐干成，郑颖人.岩土工程中屈服准则应用的研究[J].岩土工程学报，1 990，(2)：93-99.
[6] Menétrey Ph, Willam K. A triaxial failure criterion for concrete and its generalization [J]. ACI Journal, 1995,92(3):311-318.
[7] 郑颖人，孔亮.塑性力学中的分量理论——广义塑性力学[J].岩土工程学报，2000，22(3)：269-274.
[8] Zienkiewicz O C, Humpheson C, Lewis R W. Associated and non-associated visco-plasticity and plasticity in soil mechanics [J]. Geotechnique, 1975,25(4):671-689.
[9] Lechman J B, Griffiths D V. Analysis of the progression of failu re of earth slopes by finite elements [C]. Int Anal conference of slope stabil ity 2000,250-265.
[10] Wong F S. Uncertainties in FE modeling of slope stability [J]. Computer & Structures, 1984,19：771-791.
[11] Dawson E M, Roth W H, Drescher A. Slope stability analysis by strength reduction [J]. Geotechnique, 1999,49(6):835-840.
[12] Z_SOIL.PC 2001 User manual [Z]. Zace Services Ltd Report 1985-2001.Lausanne:Elmepress International.
收稿日期：2001-11-14
作者简介：张鲁渝(1974-)，男，四川万州人，现于重庆后勤工程学院攻读博士学位，主要从事岩体本构及边坡稳
定分析方面的研究。
