边坡稳定是岩土工程中的基本问题. 由于大量自然和工程边坡的稳定安全裕度不大, 近百年来岩土力学研究者与岩土工程师们持续不断地努力发展边坡稳定分析方法, 以提高边坡稳定分析的准确性和便利性. 边坡稳定分析方法可以分成刚体极限平衡法、有限元极限分析法、基于有限元应力计算的改进极限平衡法和极限分析法等.

极限平衡法通常也称作条分法, 将沿某一具有滑动可能的面上部的土体, 切成若干有限宽度的竖条或斜条, 把土条当成刚体, 根据静力平衡条件求得滑动面上的滑动力和抗滑力, 并以此为基础确定边坡的稳定安全系数. 条分法是评估给定滑动面安全性的有力方法. 最初出现的极限平衡法一般仅满足土条的力和力矩平衡中的部分条件. 如最早出现的瑞典条分法^[1]仅满足力矩平衡但不满足力平衡、简化Bishop法^[2]不满足水平力的平衡、Janbu简化法^[3]不满足力矩平衡. 这些方法的滑动面为圆弧, 无需迭代就能求得安全系数. 因其简便, 目前仍被岩土工程师所采用. 滑动面可以为任意形状且同时满足所有土条的力和力矩平衡的方法称为严格的极限平衡法. 严格的极限平衡法中, 基于摩尔库伦抗剪强度理论, 对所有土体的强度除以同一个折减系数, 使每个土条底部的剪应力与抗剪强度达到相等的临界状态, 折减系数即为给定滑动面的安全系数. 严格的极限平衡法求解的是静不定问题^[4], 需要增加额外的假定. 一般假定土条底部的法向力位于土条底部中点、且假定土条间法向力的作用位置或土条间作用力的方向后, 变成方程比未知量多1个的超静定问题, 安全系数需要迭代求解. 基于不同的条间力的假定, 严格的极限平衡法有十多个, 如Spencer法^[5]、Morgenstern-Price法^[6]、一般极限平衡法(GLE)、Sarma法^[7], Janbu精确法. 由于条间力的假定不同, 不同的严格极限平衡法的计算结果存在差异, 差异的大小与边坡情况和所给定的滑动面位置有关^[8]. 严格的极限平衡法, 尚存在以下问题: (1)复杂的情况下难以找到比较符合实际的最危险滑动面; (2)求解的土条间的作用力为合力, 渗流影响仅能通过浸润线位置反映, 土条间的作用力一般也不满足摩尔库伦强度准则, 因而, 对于复杂渗流场和多层土边坡情况, 边坡稳定安全性评估结果往往与实际差异较大.

有限元极限分析法分为载荷增加法^[8]和强度折减法^[9-13], 基于材料的弹塑性本构关系, 直接采用有限元求解边坡处于塑性极限状态时的载荷或材料强度折减系数. 此两方法并无实质差别. 在边坡稳定分析中因强度折减法在已有的程序上实施更简单而普遍采用. 1975年, Zienkiewicz等^[10]在研究土力学中的相关联流动法则与非相关联流动法则的文章中, 在算例部分用有限元法分析了一个均匀边坡的稳定性, 发现强度折减系数的倒数与极限平衡法计算的安全系数非常接近. Matsui等^[11]采用Zienkiewicz等^[10]的方法分析多个边坡的稳定性, 并把该方法正式命名为强度折减法. Griffiths等^[12]详细论述了如何将强度折减法与理想弹塑性有限元法结合分析边坡的稳定性, 并通过大量算例分析及与极限平衡法结果比较, 说明有限元强度折减法的有效性. 有限元强度折减法中安全系数大小与采用的岩土材料屈服准则密切相关^[13]. 传统的极限平衡法采用摩尔-库仑准则, 但是由于摩尔-库仑准则的屈服面是截面为不规则的六角形的角锥体表面, 存在尖顶和菱角, 给数值计算带来困难. 一些商业软件采用广义米赛斯准则, 或摩尔-库伦的外接圆屈服面. 这些与传统的摩尔-库伦准则的计算结果存在较大差异^[13]. 徐干成和郑颖人^[14]提出的等面积圆屈服面用于强度折减法计算边坡的稳定安全系数, 与Spencer法计算的误差在5%左右^[13], 表明采用合适的本构模型, 强度折减法不但可以找到最危险滑动面位置, 还能定量地计算稳定安全系数. 强度折减法中一般采用计算不收敛作为边坡临界的判据^[15], 由于有限元弹塑性迭代计算的复杂性, 计算结果与软件本身的弹塑性迭代算法关系较大. 笔者发现利用ABAQUS(版本6.14-1)软件采用摩尔-库伦模型对复杂的边坡进行计算时, 甚至得不到收敛的安全系数结果; 对于简单边坡浅表层滑动机制, 也得不到合理的安全系数, 如找不到坡角为30°单一砂土(无凝聚力) 边坡正确的塑性极限状态. 虽然强度折减法是边坡稳定分析的强有力的工具, 在国内土质边坡和岩质边坡的稳定性中也有较多的应用, 但强度折减法所得的安全系数至今仍尚未被国内相关规范所采用. 另外, 该方法对最危险滑动面以外的其他潜在滑动面的安全系数也难以分析.

极限分析用于边坡稳定分析, 是Drucker和Prager (1952)首先完成的. 他们用满足平衡条件且在任意一点都不违背屈服准则的假想应力场, 确定了破坏载荷的下限. 上限是通过一个与流动法则相容的速度场, 利用虚功原理, 即外力的虚功率等于或大于内部的能量耗损率^[16]. 陈惠发介绍了利用上限法求解直线、圆弧和对数螺线滑动面的边坡稳定极限分析方法^[16]. 极限分析法与条分法(包括垂直条分法和斜条分法)相结合, 通过求解虚功方程, 可获得给定组合刚体滑动机制的边坡稳定安全系数. 并通过优化方法, 求得临界的刚体组合滑动机制和最小上限安全系数^[17]. 条分法由于其固有的弱点, 难以适应复杂的地质情况和复杂的渗流与变形耦合作用问题. 极限分析与刚体有限元法结合的方法, 也得到了很大发展. 系统的塑性变形内能耗散仅发生在单元间的界面上, 故此类方法的计算网格, 需要适应滑动面的形状, 在优化过程中需要进行自适应网格重构^[18], 大大增加了程序算法的难度. 极限分析法还存在难以考虑复杂的渗流应力耦合相互作用的局限性.

有限元方法求解边坡的应力场, 可以考虑任意复杂的地质组成和渗流变形耦合情况. 将有限元应力结果与极限平衡法分析的概念结合起来的方法, 称为改进的极限平衡法^[19]. Kulhawy^[20]1969年提出的改进极限平衡法用一个独立的程序求出应力场后, 通过插值的方法获得滑动面上点的应力状态, 再用滑面上抗剪强度的代数和与剪切力的代数和之商计算稳定安全系数, 能求解复杂地质和复杂渗流变形耦合情况下边坡的稳定性. 这一方法对于圆弧滑动, 由于滑弧上每一点切向力对于圆心的力臂均相同, 实际上这抗滑力与滑动力在滑面上的代数和之比等于抗滑力矩与滑动力矩之比, 因而是合理的. 而对于非圆弧和非直线滑动面, 力的代数和概念不清, 存在问题. 对于任意滑动面, 刘艳章等^[21]将滑面上法向力和抗剪强度的合力定义为抗滑力矢, 并将整个滑面上抗滑力矢的和方向定义为主滑动方向, 提出采用给定滑动面上总抗滑力矢在主滑动方向上投影的代数和与总下滑力矢在此方向上投影代数和的比值定义安全系数, 进行了若干工程应用研究. 然而, 该方法对于圆弧滑动机制, 其安全系数与常规极限平衡方法显然不一致. 对于实际发生的圆弧与折线组合式滑动或多段折线组合式滑动, 其主滑面方向的物理含义也不明确, 因而, 该方法所求出的安全系数的物理意义不明确.

本文假定边坡的滑动机制为组合刚体滑动, 基于有限元应力场与运动许可的组合刚体速度场, 依据虚功原理, 提出了边坡稳定安全系数计算的新方法, 探讨了该方法所基于的有效应力场是否需要满足静力许可的问题, 并在典型算例中与其他极限平衡方法分析结果进行了对比.

极限分析法通过计算静力许可的应力场或机动容许的速度场, 根据虚功原理来计算极限载荷. 上限定理可表述为: 与机动容许速度场对应的外载荷不小于真实的极限载荷. 物体所能承受的最大外载荷, 与极限状态下虚速度作用下物体内部的虚功率成正比. 虚应变能可以作为表征物体承载力的单一指标. 下面根据边坡在机动容许的虚速度场下的虚功率来定义和计算安全系数.

非单一圆弧或非单一直线滑动面, 包括折线、非圆弧曲线、线段与圆弧的组合等, 统称为组合式滑动面. 如果将滑体作为刚体, 则不存在沿着组合式滑动面的机动容许的速度场的单一刚体, 因此, 组合式滑动面上的边坡体是刚体的组合.图1是与斜条分极限分析法^[17]中相同的滑动机制, 边坡被几个速度间断面分割成若干刚体. 速度间断面可以为圆弧、直线或对数螺线. 图中 $ABCDEF$折线称为底滑面或主滑面, $A$, $B$, $C$, $D$, $E$, $F$为主滑点, 2个端点$A$, $F$称为坡面主滑点, 其他主滑点称为内部主滑点, $BB_1$和$CC_1$等称为内部错动面或次滑面, $B_1$和$C_1$等称为次滑点. 下文所称的滑动面, 包括主滑面和次滑面.

边坡滑动时, 滑面上存在相对速度. 滑动体总的功率等于滑动面上的功率和, 可以表示为

其中, $P$为功率, $\tau $和$\sigma _n $分别为滑动面上的剪应力和正应力; $v$为滑动面上的相对速度; $\phi $为运动速度与间断面之间的夹角, 因剪胀而引起, 又称之为剪胀角, 不同滑动面上的剪胀角可以不相同; $\varGamma $为所有速度间断面.

边坡的滑动临界状态通过抗剪强度的折减达到. 刚体极限平衡法和有限元强度折减法都用临界状态时的强度折减系数的倒数作为安全系数. 假定边坡速度间断面上的正应力不受土体抗剪强度折减的影响, 临界状态时滑动面上的虚功率与强度未折减时依据有限元计算所得应力场的间断面上的总虚功率相等, 设安全系数为$F$, 则

其中$\tau _{f} $为滑面上的抗剪强度, 采用摩尔库伦强度准则

其中$c$为凝聚力, $\varphi$为内摩擦角.

由式(2)可以得到给定组合刚体边坡的滑动安全系数$F$的显式表达为

上式即安全系数等于滑动面上总的抗剪强度功率与总的剪应力功率之比. 分子和分母均是标量求和, 物理含义明确.

虚功率法基于机动容许的速度场和有限元应力场, 所得的安全系数是给定滑动机制的上限解. 边坡稳定分析中边坡的应力计算是一个独立的步骤, 因而可以不受限制地采用独立的有限元分析法获得任意复杂地质情况与耦合情况下的有效应力应力场.

速度间断面的位置可以通过分步优化的方法来确定, 以获得最危险滑动机制. 求边坡给定底滑面的稳定安全系数, 是内部错动面的位置经过优化后所得的最小安全系数. 而整个边坡的安全系数, 则是所有可能的主滑面中, 安全系数最小的那个.

二维有限元网格将滑动面分割为若干线段(或弧段), 可以首先计算滑动面上与有限元网格的交点, 找出这些线段及其所在的单元. 依据有限元计算所得的每个单元的高斯点应力, 通过等参逆变换插值法求出每个单元的节点应力. 滑动面上任意一点的应力可以通过单元的节点应力插值计算, 即

式中, $\sigma _{ij} $为所求点处的应力张量, 下标$i,j$表示笛卡尔坐标系空间坐标轴; $N_J $为计算点关于节点$J$的形函数; $\sigma _{ijJ}^{node} $为单元内节点编号为$J$处的应力张量.

滑面积分点上的应力向量$T_{i} $可由下式计算

式中, $n_j $表示滑面积分点处法线方向在整体坐标中的方向向量. 法向应力和剪应力分别为

其中$l_i $为切向方向向量. 滑动面高斯积分点的抗剪强度采用摩尔库伦公式(3)计算.

在有限元计算获得边坡的应力场以后, 通过构造组合刚体机动容许的速度场, 就可以计算得到给定滑动机制下的安全系数. 如图1所示的组合刚体滑动机制, 主滑面是$A$, $B$, $C$, $D$, $E$, $F$点构成的5段折线. 主滑面也可以由折线、圆弧(剪胀角为0)和对数螺线(剪胀角不为0)构成. 图中$AB$, $BC$, $CD$, $DE$, $EF$, $BB_{1}$, $CC_{1}$, $DD_{1}$, $EE_{1}$为速度间断面. 速度间断面上的运动方向与间断面成$\phi $角(剪胀角), 不考虑剪胀时$\phi =0$. 设第一个滑入段的速度为单位速度1, 根据速度三角形可以求出各间断面上的速度. 如图1(b)所示, 速度向量构成的速度三角形$ABC$中, 滑面上的相对速度的方向与$x$轴正方向的夹角等于滑动面的方向角加上(或减去)一个剪胀角. 可以计算出速度三角形的3个内角, $\angle A$, $\angle B$和$\angle C$. 根据三角形边角关系定理, 可得

根据滑面上的力和速度, 求得每一条速度间断面上的滑动力功率和抗滑力功率后, 依据式(4)即可计算得到给定滑动机制的边坡稳定安全系数.

寻找临界滑动面的优化方法是多变量优化问题. 虚功率法中的优化变量包括所有主滑点和次滑点的位置. 多变量优化中变量越多, 优化难度越大. 采取分步优化的方法可降低难度. 分成3层嵌套: 最内层是主滑点不动次滑点优化; 次内层是坡面主滑点不动内部主滑点优化.

最内层次滑面位置优化也是一个多变量优化问题. 对次滑点位置可以采用多轮单变量顺序优化, 将多变量优化问题, 化简为多次单变量优化问题. 次内层内部主滑点位置的优化也是多变量优化问题. 可以按照文献^[19]提出的下降差分方法和收敛准则, 采用梯度法准确找到极小值位置. 对于位于软弱结构面上的内部点, 可以规定只沿着结构面移动, 对于其他情况, 可以沿着次滑面方向移动. 最外层坡面主滑点的位置沿着坡面线移动, 其优化是2变量优化问题, 与内部主滑点优化方法相同.

作者开发的LinkFEA-Slope程序将以上理论方法纳入其中, 程序具有基于有限元有效应力计算给定圆弧或折线组合式滑动面安全系数的功能, 并采用梯度法优化寻找最危险滑动面.

1987年Donald 教授和Giam博士主持了澳大利亚计算机协会(ACADS)对澳大利亚所使用的边坡稳定分析程序进行了调查研究. 共设计了5个考题总计10个小题, 向120个单位发出测试邀请, 28个单位发回了计算答案. Donald教授也使用自编的GWEDWEM和EMU程序分析了考题, 并请以色列Baker教授用SSA程序、中国的陈祖煜教授用STAB程序和加拿大的Fredlund教授提供裁判答案. 最终综合各单位提供的答案、裁判答案和Donald教授的计算结果给出了推荐答案. 其中的第3题(EX3)的第2小题和第4题(EX4)是测试非圆弧滑动的, 用来检验本文理论方法的合理性.

EX3考题的边坡几何剖面和有限元网格见图2. 材料参数见表1. 地下水位位于软弱夹层(2#土)底部, 地下水位以上孔压为零. 第1问为圆弧或折线滑动面的最小安全系数(EX3-1), 第二问为计算指定的折线滑动面ABCD(控制点坐标如表2)的安全系数(EX3-2).

考题EX4的边坡如图3所示, 坡面上作用有竖向载荷、坡内有软弱夹层和地下水. 边坡材料参数见表3, 浸润线描述见表4. 要求确定临界滑动面和计算相应最小安全系数.

应力场采用ABAQUS(版本6.14)计算. 采用四边形孔压-应力耦合单元(CPE4P). 首先约束所有节点的位移, 求出稳定渗流场, 然后将负孔压节点的孔压修改为0后, 节点孔压作为已知条件, 采用线弹性本构模型求解有效应力场.

EX3的渗流场比较简单, 是浸润线位于弱夹层底部, 该部位以下孔隙水压力为静水压力, 其上为0孔隙水压力. EX4原题只提供了浸润线, 而有限元应力计算必须要有节点孔隙水压力. 有限元渗流计算中夹层的渗透系数取1#土渗透系数的50倍, 右边界取水面高程38.4 m的等水头边界条件, 左边界及台地地表取27.75 m水头边界, 斜坡取透水边界条件, 所得边坡渗流场浸润线与等孔压线如图4所示, 图中可见浸润线位置与原题所给有一定差异.

应力场基于线弹性参数计算结果的EX3-2与EX4的安全系数结果与文献中给出的裁判答案列于表5.

文献中29个提交的答案(同一研究者依据不同极限平衡方法提交多个答案)EX3-2安全系数平均值为1.292, 标准差(均方差)为0.064, 推荐答案为1.34$^{0}$. 本文EX3-2中给定折线滑动面的安全系数1.360, 比推荐答案大1.5${\%}$, 介于裁判答案之间. EX4文献中19个提交答案的平均安全系数为0.817, 标准差0.223, 推荐的裁判答案为0.78^[23]. 本文EX3给定折线主滑动面、优化次滑面位置的安全系数如图5所示; EX4的初始和优化后的折线主滑面折点坐标如表6. EX4折线滑动面如图6, $B$点次滑面方位角127.0$^\circ$, $C$点次滑面方位角75.4$^\circ$. 安全系数0.703, 介于裁判答案之间, 比推荐答案小9.9${\%}$. 2个含有软弱夹层的算例计算结果都介于裁判答案之间, 与推荐答案接近, 表明虚功率法对于折线组合式滑动方案的研究结果合理.

采用ABAQUS(版本6.14)软件中的强度折减法, 本构模型采用M-C模型, 单元类型采用2次三角形单元, 取不同的网格密度分别计算EX3与EX4的稳定安全系数, 发现塑性区图所显示的滑动机制基本一致, 而网格越密, 安全系数越小. 0.1 m边长网格模型计算所得的塑性区与稳定安全系数如图7所示, 安全系数分别为1.261和0.785(网格边长1.0 m的模型安全系数分别为1.285和0.810). 虚功率法应力计算的网格尺度大于1.0 m. 对比来看, 虚功率法的速度间断面位置与有限元强度折减法所得的结果基本一致. 而安全系数结果算例EX3基本一致, EX4小约10${\%}$.

依据极限分析的下限定理, 有限元强度折减法的应力场是静力许可的, 极限状态所对应的外载荷小于真实的极限承载力, 其安全系数结果理论上是一个下限值. 由于非线性迭代的收敛性和网格密度的影响, 计算所得的极限状态时的应力场并不能处处满足静力许可. 且由于塑性应变是局部化的, 因而一般需要采用细密网格和二次单元才能得到比较准确的结果(一次单元对网格密度要求更高). 网格越细密, 程序的非线性功能越好, 计算结果越接近这个下限值. 其缺点之一是要求程序具有强大的非线性计算能力, 这种能力往往在地质和渗流变形耦合情况下难以达到; 其二是需要细密的网格, 对于复杂的工程问题计算规模往往过大而成本过高, 甚至计算任务难以完成. 笔者还未见采用商业软件用强度折减法求渗流与变形双向耦合情况下的稳定安全系数的工程案例介绍.

极限平衡法对于同一土条内包含多种土层的情况误差较大, 且难以考虑复杂的渗流性状, 尤其是难以考虑渗流应力强耦合相互作用的情况. 而基于有限元应力计算的改进极限平衡法对于复杂滑动面安全系数计算又存在明显的理论缺陷.

虚功率法应力场计算独立于后续的稳定性分析, 采用有限元法完成, 复杂的地质和渗流应力耦合作用由有限元完成, 因而不受限制. 应力计算一次完成, 因而速度间断面位置优化的效率也远远高于极限平衡法. 复杂滑动机制又继承了极限分析上限法的优点, 因而对于解决地质条件和渗流应力耦合作用问题有显著的优点. 与有限元强度折减法相比, 一方面应力计算时对于材料非线性迭代计算的性能要求较低, 另一方面对于网格密度要求较低(不需要依靠计算形成塑性应变连通区得到稳定安全系数). 计算的效率较高.

虚功率法安全系数计算公式来源于极限分析的上限法, 其结果是一个上限解. 结果与实际的接近程度, 取决于假定的滑动机制的合理性及滑动机制的优化效果. 当然, 对于边坡稳定安全系数接近或小于1.0的情况, 线弹性模型计算所得的应力场仅满足静力平衡, 违反静力许可的区域较大, 安全系数计算的结果与实际值差异相对较大. 对于实际安全系数小于1.0的情况, 是不存在静力许可的应力场的. 而对于安全裕度不大的边坡, 通过采用弹塑性模型, 或对违反强度准则的单元, 进行应力迁移等非线性迭代计算, 减小违反静力许可的区域、降低单元区域应力超过破坏面的程度, 则安全系数结果会更接近于实际.

边坡的有限元应力场计算结果, 受土的本构模型及其参数影响. 基于以下条件的计算获得EX3算例4个应力场结果: (1)线弹性模型和如表1中的参数($A$); (2)线弹性模型且2#土也取表1中1#土的参数($B$); (3)摩尔库伦弹塑性模型和表1中的参数($C$); (4) 摩尔库伦模型且强度除以1.25折减($D$).

基于不同的应力场, 计算了EX3给定滑动面位置的安全系数、也计算了最危险圆弧滑动面的安全系数, 滑动面位置如图5所示. 对于给定的主滑面, 次滑面位置优化后方位角如表7. $BB_1$方位角最大差异1.5$^\circ$, $CC_1$方位角最大差异4.2$^\circ$, 优化后不同应力场的次滑面位置存在一定差异但差异不太大.

给定折线主滑面的安全系数和圆弧滑动面最小安全系数结果比较分别如表8和表9. 折线滑动面安全系数介于1.333$\sim $1.367之间, 最大差异0.034, 相对于A条件的幅度为$-$2.0${\%}$. 本例对于圆弧滑动, A应力条件下, 最危险圆弧位置如图6, 滑入、滑出和圆心点的坐标分别为(72.45,40)、 (42.54, 27.75)、 (50.35, 51.32). 不同应力条件下圆弧滑动安全系数最大差异幅度为0.9${\%}$. 可见本例中本构模型和参数造成的应力场差异对安全系数的影响很小.

弹性模型计算得到的应力场均满足静力平衡, 但不一定满足静力许可; 弹塑性模型计算所得的应力场, 既满足静力平衡, 也基本满足静力许可. 应力条件不同, 虚功率法计算得到的安全系数结果存在差异. 根据最小值原理, 一个机构的极限承载力, 不低于任意一个同时满足静力平衡和静力许可的应力极限状态所对应的外载荷. 因此对于应力场既满足静力平衡, 又满足静力许可的边坡, 无疑是稳定的, 其安全系数不低于1. 即使是稳定的边坡, 当安全裕度不高时, 无论是采用什么样的本构模型, 有限元应力计算时常常存在高斯点应力超过其抗拉和抗剪强度的情况, 而要通过弹塑性迭代或应力迁移迭代使应力结果满足强度准则, 则往往是十分困难的. 边坡稳定安全系数小于1.0的情况下, 根本就不存在静力许可的应力场. 要进行强度折减(提高)后才能获取静力许可的应力场. 因此, 对于安全系数计算所依据的应力场, 仅要求静力平衡而不要求静力许可一般还是合理的. 采用虚功率法, 依据静力平衡的应力场计算所得的稳定安全系数是边坡稳定性较好的度量.

本文提出了边坡稳定安全性计算的虚功率法理论和实施方法, 在2个典型折线滑动边坡算例中进行了检验. 并探讨了应力场静力平衡与静力许可条件对安全系数结果的影响. 得出以下结论:

(1) 虚功率法用边坡组合刚体滑动机动许可的速度间断面上抗滑力功率与滑动力功率之比计算安全系数, 解决了改进的极限平衡法中对于复杂滑动面安全系数计算公式物理含义不清晰的问题;

(2) 采用分步优化方法对组合滑动机制的速度间断面的位置优化, 能较容易准确确定最危险滑动机制;

(3) 虚功率法的应力计算独立于安全系数的计算与优化过程, 因而能考虑复杂的地质和渗流与应力耦合作用情况;

(4)采用静力许可的应力场, 安全系数结果更接近实际, 静力平衡的应力场, 其稳定安全性计算结果, 也是边坡稳定性的一个不错的度量.

虚功率法是为适应水电站库岸堆积体高边坡和复杂软弱水工构筑物地基的稳定性分析提出的新的稳定分析方法. 包含该方法的分析软件LinkFEA-Slope软件已经用于澜沧江如美和金沙江拉哇2个重大水电站工程的库岸堆积体边坡、碎裂卸荷岩体边坡或深厚软弱覆盖层高围堰的稳定性分析, 正在用于长江沿岸某物流码头软土地基上超高堆场基础稳定性分析. 采用虚功率法, 基于有限元应力场结果分析圆弧与非圆弧滑动机制的边坡、地基稳定安全系数, 得到了工程设计部门和业主的欢迎与鼓励. 目前作者研究团队一方面加紧改进复杂的渗流与变形耦合相互作用的计算程序LinkFEA-Stress, 意图应力计算不但包含渗流与应力的强耦合作用, 也包含拉应力迁移和塑性滑移在内的材料非线性模拟, 力图使应力场更接近静力许可; 另一方面加紧开发包括圆弧与折线和对数螺线管等复杂滑弧组合式底滑面滑动机制的安全系数计算与优化, 将在工程应用中不断改进与发展.