有限元强度折减系数法计算土坡稳定安全系数的精度研究

摘要：有限元强度折减系数法在边坡稳定分析中的应用正逐渐受到人们的重视。本文较为全面地分析了土体屈服准那么的种类、有限元法自身计算精度以及H(坡高)、β(坡角)、C(粘聚力)、Φ(摩擦角)对折减系数法计算精度的影响，并给出了提高计算精度的具体措施。通过对106个算例的比拟分析，说明折减系数法所得稳定平安系数比简化Bishop法平均高出约5.7%，且离散度极小，这不仅验证了文中所提措施的有效性，也说明了将折减系数法用于分析土质边坡稳定问题是可行的。关键词：强度折减系数边坡稳定屈服准那么误差分析自弗伦纽期于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各准那么参数换算编号准那么种类αφkDP1外角点外接D-P圆DP2内角点外接D-P圆DP3内切D-P圆DP4等面积D-P圆注：表中αφ、k是与D-P有关的材料参数。表2不同屈服准那么所得最小平安系数

φ/°

0.110253545DP10.5251.0441.7692.2543.051DP20.5250.9301.3321.5301.887DP30.4540.8481.2791.4991.870DP40.4770.8961.3961.6892.182简化Bishop法0.4940.8461.3161.6232.073(DP1-Bishop)/Bishop0.0630.2340.3440.3550.472(DP2-Bishop)/Bishop0.0630.0990.012-0.080-0.090(DP3-Bishop)/Bishop-0.0810.002-0.028-0.099-0.098(DP4-Bishop)/Bishop-0.0340.0590.0610.0410.053注：H=20mm；β=45°；C=42kPa。算例分析说明(表2、图2)：DP4准那么与简化Bishop法所得稳定平安系数最为接近。对有效算例(Φ≠0)的误差进行统计分析可知，中选用DP4准那么时，误差的平均值为5.7%，且离散度很小(图3)。而DP1的平均误差为29.5%，同时采用DP2、DP3准那么所得计算结果的离散度非常大，均不可用。因此在数值分析中可用DP4准那么代替摩尔-库仑准那么。

图2Φ～折减系数曲线

图3DP4准那么的计算误差3不同流动法那么的影响

有限元计算中，采用关联还是非关联流动法那么，取决于ψ值(剪胀角)：ψ=φ，为关联流动法那么；ψ≠0，为非关联流动法那么。总体说来，采用非关联流动法那么所得破坏荷载比同一类型材料而采用关联流动法那么所得破坏荷载小，如忽略剪胀角(ψ=0)，将会得到较为保守的结果。值得注意的是：当ψ=0时，正好与郑颖人等提出的广义塑性力学理论相符[7]，这时对应的塑性势面与q轴垂直。表3不同流动法那么的影响

φ=10°φ=17°φ=25°非关联0.8711.1051.363关联0.8871.1371.425相对误差0.0180.0290.045

β=45°；C=40kPa；H=20m；DP4准那么。表4网格疏密对计算结果的影响

节点数

57711112250DP40.6610.6180.593简化Dishop法0.5830.5830.583(DP4-Bishop)/Bishop0.1340.0600.017

注：H=20m；β=45°；φ=45°；c=10000Pa。笔者对采用不同流动法那么的算例进行了初步分析，表3的计算结果说明：对同一边坡，不管采用关联流动法那么还是非关联流动法那么，计算结果相差不大。这是因为它们只与坡体的体积变形有关，而在边坡稳定分析中，坡体常常为无约束天然坡体，体积变形对坡体稳定影响并不明显。然而，从破坏时位移大小及塑性区的分布来看，还是会有一些差异，有时并不能简单的忽略这种差异[8]。文中所有的算例均取ψ=0，即满足非关联流动法那么，算例结果显示出较好的精度。

4有限元法引入的误差

如前所述，本构模型的选择合理与否会对有限元折减系数法的计算精度造成较大影响，除此之外，有限元法本身也是误差的主要来源之一。

4.1网格的疏密网格疏密对单元精度的影响甚至大于单元类型的影响，对于精度较低的单元，可通过加密网格来到达较高的精度。表4列出了不同疏密的网格对计算结果的影响，由表4可知，对于折减系数法，有限元网格不能太稀，否那么结果将不可用。通过大量算例证实，对于4节点矩形单元，当单元密度到达每10m2不少于3个节点时，计算精度较为理想，如果再增加节点，计算精度应还能提高，但此时消耗的机时也将成倍增长。

4.2边界范围边界范围的大小在有限元法中对计算结果的影响比在传统极限平衡法中表现的更为敏感，在极限平衡法中只要所求滑移面在边界之内就不会对计算结果有影响，平安系数只与划分的土条有关，而与土条外的区域无关，有限元法那么不然，边界的大小直接影响到应力-应变的分布。表5边界条件对折减系数的影响

相对边距比

00.51.01.52.02.53.0L/H

R/H

B/H1.129

1.097

1.1061.124

1.078

1.1171.124

1.121

1.1201.120

1.122

1.1311.122

1.122

1.1241.121

1.120

1.1321.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/kPa1000帕松比ν0.3土体密度γ/(kN/m3)25磨擦角ø/°变量粘聚力C/Pa变量表7c为变量时的最小平安系数(节点数1111个)

C/kPa

0.120406090DP40.3040.7931.1011.3791.781简化Bishop法0.2540.7521.0361.3021.685(DP4-Bishop)/Bishop0.1970.0550.0630.0590.057

注：H=20m；β=45°；ø=17°。

表8H为变量时的最小平安系数(节点数≥1190个)

H/m

1020304050DP41.7331.1280.9230.8200.735简化Bishop法1.6121.0640.8670.7640.698(DP4-Bishop)/Bishop0.0750.0600.0650.0730.053

注：β=45°；c=42kPa；ø=17°。表9β为变量时的最小平安系数(节点数≥1210个)

坡角β/°

3035404550DP41.4551.3231.2141.1281.044简化Bishop法1.3981.2691.1561.0640.987(DP4-Bishop)/Bishop0.0410.0430.0500.0600.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]JiangGL，MagnanJP.Stabilityanalysisofembankments:comparisonoflimitanalysiswithmethodsofslices[J].Geotechnique1997,47(4):857-872.

[2]GriffithsDV，LanePA.Slopestabilityanalysisbyfiniteelements[J].Geotechnique,1999,49(3):387-403.

[3]连镇营，韩国城，孔宪京.强度折减有限元法开挖边坡的稳定性[J].岩土工程学报，2001，23(4)：407-411.

[4]龚晓南.土工计算机分析[M].北京：中国建筑工业出版社，2000.

[5]徐干成，郑颖人.岩土工程中屈服准那么应用的研究[J].岩土工程学报，1990，(2)：93-99.

[6]MenétreyPh,WillamK.Atriaxialfailurecriterionforconcreteanditsge