基于数值强度折减法的边坡失稳判定方法

彭晓钢 李 嘉 李有志 江 建

深圳市天健（集团）股份有限公司 广东 深圳 518034

准确判断边坡的稳定性对路基、边坡、基坑等工程非常重要，但是常规有限元分析方法难以直接评价边坡稳定性，Zienkiewicz等[1]提出了强度折减法，该方法定义土体所能提供的最大抗剪强度和土体内实际剪应力之比为抗剪强度折减系数，计算时先取较小的折减系数Ftrial，保证开始时边坡处于稳定状态，然后逐渐增大折减系数Ftrial，直至Ftrial增大到某一值时边坡处于极限状态，即荷载引起的土体中实际剪应力恰好等于按照实际强度指标折减后的抗剪强度指标。Duncan[2]将边坡安全系数定义为使边坡刚好达到临界破坏状态时对土的剪切强度进行折减的程度，这一概念也就是传统意义上的边坡稳定安全系数。

强度折减法提出伊始受到计算条件的限制普及难度较大，Griffiths等[3]利用Fortran语言编写有限元程序首次实现了边坡强度折减法分析。随着计算机硬件和数值软件的发展，强度折减法愈加受到重视，很多岩土工程商业软件中已经可以较为方便地应用强度折减法，如美国Itasca公司开发的FLAC、荷兰Technical University of Delft开发的Plaxis、瑞士Zace Services公司开发的Z-SOIL.PC、韩国POSCO公司开发的Midas/GTS模块等。

1 强度折减法在Abaqus软件中的实现

Abaqus是当前应用最为广泛的非线性有限元软件之一，该软件没有提供强度折减法，但可以通过场变量的设置简便地实现[4-6]。

强度折减法的基本原理实质就是强度参数c（黏聚力）、φ（内摩擦角）的逐渐降低导致土体强度不足以承担当前应力，在有限元中即反映为强度参数c、φ的逐渐降低导致某单元应力超过了屈服面产生塑性变形，超出屈服面的应力将通过塑性变形传递到周围土体单元，当应力超过屈服面的单元形成贯通的曲面时土体将失去稳定。

在Abaqus软件中，可以通过设置强度参数随场变量变化来实现强度折减，场变量没有具体的物理意义，使计算分析更少受到人为假设的影响。在Abaqus软件中，场变量Ffield随时间增量步线性变化，材料参数和场变量之间有如下关系：

其中，c、φ是土体所能提供的抗剪强度指标，c'、φ'是某时间增量步下强度折减后土体的抗剪强度指标。

对比式（1）和式（2），土体到达极限状态时的抗剪强度折减系数Ftrial和场变量Ffield有如下关系：

在Abaqus软件计算过程中，场变量随着增量步时间t而变化，通过设置材料参数随场变量变化，实现材料强度折减过程。在Abaqus软件中可实现增量步时间t的自动变化，进而实现强度的自动折减，并根据设定的破坏判据找到所对应的增量步时间，然后可得安全系数，并可得到相应的滑动面。

2 Abaqus软件中强度折减法的失稳判定

正确判断失稳状态影响到强度折减有限元法得到的折减系数，进而影响安全系数的计算结果，强度折减有限元法的失稳判据主要有3类[7-8]，即：数值迭代不收敛判据、特征部位位移突变判据、广义塑性应变或等效塑性应变贯通判据。

2.1 分析模型

为便于讨论比较，以文献[3]中的案例为对象，用Abaqus软件进行强度折减有限元分析，对上述失稳判据进行对比。

Example1为均质边坡，坡高H＝10 m，坡角β＝30°，土体重度γ＝20 kN/m3，黏聚力c＝5 kPa，内摩擦角φ＝30°。

建立边坡有限元模型，几何尺寸和文献[3]中一致，坡顶和坡脚到模型边界均取为12 m，坡脚以下土体厚度取10 m，假定弹性模量E＝100 MPa，为消除因网格划分造成的计算结果误差，采用和文献[3]一致的网格划分，采用缩减积分平面应变单元CPE4R。约束模型左右方向边界的水平方向位移，底端边界的水平方向和竖直方向位移。有限元模型如图1所示。

图1 土坡强度折减有限元模型

对于实际计算边坡安全系数的工程问题，强度折减系数Ftrial≥1，即根据土体实测强度指标进行强度折减，计算得到实际边坡的安全系数。但对于有限元计算，不考虑实际意义，强度折减系数可以小于1，即折减系数起到扩大强度指标的作用，折减后强度指标反而更趋于安全。

计算分2步进行：第1步对边坡施加竖直向下的重力，为使重力加载步顺利完成，重力加载后多数区域处于弹性状态，并使强度折减计算平稳过渡，采用Ftrial＝0.5时的参数，即黏聚力c＝10.0 kPa，内摩擦角φ＝49.0°；第2步对强度参数进行折减，设置最终状态为Ftrial＝2.0，此时对应的土体强度参数为黏聚力c＝2.5 kPa，内摩擦角φ＝16.7°。

经过计算，可得不同泊松比v时的塑性区分布，如图2所示。

图2 不同泊松比v时的塑性区分布

由图2可知，泊松比取值对塑性区的分布有较大影响，泊松比取较小值时比取较大值时的塑性区范围大，泊松比取较大值，直至计算因不收敛而终止时的塑性区范围依然很小，但当泊松比取更小值时甚至会产生塑性区贯通的现象，因此采用广义塑性应变或等效塑性应变贯通判据不能准确判断土坡失稳。

有限元最后的总位移等值线包含了第1步中施加重力荷载引起的位移，无法判断滑动面的位置，通过增量步之间的位移可反映边坡失稳趋势，计算采用最后2个增量步之间的位移进行判断，计算得到最后一步的位移增量等值线云图如图3所示，根据图3可判断滑动面的位置，和图4所示的文献[3]中滑动面形状对比也极为接近。

2.2 特征部位位移突变对边坡失稳判定的影响

分别选取坡顶点和坡脚点作为特征点，将计算得到的强度折减系数Ftrial同特征点竖向位移U之间的关系绘制成曲线，如图5所示。

图3 最后一步的位移增量等值线云图

图4 文献[3]中的滑动面形状

图5 坡顶和坡脚特征点Ftrial-U关系曲线

由图5可知，随着强度折减系数增大，土体强度指标降低，直至计算因不收敛而终止，坡脚和坡顶特征点竖向位移均呈现较小泊松比时位移较大的规律，符合弹性理论。其中坡脚特征点竖向位移没有明显突变，较大泊松比时竖向位移还有略微减小的趋势；坡顶特征点位移突变较明显，然而位移突变点的强度折减系数差异较大，选取曲线水平段结束点作为突变点，强度折减系数Ftrial如表1所示。

表1 坡顶特征点位移突变点强度折减系数Ftrial

由表1可知，塑性区的分布会影响到坡顶特征点位移突变点的强度折减系数Ftrial，随泊松比增大，塑性区范围减小，由此判断的强度折减系数Ftrial呈增大趋势，且该方法需根据经验判定位移突变点，没有明确的指标作为判别依据，有一定的经验误差。

2.3 数值迭代不收敛对边坡失稳判定的影响

将数值计算不收敛导致计算终止时的强度折减系数Ftrial汇总如表2所示。

表2 不收敛终止时的强度折减系数Ftrial

由表2可知，随泊松比的增大，虽然塑性区的分布发生了相应的变化，但计算因不收敛而终止时的强度折减系数Ftrial的差异很小，而且该方法无需根据经验判定位移的突变点。

有限元中引起计算不收敛的因素有很多，有些学者质疑以数值计算不收敛作为边坡失稳破坏依据是否存在其他人为导致的结果偏差，可通过对相同模型进行不同的设置并进行对比来说明这一问题。

该问题讨论前首先要进行假设，即计算采用的有限元模型是正确的，如果模型建得不对或采用的有限元程序本身有缺陷，由此而引起的有限元数值计算不收敛本身就是毫无意义的。

上述模型计算采用的计算步长为1，初始增量步为0.1，最小增量步为1×10－5，其中最小增量步是判断收敛与否的关键，将最小增量步进一步减小，分别设置为1×10－6和1×10－7，将计算结果进行对比如表3所示。

表3 收敛标准对计算结果的影响

Abaqus软件中默认的最小增量步为1×10－5，在将最小增量步进一步调小至1×10－6和1×10－7后，计算因不收敛而终止时的强度折减系数变化不大，说明默认的增量步已可以保证足够的计算精度，计算结果也表明有限元收敛条件的设置对结果影响不大，因此可以把有限元计算是否收敛作为边坡失稳的判据。

3 结语

利用Abaqus软件进行有限元强度折减计算的程序可靠，经过对数值迭代不收敛判据、特征部位位移突变判据、广义塑性应变或等效塑性应变贯通3类判据的计算分析，得到如下结论：

1）泊松比取值对塑性区的分布有较大影响，泊松比取较小值时比取较大值时的塑性区范围大，泊松比取较大值，直至计算因不收敛而终止时的塑性区范围依然很小，但当泊松比取更小值时甚至会产生塑性区贯通的现象，因此采用广义塑性应变或等效塑性应变贯通判据不能准确判断土坡失稳。

2）塑性区的分布会影响到坡顶特征点位移突变点的强度折减系数Ftrial，随泊松比增大，塑性区范围减小，由此判断的强度折减系数Ftrial呈增大趋势，且该方法需根据经验判定位移突变点，没有明确的指标作为判别依据，有一定的经验误差。

3）数值计算是否收敛与边坡是否失稳存在对应关系，计算因不收敛而终止时的强度折减系数变化不大，说明默认的增量步已可以保证足够的计算精度，计算结果也表明有限元收敛条件的设置对结果影响不大，建议利用有限元计算收敛作为失稳判定的依据。