库水位变化速率对非饱和边坡稳定性影响分析研究

翁新海，王家鹏，翁郁馨

（1.浙江省水利水电勘测设计院，浙江 杭州 310002；2.浙江大地岩土勘察有限责任公司，浙江 杭州 310002；3.东南大学研究生院，江苏 南京 210096）

1 概述

库水位变化是水库区滑坡失稳关键因素，引起破坏性极大，造成的损失巨大，如著名瓦依昂水库库区滑坡造成生命财产巨大损失。因此，研究库水位变化对非饱和边坡的稳定性具有重要的现实意义。以西北部库区边坡为例，使用Geo-slope 有限元计算软件，进行渗流分析及稳定性计算。分析不同库水位上升和下降幅度对非饱和边坡渗流场和稳定性的影响规律。

2 饱和—非饱和渗流研究现状

土中水的渗流研究历史悠久。早在1856 年，法国工程师达西Darcy 提出了线性渗流的达西定律。1889年，俄国的茹可夫斯基首先推导了渗流的微分方程。1922 年，俄国学者巴普洛夫斯基提出了求解渗流场的电模拟法。由于渗流的微分方程在复杂边界条件下很难得到其解析解，所以人们力图用数值算法解决它。1910 年，英国学者理查森首先提出了有限元法。在20世纪70年代之后，由于计算机及计算技术的迅速发展，人们开始广泛应用有限元法、边界元法和许多其他计算方法计算解决渗流问题。

Frediund 对饱和—非饱和土坡渗流问题做了较完整的论述，把非饱和土壤水运动理论与非饱和土固结理论相结合，得到了饱和—非饱和渗流控制方程，并运用二维有限元方法对复杂地下水流动系统的几个暂态渗流实例问题进行了数值模拟等[1]。

吴梦喜、高莲士[2]对饱和—非饱和土体非稳定渗流做了数值分析，对一般的非饱和渗流有限元计算方法加以改进，以消除非饱和渗流数值计算存在的数值弥散现象。同时还提出了一种逸出面处理的新方法，并给出了非饱和非稳定渗流计算的实例。

针对水位降落所引起的边坡内渗流场的研究，目前采用最多的是以非饱和土理论为基础，采用有限元法进行分析。时卫民等[3]以布西涅斯克非稳定渗流基本微分方程为基础，给出了坡前水位降落情况下的坡内浸润线的计算方程,并进一步利用多项式拟合原理对其进行了简化，使其变的易于应用。

朱军等[4]对常用的饱和—非饱和三维多孔介质非稳定渗流的有限元算法进行了改进，提出以高斯点为研究对象。在单元中不同的高斯点有不同的压力水头，进而可得到不同的渗透系数k及水体积变化系数，这样修正得的渗透矩阵和贮水降雨条件下饱和—非饱和土坡的稳定分析矩阵比一般方法更合理。

张培文等[5]在进行饱和—非饱和渗流模拟时，克服了传统降雨边界处理方法的缺陷，考虑了降雨过程中地表边界的条件转化，为以后的径流渗流耦合的模拟打下了良好的基础，为降雨诱发滑坡研究提供了定量的分析手段。

3 饱和—非饱和渗流问题的数学描述及边界条件

求解饱和—非饱和渗流问题，需要建立其控制方程和边界条件。在饱和—非饱和渗流分析中，定解条件方程表达式为：

式中：Γ1——水头分布规律已知的边界；

H1——边界水头，式（1）称为第一类边界条件；

Γ2——流量清况已知的边界；

qn——单位时间边界法向流量；

Γ3——渗出面边界；

z——渗出面节点坐标，式（4）可以归为第一类边界条件；

H0——初始渗流场水头，式（5）为初始条件。

在边界节点处规定水头为第一类边界条件，而规定通过边界的流量称为第二类边界条件。正的节点流量表示节点处有入渗，负节点流量表示该节点处有蒸发或蒸腾，当通过边界的流量为零时，表示不透水边界。

4 工程实例分析

4.1 边坡概况

边坡分布于南侧分水岭一带，边坡分布高程1080～1450m，在平面上似椭圆形，其长轴方向为NW向，长约900m，中上部最宽处约694m。边坡后缘地形较陡峻，坡度40°～60°，中部地形较平缓，坡度25°～30°，前缘地形较陡峻，坡度60°，总体形态上呈上陡中缓下陡的折线状，总体方量约3100×104m3。

边坡表部，上部树木较粗较大，有的树龄已达数百年，但未见马刀树和醉汉林，表部未见裂缝，说明该边坡形成时间较久远，现状处于稳定状态。

边坡上游边界发育有干沟近代冲沟，顶部宽，下部窄，沟底纵坡覆盖层薄，基岩面有起伏，而成台阶状分布。下游边界也为干沟，沟底无水，仅长时间降雨时才有水，此外边坡内有一条顺坡向切割深度3～15m 宽，宽20～30m小冲沟，无水，长时间降雨都无水。

4.2 数值模拟

对边坡根据地质勘探成果进行有限单元格划分，综合考虑计算精度和计算时长。边坡岩土体对库水位较敏感，同时边坡底部与岩体接触带比较敏感，这些部位也进行了加密。有限元网格模型如图1所示。

图1 有限元网格模型

4.3 岩土体物理力学参数

经现场和室内试验及工程地质类比法，边坡体岩土体物理力学参数见表1。

表1 岩（土）体物理力学参数

4.4 计算工况与结果

水库边坡应该考虑水库蓄水对边坡地下水的影响，地下水水位线上的部分设置容重时应采用天然重度进行计算，而地下水位线以下的部分采用饱和重度进行计算。根据库区水位运行调度资料，水位缓慢上升速度为每天0.5m，经过30d 从枯水位1155m 上升到1170m：水位上升为每天2m，经过7.5d从枯水位1155m上升到1170m；水位缓慢下降速度为每天0.5m，经过30d 由蓄水位1170m 降低到1155m。水位快速下降时为每天2m，经过7.5d 从蓄水位1170m 降低到枯水位1155m。故分别设置以下4种工况进行模拟：工况1：水位上升速率为0.5m/d；工况2：水位上升速率为2m/d；工况3：水位下降速率为0.5m/d；工况4：水位下降速率为2m/d。

4.4.1 库水位上升阶段

首先对边坡在枯水位1155m 处进行稳态渗流分析。枯水位1155m处孔隙水压力分布图随着库区水位上升，地下水位线也逐渐上升，坡内孔隙水压力也随之上升。

经过枯水位1155m上升到蓄水位1170m进行渗流模拟，可以了解到每天滑坡体的渗流状况以及地下水位线，如图可以看到，库区水位上升速度高于坡内水位上升速度。可以看出，水位线呈凹形，即坡内水位线低于库区水位线。主要是因为滑体渗透性小于水位变动速率，滑坡体中水位上升有一个滞后效应，坡内水位变动滞后于坡外水位变动。

通过对边坡进行稳定性分析，以整个边坡为研究对象，指定滑动面,经计算可以得到边坡在水位变动工况下各时段的安全系数，如图2所示。可以看到随着水位上升，边坡稳定性先增强后降低，水位上升速度越快，边坡稳定性增强越快。水位上升到蓄水后，2m/d时边坡的稳定性仍然大于0.5m/d 时的稳定性。这是由于坡外水位高于坡内水位，产生指向滑坡体内的动水压力，有利于边坡稳定。同时到达蓄水位后，坡内水位与坡外水位慢慢持平，此时指向坡内的动水压力逐渐降低，以及上层土体基质吸力降低，边坡稳定性逐渐降低。水位上升0.5m/d 和2m/d 情况下，均是在最高水位处稳定性安全系数最大。由于整个边坡体型较大，水位上升对整个边坡稳定性影响较小。

图2 水位上升整个坡体安全系数图

4.4.2 库水位下降阶段

首先对边坡在蓄水位1170m处进行稳态渗流分析（图3）。随着库区水位下降，地下水位线也逐渐下降，坡内孔隙水压力也随之下降。水位下降到枯水位1155m时孔隙水压力分布如图4所示。

图3 稳态下处于蓄水位1170m坡体孔隙水压力分布图

图4 水位下降到枯水位1155m后孔隙水压力分布图

经过蓄水位1170m下降到枯水位1155m进行渗流模拟，可以了解到每天边坡体的渗流状况以及地下水位线，可以看到，库区水位下降速度高于坡内水位下降速度。可以看出，水位线呈凸形，即坡内水位线高于库区水位线。主要是因为滑体渗透性小于水位变动速率，滑坡体中水位下降有一个滞后效应，坡内水位变动滞后于坡外水位变动。

通过对边坡进行稳定性分析，以整个边坡为研究对象，经计算可以得到边坡在水位变动工况下各时段的安全系数，如图所示。可以看到随着水位下降，边坡稳定性先降低后上升，水位下降速度越快，边坡稳定性降低越快。水位下降到枯水后，水位2m/d下降时边坡的稳定性仍然小于水位0.5m/d 下降时的稳定性。这是由于坡外水位低于坡内水位，产生指向滑坡体外的动水压力，不利于边坡稳定。同时到达枯水位后，坡内水位与坡外水位慢慢持平，此时指向坡外的动水压力逐渐降低，以及上层土体基质吸力的升高，边坡稳定性逐渐升高。水位下降0.5m/d 和2m/d 情况下，均是在最低水位处稳定性安全系数最小。由于整个边坡体型较大，水位下降对整个边坡稳定性影响相对较小。

通过对边坡进行稳定性分析，以局部涉水边坡为研究对象，得到边坡前缘涉水段在水位变动工况下各时段的安全系数，如图5所示。可以看到边坡安全系数与水位变动同边坡整体考虑有相同的规律。局部涉水边坡的特点是在水位变动下，稳定性变化幅度相当大。

图5 水位下降边坡整体安全系数图

5 结论与建议

本水电站堆石坝为2 级建筑物，依据《水电水利工程边坡设计规范》（DL/T5353-2006），枢纽工程区边坡为Ⅱ级边坡，安全系数需满足表2的要求。

表2 水电水利工程边坡涉及安全系数

对该水电站库前非饱和边坡在水位变动及降雨工况下进行模拟分析，对边坡整体考虑，具体安全系数如表3所示。可以得出以下结论，边坡在水位变动均处于稳定状态，水位变动对边坡整体稳定性产生了一定的影响，其中水位骤降对边坡稳定性影响较大，边坡稳定性下降较多，但边坡安全系数均大于1.25，边坡仍处于稳定状态。

表3 边坡在不同工况下的最小安全系数

根据有限元数值模拟分析，库水位按2m/d上升速率，预计8d 左右安全系数达到最大，库水位按0.5m/d上升速率，预计30d左右安全系数达到最大，对边坡稳定有利，因此建议在水库蓄水时可以采用先快后慢运行方式；当水库泄洪时，按0.5m/d 下降速率，预计30d左右安全系数达到最小，按2m/d 下降速率，预计7d左右安全系数达到最小，因此采用先慢后快运行方式，对边坡的稳定最为有利。