有限差分强度折减法求解边坡稳定性

张 超， 杨春和

(中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071)

边坡稳定分析是经典土力学最早试图解决而至今仍未完全解决的课题，各种稳定分析方法在国内外水平大致相当。目前，研究边坡稳定性的传统方法主要有:极限平衡法，极限分析法，滑移线场法等。这些建立在极限平衡理论基础上的各种稳定性分析方法没有考虑土体内部的应力应变关系，无法分析边坡破坏的发生和发展过程，没有考虑土体与支挡结构的共同作用及其变形协调，在求安全系数时通常需要假定滑裂面形状为折线、圆弧、对数螺旋线等［1］。对于均质土坡，可以通过各种优化方法来搜索危险滑动面。但对于非均质土坡，利用圆弧搜索滑面就不是很合适。而岩质边坡由于实际岩体中含有大量不同构造、产状和特性等不连续结构面，因此传统分析方法在岩质边坡的稳定分析方面就遇到了一定的困难。而有限单元法不但满足力的平衡条件，而且考虑了材料的应力应变关系，使得计算结果更加精确合理。在有限元法中通过强度折减，使系统达到不稳定状态，有限元计算不收敛，此时的折减系数就是安全系数。有限元强度折减法在20世纪70年代提出，在国外80年就有采用，但由于力学概念不十分明确，且计算力学还在起步阶段，缺少严密可靠的大型商用程序，有限元前后处理技术水平较低，采用有限元法对岩土工程建模分析需要很长时间，计算成本高，从而阻止了有限元强度折减法的应用。随着计算机技术的发展，尤其是岩土材料的非线性弹塑性有限元计算技术的发展，出现了许多适合于岩土材料的大型通用有限元软件，其前、后处理的功能越来越强大，为利用有限元法进行边坡稳定分析创造了条件。1999年美国科罗拉多矿业学院的Griffith等人［2］采用有限元强度折减法与传统方法得到的稳定性系数比较接近，再次引起了国内外学者广泛关注。而国内的郑颖人等［3～5］在有限元强度折减法方面做了大量的工作，也取得了很多的研究成果。特别需要提到的是，传统的有限元法采用莫尔－库仑准则，但因莫尔－库仑准则的屈服面为不规则的六角形截面的角锥体表面，存在尖顶和棱角，给数值计算带来困难。为了解决这个问题，徐干成、郑颖人提出的莫尔－库仑等面积圆屈服准则代替莫尔－库仑准则［5］，并导出各准则间的换算关系，由此可将求得的安全系数折算成莫尔－库仑等面积圆屈服准则下的安全系数。显然经过这种近似处理后能够解决计算中的困难，也将有限元强度折减法在分析边坡稳定性方面的应用向前推进了一大步。但是这种莫尔－库仑等面积圆屈服准则对计算结果是否有影响，有多大影响，本文将对该问题进行一定的分析。

近年来，随着美国Itasca咨询公司根据Cundall等人提出的显式有限差分法而开发的有限差分软件FLAC的问世及逐步完善，有限差分方法在岩土工程中获得了广泛的应用。FLAC法是一种分析边坡稳定的新的数值方法。对于大部分边坡，边坡自身的弹性变形、塑性变形或流变变形较大，使变形体在后缘张开，前缘压缩隆起，滑动面也不再是滑动前的形状，在这种情况下，必须考虑用与变形介质相应的分析方法。近年来发展起来的FLAC法是一种比较有效的方法。FLAC法能处理一般的大变形问题，而且能模拟岩体沿某一软弱面产生滑动变形。FLAC法能针对不同的材料特性，使用相应的本构方程来比较真实地反映实际材料的动态行为。用FLAC法分析边坡稳定问题可以考虑岩土体的非线性本构关系，适用于复杂边坡稳定性分析。而且，Itasca咨询公司推出和完善了FLAC/Slope程序，该程序是FLAC程序的微型版本，专门设计用来求解边坡稳定安全系数，采用的是有限差分强度折减法对边坡的稳定性进行计算分析，该程序计算建模方便，易于使用。基于有限差分强度折减法，杨光华［6］等用局部强度折减法进行了边坡稳定性分析，冯磊［7］等采用有限差分强度折减法对黄延高速公路边坡进了稳定性研究。可以看出有限差分法求解边坡稳定性已经开始在一定范围内应用，但是对于其可靠性与适用性却没有相关研究，本文拟以FLAC有限差分商用程序为基础，对有限差分方法进行边坡稳定性分析进行详细介绍，并与传统的有限元强度折减法进行对比分析，探讨有限差分强度折减法的适用性与可靠性。

1 有限差分强度折减法基本原理

FLAC是在较好地吸收了其它数值方法的优点并克服其缺点的基础上而形成的一种新型显式差分程序，该程序的基本原理和算法与离散元法相似，但它应用了节点位移连续条件，可以对连续介质进行大变形分析。FLAC程序可以模拟多种模型的材料，还可以自定义材料本构模型，可以进行边坡、基础、坝体、隧道、地下采场、洞室等问题的模拟分析。FLAC还具有很强的前处理功能和后处理功能。由于该程序特别适合于岩土工程问题分析，推出以后立即受到重视，目前在国际岩土界相当流行。

FLAC用有限差分方法求解边坡稳定性问题的解是耦合应力与位移的全数值解，FLAC在计算边坡稳定性的过程中考虑了边坡的平衡方程和本构方程。只要给出一系列参数，FLAC即可通过计算发现系统是否稳定。在给出不同参数后，FLAC将进行一系列模拟运算，计算出安全系数，该系数表示边坡稳定程度。与极限平衡方法相比，该方法计算速度慢，但对问题的考虑更加全面。近几年由于计算机速度的大大提高，该方法开始成为极限平衡方法的一种可行的代替方法。

有限差分法求得的全数值解(FLAC的解)与极限平衡方法的计算有很大的区别，具体差异见表1。从表1可以看出，有限差分的数值解考虑的问题非常全面，而采用极限平衡方法的时候有很多问题都没有考虑，很多条件也不满足，且极限平衡法只能给出一个安全系数，对于应力和变形则是无能为力。因此可以说在求解安全系数方面，除了计算时间方面有一定的劣势以外，有限差分方法比极限平衡方法在其它各方面均有明显的优势。

表1 FLAC数值方法与极限平衡法的差别

有限差分方法是用给定的参数值进行单个数值模型模拟，即可显示是否破坏。有限差分强度折减系数法的基本原理是将坡体强度参数:粘聚力c和内摩擦角φ值同时除以一个折减系数Ftrial，得到一组新的 Ctrial、φtrial值，然后作为新的参数输入，再进行试算，直到计算出坡体的破坏和不破坏的界限值，对应的Ftrial被称为坡体的最小稳定安全系数，此时坡体达到极限状态，发生剪切破坏，同时可得到坡体的破坏滑动面。有限差分强度折减法的思想和有限元强度折减法的思想完全一样，其安全系数的概念也是相同的，即为强度折减到极限平衡状态条件下的强度折减程度。

强度参数(内聚力和内摩擦角)除以安全系数的“试验值”，如下:式中:C为边坡体材料的内聚力，φ为材料的内摩擦角，Ctrial为内聚力试算值，φtrial为内摩擦角试算值，Ftrial为安全系数的试算值。

通过对一系列的Ctrial和φtrial值进行试算，然后找到其临界状态下的Ftrial值，即为所求的安全系数。如何通过最少数量的模拟，快速得到安全系数的准确值，这个问题是解决有限差分强度折减法的关键问题之一，直接决定有限差分强度折减法的效率和实用性。Dawson等人［8］给出数值计算安全系数的基本依据，根本上讲，它是不断地等分两个分别稳定和破坏的Ftrial，这个过程可以很快收敛，当这两个Ftrial之间的间隔很小(比如小于0．005)时，即停止计算。FLAC/Slope也正是利用这个方法对边坡的稳定性进行强度折减分析，从而可以在尽可能短的时间内找出安全系数。

2 有限差分强度折减法步骤

采用FLAC实现有限差分强度折减法的求解步骤如下:

(1)根据系统达到平衡需要的计算步数来确定“特征反应时间”，称之为NC;

(2)设F=1．0，不断平分该值直至下限(稳定首次出现)，这时的F值称为FS;

(3)不断翻倍增加F直到到达上限(首次出现不稳定)，这时的F值称为FU;

(4)设F=(FS+FU)/2，计算是否稳定，如果稳定，则设FS=F，如果不稳定，则设FS=F，继续进行下一步计算;

(5)如果 FU－FS＜0．005，停止计算，否则回到第(4)步。

边坡的稳定或不稳定由程序中给定的方法确定。例如，在采用隐式矩阵计算的有限元程序中，不稳定条件通常是根据方程的不收敛所决定［6］。在FLAC中，不稳定状态由检测模型的运动能量确定，运动能量的变化则通过不平衡力率来衡量。

在FLAC计算中，检测稳定或是不稳定的步骤:

(1)运行直至NC步，记录不平衡力率RU;

(2)运算过程中，如果 RU降低至0．001以下，认为是稳定并退出计算。

(3)如果 RU－RU(old)/RU＜0．1，认为不稳定并退出计算。

(4)如果总迭代数(即第(1)至第(3))大于6，认为不稳定并退出计算。

(5)返回(1)。

在整个过程中，下列信息是被展示出来的:

(1)上面第(1)步中完成的计算步骤数，以NC的百分数表示;

(2)完成的求解循环数(第(1)至第(3)步);

(3)当前的FS和FU值，以区间的形式给出。

3 有限差分强度折减法可靠性验证

FLAC强度折减法的可靠性到底怎么样，下面就以一个典型的算例［8］进行验证。这个算例是1987年澳大利亚计算机应用协会(ACADS)对澳大利亚所使用的边坡稳定分析程序进行的一次调查中的第一考核题。该考核题的剖面如图1所示(单位:m)，材料性质见表2。本题推荐的裁判答案为Donald给出的，结果为1．00。

图1 ACADS考核题－EX1

表2 考核题EX1材料性质

下面采用FLAC计算该考核题，不同密度的网格计算得到的安全系数见表3，计算得到各种网格的滑动面如图2～图4所示。从计算结果可以看出，用FLAC进行有限差分的强度折减法计算得到的结果和裁判结果非常接近，几乎相等。因此，只要网格划分满足一定的要求，其误差就会相当的小，对于数值分析，这样的误差是完全可以接受的。

表3 FLAC计算的安全系数

图2 计算得到的破坏面(粗糙网格)

以下再看一验证算例，剖面图及计算参数如图5所示，该问题的极限分析安全系数精确解为1．0［8］。在这个例子中，重力加速度为 10 m/s2，坝体材料的体积应变服从关联流动法则。按照以上条件采用FLAC对该算例进行计算，计算结果如下:采用中密网格计算时，安全系数为1．01;采用精细网格时，安全系数为0．99。这个算例再一次证明了有限差分强度折减法和FLAC程序的可靠性。

图3 计算得到的破坏面(中密网格)

图4 计算得到的破坏面(精细网格)

图5 均质土坝剖面图

4 两种强度折减法比较

以图6所示的一个均质边坡为例，分别采用有限元强度折减法和有限差分强度折减法进行计算，并对计算结果进行比较。该边坡基本参数如下:坡高H=20 m ，材料容重γ=25 kN/m3，粘聚力C=42 kPa，内摩擦角φ=17°。采用有限元强度折减法和有限差分强度折减法求坡角α分别为 30°，35°，40°，45°，50°时边坡的安全系数。计算结果见表4。表4中给出的有限元强度折减法的安全系数参考文献［3］的计算结果，该计算结果为精细网格条件下得到的。FLAC给出了三种不同网格密度情况下的安全系数值。

图6 边坡剖面图

表4 安全系数计算结果

从表4的计算结果可以看出，用有限差分强度折减法和有限元强度折减法得到的安全系数比较接近，但是有限差分强度折减方法得到的安全系数要小于有限元强度折减法，究其原因，其差异应该来自于有限元强度折减法中采用莫尔－库仑等面积圆屈服准则代替莫尔－库仑准则，该莫尔－库仑等面积圆屈服准则满足在π平面上各屈服面面积相等，严格来讲这种等效只是一种近似，而并非对计算过程和结果的完全等效。而采用有限差分方法计算岩土工程体的稳定性时，可直接采用莫尔－库仑准则进行计算，不需要对屈服准则进行处理，从而避免了误差。因为有限差分强度折减方法不需要进行屈服条件的转换，因此可以给计算带来一定的方便。从表4还可以知道对于用FLAC进行的强度折减计算，其计算结果受网格划分的精度影响，网格越细，结果越理想，如果网格划分太粗，将会造成较大的误差，故计算时必须考虑适当的网格密度。

除以上差别之外，有限差分强度折减法与传统的边坡稳定性计算方法相比，基本上具有有限元强度折减法所有的优点。如果使有限差分法保持足够的计算精度，那么有限差分法较传统的方法具有如下优点:能够对具有复杂地貌、地质的边坡进行计算;可以考虑土体的非线性弹塑性本构关系，以及变形对应力的影响;能够模拟土坡的失稳过程及其滑移面形状;能够模拟土体与支护的共同作用;求解安全系数时，可以不需要假定滑移面的形状，也无需进行条分。

5 结论

(1)有限差分强度折减法不需要对滑动面形状和位置做假定，通过强度折减使边坡达到极限平衡的临界状态，此时的折减系数就是稳定安全系数，同时可得到边坡破坏时的滑动面。

(2)有限差分强度折减系数法由程序自动求得边坡的危险滑动面以及相应的稳定安全系数。通过算例分析表明了有限差分强度折减法和FLAC的可靠性。只要保证一定的网格密度，则得到的计算结果误差就很小，一般在1%左右，足够满足数值计算和工程分析的要求。

(3)如果使有限差分计算保持足够的计算精度，那么有限差分法较传统的方法具有如下优点:能够对具有复杂地貌、地质的边坡进行计算，包括非均质土坡和岩质边坡;考虑了土体的非线性弹塑性本构关系，以及变形对应力的影响;能够模拟土坡的失稳过程及其滑移面形状;能够模拟土体与支护的共同作用。

(4)与有限元强度折减法相比，有限差分强度折减法具有以下优点:不用转化屈服准则，可以直接采用摩尔－库伦屈服准则;适用于可能产生大变形的边坡。

［1］ 陈祖煜．土质边坡稳定性分析:原理·方法·程序［M］．北京:中国水利水电出版社，2003．

［2］ Griffiths D V，Lane P A．Slope stability analysis by finite elements［J］．Geotechnique，1999，49(3):387-403．

［3］ 赵尚毅，郑颖人，时卫民，等．用有限元强度折减法求边坡稳定安全系数［J］．岩土工程学报，2002，24(3):343-346．

［4］ 郑颖人，赵尚毅，宋雅坤．有限元强度折减法研究进展［J］．后勤工程学院学报，2005，(3):1-6．

［5］ 徐干成，郑颖人．岩土工程中屈服准则应用的研究［J］．岩土工程学报，1990，12(2):93-99．

［6］ 杨光华，钟志辉，张玉成，等．用局部强度折减法进行边坡稳定性分析［J］．岩土力学，2010，31(S2):53-58．

［7］ 冯 磊，韩 飞，梁 栋．基于有限差分强度折减法的黄延高速公路边坡稳定性分析［J］．华北水利水电学院学报，2009，30(3):86-88．

［8］ Dawson E M，Roch W H，Drescher A．Slope stability analysis by strength reduction［J］．Geotechnique，1999，49(6):835-840．