库水位上升对茨哈峡4#倾倒体的稳定性研究

杨 敏，李宏儒，李 宁，李国锋，安晓凡

(西安理工大学 岩土工程研究所，陕西 西安 710048)

库岸滑坡是伴随着水电工程的建设而加剧呈现的一种重大工程问题，库岸滑坡失稳，滑坡体高速滑入水库，会造成巨大的涌浪，不仅危害水库本身，影响水库的正常运营，而且还会带来严重的次生灾害，甚至直接危及大坝的安全，给人类的生命财产带来重大威胁[1].许多专家学者结合各自参与的具体工程开展了大量有针对性的研究工作.刘才华等[2]以湖北省秭归县三峡库区的千将坪高速滑坡为例，论述了库水位上升有可能诱发边坡失稳破坏机理.莫伟伟等[3]根据水岩相互作用机理，分析了库水涨落对滑坡岩土体的影响.王新刚等[4]以三峡库区马家沟I号滑坡为例，建立真三维模型，研究了应力-渗流耦合作用下抗滑桩加固库区滑坡位移和受力特征，探讨了滑坡-抗滑桩相互作用体系的防治效果.马淑芝等[5]以赵树岭滑坡为例，运用水-岩耦合三维有限元数值模拟对不同库水位变动状态下的地下水水位变化和滑坡稳定性进行了研究.周永强等[6]以付家坪子高陡滑坡为例，分析了库水位不同升速和降速、降雨不同强度以及库水位变化与降雨联合作用工况下滑坡的稳定性.孙冬梅等[7]通过建立水一气二相流模型，研究了库水位下降的岸坡非稳定渗流场.刘钊等[8]利用有限元法研究了水位骤降时影响上游坝坡稳定性的因素.杨超[9]以寨沟村滑坡为例，进行水库型滑坡的稳定性分析的研究.廖红建[10]以三峡水利工程为背景，数值模拟分析出库区降水速度和渗透系数与边坡稳定性之间的变化规律.李壮壮[11]以三峡库区奉节县龙潭滑坡为研究对象，运用Geo-studio 软件模拟库水位从175 m 降至145 m 的滑坡渗流场.上述研究成果在解决具体工程问题方面起到了显著的作用，并对类似边坡问题的解决具有一定的指导意义，但针对茨哈峡4#倾倒体的实际情况，仍需要开展专门的研究工作.

本文研究的茨哈峡4#倾倒体是典型的反倾层状岩质边坡，位于库区内部，整体规模较大，其所在斜坡浅表部岩体发生严重的倾倒变形，倾倒变形主要发生在2 900 m高程以上，有向上游侧倾倒的趋势，岩体拉裂、破碎、解体，稳定性较差，形成危岩，局部发生小规模的垮塌变形.倾倒体在后期正常蓄水时，倾倒体前缘岩体在浸泡作用下处于饱和状态，一方面岸坡的受力条件将发生变化，库水的悬浮效应将减小倾倒体前缘的抗滑能力，另一方面倾倒体前缘岩土体材料的力学参数降低，抗滑能力进一步降低.则倾倒体前缘可能发生拉裂滑移现象，从而牵引中上部高程层状岩体发生倾倒变形.因此可知，4#倾倒体的失稳模式为牵引式倾倒破坏.因此，对该倾倒体的稳定性进行分析评价是十分必要的.

1 工程概况

1.1 试验模型方案及力学参数标定

茨哈峡水电站4#倾倒体位于大坝上游、黄河左岸，倾倒体前缘位于黄河水位以上，高程在2 785 m左右.倾倒体后缘位于左岸边坡顶部，高程在3 050～3 100 m之间，整个倾倒体边坡垂直高差在300 m左右，顺黄河方向长近千米.在倾倒体中部冲沟较发育，发育多条规模不等的冲沟[12]，如图1～图2所示.

图1 4#倾倒体地形地貌特征Fig.1 4#Dumping landform characteristic

图2 4#倾倒体表部冲沟Fig.2 The surface gullies of 4# toppling slope

总体来说，斜坡岩体的倾倒水平深度在70～160 m之间.并且随着斜坡高程的增加，岩体的倾倒水平深度也逐渐增大.某设计院提供了4#倾倒体的地质平面横1至横5共五个剖面图.其中三个剖面的钻孔地质具体如下：

①横2剖面，PD17号平洞岩体无倾倒变形，PD65号平洞岩体倾倒深度114 m.整个剖面上倾倒体呈下窄上宽度特征，随着高程的增大，斜坡岩体倾倒深度增大，在斜坡顶部，倾倒水平深度接近150 m；

②横3剖面，PD42和PD66号平洞岩体倾倒深度分别是90 m和80 m.受地形影响，在高程2 860～2 920 m之间为一个突出的山梁，导致这一段斜坡岩体倾倒深度较大，水平深度超过100 m.在斜坡上部，高程2 950 m以上为冲沟地形，导致边坡岩体倾倒深度相应的变小，倾倒深度在80～90 m之间.

③横5剖面，PD52平洞倾倒深度135 m，ZK11钻孔倾倒深度55 m，PD53平洞岩体无倾倒变形.在高程2 920m 以下，随着高程的增加，斜坡岩体倾倒深度逐渐增大；在高程2 920 m以上，岩体倾倒深度基本上在150～160 m之间.

根据《黄河茨哈峡水电站可行性研究——枢纽布置选择专题报告》，以平洞和钻孔中岩体倾倒深度为基础，边坡岩体的倾倒水平深度在70～160 m之间.并且随着斜坡高程的增加，岩体的倾倒水平深度也逐渐增大.

4#倾倒体位于库区内部，其所在斜坡浅表部岩体发生严重的倾倒变形，岩体拉裂、破碎、解体，稳定性较差，形成危岩；边坡中上部，倾倒后岩层倾角普遍在30～40°，边坡下部，倾倒后岩层倾角一般在50°左右，变形程度相对稍低.边坡在后期库水位上升等作用下，可能存在倾倒失稳的可能.倾倒方向为河谷偏上游.其稳定性对后期的施工与运营影响大，故需对该边坡进行三维仿真分析计算，研究其稳定性及软弱区，为后期边坡开挖与加固分析做好铺垫，确保施工运营期的安全与稳定.

2 数值模型和计算参数

2.1 数值模型

根据地质分析，边坡地质模型可概化为：强风化、弱风化、下部基岩.某设计院提供的4#倾倒体的地质平面横1至横5五个剖面图.在此基础上依据地质概况及岩层岩体及结构面分区分级，建立三维地质概化模型如图3所示.

图3 4#倾倒体三维地质概化模型Fig.3 3D generalization geological model of 4# dump

三维地质模型范围平行于滑体主轴方向为750 m(X轴方向)，边坡后缘垂直于滑体主轴方向为390 m(Y轴方向)，边坡前缘垂直于滑体主轴方向为645 m(Y轴方向)，垂直方向最高为403～508 m(Z轴方向，高程为0～508 m)，模型4个侧壁施加法向位移约束，底部边界施加位移全约束，该区域构造应力场较小，不予考虑，只计算重力场.坡面在自然状态下为自由边界，水库蓄水后施加库水压力.采用国际通用有限元软件ANSYS建立计算模型，计算模型如图5 所示，共划分44 834个单元，8 880个节点，最后导入FLAC3D三维有限差分程序软件进行分析.

图4 三维计算模型及监测点示意图Fig.4 Three-dimensional calculation modeland key schematic

图5 三维计算模型分离示意图Fig.5 Three dimensional calculation separation model

2.2 计算参数

根据地质报告，划分坡表强风化强卸荷土体为Ⅳ2级岩体；在断层影响带内岩体局部拉裂变形属弱风化弱卸荷岩体为Ⅲ2级；在影响带以里岩体属微未风化未卸荷基岩体完整性差或者较完整岩体为Ⅱ级；对于Ⅱ级岩体考虑到地下水的作用又分为(1)和(2)两区.对于结构面，倾倒断层结构面为Ⅱ级，底部剪切面为Ⅳ级，其它小断层裂隙可忽略不计.

参数选取时，水位以下的岩体参数采用饱水状态下岩土体参数，水位以上的岩体参数采用岩土体相应的天然强度参数.其中地下水位下按照饱和容重计算，并考虑岩体软化，根据地质资料软化系数平均取0.72.库水位抬升阶段同时包括湿化软化变形和渗流变形，对岩体的弱化效应通过软化系数来考虑，即：当库水位升高至某一高程时，对该高程以下岩体进行弱化，弱化程度和范围随时间的变化而增加，本次计算所取时间为6个月，计算边坡的位移及塑性区；通过耦合求解建立孔压场的分布，来计算边坡的渗流场和应力场，参照相关论文和地质资料软化系数的范围取0.5～0.7折减[13].

有限差分程序软件的网格模型按不同的风化程度、地下水、滑体及结构面等因素共同作用进行分体，总共分为13个体积，并在4个结构面的分界面处生成界面单元，水面上岩土体赋予真实岩土体物理参数(见表1)，岩土体采用莫尔-库仑模型.法向刚度和剪切刚度，取周围“最硬”相邻区域的等效刚度10倍：Kn和Ks可采用如下公式计算[14]：

Kn=Ks=10max[(K+4G/3)/Δzmin]

(1)

式中：Kn和Ks分别为法向刚度和切向刚度(MPa)；Δzmin为接触面法向上连接区域的最小尺寸，无量纲.

在FLAC3D中，岩体的变形参数采用的是体积模量K和剪切模量G，其与岩体的弹性模量E和泊松比μ的关系如下：

(2)

(3)

表1 结构面参数选取表

表2 边坡岩体参数选取表

2.3 计算方法

边坡的稳定性分析与评价是边坡研究的核心问题，不同的变形破坏模式存在不同的破坏类型，所以选择合理的边坡稳定性进行分析与评价方法是非常必要的.

目前在采用强度折减弹塑性有限元数值方法进行边坡稳定性分析时，尚缺乏统一的失稳评判标准，因此由不同的失稳判据所估算的边坡总体安全系数可能有所不同.通常采用的边坡失稳判据有：边坡特征点的位移随折减系数的变化关系曲线出现突变，有限元数值无法在给定的迭代次数内达到收敛，或大于某一量值的广义剪应变区域，或广义塑性应变区域在边坡中上下贯通.研究表明，对不同强度折减系数情况下的位移特性或塑性区特性进行对比分析，可以减小非确定性因素对安全系数的影响[15]，因此，本文所研究的茨哈峡4#倾倒体地边坡三维稳定性分析是在ANSYS中建立三维模型后导入FLAC3D中进行计算，联合采用关键点处的位移是否突变和塑性区是否贯通作为边坡的失稳判据.

本文采用FISH自编的强度折减法，其中折减参数为C0和tanφ0同比折减，即折减系数Ki可表示为

(4)

本文根据某设计院提供的库水位资料，拟定了四种计算工况，分别为自然库水位2 760 m (125 m)工况、水位升至2 840 m (205 m)工况、水位升至2 935 m (300 m)工况、水位升至2 990 m (355 m)，其中，括号里的高度是水位距离模型底部的高度.

3 4#倾倒体三维结果分析

3.1 自然库水位至2 760 m(125 m)计算结果

为了节省篇幅，本文仅给出不同折减系数下各关键点的总位移曲线如下图6中图(a)所示，从左到右依次为上游、中游和下游，图6中图(b)是中游各关键点总位移随折减系数变化曲线.计算云图仅给出折减系数为1.00时对应的塑性区云图和剪应变增量与速度云图如下图7(a)～图7(b)所示.

由图6(a)～图6(b)可知：中下游位移、速度变化相对较大，随着折减系数的增大，边坡呈现出向上游滑移.中游关键点的位移、速度均为中高高程变化相对较大，其中位移随着折减系数的增大而增大，折减系数在1.78时位移发生较为明显的转折，随后急剧增大，而中游关键点的速度在折减系数1.83时有明显的转折.

图6 各关键点位移速度与折减系数关系曲线Fig.3 Displacement speed of the key point and thereduction factor relationship curve

由图7(a)～图7(b)可知：当折减系数为1.0时，塑性区范围在结构面处比较集中，剪应变增量在下游和结构面处比较明显，变形速度主要表现在边坡下游处.随着折减系数的增大，塑性区范围逐渐增大，剪应变增量与速度逐渐在中游及上游位置增大，在结构面附近增大的比较明显，中游在结构面附近的拉应力区增大，在下游的高高程位移区相对逐渐增大.

根据计算不收敛判据确定安全系数为1.83，根据位移突变判据确定安全系数为1.78.因此，该方案安全系数确定为1.78.

3.2 库水位升至2 840 m (205 m)计算结果

图7 塑性区及剪应变增量与速度云图Fig.7 The plastic zone, the shear strain incrementand velocity contours

图8 各关键点位移速度与折减系数关系曲线Fig.8 Displacement speed of the keypoint andthe reduction factor relationship curve

为了节省篇幅，本文仅给出不同折减系数下各关键点的总位移曲线如下图8中图(a)所示，从左到右依次为上游、中游和下游，图8中图(b)是中游各关键点总位移随折减系数变化曲线.计算云图仅给出折减系数为1.00时对应的塑性区云图和剪应变增量与速度云图如下图9(a)～图9(b)所示.

由图8(a)～图8(b)可知：上中游位移、速度变化相对较大，上游整体向下游变形而中下游向上游变形.位移和速度都是从上部监测点到下部监测点依次减小，可以说明此时坡体有发生倾倒的可能.位移和速度均随着折减系数的增大而增大，且折减系数在1.51时位移发生较为明显的转折，随后急剧增大.

图9 塑性区及剪应变增量与速度云图Fig.9 The plastic zone, the shear strain incrementand velocity contours

由图9(a)～图9(b)可知：当折减系数为1.0时，塑性区范围也是在结构面处比较集中，相对于自然库水位下的塑性区，在该工况下的塑性区范围增大，并且主要增大在结构面和边坡前缘附近.剪应变增量与在下游和结构面处比较大.随着折减系数的增大，塑性区范围逐渐增大，剪应变增量与速度也逐渐增大，在结构面附近增大的比较明显.在结构面附近出现了拉应力区，并且X方向的位移值较大.坡体后缘部分的总位移值较大.

根据计算该方案根据位移突变判据确定安全系数为1.51.

3.3 库水位升至2 935 m(300 m)计算结果

为了节省篇幅，本文仅给出不同折减系数下各关键点的总位移曲线如下图10(a)所示，从左到右依次为上游、中游和下游，图10(b)是中游各关键点总位移随折减系数变化曲线.计算云图仅给出折减系数为1.00时对应的塑性区云图和剪应变增量与速度云图如下图11(a)～图11(b)所示.

由图10(a)～图10(b)可知：上中游位移、速度变化相对较大，上游整体向下游变形而中下游向上游变形.位移和速度都是从上部监测点到下部监测点依次减小，可以说明此时坡体有发生倾倒的可能.位移和速度均随着折减系数的增大而增大，其中位移在折减系数为1.27时发生较为明显的转折，随后急剧增大，而速度在折减系数为1.24时发生了明显转折.

图10 各关键点位移速度与折减系数关系曲线Fig.10 Displacement speed of the key point andthe reduction factor relationship curve

由图11(a)～图11(b)可知：当折减系数为1.0时，塑性区范围明显比库水位升至2 840 m时的范围增加了，主要在结构面处比较集中，剪应变增量与速度也增大了.随着折减系数的逐渐增大，塑性区范围和剪应变增量与速度也分别逐渐增大，而且在结构面附近增大的比较明显.出现了拉应力区，并且X方向的位移值较大.坡体后缘部分和上游的总位移较大.

根据计算该方案根据位移突变判据确定安全系数为1.27.

3.4 库水位升至2 990 m(355 m)计算结果

为了节省篇幅，本文仅给出不同折减系数下各关键点的总位移曲线如下图12(a)所示，从左到右依次为上游、中游和下游，图12(b)是中游各关键点总位移随折减系数变化曲线.计算云图仅给出折减系数为1.00时对应的塑性区云图和剪应变增量与速度云图，见图13(a)～图13(b).

图11 塑性区及剪应变增量与速度云图Fig.11 The plastic zone, the shear strain incrementand velocity contours

图12 各关键点位移速度与折减系数关系曲线Fig.12 Displacement speed of the key doint and thereduction factor relationship curve

由图12(a)～图12(b)可知：上中游位移、速度变化相对较大，上游整体向下游变形而中下游向上游变形.位移和速度都是从上部监测点到下部监测点依次减小，可以说明此时坡体有发生倾倒的可能.位移和速度均随着折减系数的增大而增大，其中位移在折减系数为1.时发生较为明显的转折，随后急剧增大.

图13 塑性区及剪应变增量与速度云图Fig.13 the plastic zone， the shear strain incrementand velocity contours

由图13(a)～图13(b)可知：当折减系数为1.0时，整个坡体都发生塑性变形，剪应变增量与速度相对于库水位升至2 925 m时增大了.随着折减系数的增大，塑性区范围逐渐增大，剪应变增量与速度也逐渐增大，在结构面附近增大的比较明显.在结构面附近出现了拉应力区，并且X方向的位移值较大.坡体后缘部分和上游的总位移较大.在坡体的下表面孔压增大.

根据计算该方案根据位移突变判据确定安全系数为1.12.

综上所述，本文针对4#倾倒体边坡三维稳定性分析，采用强度折减法并采用结构面和结构面上部岩体折减的方案进行计算分析.结合位移突变及塑性区贯通等判据确定综合安全系数，对于水位上升工况采用逐步增加折减系数的方法，根据位移突变判据确定安全系数.各工况下的安全系数如下表3所示.

表3 各工况下安全系数汇总表

4 结论

茨哈峡4#倾倒位于茨哈峡水电站上游，整体规模较大，其所在斜坡浅表部岩体发生严重的倾倒变形.本文采用FLAC3D有限差分软件对该倾倒体进行数值模拟，取得了以下结论:

(1)4#倾倒体稳定性受库水位升降影响较为明显，在库水位上升过程中，倾倒稳定性逐渐降低，水位处于2 990 m稳态时，倾倒稳定系数最小，其值为1.12.

(2)边坡各关键点位移和速度随着高程的降低依次减小，中高部高程岩体有明显的向坡外倾倒变形趋势，拉应力和塑性区的均沿结构面分布.

(3)随着折减系数的增大，塑性区范围逐渐增大，剪应变增量与速度也逐渐增大，在结构面附近增大的比较明显.

(4)鉴于在库水位升降作用条件下的倾倒稳定性最低值为1.12，小于设计安全系数，建议对倾倒体采取监测控制措施，为倾倒的防治提供预警机制.