有限元强度折减法计算边坡稳定的对比分析

程灿宇，罗富荣，戚承志，王 霆

（1. 北京建筑工程学院 工程结构与新材料北京市高校工程研究中心，北京 100044；2. 北京市轨道交通管理公司，北京 100037）

1 引 言

我国是一个地质灾害多发的国家，据不完全统计，2010年仅四川省就因地质灾害造成直接经济损失90多亿元[1]，而且滑坡泥石流等地质灾害每年都给人民群众的生命和财产安全带来巨大的损失。我国西高东低的地势导致我国西部地区地质灾害问题更加严重，伴随着西部大开发的推进和大型工程项目的快速上马，在我国西南和黄土高原地区以及西北部新疆地区地质灾害问题变得更加严重，涌现出一系列亟待解决的边坡稳定问题。因此，在我国边坡稳定分析现阶段仍是一个很大的课题。

在边坡稳定分析领域一种常用的方法是有限元强度折减法，早在1975年该方法就被Zienkiewice[2]用来求解边坡稳定问题，只是当时由于需要花费大量的机时而在具体应用领域受到限制。Wong[3]给出了用有限元方法分析边坡稳定时误差产生的原因。现在，随着计算机硬件技术和有限元软件技术的飞速发展，运用有限元强度折减法分析边坡稳定已经成为新的趋势[4]。Ugai[5]、Matsui和 San[6]、Ugai和Leshchinsky[7]、Griffiths 和 Lane[8－9]、Manzari 和Nour[10]等都对此做了进一步的研究，研究结果均表明：有限元强度折减法能得到比较稳定的与极限平衡法接近的安全系数和临界滑动面。应用于岩土工程的分析软件也都可以采用强度折减法计算边坡稳定问题，但目前还没人研究针对具体的工况条件各软件计算边坡稳定时的偏差以及各软件计算特点问题。

本文研究了FLAC3D、MIDAS/GTS和ANSYS在分析边坡稳定时各自的特点；配合D-P屈服准则和M-C屈服准则，计算了软黏土、弱膨胀土、硬黏土3种工况下不同坡度边坡的安全系数；研究了3种工况条件下M-C屈服准则和D-P屈服准则计算结果的偏差；分析了3种工况下滑动面的规律以及各软件得到的滑动面特征；总结了软件计算时的注意事项。本文的研究可以为进行边坡稳定分析的研究人员提供参考。

2 有限元强度折减法计算原理

Duncan[11]指出，边坡安全系数可以定义为使边坡刚好达到临界破坏状态时对土的剪切强度进行折减的程度。通过逐步减小抗剪强度指标，将黏聚力c、内摩擦角φ值同时除以折减系数K，从而得到一组新的强度指标c′、φ′，反复计算直至边坡达到临界破坏状态，此时的强度折减系数即为该边坡的安全系数。

用这种方法分析边坡稳定性时，在弹塑性有限元数值分析中，先将折减系数K取的小一些，以保证开始时是一个近乎弹性问题，然后不断增加折减系数K，反复分析边坡直至K增加至某一值时边坡达到临界状态。有限元强度折减法的优点是安全系数可以直接求出，滑裂面的形状和位置都不需要事先求出，而且从计算过程中可以看到土坡逐步破坏的过程。

3 强度准则及失稳判据

岩土材料具有复杂的本构特性，而边坡的稳定性分析主要关心的是力和强度问题，而不是位移和变形问题，因此，可在有限元强度折减法中采用理想弹塑性本构模型[12]。

对于岩石、混凝土和土壤等颗粒状材料其受压屈服强度远大于受拉屈服强度，且材料受剪时颗粒会膨胀，常用的Mises屈服条件不适用于这类材料，岩土工程中常采用M-C屈服准则和D-P准则[13]。其表达式如下：

式（4）中的α、k与（3）中的常数c、φ关系为

此时D-P准则屈服面在π平面上的投影为M-C屈服准则在π平面上投影的六边形的外角点外接圆[14]。

边坡失稳判据目前有3种比较主流的观点：①有限元数值迭代不收敛判据[15－17]；②特征部位位移突变判据[18－20]；③等效塑性应变区贯通判据[21－22]。前人大量的研究成果表明，这3种有限元强度折减法边坡失稳判据具有统一性[23]，因此，本文根据实际情况综合考虑判据1和判据3作为本文计算中判定边坡失稳的依据。

4 不同软件计算结果的对比分析

FLAC3D、MIDAS/GTS和ANSYS这3款软件中除了ANSYS都既提供了D-P模型又提供了M-C模型，所以在对比不同软件的计算结果的同时，本文同时对比了不同软件不同模型下计算结果的差异。严格地说，这4款软件中FLAC3D应用的是有限差分程序，不是严格意义上的有限元程序，在FLAC3D中对控制偏微分方程进行的是直接逼近来推导方程[24]，这与有限元法不同。因此，通过本文对FLAC3D计算结果和其他软件计算结果的对比，在一定程度上也可以看到在边坡分析中有限元法和有限差分法的差别。

本文计算中选用3种材料参数，分别为：工况1、软黏土；工况2、硬黏土；工况3、弱膨胀土[25]。材料参数如表1所示。

表1 材料参数Table 1 Material parameters

4.1 算例

针对研究的问题：不同软件计算结果的对比分析，本算例采用在许多论文中都得到印证的一个算例[8,23－25]作为对比的参照：坡高为20 m，边坡倾角为 45°。

4.1.1 M-C屈服准则下计算结果的对比

3种工况下利用FLAC3D和MIDAS计算，利用其中内置的强度折减计算模块调整相同的收敛参数得到的安全系数如表2所示。

表2 不同工况下求得的稳定安全系数Table 2 Safety factors by different conditions

从表2可以看出，3种工况下MIDAS的计算结果都要比FLAC3D的计算结果大，其中硬黏土工况下二者偏差最大，达到10.81% 。

为了得到进一步的规律，调整边坡坡度利用MIDAS和FLAC3D分别进行3种工况下边坡安全系数的计算，计算发现MIDAS和FLAC3D的计算结果差异不仅与工况有关而且还与边坡的坡角有关，计算结果如图1所示。

从图1中可以看出，当边坡土体为软黏土时，FLAC3D和MIDAS的计算结果吻合得比较好，经计算结果最大偏差不超过 3%，而且随着坡度的增加二者计算结果的差别在 45°坡时趋于最小。与此相对，当边坡土体为硬黏土或弱膨胀土时，FLAC3D和MIDAS的计算结果偏差相对较大，最大偏差发生在45°的硬黏土坡，最大偏差为10.81%。当坡度较小时FLAC3D的计算结果要明显大于MIDAS的计算结果，当坡度达到 30°时二者计算结果趋于吻合，且MIDAS的计算结果略大于FLAC3D的计算结果。而当坡度继续增加时二者的计算结果偏差也在增大，而且这时MIDAS的计算结果大于FLAC3D的计算结果，与软黏土工况不同的是：在工况2、3下边坡坡度达到 45°时二者偏差达到最大，随着坡度的继续增加，相对偏差趋于减小。

图1 3种工况下MIDAS和FLAC计算的安全系数对比Fig.1 Comparison of safety factors for three conditions calculated by MIDAS and FLAC

4.1.2 D-P屈服准则下计算结果的对比

选择同样的3种工况在ANSYS和MIDAS中选择D-P准则进行手动折减，分别计算15°、30°、45°、60°时边坡的稳定安全系数。手动折减时设置和上面M-C准则下相同的收敛容差值。得到典型的45°坡时的安全系数如图2所示。

图2 45°坡不同工况下安全系数对比（D-P准则）Fig.2 Safety factors for 45 degrees slope under different conditions (D-P criterion)

从图2中可以看出当坡体材料为软黏土和硬黏土时，MIDAS计算的安全系数比ANSYS计算得到的结果大，最大相差5.73%，边坡土体为弱膨胀土时ANSYS计算得到的安全系数比MIDAS计算得到的结果大1.69%。可见45°坡时MIDAS和ANSYS采用D-P准则计算时得到的安全系数具有较高的吻合度。

表3 改变坡角时3种工况下边坡的安全系数Table3 Slope safety factors under different slope angles and conditions

在研究坡角对软件计算偏差的影响时，为了减少采用塑性区贯通作为判别准则产生的人为误差，计算中设定相同的收敛容差后统一以计算不收敛作为边坡破坏的依据。从表3中可以得到比较一致的结论：采用D-P模型时，MIDAS计算的安全系数大于ANSYS计算的安全系数。对硬黏土进行分析时偏差最明显，最大相差14.58%，最小相差1.83%，且呈现出偏差随边坡坡度的增加而增大的趋势。在对软黏土边坡进行分析时也呈现出这种趋势。但对弱膨胀土边坡进行分析时，偏差随边坡坡度变化的规律没有硬黏土时那么明显，大约都在5%左右。

4.1.3 不同软件在D-P、M-C准则下计算结果的对比

在MIDAS和ANSYS中调用的D-P屈服准则的屈服面在π平面上的投影为 M-C屈服准则的屈服面在π平面上投影的外接圆（见图3），因此，从理论上说采用 M-C屈服准则得到的边坡稳定安全系数小于采用D-P准则得到的边坡稳定安全系数，本例通过计算印证了这一推导，并且综合各软件计算结果得到了针对具体问题时，两个准则计算结果的偏差大小（见图4）。

图3 各屈服准则在π 平面的投影Fig.3 Yield surfaces on the π plane

图4 3种工况下4种计算方式的结果对比Fig.4 Comparison of results by using four calculation methods for three conditions

从图中可以看到利用MIDAS（D-P）的计算结果和MIDAS（M-C）的计算结果有较好的一致性，同时ANSYS（D-P）的计算结果和FLAC（M-C）的计算结果有较好的一致性。因此，我们进行如下的偏差统计，如表4所示。

表4 4种计算模式的相对偏差统计Table 4 Statistics of relative deviation of four calculated modes

从表4可以看出，D-P准则和M-C准则同时用于软黏土工况和弱膨胀土工况时计算结果的偏差一致且相对较小，而当D-P准则和M-C准则分别被用来计算硬黏土边坡时，计算结果的偏差明显变大。由于D-P准则计算得到的安全系数偏大，因此，当D-P准则被用于硬黏土边坡分析时应该特别小心。

4.2 计算分析

4.2.1 关于折减系数的选取

采用手动折减安全系数进行边坡稳定分析时，最初一定要选取较小的折减系数。比如，采用ANSYS进行分析时，当边坡坡角较小时，计算中先出现塑性区贯通，然后出现计算不收敛，而随着边坡坡度的增大，计算不收敛和塑性区贯通可能会同时发生，甚至计算不收敛时塑性区仍未贯通，而且在计算不收敛的情况下会出现塑性区贯通后又不贯通的情况。图 5～7反映了 ANSYS（D-P）分析45°软黏土边坡时塑性区随安全系数变化的关系。

图5 K=1.43（计算不收敛，塑形区恰好贯通）Fig.5 K =1.43(without convergence of calculation, the plastic zone just run through)

图6 K =1.45（计算不收敛，塑性区未贯通）Fig.6 K =1.45(without convergence of calculation,and without breakthrough of plastic zone)

图7 K =1.46（计算不收敛，塑性区贯通）Fig.7 K =1.46 (without convergence of calculation, but with breakthrough of plastic zone)

4.2.2 工况对滑动面的影响

上面算例定义了3种典型工况：软黏土、硬黏土、弱膨胀土。以45°坡为例，3种工况下ANSYS（D-P）计算的滑动面如图8所示。

图8 45°坡角的各边坡滑动面Fig.8 Sliding surfaces of slopes with 45 degrees

从图中可以看出，硬黏土的滑动面最浅，软黏土和弱膨胀土的滑动面都较深，而且模拟结果显示边坡破坏时软黏土的位移量最大，弱膨胀土次之，硬黏土最小。改变坡角可以得到同样的规律。

不同软件得到的关于工况影响的规律是相同的，但是不同软件计算得到的滑动面的位置却是有差异的，以30°软黏土边坡为例，ANSYS和MIDAS计算的安全系数分别为1.70和1.78。而且ANSYS计算得到的滑动面较深切，变换坡度和工况后仍能发现这一规律。

图9 30°软黏土坡的屈服面Fig.9 The yield faces of 30 degrees slope

5 结 论

（1）采用 ANSYS(D-P)、MIDAS(D-P)、MIDAS(M-C)、FLAC(M-C)共 4种分析方法，在软黏土、硬黏土和弱膨胀土 3种工况下的计算表明：FLAC(M-C)最保守，MIDAS(D-P)的结果最大。

（2）采用M-C屈服准则时当边坡坡度较小（15°左右）时MIDAS/GTS的计算结果小于FLAC，计算中显示30°及以上边坡FLAC得到的安全系数均小于由MIDAS/GTS得到的。

（3）边坡的滑动面位置与土体材料和坡角大小有关，软黏土边坡的滑动面较深，硬黏土滑动面较浅。

（4）计算软件的选取也会影响对滑动面深浅的判断，MIDAS计算得到的滑动面比ANSYS计算得到的滑动面浅。

（5）对软黏土边坡分析时，D-P准则的计算结果和M-C准则的计算结果偏差较小（12%～25%）；对硬黏土边坡进行分析时，D-P准则的计算结果与M-C准则的计算结果偏差较大（12%～42%）；弱膨胀土工况的情况和软黏土工况类似。工况条件对分析结果的影响规律通过不同软件分析印证是一致的。

综上所述，当边坡土体为软黏土或弱膨胀土时采用M-C屈服准则和D-P屈服准则分析均可；当边坡土体为硬黏土时D-P准则计算结果和M-C准则计算结果偏大较大，从安全考虑本文推荐硬黏土边坡尽量采用M-C屈服准则进行分析；当边坡支护时，如果参照 MIDAS的分析结果，设计者应考虑MIDAS计算得到的滑动面较其他软件得到的浅。

[1]四川省国土资源厅. 四川省地质环境公报[J]. 资源与人居环境, 2011, 8: 53－55.

[2]ZIENKIEWICZ O C, HUMPHESON C, LEWIS R W.Associated and non-associated visco-plasticity in soil mechanics[J]. Geotechnique, 1975, 25(4): 671－689.

[3]WONG F S. Uncertainties in FE modeling of slope stability[J]. Computers & Structures, 1984, 19: 777－791.

[4]连镇营, 韩国城, 孔宪京. 强度折减有限元法研究开挖边坡的稳定性[J]. 岩土工程学报, 2001, 23(4): 407－411.LIAN Zhen-ying, HAN Guo-cheng, KONG Xian-jing.Stability analysis of excavation by strength reduction FEM[J]. Chinese Journal of Geotechnical Engineering,2001, 23(4): 407－411.

[5]UGAI K. A method of calculation of total factor of safety of slopes by elasto-plastic FEM[J]. Soils and Foundations, 1989, 29(2): 190－195.

[6]MATSUI T, SAN K C. Finite element slope stability analysis by shear strength reduction technique[J]. Soils and Foundations, 1992, 32(1): 59－70.

[7]UGAI K, LESHCHINSHY D. Three-dimensional limit equilibrium and finite element analysis: A comparison of results[J]. Soils and Foundations, 1995, 35(4): 1－7.

[8]GRIFFITHS D V, LANE P A. Slope stability analysis by finite elements[J]. Geotechnique, 1999, 49(3): 387－403.

[9]GRIFFITHS D V, LANE P A. Slope stability analysis by finite elements[J]. Geotechnique, 1999, 49(3): 387－403.

[10]MANZARI M T, NOUR M A. Significance of soil dilatancy in slope stability analysis[J]. Journal of Geotechnical and Geoenvironmental Engineering,ASCE, 2000, 126(1): 75－80.

[11]DUNCAN J M. State of the art: Limit equilibrium and finite-element analysis of slope[J]. Journal of Geotechnical Engineering, ASCE, 1996, 122(7): 577－596.

[12]水利部水利水电规划设计总院, 黄河勘测规划设计有限公司. SL386－2007水利水电工程边坡设计规范实施指南[S]. 北京: 中国水利水电出版社, 2009.

[13]谭长建, 张娟, 董城. ANSYS高级工程应用实例分析与二次开发[M]. 北京: 电子工业出版社, 2006.

[14]余天庆, 王勋文, 刘再华. 弹性与塑性力学[M]. 北京:中国建筑工业出版社, 2009.

[15]DAWSON E M, ROTH W H, DRESCHER A. Slope stability analysis by strength reduction[J]. Geotechnique,1999, 49(6): 835－840.

[16]赵尚毅, 郑颖人, 孙玉芳. 有限元强度折减法中边坡失稳的判据探讨[J]. 岩土力学, 2005, 26(2): 332－336.ZHAO Shang-yi, ZHENG Ying-ren, SUN Yu-fang. Study of slope failure criterion in strength reduction finite element method[J]. Rock and Soil Mechanics, 2005,26(2): 332－336.

[17]王栋, 年廷凯, 陈煜淼. 边坡稳定有限元分析中的三个问题[J]. 岩土力学, 2007, 28(11): 2310－2313, 2318.WANG Dong, NIAN Ting-kai, CHEN Yu-miao. Three problems in slope stability analyses with finite element method[J]. Rock and Soil Mechanics, 2007, 28(11):2310－2313, 2318.

[18]李红, 宫必宁, 陈琰. 有限元强度折减法边坡失稳判据[J]. 水利与建筑工程学报, 2007, 5(1): 79－82.LI Hong, GONG Bi-ning, CHEN Yan. Study of criteria for evaluating stability of slope with FEM based on shear strength reduction methods[J]. Journal of Water Resources and Architectural Engineering, 2007, 5(1):79－82.

[19]MANZARI M T, NOUR M A. Significance of soil dilatancy in slope stability analysis[J]. Journal of Geotechnical and Geoenvironmental Engineering,ACSE, 2000, 123(1): 75－80.

[20]宋二祥. 土工结构安全系数的有限元计算[J]. 岩土工程学报, 1997, 19(2): 1－7.SONG Er-xiang. Finite element analysis of safety factor for soil structures[J]. Chinese Journal of Geotechnical Engineering, 1997, 19(2): 1－7.

[21]栾茂田, 武亚军, 年廷凯. 强度折减有限元中边坡失稳的塑形区判据及其应用[J]. 防灾减灾学报, 2003, 23(3):1－8.LUAN Mao-tian, WU Ya-jun, NIAN Ting-kai. A criterion for evaluating slope stability based on development of plastic zone by shear strength reduction FEM[J]. Journal of Disaster Prevention and Mitigation Engineering,2003, 23(3): 1－8.

[22]周翠英, 刘祚秋, 董立国. 边坡变形破坏过程的大变形有限元分析[J]. 岩土力学, 2003, 24(4): 644－647, 652.ZHOU Cui-ying, LIU Zuo-qiu, DONG Li-guo. Large deformation FEM analysis of slopes failure[J]. Rock and Soil Mechanics, 2003, 24(4): 644－647, 652.

[23]裴利剑, 屈本宁, 钱闪光. 有限元强度折减法边坡失稳判据的统一性[J]. 岩土力学, 2010, 31(10): 3337－3341.PEI Li-jian, QU Ben-ning, QIAN Shan-guang.Uniformity of slope instability criteria of strength reduction with FEM[J]. Rock and Soil Mechanics, 2010,31(10): 3337－3341.

[24]刘波, 韩彦辉. FLAC原理、实例与应用指南[M]. 北京:人民交通出版社, 2005.

[25]闫宇, 赵旻. 渠道膨胀土边坡土体强度参数取值原则初探[J]. 水利水电工程设计, 2008, 27(2): 37－39, 46.YAN Yu, ZHAO Min. The preliminary discussion on the principle of taking value for strength parameters of canal expansive soil side slope[J]. Design of Water Resources and Hydroelectric Engineering, 2008, 27(2): 37－39,46.