吴敏强，蔡宜洲，周 豪

(三峡大学 水利与环境学院，湖北 宜昌 443002)

1 概 述

极限平衡法于1927年由Fellenius首次提出，是边坡稳定性分析中最常用的方法[1]。当今应用较为广泛的极限平衡法有Fellenius法[2]、简化Bishop法[3]、Janbu法[4]、Sarma法[5]、Morgenstern-Price法[6]等。极限平衡法的研究领域分布广泛，条间力的假设与滑动面的搜索是其主要应用方向，而且在三维扩展方面也成果显著。

关于刚体极限平衡法的研究，国内学界大致有如下几种观点。朱大勇等[7]通过假设一组合式为零的条间剪力分布，使整个滑体所有平衡条件得到满足，验证了简化Bishop法计算公式中没有条块间剪力，是因某种组合式为零，而不是条间剪力实际为零的设想。邓东平等[8]找到了一种用随机角来搜索随机滑动面、简化Janbu计算安全系数的新方法，通过随机模拟曲线，发现随机角物理意义明确，优势独特。郑颖人等[9]针对Sarma法条间剪力方程存在的两个问题：不满足边坡稳定性分析的合理性要求与不能正确表示任意条块的条间剪力，给出了新的条间剪力方程。

在整个发展过程中，学者们一直将多余未知量的假设选在条间作用力或作用位置上。对于严格的二维极限平衡法而言，这种处理方法尚能够得到令人比较满意的安全系数结果，但有时也存在收敛性问题。对于三维极限平衡法而言，由于假设较多，例如假设合理性、力学严格性，这种研究思路已经渐显不足。

Bell、朱大勇和郑宏等学者将研究视线转到滑面正应力的分布模式上来。这种极限平衡分析方法具有力学严格性，与其他严格的条分法相比，求解过程更便捷、易于编程实现且不存在收敛性问题。

严格极限平衡法发展潜力巨大，如采用此种方法对地震力作用下的边坡稳定性进行研究，将会产生可靠有效的理论成果，也有望取得不错的实际效益。如若研究成果的合理性与适用性被证明，并被广泛接受，大量应用，便可在工程边坡的稳定性分析中，获得省时省力、节约成本、提升速率的良好效果。

2 基本原理

在滑面正应力修正法计算分析中，若能准确构建滑面上的正应力分布模式，就可以求解出准确的安全系数。对于实际计算过程，滑面正应力分布基本无法直接获得。但若是假设一个合理的分布模式模拟真实状态，计算精度达到要求也是允许的。对于滑面正应力的分布模式，研究工作主要存在于两个部分：①滑面正应力假设模式；②求解安全系数。

2.1 滑面正应力假设模式

Bell提出将滑面正应力函数假设成含有两个待定参数的方程，这样就与安全系数这个未知数一起构成3个未知数，3个未知数3个方程理论上一定能求出数值解，这样安全系数在该假设模式下一定可以求出。假设形式见式(1)。

(1)

式中：λ1，λ2为假设的未知参数；σ(x)为滑面正应力分布函数；σ0(x)为瑞典法的分布函数；a，b为滑体两端点的横坐标。

郑宏教授提议将滑体看作积分域，通过数学原理将其转变为边界积分的研究方法不再局限于对滑体进行条分，可直接获得其安全系数。他建立了条块的平衡条件后发现，滑面正应力由两部分构成，其中滑体体积力和条间作用力各占一部分。通过进一步的研究使滑面正应力的具体力学组成结构关系更加清晰明朗[10]。同时，在他的研究中提出了一种新的滑面正应力分布的假设模式：

(2)

2.2 求解安全系数

图1为具有一般形状滑动面的边坡，y=s(x)为任意形状滑动面；y=g(x)为坡表地形线；单位宽度条块上重力w(x)；单位宽度条块上水平地震力为kcw(x)；表面水平方向上的分布荷载为qx(x)、竖直方向上的分布荷载为qy(x)；滑面上孔隙水压力为u(x)；总正应力为σ(x)；剪应力为τ(x)；y=yt(x)为边坡体推力线；滑面形心点坐标为(xc,yc)。向坡体设置预应力锚索并视为一个外部集中力对边坡整体进行极限平衡分析，其中P为每根锚杆(预应力锚索)所提供的锚固力，n为锚杆(预应力锚索)的数量，m为锚杆(预应力锚索)编号，i为锚杆(预应力锚索)在锚固过程中与水平线的夹角，(xpm,ypm)为锚杆形心点坐标。

图1 边坡滑体及受力分布

设边坡安全系数为Fs，滑面摩擦系数为tanψ(x)，有效黏聚力为c(x)，根据摩尔-库仑准则，可得：

(3)

安全系数的平衡方程如下：

(4)

(5)

(6)

式中：s′=ds/dx=tanα;α为计算点滑面的倾角。若能确定滑面正应力的分布函数，就可得到安全系数。

对边坡滑体的任一条块进行受力分析，见图2。

图2 边坡滑体受力模型

图2中，E为条间水平向作用力；T为条间竖向作用力；h为条块高度；ht为条间力的作用点高度；Δx为条块宽度。

(7)

(8)

根据摩尔-库仑强度准则：

(9)

由式(9)可见，滑面正应力由滑体体积力和坡面外力对滑面正应力的贡献(记为σ0(x))和土条间作用力对滑面正应力的贡献(记为σ1(x))两个部分构成，则可简记为：

σ(x)=σ0(x)+σ1(x)

(10)

(11)

(12)

在滑面已知前提下，σ0(x)属于已知函数。为了方程组能够求解，可以采用以下形式来对滑面正应力进行假设：

σ(x)=σ0(x)+f(x)

(13)

式中：f(x)为逼近函数，f(x)构造为含有两个待定参数的函数才能保证平衡方程能够顺利求解。逼近函数对滑面正应力的贡献是非常有限的，不需过分追求与真实分布之间的吻合度，可以选用线性函数逼近：

(14)

将式(11)、式(13)、式(14)和摩尔-库仑准则带入平衡方程，经整理可得：

(15)

(16)

(17)

解平衡方程，可以得到：

(18)

其中，

(19)

(20)

(21)

其中，

T0=A3B2-A2B3

(22)

S0=A1B3-A3B1

(23)

G0=A1B2-A2B1

(24)

则式(18)可以写成：

(25)

(26)

解式(25)，取实根得：

(27)

方程具有解析意义，是边坡严格极限平衡法显式解答，适用一般形状滑裂面。

3 地震力作用下边坡稳定性计算

3.1 均质边坡稳定性计算

取某均质土质边坡，边坡长100 m，边坡高20 m，边坡坡角为45°，边坡的断面尺寸见图3。该坡体重度为19 kN/m3，黏聚力为30 kPa，内摩擦角为35°。

图3 土质边坡断面

均质边坡的稳定性分别采用Slide软件和Geostudio软件来计算。由于边坡的坡比为1∶1，可知该边坡最大水平动力系数为2，最大竖向动力系数为1，地震烈度为Ⅶ度，参照《水利水电工程边坡设计规范》(SL 386-2007)取地震作用综合系数为 0.25，取地震作用重要性系数为 1.0。Slide软件Bishop法计算得到安全系数为1.571，Janbu法计算得到安全系数为1.488，见图4。

3.1.1 数值应力场获取

通过对均质边坡的二维建模以及材料参数的赋值，进行分析计算，以便从计算软件中提取出滑动面上各节点的有用应力信息。土坡的稳定性计算分析采用理正岩土计算软件，采用软件中内置的地震烈度考虑地震效应。

图4 均质边坡稳定性计算

3.1.2 有限元法与滑面正应力修正法结合

初始正应力分布函数σ0(x)通过软件得取可按照以下步骤操作来进行：

1) 条分整个滑体，读取滑面穿过的条块相对应位置的坐标。

2) 读取相应坐标点的应力信息(σx,σy,τxy)，参照式(28)可以换算成对应点的初始正应力信息：

(28)

式中：条块的底面中点与数值方向的夹角为β。单元的应力信息可直接由计算软件导出,当正应力初始分布形式通过σ0(x)数值方法得到后，再结合线性修正方法，对任意边坡安全系数的求解会更加方便、准确。

在Slide软件中，根据断面图所提供的信息建立计算模型，由提供的已知条件赋予相应的荷载，地震烈度采用7级，可由计算软件根据极限平衡理论自动搜索安全系数最小的滑面。在ANDIA软件中提取滑面上的应力信息，经过公式转换拟合出对应的初始滑面正应力分布函数，见图5。

图5 初始滑面正应力函数分布曲线

将得到的初始滑面正应力分布函数进行修正，计算得出边坡在地震力作用下严格极限平衡法的安全系数。

经过曲线拟合，边坡圆弧滑面正应力分布函数为：

y=0.0165x3-1.0566x2+17.38x+5.1878

计算在地震工况下边坡的稳定性，将得到的滑面正应力的初始正应力分布函数σ0(x)，根据式(29)对滑面正应力进行修正，其中λ和ξ为假定未知参数。

σ(x)=σ0(x)+λζ1(x)+ξζ2(x)

(29)

(30)

(31)

修正后的滑面正应力分布函数见图6。

图6 修正后滑面正应力函数分布曲线

经过曲线拟合，修正后边坡圆弧滑面正应力分布函数为：

y=0.0177x3-1.1266x2+15.366x+2.405

根据式(29)对滑面正应力进行修正，根据式(3)可以得出安全系数的显示解。选取Geostudio软件与严格极限平衡法结果进行比较，具体结果见表1。

表1 均质边坡安全系数汇总

根据计算结果可以看出，基于极限平衡理论的滑面正应力修正法与数值方法相结合的稳定性计算，可以直接从数值模拟软件中了解滑面底部的正应力分布情况，通过简单的数值计算即可获得比较准确的初始正应力，再对初始正应力进行修正，所得滑面正应力函数比较接近实际情况，所得安全系数的值与传统极限平衡法所得结果的相对误差非常小，从而验证了严格极限平衡法的合理性，便于实际应用。

3.2 非均质边坡稳定性计算

选取具有一般形状典型的非均质边坡，其剖面形状见图7。边坡坡高15 m，边坡比为1∶2，共有4个土层，各土层物理力学参数见表2。取地震峰值加速度为1.63 m/s2，对照中国地震烈度表取对应的地震烈度为Ⅶ度；参照《水利水电工程边坡设计规范》(SL 386-2007)取地震作用综合系数为 0.25，取地震作用重要性系数为1.0。

图7 非均质边坡剖面图

表2 土层物理力学参数

基于数值应力场的非均质边坡稳定性计算，可参照以下步骤进行：

1) 在Geostudio软件中，根据剖面图所提供的信息建立计算模型。

2) 可由计算软件根据极限平衡理论自动搜索出安全系数最小的圆弧滑面。

3) 在大型有限元计算软件ANSYS中，建立该剖面的二维有限元计算模型，其二维模型见图8，采用D-P弹塑性本构模型。

4) 在有限元软件中提取两种滑动面上的应力信息，经过公式的转换拟合出对应的初始滑面正应力分布函数。

5) 将得到的初始滑面正应力分布函数进行修正，计算得出安全系数。

6) 对比严格极限平衡法与传统极限平衡法所得结果。

图8 非均质边坡二维有限元模型

在进行材料参数赋值后对该边坡进行稳定性分析，对应极限平衡法搜索出的滑面位置在相应位置导出其应力信息，通过式(28)的转换，可得到选取单元上的初始正应力值。经过正应力修正后便可以拟合出修正后滑面正应力分布函数，见图9。

经过曲线拟合，修正后边坡圆弧滑面正应力分布函数为：

y=-0.0071x3+0.2277x2+8.5385x+41.095

按照均质边坡的计算过程，可以得出非均质边坡在地震作用下严格极限平衡法的安全系数，其结果汇总见表3。

图9 修正后滑面正应力函数分布曲线

表3 非均质边坡安全系数汇总

根据计算结果可以看出，在地震作用下严格极限平衡法计算得出的结果小于Bishop法、M-P法和瑞典圆弧法的安全系数，但是大于Janbu法和Bishop法折线滑面的计算结果。与Bishop法相比，结果相对误差为1.9%；与Bishop法折线滑面相比，结果相对误差为2%；但是与Janbu法相比，相对误差为11.1%。总体来说，误差范围在1.9%～11.1%之间。造成误差的原因主要在于正应力提取的不完整和正应力修正造成一定的误差，但在容许的范围内，说明严格极限平衡法在非均质边坡应用的合理性与适用性。

4 工程应用

4.1 工程概况

牙根水电站在雅砻江干流上，位于四川省甘孜州雅江县境内，是我国重要水电能源基地雅砻江流域的大型电站工程。勘探揭露坝址区右岸Ⅰ-Ⅴ勘探线之间发育拉裂松动岩体，其地表后缘边界高程自Ⅰ-Ⅴ线从2 640 m→2 667 m→2 690 m→2 697 m→2 597 m呈弧形变化，上游从PD17硐下部基岩突出处延伸至江边，下游侧从PD15硐延伸至江边。此外，在Ⅶ线高高程PD31平硐(硐口高程2 612.67 m)硐深0～57 m，Ⅸ线低高程PD33平硐(硐口高程2 513.74 m)硐深0～129.5 m，中高程PD35平硐(硐口高程2 564.66m)等处亦发现有拉裂松动岩体分布。通过勘探，坝址区右岸存在的拉裂松动岩体对大坝的整体布局、溢洪道等水工建筑物全面布置，以及对施工和运行期的安全都存在隐患，评价其在不同工况下的稳定性，提出合适的加固处理措施。选取2号拉裂松动岩体为对象进行稳定性研究。

2号拉裂松动岩体分布于Ⅲ区，即Ⅶ-Ⅸ线下游侧，拉裂松动岩体自Ⅶ线高高程斜向下延伸至Ⅸ线下游低高程江边，呈现出从上游至下游自高向低的斜河向形态分布，横向上受NE向陡倾角断层的控制，其地表出露面积约6.3×104m2，体积约为300×104m3。该拉裂松动岩体内地表发育的断层有：N20-40°E/SE∠80-85°的近横斜河陡倾的f11-f14共4条Ⅲ级结构断面(图10)。

图10 2号拉裂松动岩体平面地质

选取横Ⅶ剖面对2号拉裂松动岩体的稳定性作计算。其计算剖面简图见图11。

各地层物理力学参数见表4。

图11 2号拉裂松动岩体横Ⅶ剖面边坡计算简图

表4 岩体物理力学参数表

4.2 计算成果分析

利用严格极限平衡法，计算边坡在天然、暴雨、地震及地震+暴雨4种工况下的安全系数。暴雨工况表层参数c，f折减至90%，容重取饱和容重，地震工况的加速度为0.173 g，采用拟静力法来计算，利用Slide软件自带的Bishop法和Janbu法计算不同工况下滑面的安全系数，结果见表5。

表5 计算结果汇总

注：表5中数值表示差值相对于Janbu法结果的百分比。

按照前述的方法计算不同工况下滑面的应力分布函数，按照式(29)进行修正并进行拟合，结果见图12-图15。

天然工况下圆弧滑面正应力分布见图12。

则天然状态圆弧滑面正应力分布函数为：

y=0.00008x3-0.051x2+8.1193x+187.36

带入平衡方程中，解出该工况下的安全系数为1.193。

暴雨工况下圆弧滑面正应力分布见图13。

则暴雨工况下圆弧滑面正应力分布函数为：

y=0.00007x3-0.0528x2+9.3091x+64.77

带入平衡方程中，解出该工况下的安全系数为1.126。

地震工况下圆弧滑面正应力分布见图14。

则地震工况下圆弧滑面正应力分布函数为：

y=0.0001x3-0.0678x2+10.537x+45.782

带入平衡方程中，解出该工况下的安全系数为0.972。

地震+暴雨工况下圆弧滑面正应力分布见图15。

则地震+暴雨工况下圆弧滑面正应力分布函数为：

y=0.00006x3-0.0481x2+8.9016x+53.2541

带入平衡方程中，解出该工况下的安全系数为0.897。

图12 天然工况下圆弧滑面正应力分布

图13 暴雨工况下圆弧滑面正应力分布

图14 地震工况下圆弧滑面正应力分布

图15 地震+暴雨工况下圆弧滑面正应力分布

不同方法、不同工况计算结果汇总见表5。

通过对牙根水电站2号拉裂松动岩体采用基于数值应力场的滑面正应力修正法，分析边坡在不同工况下的稳定性以及得出相对应的安全系数，并与Bishop法、Janbu法的安全系数作对比，可得出如下结论：3种极限平衡方法的计算结果相差在3%以内，三者得出的结论基本吻合。将滑面正应力修正法与有限元法相结合，理论明确，操作简单，计算结果精度较高，适合大规模推广。

5 结 语

本文采用滑面正应力修正法与数值应力场的分析方法对边坡进行稳定性分析，其结果与传统刚体极限平衡法的分析结果相比较相对误差较小，进一步证明该方法可行性。通过对不同工况下的边坡进行稳定性计算发现，边坡在天然状态下安全系数较小，当遇到暴雨或地震时安全系数下降幅度大，需要及时治理。通过对牙根水电站2号拉裂松动岩体采用基于数值应力场的滑面正应力修正法，分析边坡在不同工况下的稳定性以及得出相对应的安全系数，并与Bishop法、Janbu法的安全系数作对比，可得出结论：3种极限平衡方法的计算结果相差在3%以内，三者得出的结论基本吻合。将滑面正应力修正法与有限元法相结合，理论明确，操作简单，计算结果精度较高，适合大规模推广。