基于ABAQUS强度折减法分析库水位下降对边坡稳定性影响

郭 飞 付调金 阮荣乾 管宏飞

（1.三峡大学 三峡库区地质灾害教育部重点实验室，湖北 宜昌 443002；2.长江三峡勘测研究院有限公司，武汉 430074；3.水利部海委漳卫南运河乐陵河务局，山东 德州 253600）

水是诱发滑坡的主要因素.据不完全统计，三峡库区在175m水位范围内共有大小滑坡2 000余个.在日本，大约60%的水库滑坡发生在库水位骤降时期，40%发生在水位上升时期，包括蓄水初期［1].库水位骤然下降时，由于坡体内地下水位下降相对滞后，导致坡体内产生超孔隙水压力，对滑坡的稳定性会产生不利影响，因此，三峡库区水位下降对滑坡稳定性的影响是当前滑坡研究领域十分重要的研究课题.

目前对三峡库区滑坡稳定性评价和治理都是建立在饱和土假设基础上，对暴雨、库水位下降产生滑坡的机理研究较少，已有的工程设计均采用对暂态饱和区及暂态水压力进行假定的方法［2].同时，传统的极限平衡法难以有效地考虑水与岩土体的相互作用机理以及渗流场和应力场的耦合效应.基于饱和与非饱和的流固耦合理论的数值方法虽然能解决这些问题，但它是根据边坡的位移场、应力场、塑性区来间接地评价，不能直接给出安全系数以直观评价边坡稳定性.本文通过ABAQUS软件，考虑饱和与非饱和流固耦合效应的强度折减分析，可以直接获得一个安全系数，此方法不仅保持了有限元在模拟复杂问题上的优点，而且概念明确，结果直观，对实际边坡治理工程具有一定的参考意义.

1 树坪滑坡概况

树坪滑坡［3]位于湖北省秭归县沙镇溪镇长江干流南岸斜坡地段的古滑坡.上距沙镇溪镇千将坪滑坡约3km处，下距三峡工程三斗坪坝址47km.树坪滑坡属多期性巨型滑坡，它是由两个区块（见图1中区块1和区块2）组成，物质组成较复杂（见图2），树坪滑坡区块1大致可以分为滑体、滑带、滑床3部分，其中滑体有两种组成物质.

图1 树坪滑坡全貌

滑体1为坡积物（Qdl），主要为粉质粘土层夹碎块石，碎块石多呈次棱角状，以粉砂岩、泥质粉砂岩为主，结构松散～稍密，透水性较差.

图2 树坪滑坡工程地质剖面图

滑体2为滑坡堆积物（Qdel），主要为碎块石粘土层，碎块石成分主要为泥质粉砂岩、泥灰岩和灰岩.土为壤土、粉质粘土，填充于碎块石中，结构不均匀，稍密，滑体堆积物结构松散，透水性较好.

滑带为紫红色角砾石土层，较湿，结构紧密，土可塑.碎石呈次棱角状～次园状，碎石上可见擦痕，土层中可见明显揉皱、光滑镜面.滑带埋深较大，层厚一般在10～20cm.

自2003年6月三峡水库第一次蓄水至高程135 m后，库区发生大量的滑坡复活及失稳现象.2004年1月，树坪滑坡出现多处变形，宏观变形以地面裂缝和房屋开裂为主.根据地质勘察结果，坡体侧缘剪切带内出现新鲜错位，江面以下剪出口处，特别是在降雨时节，江面出现大面积浑水现象.但根据自动位移计的监测资料表明，目前树坪滑坡处于较稳定状态.

2 库水下降饱和与非饱和流固耦合数值模拟

2.1 滑坡饱和与非饱和渗流场与应力场耦合模型［4]

由应力场平衡方程（1）和非饱和渗流场连续方程（2）共同组成了流固耦合问题的控制方程组：

其中，方程（1）中IN是内力矩阵，IP是外力矩阵，此方程由虚功原理推导简化而来.

流固耦合数学模型的求解，还需要有相应的定解条件，即确定模型的边界条件和初始条件才能求解.本文研究库水下降对滑坡的影响作用，其边界条件需满足方程（3）和方程（4）：

其中，方程（3）为水头边界条件，方程（4）为初始孔隙水压力的分布函数.

该方程组还需同时满足以下初始条件：位移边界条件，即U＝u（x，y，z）；初始时刻（t＝0），位移或质点速度的初始值.

将流固耦合控制方程组（1）、（2）与在库水下降时相应的渗流场与应力场边界条件和初始条件相结合，就可得到解决非饱和流固耦合问题的控制方程，然后便可求解方程组.

2.2 计算模型确定

根据树坪滑坡的野外地质勘察和钻孔资料，选择计算剖面如图3所示，该剖面为树坪滑坡1号块体的主滑方向，沿该方向上布置有日本坂田电机株式会社生产的自动位移计，计算模型尺寸为：X方向上最大长度为946m（0＜X＜946m），模型后缘Y方向上最大高程为468m（0＜Y＜468m），前缘Y方向上45m.采用ABAQUS数值计算软件进行工况模拟，计算域包含滑体、滑带和基岩，整个计算域剖分网络单元12 390个，结点共25 330个.

图3 树坪滑坡计算模型

约束边界条件：左右边界约束水平位移，下边界约束水平和竖直位移，上边界自由.

渗流边界条件：右侧边界库水位以下为水头边界，库水位以上为零流量边界；左侧边界为水头边界.初始地下水位线根据地质钻孔资料得到的实际地下水位线，并通过两侧水头计算得出.

2.3 计算工况

根据三峡水库运行曲线图（见图4），每年5月至9月为汛期，库水位保持145m低水位运行；9月底至11月初，汛期过后为满足发电需要，水位从145m大幅上升至175m；11月至次年1月为枯水期，水位维持在175m高水位运行；1月至5月为腾出防洪库容，水位缓慢降落至145m.目前，较多学者研究库水位上升及骤降对库区边坡的影响，未充分考虑到从1～5月，水位缓慢降落情况，从水位调度图中可以看出三峡库区在此时间段内库水位下降平均速度约为0.21m／d，而大多数学者考虑水位骤降的影响，对于揭示古滑坡在库水位变动作用下的复活机理有一定意义，但不能反映水库正常运行下对古滑坡的影响.本文考虑水位正常调度情况下，库水位下降（下降速度约为0.21m／d）对树坪滑坡的影响.

图4 三峡水库运行水位调度图

2.4 计算参数选取

2.4.1 非饱和土渗透性参数

非饱和渗流分析涉及的参数除了饱和渗透系数外，还需要材料的初始含水率、土水特征曲线及非饱和渗透性函数.由于树坪滑坡非饱和土性状还未进行实验研究，则其滑体和滑带的非饱和参数采用工程类比法，同时根据地质勘察资料确定树坪滑坡3种材料的渗透性参数分别如下：

1）1号滑体：渗透系数取为2.8×10－5m／s，即2.419m／d.饱和度（Sr）与孔隙水压力（Uw）及渗透系数折减系数（Ks）关系见表1～2.

表1 树坪滑坡1号滑体Sr－Uw对应关系 （孔压单位：kPa）

表2 树坪滑坡1号滑体Sr－Ks对应关系

2）2号滑体：渗透系数取为6.0×10－5m／s，即5.184m／d.饱和度（Sr）与孔隙水压力（Uw）及渗透系数折减系数（Ks）关系见表3～4.

表3 树坪滑坡2号滑体Sr－Uw对应关系 （孔压单位：kPa）

表4 树坪滑坡2号滑体Sr－Ks对应关系

3）滑带：渗透系数取为5.0×10－7m／s，即0.043 m／d.饱和度（Sr）与孔隙水压力（Uw）及渗透系数折减系数（Ks）关系见表5～6.

表5 树坪滑坡滑带Sr－Uw对应关系 （孔压单位：kPa）

表6 树坪滑坡滑带Sr－Ks对应关系

4）滑床：渗透系数取为2.0×10－10m／s，即1.728×10－5m／d.因不考虑基岩的非饱和特性，所以，计算参数不含饱和度（Sr）与孔隙水压力（Uw）及渗透系数折减系数（Ks）关系.

2.4.2 材料力学参数

应力应变分析采用 Mohr－Coulomb非饱和土本构模型，计算参数主要有岩土体变形模量（E）、泊松比（μ）、容重（γ）、凝聚力（c）、内摩擦角（φ），主要由地质勘察资料所建议.滑坡计算参数见表7.

表7 树坪滑坡岩土体物理学力学参数

3 ABAQUS强度折减法原理

3.1 强度折减法原理［5]

强度折减法最早由Zienkiewicz等提出，后被许多学者广泛采用.他们提出了一个抗剪强度折减系数（SSRF：Shear Strength Reduction Factor）的概念，其定义为：在外载荷作用保持不变的情况下，外载荷所产生的实际剪应力与抵御外载荷所发挥的最低抗剪强度是按照实际强度指标折减后所确定的、实际中得以发挥的抗剪强度相等.当假定边坡内所有土体抗剪强度的发挥程度相同时，这种抗剪强度折减系数相当于传统意义上的边坡整体稳定安全系数Fs，又称为强度储备安全系数，与极限平衡法中所给出的稳定安全系数在概念上是一致的.

折减后的抗剪强度参数可分别表示为

式中，c和φ是土体所能够提供的抗剪强度参数；cm和φm是维持平衡所需要的或土体实际发挥的抗剪强度；Fr是强度折减系数.

3.2 失稳判据

1）以数值计算收敛与否作为评价标准.即根据有限元解的收敛性确定失稳状态.该失稳判据与一定的计算方法相关.因此，以数值计算收敛与否为失稳判据是不合理的，适用性差［6－7].

2）以特征部位的位移拐点作为评价标准.根据计算域内某一部位处结点的位移与折减系数之间关系密切曲线的变化特征确定失稳状态，当某处一结点位移突然迅速增大，则认为边坡发生失稳［8].

3）以是否形成连续的贯通区作为评价标准.理论上，边坡的变形过程总是伴随着一些物理量的出现和发展，如塑性区、塑性应变、广义剪应变和应力水平等，当这些物理量达到一定的值时，边坡失稳［9].

本文采用强度折减法分析库水下降对树坪滑坡影响时，主要采用第2）和3）两种判据，分析在库水下降过程中，滑坡的位移场、塑性区及稳定系数，具有较为明确的物理意义.

4 结果分析

通过ABAQUS软件，采用弹塑性有限元强度折减法对树坪滑坡进行饱和－非饱和流固耦合计算，其计算结果如图5～6所示.

图5 树坪滑坡不同时刻滑体内塑性区分布图

图6 树坪滑坡位移增量云图（单位：m）

由图5可以看出，库水由175m以0.21m／d的速度下降至145m过程中，塑性应变主要发生在滑带的前缘，滑带的后缘和中部基本上没有产生塑性应变.最大的塑性应变发生在滑带的前缘，对应着滑坡的剪出口.而且，随着库水不断下降，塑性应变有扩大的趋势.根据地质资料分析，在库水作用下，树坪滑坡变形破坏方式主要是牵引式解体，慢速滑移.从塑性区的分布来看，计算结果基本上和地质上的认识和判断是较符合的.树坪滑坡之所以未产生整体性滑移，主要可能由于滑带中部提供强有力的抗滑力，如果滑带中部产生较大的塑性变形，则会导致树坪滑坡滑带塑性区全部贯通，进而导致其整体性失稳.

从位移云图分析得出，树坪滑坡的最大位移发生在滑体1和滑体2前缘位置.从滑体前缘至后缘，位移逐渐减小，后缘和中部基本无位移，与塑性区的分布吻合，此时滑坡的位移主要由前缘带动.当采用特征部位的位移拐点判据得到的安全系数为1.112，此时滑坡处于较稳定状态.

5 结 论

1）通过数值计算，在库水以0.21m／d的速度由175m降至145m过程中，树坪滑体1和滑体2前缘出现较大位移，滑体后缘和中部基本上未产生位移，计算结果表明，树坪滑坡由于前缘的牵引，慢速滑移.这与地质资料较符合.

2）通过数值模拟，三峡水库在正常泄水的过程中，树坪滑坡塑性区和位移出现的位置，表明在该工况下，树坪滑坡的稳定性受到了一定的影响，引起滑体前缘变形，但计算稳定性系数为1.112，与监测资料较吻合.该计算结果为三峡水库的调度及灾害防治提供一定的参考.

［1]朱冬林，任光明，聂德新，等.库水位变化下对水库滑坡稳定性影响的预测［J].水文地质工程地质，2002，3（2）：6－9.

［2]黄润秋，戚国庆.非饱和渗流基质吸力对边坡稳定性的影响［J].工程地质学报，2002，10（4）：343－348.

［3]Wang Fawu，Zhang Yemin，et al.Deformation Features of Shuping Landslide Caused by Water Level Changes［J].Chinese Journal of Rock Mechanics and Engineering，2007，26（3）：509－517.

［4]张 欣.基于ABAQUS流固耦合理论的库岸滑坡稳定性分析［D].济南：山东大学，2005.

［5]费 康，张建伟.ABAQUS在岩土工程中的应用［M].北京：中国水利水电出版社，2010.

［6]Dawson E M，Roth W H，Drescher A.Slope Stability Analysis by Strength Reduction［J].Geotechnique，1999，49（3）：835－840.

［7]吕擎峰.土坡稳定分析方法研究［D].南京：河海大学，2005.

［8]宋二祥.土工结构安全系数的有限元计算［J].岩土工程学，1997，19（2）：1－7.

［9]栾茂田，武亚军，年廷凯.强度折减有限元方法中边坡失稳的区判据及其应用［J].防灾减灾工程学报，2003，23（3）：1－8.