0 引言

边坡稳定性分析是岩土工程领域一个极为重要的课题和任务。目前岩土工程设计领域主要以简化的二维模型进行边坡稳定性计算^[1-4]，但是随着勘测水平和理论方法的提高以及三维工程设计的推广，三维边坡稳定性计算逐渐普及^[5-8]。三维边坡稳定性计算可以整合三维地形和岩土体空间分布，可以为工程设计提供更为全面而真实的岩土工程信息^[9-11]。三维边坡稳定性计算方法主要分为极限平衡法、极限分析法和数值分析方法，其中以强度折减法(SRM)为理论依据的数值分析方法(如有限元法、有限差分法)最为常用^{[9, 12-20]}。虽然三维边坡稳定性计算具备诸多优势，但是相对于传统的二维计算方法，计算耗时是其最大的不足，尤其对于复杂边坡，三维边坡稳定性计算往往耗时几个小时，甚至更长，因此一定程度上限制了其在工程实际中的应用。

FLAC^3D是基于快速拉格朗日分析方法的岩土有限差分软件，能够准确描述岩土体材料的非线性破坏和流动特征，数值计算收敛效果好，广泛应用于各种三维岩土工程研究与设计^[21-25]。FLAC^3D内置三维边坡强度折减法计算方案(factor-of-safety)^[21]，但是对于复杂的边坡模型，该内置方案计算耗时较长。陈育民等^[25]采用“二分法”对该内置方案进行了改进，有效缩短了强度折减系数(SRF)的计算时间。虽然上述二分法强度折减法有效地提高了SRF计算效率，但每次迭代循环都重复计算模型的初始应力状态，这在一定程度上影响了计算效率。苏飞等^[26]采用FLAC^3D内置的Fish语言在每次循环之初载入初始应力场矩阵，省去重复计算的过程，所以计算时间明显缩短。但此方法每次载入初始应力场矩阵时，都需将矩阵中所有数值元素一一对应到模型单元中，操作复杂且影响计算效率。此外，上述两种改进方法并没有对迭代计算的时步上限进行合理化分析，而时步上限值将直接影响计算结果的准确性和计算耗时。因此上述两种改进方法仍存在较大的优化空间。

本文基于FLAC^3D 6.0集成的Python编程语言，对二分法强度折减法进行了初始应力场载入和计算时步上限的优化，从而提高三维边坡稳定性计算效率，有效缩短计算耗时。利用该优化方法对某公路高填方路基边坡进行了三维稳定性计算，结果表明该优化方法可以快速而准确地计算三维边坡稳定性，对类似填方路基边坡稳定性分析和支护设计具有重要的参考价值。

1 FLAC^3D内置强度折减法及已有的改进方法 1.1 内置折减计算方案

FLAC^3D采用数值计算的收敛性作为边坡失稳判据。内置折减计算方案：首先将所有岩土体强度参数设置为一个大值；然后计算体系达到平衡的典型时步N_r，并设定折减系数上下限分别为64和0；最后在此区间内选定初始迭代值F_tri进行体系平衡计算。若体系不平衡力比小于10^-5，则认为计算收敛，F_tri向上取值；反之计算不收敛，F_tri向下取值。除采用不平衡力比小于10^-5作为收敛依据外，当某次迭代不收敛时，程序将保持折减系数不变，并继续计算一个N_r，直到两次计算的不平衡力比差值小于10%或同一折减系数下完成6次N_r计算，判定体系不收敛；反之则判定收敛。折减系数迭代值容差为前后两次试算的折减系数平均值的0.5%，当不大于该容差值时，整个强度折减法迭代运算结束。软件还提供了折减系数上下限、不平衡力比及迭代值容差的自定义功能，以满足用户的特殊需求。

1.2 强度折减法改进

针对软件内置的强度折减法存在的问题，陈育民等^[25]对该方法进行了改进，开发了“二分法”强度折减法：首先设定折减系数上下限和容差值(默认0.02)；然后计算折减后的强度参数，并给模型赋值，设定每次迭代计算时步(默认10 000步)，进行应力场计算，根据设定的不平衡力比判定体系是否收敛；最后根据设定的折减系数容差，判定循环是否继续。

该方法每次迭代计算都要清空模型并重新计算体系的应力场，所以额外增加了较多的计算耗时。鉴于此，苏飞等^[26]利用FLAC^3D自带的Fish语言，在每次循环计算之初将初始应力场矩阵映射到模型中，从而省去了重复计算应力场的累赘。但是该改进方法的应力场映射过程操作复杂，仍然存在明显的计算耗时，并且默认每次迭代时步上限为60 000步，当某次迭代无法收敛时，将极大延长计算耗时。

2 强度折减法优化开发 2.1 实现过程

Itasca开发了一个专门用于Python扩展编程的程序库，可以实现在Python环境下利用FLAC^3D求解器完成扩展操作。本文对强度折减法的优化开发正是通过Python编程语言的灵活性实现的。初始应力场可以直接载入到新模型中，从而使得操作变得极为简单，且几乎不新增计算量。具体实现及开发过程如下。

1) 导入Itasca模块，并设置Python计算状态：

import itasca as it

it.command (‘python-reset-state false’)。

2) 设定折减系数容差δ、上下限(SRF_upper、SRF_lower)以及初始折减系数(SRF_trial)。

3) 设置循环条件(折减系数上下限>容差值)，开始迭代运算：

① 折减强度参数；

② 用“it.command()”命令载入初始应力场，并进行应力-应变计算；

③ 获得体系单元的不平衡力比R(it.zone.mech_ ratio())，判断在设定的时步上限内体系是否收敛。如果

a. R=it.zone.mech_ratio()＜R_set=10^-5，体系收敛，增大折减系数。

b. R=it.zone.mech_ratio()>R_set =10^-5，体系不收敛，减小折减系数。

4) 当折减系数上下限＜容差值时，迭代结束，强度折减完成。

整个优化开发实现流程见图 1。

2.2 时步上限优化

时步上限指的是一次迭代计算所能达到的最大计算时步数，对计算结果影响明显：如果时步上限太低，原本可以收敛获得更低不平衡力比的体系会被约束在较大的不平衡力比，因此SRF会朝着小值侧迭代，从而使最后结果偏低；如果时步上限太高，当体系无法达到设置的最小不平衡力比时，迭代计算必须在达到时步上限时才结束，因此浪费计算耗时；而当选择合适的时步上限时，既可以获得准确的SRF值，也可以减少计算耗时(图 2，其中N_u1、N_u2和N_u3分别表示较低的时步上限、适中的时步上限和较大的时步上限)。

软件初学者一般很难准确选择合适的时步上限，因此计算结果的准确性难以保证。为了提供一种方便快捷的确定时步上限的方法，本文对此进行了优化研究。

为了考察时步上限对复杂程度各异的边坡强度折减法计算结果的影响，本文精选了4个边坡模型算例，所有模型网格划分采用六面体网格为主、四面体网格为辅的方式(图 3)，模型岩土体物理力学参数见表 1。所有初始应力场计算采用弹性本构方法，初始应力场计算时步为N，然后将时步上限设置为1.00N、1.25N、1.50N、1.75N、2.00N、2.25N、2.50N、2.75N和3.00N，分别进行强度折减法计算，并记录计算得到的SRF。4个边坡模型相关内容介绍如下：

模型1    由Cai等^[27]设计的理想边坡模型，高10 m，坡比为1∶1.5，由单一土层组成(图 3a)。Cai等^[27-28]采用强度折减法计算得到此边坡SRF分别为1.14和1.20，Cheng等^[29-30]采用三维Spencer极限平衡法计算得到该边坡的安全系数(FOS)为1.18。模型网格单元数为12 572，初始应力场计算时步为3 384。

模型2    为Arai等^[31]提供的第2个边坡模型，由上中下3层土组成，其中中层土相对较为软弱(图 3b)。Arai等^[31]采用简化Janbu法计算得到的安全系数为0.405，Greco^[32]采用Spencers法结合蒙特卡洛随机搜索法计算得到安全系数为0.39，Kim等^[33]采用Spencers法结合随机特征搜索法计算得到安全系数为0.44。模型网格单元数为7 221，初始应力场计算时步为10 558。

模型3    根据Tandjiria等^[34]的填方路堤模型进行修改，在黏土地基上填筑路基，并在坡顶施加10 kN/m的交通荷载(图 3c)。采用简化Bishop和GLE/Morgenstern-Price法计算获得的安全系数分别为1.564和1.556。模型网格单元数为3 477，初始应力场计算时步为8 595。

模型4    根据澳大利亚计算机应用协会(ACADS)边坡分析第4考题进行修改，由土层①和软弱饱和夹层②组成，土层中含地下水，并在坡脚外形成静水位面，水位高出地表 3 m(图 3d)。采用简化Bishop和GLE/Morgenstern-Price法计算获得的安全系数分别为0.739和0.765。模型网格单元数为18 849，初始应力场计算时步为10 619。

4个边坡模型不同时步上限的强度折减系数计算结果见表 2和图 4。当时步上限较低时，获得的SRF值明显小于最终的SRF值，并且随着时步上限值的升高，SRF值趋于收敛，收敛获得的SRF值与已有的极限平衡法计算结果基本一致。根据计算结果分析，当时步上限设定在1.50N时(图 4红色虚线所示)，既能保证计算结果的准确性，又能提高计算效率，并且从剪切应变率云图上也能清晰反映出屈服区贯通的特征(图 5)，这与其他有限元软件整体失稳判据吻合^{[15, 21]}。

3 工程实例应用与分析 3.1 工程背景介绍

此工程实例为某库区复建二级公路一高填方路基，填方边坡为7级边坡，第1级边坡高8 m，坡率为1∶1.50，第2-7级边坡高度均为10 m，坡率均为1∶1.75(图 6)。填土弃渣来自复建公路沿线开挖土石料，其强度参数见表 3。

填方边坡进行稳定性分析时计算工况分为：1)正常工况，边坡处于天然状态下的工况；2)非正常工况Ⅰ，边坡处于暴雨或连续降雨状态下的工况；3)非正常工况Ⅱ，边坡处于地震等荷载作用状态下的工况。

其中，非正常工况Ⅰ中填土全坡体饱和，采用饱和强度参数进行稳定性计算；非正常工况Ⅱ计算时，考虑本项目位于地震烈度Ⅷ度区，设计地震动峰值加速度为0.2g，该工程为二级公路抗震重点工程，因此抗震重要性修正系数C_i=1.3，水平和竖向地震影响系数分别为：

式中：α_h为水平影响系数；α_v为竖向影响系数；C_z=0.25, 为综合影响系数；A_h=0.2, 为水平向设计基本地震动峰值加速度；A_v为竖向设计基本地震动峰值加速度(A_v=0.5A_h)；采用拟静力法进行该工况边坡稳定性计算。

3.2 边坡稳定性计算

为了对比本文优化的强度折减法在计算效率和准确性上与FLAC^3D内置方法、陈育民等“二分法”强度折减法和苏飞等改进的强度折减法的差别，分别用上述4种方法对该填方边坡进行三维强度折减法稳定性计算。为了便于横向对比，所有模型均在相同的计算机配置和工作状态下完成，计算机配置为：Intel酷睿i7 5500U(2.4 GHz)处理器，8 G DDR3L内存。

模型采用六面体为主的网格划分方法，共划分网格130 875个，初始应力场计算时步为4 908，因此在优化方法中，时步上限设定为7 362。陈育民等方法时步上限为默认的10 000，苏飞等方法并没有提供时步上限取值的信息，故而按照陈育民等方法，设置为10 000。模型边界条件设置方法：底面为全约束，4个侧面采用法向速度约束。

内置方法无法处理弹性本构模型，所以在使用该方法进行强度折减计算时，可以将基岩强度参数设置一个大值(c=1 MPa，φ=50°)；在其他方法中，基岩只参与初始应力场计算，不参与强度折减计算。数值计算过程中采用非关联流动法则，材料抗拉强度为0，初始应力场为自重应力场。计算结果见表 4。

计算结果表明，3种不同工况下，4种方法获得的SRF值基本一致。从剪应变率云图(图 7)上也可以看出，4种方法强度折减后边坡的屈服特征也保持一致。从计算耗时方面进行比较，本文优化方法明显短于其他3种方法，平均计算效率提升率分别为331%、103%和53%，随着模型复杂程度的增加，计算效率差异会更加明显。

本工程实例研究表明，在陈育民“二分法”强度折减法基础上经过优化以后，可以在保证计算结果准确性的同时极大地提高三维复杂边坡稳定性计算效率，从而克服三维稳定性计算耗时的缺点，便于在实际工程中推广应用。

4 结论

本研究在现有“二分法”强度折减法的基础上，利用FLAC^3D 6.0版本自带的Python编程语言，对上述方法进行了优化，主要获得以下结论：

1) 通过Python编程语言，直接将初始应力场载入到强度折减流程中，简化了初始应力场载入过程。

2) 通过多个模型的对比计算研究，认为将时步上限设定在1.50倍的初始应力场计算时步时，既可以准确获得边坡SRF，同时减少迭代计算时间。

3) 通过工程实例分析，表明本文优化方法可以明显提高三维边坡强度折减法计算效率，且与其他方法计算结果一致。

4) 本文优化方法操作简单，有利于三维强度折减法在边坡工程中的推广应用。