边坡稳定性分析的离散元极限平衡法研究

2018-07-02曾亚武

水利与建筑工程学报 2018年3期

陈颉，曾亚武

(1.武汉大学土木建筑工程学院，湖北武汉 430072； 2.湖北工业大学计算机学院，湖北武汉 430065)

极限平衡法是最古老的边坡稳定性分析方法，其基本理论是假定边坡为刚体，可能沿着刚体中的某一滑动面滑动，计算时将可能滑动部分切分成条块，并对条块之间的作用力进行合理的假设，通过迭代试算获得边坡的安全系数。根据条间力的假定不同，可将极限平衡法分为瑞典圆弧法、毕肖普法、简布法、斯宾塞法、摩根斯坦-普莱斯法、莎尔玛法、不平衡推力法，以及通用条分法等[1]。这些方法都是对滑动面和条间力的三要素(力的大小、方向、作用点)进行了合理假设，降低了问题复杂度并简化了计算。但是，几乎每种方法都有遇到失效或者迭代计算不收敛的情况[2]，而且随着计算模型精密程度提高，计算量也大为提高。

有限元法[3]是最早引入边坡稳定分析的一种数值方法，其中以基于强度折减的有限元法[4-5]使用最为广泛，其基本思想也与极限平衡法中的强度折减相类似，不同的是该方法考虑了边坡岩土体的受力变形特征。实际的使用过程中，如何判定边坡的极限状态、确定准确的安全系数仍然是值得研究的问题。此外，在计算精度要求较高的情况下，有限元法的计算量较大，尤其在进行可靠度分析时，计算量巨大。为此李典庆等[6]将随机模型引入到有限元法进行边坡稳定可靠度分析时，结合了极限平衡法思想，有效的减少了计算时间，是一种有益的尝试。

离散元法引入边坡稳定分析中以后[7]，很快显示出对传统的极限平衡法以及有限元法的互补性优势。由于离散元法的物理模型简单，数学理论较为完备，在将边坡的宏观土力学物理特征参数合理的转化为细观颗粒物理力学参数后，能比较准确的预测边坡在各种工况条件下的运动趋势。但是，离散元法计算边坡稳定性一般都需要消耗大量的计算时间，而且计算时间随着颗粒数量呈指数级增长。即使在计算机硬件性能不断提高的情况下，离散元软件的计算时间并没有显著降低，其本质原因是离散元模型是一种在时间和空间上串行的计算模型。针对离散元模型计算效率不高的缺陷，结合计算机硬件和软件的优化方案也是一个当前研究的热点。常新正[8]将GPU并行计算引入离散元计算中，取得了一定的效果，但是该方案主要是针对颗粒在空间的碰撞检测过程进行优化，计算机效率提高程度受到了制约。田健[9]将Open-CL并行计算引入离散元计算中，使得加速过程与计算机硬件的相关性大大降低，适用性大为提高，但优化的主要方向依然是加速颗粒在空间碰撞的检测过程。

考虑到离散元方法和极限平衡方法各自的特点，本文提出了一种新的边坡稳定可靠度分析方法——离散元极限分析法。该方法结合了离散元法和极限平衡法的优点，模型离散化方面采用了离散元的思想，而在计算过程以及边坡稳定性判定方案的设计上运用了极限平衡法的主要思想，因而大大降低了计算过程的复杂程度，节约计算时间，且保证计算精度基本不变，可以为类似的工程案例提供一种参考。

1 边坡稳定性分析的离散元极限平衡法

离散元极限平衡法主要分为四步实现：第一步是将实体边坡截面离散化为许多个直径较大的圆盘，同时保证边坡的土力学性质，比如土密度，力学边界的刚度，空隙率等，圆盘与周边圆盘以及圆盘与刚性边界保持接触关系；第二步是利用极限平衡理论，计算出各个圆盘接触点上的法向压力与切向静摩擦力的大小，同时计算出切向摩擦力与法向压力的比值——将该比值定义为相对摩擦系数，并按照相对摩擦系数的大小对圆盘接触点进行排序，进而勾勒出潜在的可能滑动带；第三步是将可能滑动带以外的区域设置成刚性块体，仅将潜在滑动带区域用半径更小的圆盘进行覆盖，以减少圆盘的数量，利用刚性块体和小圆盘的静力平衡，再次计算出相对摩擦系数较大的新接触点集合，勾勒出最危险滑动线(面)；第四步是按照两次计算出的滑动带(面)中各个圆盘接触点的相对摩擦系数的分布律，综合考虑这些接触点摩擦系数在整体安全系数计算中的权重值，按照正态分布模型迅速的计算出在极限平衡滑动带区域内在统计意义上的最小边坡稳定系数。针对土质边坡的计算结果表明，该方法计算的精度是可以保证的，与传统的极限平衡法相比误差很小，比单独采用离散元法或者有限元法进行计算时间更短，效率更高。

1.1 实体边坡的离散化模型

边坡离散化的主要目的是将边坡按照有利于迅速找到潜在滑动区域带，有利于加速后续的稳定性计算过程的原则，划分为规则的若干区域，在保证土力学属性以及力学模型的前提下降低力学计算的复杂程度。在综合比对了诸多极限平衡法和离散元法的计算实例后[10]，提出一种综合两者优势的离散化模型，分两步计算出潜在滑动线(面)。第一步用直径较大的相互接触良好的圆盘覆盖整个计算剖面区域，如图1所示，并按照静力计算结果初步确定潜在滑移带区域；第二步在第一步计算结果的基础上采用较小直径的相互接触良好的圆盘覆盖潜在滑动带区域，而带外则采用刚性块体替代，如图2所示，再次计算并细化潜在滑移区域带，确定最危险滑动线(面)。

计算时，首先采用如图1所示的离散模型，计算出各圆盘接触点的法向压力和切向摩擦力，以及相对摩擦系数(切向摩擦力和法向压力之间的比值)，按照相对摩擦系数大小排序，初步定位出潜在滑动带区域；然后将边坡计算剖面划分为潜在滑动带和潜在滑动带以外的区域，滑动带以外的区域按照极限平衡的思想整体化为一个刚性块体，而潜在滑动带区域则用直径较小的圆盘覆盖，整个过程中要保持圆盘与圆盘、圆盘与刚体之间的相切接触。根据需要，可以设置自由边界(没有约束的边界)，如图1所示。

图1 边坡计算剖面的第一次离散化(整体覆盖)

图2边坡计算剖面的第二次离散化(分区覆盖)

这种计算方法可以减少圆盘的数量，简化计算模型进而加速计算。需要指出的是，在采用圆盘进行覆盖时，应该采用稳定的三角形排列方式，而不是正方形排列方式，因为前者在力学上是稳定的，后者不稳定。

在圆盘覆盖的过程中，对于不同的边界，圆盘的放置有一些规则。对于潜在滑动带与刚体之间的分界，圆盘要尽量保持与刚体的相切接触，对于垂直边界和水平边界也应如此，这些一般都可以通过调整圆盘直径和刚体尺寸来实现。针对个别不能严格相切的区域，采取一个折中方法使得相互接触的条件能够满足：用折线边界与邻近圆盘实质性接触，此法相对于圆盘填充方法，可以在不损失计算精度的前提下大量节省计算时间[11-12]。

1.2 离散化后各单元接触点作用力计算

离散后的模型各单元(圆盘或刚体)之间没有粘结，只有接触。对于圆盘单元来说，有两种接触，一种是圆盘与圆盘之间的接触，另一种是圆盘与刚体边界或界面之间的接触。根据本文约定的圆盘排列规则，一个圆盘最多与其他六个圆盘接触，且此圆盘属于“中间”圆盘，即只与圆盘接触，如图3(a)所示；除此以外，与刚体边界或界面接触或相邻的圆盘，称作“边界”单元，其接触数应小于六个，如图3(b)所示。

图3圆盘单元的接触

对于刚体来说，可以与若干个圆盘单元接触。由于要保证边界单元尽可能与刚体接触，刚体边界可以是直线，也可以是折线，理论上也可以是曲线，本文为简化计算，不考虑曲线边界情况。如图4(a)所示为刚体直线边界与多个圆盘的接触，图4(b)所示为刚体折线边界与多个圆盘的接触。

图4刚体与圆盘的接触

所有单元除了重力外，还受到相邻单元的作用力，该作用力作用在相应的接触点上，可分解为法向作用力和切向摩擦力，且法向作用力只能为压力，不能为拉力。单元受力分析见图5。

图5单元受力示意图

假定所有单元(圆盘和刚体)在重力、接触力和边界约束力的作用下保持静力平衡状态，即假定单元的强度及单元之间的摩擦系数无穷大，单元不会破坏，也不会产生相对运动。按照静力平衡要求，写出每个单元水平方向、垂直方向的受力平衡方程和相对于单元形心的力矩平衡方程，有

(1)

(2)

(3)

式中：i为单元接触编号，i=1，2，……，6；j为单元编号，j=1，2，……，n，其中n为单元总数。

对于圆盘单元而言，其形心为圆盘的圆心，重力作用于圆心，法向力都指向圆心，所以都不产生力矩，所有力矩均由接触点的摩擦力产生，且力臂都等于圆盘半径，因此圆盘单元的力矩平衡方程式(3)可以写成摩擦力的平衡方程，见式(4)。

(4)

将所有单元的静力平衡方程联立，得到整个计算模型的代数方程组，其中方程个数为3n，未知数个数为2m，m为单元与单元之间、单元与边界之间的接触点总数。由于每个单元至少与临近单元或边界有2个接触点(否则不稳定)，即m>2n，因此未知数个数2m应大于方程总数3n，即离散模型是超静定系统，需要补充其他条件或假定才能求解。

因此，为了使得整个离散模型的代数方程组可解，参照极限平衡法思想，做出如下三个假定：

(1) 圆盘与圆盘之间在水平方向的接触点具有相离的趋势，因此假定其法向作用力和切向摩擦力为零。

(2) 假定圆盘与垂直刚性边界之间的法向压力等于土力学中的静止土压力，即

N=(1-sinφ)×∑G

(5)

式中：N为圆盘与垂直刚性边界之间的法向压力(即水平作用力)；φ为内摩擦角；G为接触圆盘及其正上方所有圆盘的重力之和；并且假定圆盘与垂直刚性边界接触点的摩擦力为零(光滑接触)。

(3) 假定圆盘与水平刚性边界之间的压力等于该圆盘及其上部相应范围内圆盘的重力之和。

通过上述三个假定，可以使整个离散模型的代数联立方程组中的未知量大大减少，变成静定问题，从而可解，求出单元与单元之间、单元与边界之间接触点的作用力。在计算单个圆盘的受力时，先从边界的单元开始计算，进而通过迭代求解得到周围的圆盘的受力情况。

1.3 边坡潜在滑动带的搜索

在获得了边坡离散体系中每个接触点的法向压力和切向摩擦力后，就可以利用接触点的作用力来确定边坡潜在的滑动带区域。

首先，将边坡离散体系中接触点的最大允许切向力与法向力的比值定义为极限摩擦系数μmax，根据文献[13]中关于离散元细观参数与土力学宏观参数的实验结果，该极限摩擦系数与土体的内摩擦角φ之间存在如下的函数关系

μmax=φ/37.292-0.1044

(6)

文献[12]的实验显示极限摩擦系数μmax与土体黏聚力c关系不大，而且在数值上极限摩擦系数大于内摩擦角的正切值，说明实验结果中黏聚力的作用已经体现在了极限摩擦系数的计算结果中。而黏聚力为零的沙质土情况下的极限摩擦系数与其自身的级配、孔隙比以及含水率等有着密切关系，不在本文讨论范围内。

其次，将计算结果中每个接触点的切向力和法向力的比值定义为相对摩擦系数μre，即

μre=f/N

(7)

式中：f为接触点的切向摩擦力绝对值；N为接触点的法向压力值。

根据整个离散体系静力平衡方程，可求解出所有接触点的作用力，进而获得相对摩擦系数在边坡二维平面的一个分布。计算出相对摩擦系数可能大于、等于或小于极限摩擦系数，显然，相对摩擦系数大于极限摩擦系数的接触点，在受力过程中会发生破坏；相对摩擦系数等于极限摩擦系数的接触点，处于极限平衡状态；相对摩擦系数小于极限摩擦系数的接触点，则处于安全状态。也就是说，接触点相对摩擦系数的大小以及它们与极限摩擦系数的接近程度，反映了接触点的安全状态，而整个离散体系接触点相对摩擦系数的分布，则反映了整个边坡的安全程度。为此，进一步定义每个接触点的相对摩擦系数与极限摩擦系数之比为极限滑移比K：

K=μre/μmax

(8)

显然极限滑移比K值越大，则该接触点发生滑动破坏的可能性越大。将所有接触点的极限滑移比按照数值大小排序，并在整个边坡离散体系中搜索一系列极限滑移比较大，相互毗邻的圆盘接触点。这些圆盘在空间分布上一般是从坡面的某一位置一直连续的延伸到坡顶。如图6所示，给出了一个由极限滑移比较大的圆盘组成的潜在滑动带。

图6通过极限滑移比大小排序搜索获得的边坡潜在滑动带

只要圆盘单元直径相对于边坡尺度足够小，连接潜在滑动带中相邻圆盘中心的折线即可以认为是极限平衡法中的潜在滑动线(面)。但由于圆盘单元直径过小，会导致整个离散体系单元数量过多，计算量过大，难以快速获得精确结果。因此，采用本文提出的计算方案，先以直径较大、数量较少的圆盘单元来建立边坡的离散体系，给出边坡潜在滑移带的初步范围；然后再在初步计算的基础上，将潜在滑移带区域用更小的圆盘单元来替代，而该区域外则采用刚性块体来替代，以减少离散单元体系的单元数量，再进行精确计算，确定边坡最危险滑动线(面)。在第一次初定潜在滑移带范围后，第二次用小圆盘加刚性块体模型计算获得的更精细的边坡潜在滑移线(面)。

1.4 潜在滑动带内安全系数的计算

获得边坡的潜在滑动线(面)后，最简单的方法就是利用传统的极限平衡法求解获得边坡的稳定安全系数。

实际上，从概率统计的角度来看，可以认为，边坡离散体系中各接触点的极限滑移比与边坡的整体安全性是相关的。首先定义各接触点的极限滑移比的倒数为该接触点的安全系数Ks

Ks=1/K

(9)

显然Ks越大则该接触点产生滑移的概率越低，反之亦然。按照离散体系中约定的圆盘排列方式，一个圆盘最多可能和周边圆盘都接触产生六个接触点，如图3(a)所示。这些接触点中并不是每个接触点都被纳入计算滑动带内安全系数的统计样本中，因为有些接触点的安全系数较大，将会导致边坡的安全系数计算值偏大，这是偏于不安全的。为此，以第二次计算得到的最危险滑动线(面)相关的圆盘为中心，再加上与之相接触的两层圆盘为整体考虑的区域，将这个区域内圆盘接触点作为一个初步统计样本，但是要去除区域内圆盘与区域外圆盘之间的接触点。这些接触点的安全系数可以构成一个集合，计算其均值、方差、中位数值、最大值、最小值等。采用那些大于中值与中位数值，且与最大值相比大于三分之二的圆盘接触点作为最终需要在计算边坡安全系数中考虑的集合，该集合应该是滑动带内圆盘接触点的一个子集合，将这个集合定义为P。

对于集合P中的所有接触点的安全系数，可以计算出它们的平均值K′以及方差值σ。这些值符合正态分布，可以将其规则化为一个标准的正态分布，这些安全系数点的分布可以计算出安全系数。

2 关键问题的讨论

2.1 建模中刚性块体的设定

按照传统的极限平衡理论分析，在边坡达到极限平衡状态即将产生滑移失稳的瞬间，边坡在滑动面区域以外的应力分布并没有同时到达极限状态，其形状的整体性以及其它土力学性质基本是保持不变的。应力达到极限状态的区域大致在滑动面附近一个不大的区域内，因此将其滑动带以外大部分区域按照极限平衡法的思想设定为若干个刚性块体是基本符合边坡的土力学特征的，这样可以缩小圆盘所需覆盖的范围，简化进一步的计算。

从土力学原理分析，在潜在滑移带以外的区域有一部分是位于滑动面曲线的上部的(滑移带的凹向)，该部分在发生边坡失稳时有一定的整体性，因此可以整体上将这一个部分设置为刚性块体，这样在第二次建模计算中可大大减少圆盘单元数量，加快运算速度。但是如果采用三角形的刚性块体，在静力学分析中，其尖角受力分析比较困难，在文献[14]中特别设置了性质复杂的力学模型，这样无疑会加大计算量。为此，本文建立的是一个梯形刚体，计算分析表明是一个比较合理的方案。

2.2 圆盘直径选择及其与边界接触的处理

按照PFC2D的试验分析，圆盘的直径越小，得到的应力分析结果越准确，但计算量也会呈指数级增加。而且，圆盘直径减小到一定尺度后，进一步的减小不会实质性的提高计算准确度[15]。显然，模型中的圆盘直径越大，覆盖区域所需的圆盘越少，计算速度越快，但是计算精度不高；反之亦然。在分析实验结果的基础上，第一次计算时，模型中采用边坡高度值的2%至5%作为圆盘直径可以保证后续计算结果的精确度。

由于边坡离散化方案中采用的是三角形稳定覆盖，因此在圆盘覆盖的过程中，底部的圆盘可以与底部的水平边界完全接触；但是对于垂直边界及倾斜刚体边界，若采用直线边界，则难以做到与边界上的所有圆盘相接触，此时圆盘受力分析可能会产生较大的误差。传统的处理方法是用更小的圆盘进行填充，使边界尽可能与边界圆盘接触，这样做会增加计算量，且不符合本文的离散方案。为此本文考虑采用折线边界来代替直线边界，或者用矩形小刚块来契合边界和邻近的圆盘(可以简单的将这些小矩形刚性块体理解为刚性边界的突出部分)。这样处理，可以在不增加计算量的情况下满足计算精度要求。

3 计算程序与流程

上述离散元极限平衡法计算边坡稳定性的方法已基于MATLAB平台编制了相应的边坡潜在滑动带安全系数计算程序Slope-zone，其基本计算流程如图7所示。

图7离散元极限平衡法计算边坡稳定性的基本流程

4 算例分析

4.1 均质土坡稳定性计算

本算例为澳大利亚边坡委员会1986年给出的一个经典算例。某均质土层边坡，边坡高度10 m，上坡长度20 m，下坡长度40 m，土重度20 kN/m3，内摩擦角19.6°，黏聚力为3 kPa，建议的安全系数标准答案为1.00，也可以按照式(10)对安全系数进行预估。

(10)

式中：φ为内摩擦角；c为黏聚力；h为边坡高度；γ为土的重度；β为边坡的倾斜角度。

采用本文方法进行分析(第一次采用了3 256个圆盘，圆盘直径为边坡高度的4%，第二次采用4 230个圆盘，圆盘直径为边坡高度的2.5%，两块刚体区域)，并与毕肖普法、式(10)、离散元法结果进行比较，结果如表1所示。

表1 均质土坡算例计算结果

由表1的计算结果可以看出，本文方法的计算精度优于离散元法，且计算时间远小于离散元法。

4.2 水平层状边坡稳定性计算

某水平层状土质边坡[16-18]，由上下两层土体组成，坡高24 m，坡角37°，两层土体的厚度分别为18 m和6 m，内摩擦角分别为30°和20°，黏聚力分别为7 kPa和2 kPa，如图8所示。

图8分层边坡示意图

利用极限平衡法中的毕肖普法按照0.618法搜索潜在滑动圆弧得到的最小安全系数为1.424，利用PFC2D的计算结果为1.472，用本文本文方法计算得到的结果为1.441。其中，第一次采用了3 565个圆盘，圆盘直径为边坡高度的4%；第二次采用4 623个圆盘，圆盘直径为边坡高度的2.5%。计算结果如表2所示。

表2 非均质土坡算例计算结果

5 结语

本文提出了一种结合离散元法和极限平衡分析法各自优点的离散元极限平衡法，给出了分析边坡稳定性的步骤，并基于MATLAB平台编制了相应的计算程序。采用本文方法分别针对均质土坡和水平分层边坡进行了稳定性分析，并与常用的极限平衡法和离散元法计算结果进行了比较，结果表明本文方法不仅能保证较高的精度，而且比离散元法效率更高。

(1) 结合离散元法和极限平衡法的核心思想，提出了一种用于边坡稳定性分析的离散元极限平衡方法，将搜索边坡最不利滑动面的复杂工作转化为离散单元间相对摩擦系数的排序，计算复杂度大大降低，且能保证较高的精度。

(2) 本文方法的核心思想之一是将边坡最危险滑移面(线)的搜索分两步计算完成。第一步是初步计算，得到潜在的滑移带区域；第二步是精细计算，将第一次计算确定的潜在滑移带区域以外的部分以刚性块体替代，而在潜在滑移带区域以更小的圆盘单元进行覆盖，从而减少圆盘单元数目，提高计算效率。

(3) 针对本文提出的离散元极限平衡法，基于MATLAB平台编制了相应的计算程序，可用于简单边坡和复杂边坡的稳定性分析。针对均质土坡和水平层状边坡的算例结果表明，相较于离散元法，本文方法不仅具有更高的计算效率，而且具有更高的计算精度。

参考文献：

[1] 陈祖煜.土质边坡稳定分析——原理·方法·程序[M].北京:中国水利水电出版社,2003.

[2] 陈祖煜.岩质边坡稳定分析——原理·方法·程序[M].北京:中国水利水电出版社,2003.

[3] 郑颖人,赵尚毅,李安红,等.有限元极限分析法及其在边坡中的应用[M].北京:人民交通出版社,2011.

[4] 马晓雨.基于强度折减有限元法的土坡稳定性分析研究[D].大连：大连理工大学,2008.

[5] 刘汉东，贾聿颉.基于FLAC3D的边坡稳定性分析自编强度折减程序的修正[J].华北水利水电大学学报(自然科学版)，2015,36(6):59-62.

[6] 李典庆，肖特，曹子君，等.基于极限平衡法和有限元法的边坡协同式可靠度分析[J].岩土工程学报，2016,38(6):1004-1012.

[7] 王泳嘉，邢纪波.离散单元法及其在岩土力学中的应用[M].沈阳:东北工学院出版社,1991.

[8] 常新正.基于GPU的颗粒离散元计算方法研究[D].大连：大连理工大学，2013.