基于离散上限机构的复杂边坡稳定性分析*

侯超群 张亭亭 孙志彬 李见飞

(合肥工业大学汽车与交通工程学院，合肥 230009)

边坡稳定性分析既是经典的土力学问题，也是实际边坡工程中的重要课题。传统的边坡稳定性分析方法是极限平衡法。该方法基于摩尔-库仑屈服准则及静力平衡条件，计算简单且具有一定的精度，在工程中得到广泛应用［1］。但该方法忽略了土体的应力-应变关系，且对土条间的作用力分布进行了假设，不能代表边坡土体的真实应力状态［2］，因此在稳定性分析中具有一定的局限性。

随着计算机的发展，数值分析方法在边坡稳定分析中也受到了工程人员的青睐。相对于极限平衡法，数值分析方法不仅可以判断边坡的稳定状态，同时可以得到土体内部位移、应力以及应变的分布规律［3］。但数值方法的结果易受网格大小、迭代次数、边界条件等因素的影响，文献［4］指出：对于某些存在软弱夹层、节理等的复杂边坡，上述参数设置的不同可能会导致计算结果的较大改变。对于多级复杂边坡，数值模拟的结果会忽略除主要滑动面外的潜在滑动面，因此其适用性受到了限制。

除上述方法外，极限分析方法也是边坡稳定性分析的重要方法。该方法建立在塑性力学的基础上，利用上、下限定理对边坡临界荷载或者稳定性指标进行求解。运用上限定理时需构建机动速度场，运用下限定理时需构建平衡应力场。由于机动速度场的构建较为方便，上限定理的应用也更为普遍，其分析过程又被称为上限分析［2］。

传统极限分析大多针对单级均质边坡［5-7］。近年来，也有学者对非均质边坡和多级边坡(统称为复杂边坡)进行分析［8-11］。但其破坏机构的构建过程和功率计算较复杂，适用性不强，因此存在一定的局限性。年廷凯基于“组合对数螺旋线”的破坏机制对分层多级复杂边坡进行上限分析［12］。该方法具有一定的普遍性，但仍需要预设滑动面形状为对数螺旋线，且功率计算时采用数值积分方法，推导过程仍不十分简便。

近年来，Mollon针对隧道的非均质掌子面稳定性问题提出了基于离散的上限分析机构［13-14］，文献［15］将该机构推广至边坡上限分析，并研究了土体参数线性变化的非均质边坡稳定性问题。该机构不需提前假设滑动面的形状，在严格满足相关联流动法则的前题下，采用“点到点”的方式构建速度间断面，计算过程简单明晰。在前人研究的基础上，本文将该机构拓展至复杂边坡的上限分析。通过对数个算例的对比，验证该方法的计算精度，为解决分层多级复杂边坡上限分析提供了新的思路。

1 上限分析离散法基本原理

边坡离散破坏机构如图1所示。与传统的对数螺旋线滑动面不同，该机构的滑动面被离散点Pi划分为若干直线单元PiPi+1。

图1 离散结构示意Fig.1 Schematic diagram of the discretization surface

机构采用“点到点”的方式确定离散点 Pi，即某离散点Pi+1的位置由上一个离散点 Pi决定。机构生成从坡趾C(点P0)开始，并满足以下条件：

1)相关联流动法则。间断面的滑移速度v→i与滑移面方向的夹角为内摩擦角φ。

2)单元长度的均匀性。点 Pi、Pi+1与旋转中心O连线的夹角固定为δθ。

3)机构假设为刚体，即间断面上下均为刚体。

结合几何条件及向量关系［15］，由点 Pi推导点Pi+1的计算式如下：

式中：xi、yi分别为点 Pi的横、纵坐标;xi+1、yi+1分别为点Pi+1的横、纵坐标;xO、yO分别为点O的横、纵坐标;θi为水平方向与OPi的夹角［15］。

当离散点Pi的纵坐标yi大于(或者等于)边坡高度H时停止计算，破坏面生成结束。若yi＞H，采用线性插值的方法调整yi的值，使yi=H。

利用上限定理进行求解时，需要计算破坏机构的外力功率及内能耗散。该机构的外力功率为破坏块体ABC的重力功率，计算时需要对该块体进行离散。文献［15］采用了三角形块体的离散方式，即将破坏面AC上的相邻两个离散点(Pi、Pi+1)分别与点B连接(图2)，将ABC划分为若干三角形离散块体BPiPi+1。边坡重力功率即为离散块体的重力功率之和：

式中：γ为边坡土体的重度;ω为边坡块体的旋转角速度;Si为离散块体的面积;RGi为离散块体的重心PGi到 O点的距离;θGi为水平方向与OPGi的夹角。

图2 三角形离散块体功率计算分析简图Fig.2 The power analysis of the triangle discrete block

由于机构ABC为刚体，因此只考虑发生在滑动面上的能量耗散。机构的内能耗散即为各直线单元PiPi+1上能量耗散之和：

式中：c、φ分别为点Pi处的黏聚力和内摩擦角;Li为PiPi+1的长度;Ri为点O到点Pi的距离。

需要指出的是：计算重力功率时块体离散方式并不唯一，除三角形离散块体外，另有四边形离散块体等形式，本文将在下节复杂边坡的离散上限分析时进行讨论。

2 复杂边坡离散机构分析

2.1 复杂边坡离散机构的生成

分层多级复杂边坡如图3所示。边坡高度为H，共分为 m级，每级高度 αkH，α1+α2+… +αm=1，坡角为 βk。边坡每级设置台阶，宽度为 AkBk。边坡土体非均质，共分为n个土块，相邻土层分界线为lj。各层土体重度γ相同，内摩擦角φj和黏聚力cj各不同。其中 k∈［1，m］，j∈［1，n］。

图3 复杂边坡离散机构Fig.3 Discretization mechanism of complex slope

针对上述分层多级复杂边坡构建离散上限机构，由于不同土层的内摩擦角不同，采用分段形式进行构建。

1)滑动面构建从坡趾A0开始，构建过程中应保证单元体PiPi+1方向与Pi点速度方向夹角为土块1的内摩擦角φ1。Pi点超过分层线l1后暂停滑动面生成，利用线性差值使得目前滑动面上的最后一点D1落在线l1上。

2)土块2中的滑动面构建从点D1开始，构建过程中PiPi+1方向与Pi点速度方向夹角为φ2，并采用与步骤1)相同的方法生成D2。其余土层中的滑动面生成方式类推。

3)当Pi点纵坐标大于H时生成结束，采用线性插值求出滑动面与坡顶交点。

2.2 功率计算

复杂边坡离散机构不考虑其他工况的影响，外力功率为土体的重力功率，内能耗散为土体滑动面上的能量耗散。

2.2.1 外力功率

外力功率的计算方法与离散块体的选择有关，本文分别采用三角形块体离散与四边形块体离散方法对外力功率进行求解。

利用三角形离散块体进行功率计算时，需要根据边坡的几何轮廓进行土体划分。如图4所示，将台阶AkBk延长直至速度间断面A0C，交点为Ek，该多级边坡划分为土块 1(A0A1E1A0)、土块 2(B1A2E2E1B1)、…、土层 k(Bk－1AkEkEk－1Bk－1)…。假 设土块k的重力功率为Wk，则该复杂边坡的外力功率为：

图4 边坡块体划分Fig.4 The block division of complex slope

注意土块1的形状与土块2—土块m的形状并不相同，计算土块的重力功率时应分开考虑。对于土块1，其形状与简单边坡相似，计算时采用以A1、Pi以及Pi+1为顶点的三角形对其进行离散，其外力功率的计算方法与前节相同，在此不再赘述。

对于土块2—土块 m，由于其形状为不规则的曲边多边形，计算其土块功率时需做进一步处理。以土块k为例进行说明。

如图5所示：延长边坡斜线AkBk－1交速度间断面于点 Fk－1，此时土块 k的重力功率为土块 AkFk－1EkAk和土 块 Bk－1Fk－1Ek－1Bk－1的 重 力 功 率 之 差，将 土 块AkFk－1EkAk简写为 Ak，土块 Bk－1Fk－1Ek－1Bk－1简写为Bk－1，则土块 k 的重力功率为：

式中：WAk为土块 AkFk－1EkAk的重力功率;WBk－1为土块 Bk－1Fk－1Ek－1Bk－1的重力功率。

WAk和WBk－1的计算方法与土块1类似。

图5 土块k外力功率计算原理Fig.5 Principle of the gravity power calculation for block k

2.2.2 四边形离散块体

在文献［15］利用三角形离散块体进行功率计算的基础上，本文利用四边形块体进行功率计算，以验证其在复杂边坡分析中的适用性。如图6所示，经点Pi和点Pi+1沿横轴负方向做水平线，与坡面交于点Qi和点Qi+1，则边坡体可以划分为若干个四边形离散块体PiPi+1Qi+1Qi。

图6 四边形离散块体PiPi+1Qi+1Qi能量计算分析简图Fig.6 Analysis of the discrete block PiPi+1Qi+1Qi

求解每个四边形离散块体的外力功率并累加，即可得到该块体的外力功率，即：

结合O点的坐标，则可求得RGi的长度。

复杂边坡机构的内能耗散计算方法与文献［15］相同，与离散块体的选择无关，因此采用两种块体划分方式计算得到的内能耗散相同。

为验证上述计算过程的正确性，以简单边坡为例，对比离散机构与对数螺旋线滑裂面解析法得到的内外功率。计算参数为黏聚力c=20 kPa，重度γ=20 kN/m3，边坡倾角 β=60°，结果如表 1所示。可以看出：两种离散方法计算得到的外力功率与内能耗散相同，与解析法计算得到的结果十分相近。本例中，重力功率的最大误差为0.49%，内能耗散最大误差为0.32%，表明该机构有足够的精度。由于四边形离散方法的计算过程不需要进行土块划分，且根据计算经验，计算时间与三角形离散方法较为接近，因此优先推荐四边形离散。

表1 两种离散方法与对数螺旋线滑裂面解析法功率对比Table 1 Comparison of powers from the twodiscretization approaches and logarithmic spiral slip surface analysis

2.3 安全系数计算

待分析边坡可能处于稳定或者失稳状态。要求解边坡的安全系数，一般采用强度折减法，通过将土体的抗剪强度折减，使边坡达到临界破坏状态。折减系数Fs即为边坡的安全系数，经过折减后的抗剪强度参数 cf和 φf为［3，16］：式中：c、φ为原始抗剪强度参数;cf、φf为折减后的抗剪强度参数。

利用传统上限法结合强度折减法计算安全系数，可以得到关于Fs的隐函数，然后对其进行优化计算。在本离散机构中，由于内、外功率通过累加求和方式得到，其安全系数的计算需结合二分法。具体计算流程如下：

步骤1。设置安全系数的搜索范围［Fs1，Fs2］，令边坡安全系数 Fs=(Fs1+Fs2)/2，将其代入式(9)，得到折减后的抗剪强度参数。设置安全系数的初始搜索范围时，若没有可参考的Fs值，一般设为［0，5］，可涵盖一般的工程情况。

步骤2。根据上限定理，若能够在搜索域中找到一组 θ0以及R0，使得min W －D ＜ ε，则Fs为安全系数的上限解。此时令 Fs2=Fs，否则，Fs1=Fs。其中，ε 为规定阈值，当 δθ=0.1°时，ε 可取 5。

步骤3。重复上述过程，直至abs(Fs2－Fs1)＜δFs，优化结束。此时得到的Fs为安全系数的最小上限解。δFs为规定阈值，一般取 0.01。

具体求解过程如图7所示。

上述求解过程的关键在于步骤2中判断Fs是否为安全系数的上限解，求解过程可利用全局法或序列二次规划法对O点位置，即θ0和R0值进行优化搜索。

3 算例对比和分析

为验证离散机构在复杂边坡中的适用性，本文分别利用分层边坡、几何形状复杂边坡以及分层多级复杂边坡3个算例进行验证。

图7 边坡安全系数计算流程Fig.7 Calculation chart of the slope safety factor

3.1 非均质边坡

图8 为一分层边坡［16］，高度为 10 m，坡角为60°，该边坡土体划分为两层，土层1的高度为6 m。上层和下层的重度相同，摩擦角和黏聚力不同，具体参数列于表2。

图8 临界滑动面对比Fig.8 Comparison of critical sliding surfaces

表2 材料强度参数Table 2 Parameters of geotechnical materials

已有学者采用分段对数螺旋线法，简化Bishop法以及 Janbu法对该边坡进行了分析［9，17］，本文采用离散上限法对该边坡的安全系数进行求解，结果如表3所示，并在图8中绘制了不同方法得到的滑动面。

表3 安全系数对比Table 3 Comparison of the safety factors

由表3可知：采用离散机构得到的安全系数与已有结果较为接近，说明本文机构具有较高的准确性。从计算结果对比可以看出：由于极限分析得到的是严格上限解，因此该边坡安全系数的真实解应等于1.0，简化 Bishop法以及 Janbu法仍有改进的空间。本文方法与组合对数螺旋线机构的误差可能来自于离散过程，通过减小参数δθ可以进行优化。

由图8可以看出：本文得到的临界滑动面与简化Bishop法、Janbu法得到的圆弧滑动面相对接近。简化Bishop法、Janbu法得到的临界滑动面是连续的圆弧曲线，没有考虑内摩擦角对滑动面的影响。本文的离散法得到临界滑动面与组合对数螺旋线有相同的变化趋势。

注意到土层2相对于土层1变缓一些，这是因为随着土层2中土体内摩擦角φ的减小，单元块体PiPi+1与水平方向的夹角减小，导致破坏面变缓。工程实例表明：边坡破坏并非完全符合圆弧，尤其对于土体性质多变的复杂边坡，其滑动面情况更加复杂［18］，本文得到的滑动面更符合工程实际。

3.2 多级边坡

图9为某多级复杂边坡的剖面，边坡高17 m，从下到上每级占比分别为5∶5∶3∶4，土体均质，内摩擦角 φ 为30°，黏聚力 c为5 kPa，重度 γ 为20 kN/m3。

由于多级边坡的失稳不一定沿整个坡体，某局部滑动面也可能导致边坡滑移。文献［4，19］介绍了利用极限平衡法与联合搜索法，获得了该边坡的数个临界滑动面，并计算了不同滑动面的安全系数。

本文利用离散机构对上述滑动面进行验算，并将结果与上述文献进行对比，安全系数如表4所示，边坡临界滑动面如图9所示。

表4表明本文结果与已有解相差不大，证明了离散上限方法对该复杂边坡的适用性。特别地，本文结果大于极限平衡法的结果而略小于联合搜索法的结果。由于极限分析法给出的是严格上限解，表明Cheng的结果［4］更接近真实解。由图9可以看出：本文得到的临界滑动面与已有文献吻合较好，表明离散上限分析得到的滑动面有较高的可信度。特别地，本文的滑动面在某些局部与联合搜索法吻合更好。

图9 临界滑动面对比Fig.9 Comparison of critical slip surfaces

表4 安全系数对比Table 4 Comparison of the safety factors

3.3 分层多级边坡

利用离散机构对分层多级边坡进行稳定性分析。边坡情况如图10所示，坡高10 m，共4级，每级坡高2.5 m，相邻两级边坡之间的平台宽2 m，坡脚分别为 45°、30°、20°、45°。边坡填土共分为三层，各层土的材料参数见表5。

图10 填土分层边坡剖面Fig.10 The geological structure of the layered strata slope

表5 材料强度参数Table 5 Parameters of geotechnical materials

由于暂无其他理论解结果，将本文结果与FLAC3D中强度折减法的数值模拟结果进行对比。经计算，离散上限机构得到的安全系数为1.60，通过FLAC3D内置安全系数求解命令FOS得到的安全系数为1.63，两者较为接近。图11绘制了 FLAC3D得到的剪应变增量图及离散机构得到的临界破坏面。由于边坡滑移为土体剪切破坏，因此土体剪应变增量较大处一定程度上可代表边坡滑动面。可见本文结果与FLAC3D结果较为接近。该算例表明离散机构适用于较为复杂的分层多级边坡。

图11 临界滑动面的对比Fig.11 Comparison of critical slip surfaces

4 结束语

将离散上限分析方法应用于复杂边坡的安全系数求解，针对性地提出了复杂边坡的机构生成方式，内、外功率计算以及安全系数求解流程，其中计算边坡土体重力功率时，可采用三角形块体离散或四边形块体离散两种方式进行，两者结果基本相同。鉴于四边形块体离散不需要进行土块划分，较三角形离散更优。

通过分层边坡、分级边坡以及复杂分级边坡等三个算例的对比，离散机构能够用于分层多级复杂边坡的上限分析，得到的安全系数与现有结果相比误差较小，得到的滑动面也与其他方法较为接近。在保证计算精度的条件下，避免了复杂的积分计算。

该机构不需要提前假设滑动面的形状，可以根据实际地质情况构建滑动面，计算简便，适用性强，为解决分层多级复杂边坡的上限分析提供了有效工具和新的思路。但本机构目前只能分析以土体转动为主的边坡失稳，对于平动+转动的破坏模式下的离散上限分析，需要进一步研究。