库水位变化和降雨作用下三门洞滑坡流固耦合分析

徐兴倩 屈新 张新启 窦思军 陈肖肖

摘要：目前对三门洞滑坡的已有研究都重点关注稳定系数和变形特征，缺乏对孔隙水压力及应力演化特征的系统性研究。采用Abaqus建立滑坡的三维计算模型，推导库水位变化和降雨入渗的耦合边界条件，提出基于莫尔一库伦和最大拉应力准则相结合的复合准则，模拟滑坡在库水位变化和降雨共同作用下的变形演化过程。系统探讨了滑坡的孔隙水压力分布、位移场、应力场和塑性区分布特征，并综合数值模拟结果分析了库水位变化和降雨入渗条件下的边坡稳定性问题。结果表明：在强降雨和库水位缓慢下降综合作用下，三门洞滑坡在滑坡中后部和前缘分别出现局部的张拉和剪切破坏，滑体内渗流场发生较大变化，导致滑体内孔隙水压力、应力、位移和等效塑性应变整体增大，但并未出现塑性破坏区，边坡整体處于稳定状态。

关键词：滑坡稳定性;流固耦合;库水位变化;降雨入渗;三门洞滑坡;三峡库区

中图法分类号：P642.22 文献标志码：A DOI：10.15974/j.cnki.slsdkb.2020.04.005

1研究背景

三峡库区的地质环境因素复杂，库水位变化大，降雨频繁，导致库区边坡极易出现不稳定的情况，给当地人民的生命和财产安全造成了威胁。已有很多学者从库水位变化或降雨角度对边坡稳定性开展了大量研究工作。Ng和Shi对暴雨条件下土坡的暂态渗流场和稳定性进行了分析，探讨了正负孔隙水压力演化规律，并发现降雨入渗导致坡体内部负孔隙水压力减小，边坡稳定系数降低。刘才华等探讨了库水位上升边坡失稳的诱发机理。刘金龙等采用有限元模拟了库水位变化对边坡稳定性的影响。岳庆河等分析了库水位回落对土石坝边坡稳定性的影响。吴海燕等研究了降雨人渗导致的黄土边坡破坏面形成过程及滑动机理。齐云龙等讨论了地震和降雨共同作用下边坡稳定性，针对性地提出了边坡防治建议。娄一青等综合有限元数值模拟分析结果得出了降雨条件下边坡地下水的渗流特征。尚敏等采用Geo-studio模拟了白水河滑坡在降雨入渗作用下的渗流场，探讨了地表变形幅度与降雨强度之间的关系。上述研究揭示了降雨或者库水位变化对边坡稳定性的影响规律，但都是单一地从库水位变化或降雨入渗的角度分析边坡稳定性。莫伟伟研究了库水位变化和降雨条件下滑坡水岩流固耦合的作用机理，并基于BP神经网络技术和遗传算法获取了库水位变化和降雨条件下岩土体的物理力学参数。帅红岩采用有限元模拟了三峡库区晒盐坝滑坡在库水位和降雨共同作用下的滑坡渗流场变化情况，并基于非饱和土抗剪强度理论和Fredlund的极限平衡法理论，分析了滑坡的稳定性。秦洪斌提出了库水和降雨共同作用下复活型滑坡的复活判据，并基于3个具体滑坡验证了该判据的有效性。上述研究探讨了库水位变化和降雨共同作用下的边坡变形稳定性问题，为本文的研究工作提供了理论支撑。

已有很多学者对三峡库区的典型滑坡——三门洞滑坡进行了研究。薛聪聪等综合考虑库水位和降雨渗流的作用，分析了不同工况下三门洞滑坡的稳定系数，并得出了定性结论。易庆林等n。基于三门洞滑坡的地质勘探数据，综合分析了该滑坡体的动态变形特征和主要影响因素。上述研究均探讨了三门洞滑坡的变形及稳定系数，但很少对滑坡体的孔隙水压力分布、位移场、应力场和塑性区等进行系统性研究。

基于此，本文拟以三门洞滑坡为研究对象，采用有限元软件Abaqus模拟其在水库水位变化和降雨共同作用下的变形演化过程，获取该滑坡的孔隙水压力分布、位移场、应力场和塑性区，并综合数值模拟结果分析库水位变化和降雨对边坡稳定性的影响规律。

2三门洞滑坡地貌及有限元网格模型

三门洞滑坡位于长江支流青干河右岸，距河口约8.5km。如图1所示，滑坡前缘较缓，中部和后缘稍陡，西高东低。滑坡呈凹形坡，滑体前缘高程125m，后缘高程350m，呈圆弧状，以基岩为界;左侧以基岩山脊为界，右侧以陡坎临空面为界，平均坡度15°，均宽300m。滑体长830m，面积24.9万m²，平均厚度22m，体积548万m³，主滑方向60°。

根据滑坡的工程地质勘察资料，分析了滑坡地质结构特征，建立了三维有限元计算网格模型，如图2所示。该模型沿青干河水流方向为1080m，垂直青干河水流方向为1468m，底面高程为-142m，计算区域包含滑体、滑带和基岩，对其划分网格后，共生成94090个六面体单元，105644个节点。

在计算模型中考虑了滑坡体的3种典型土体：滑体、滑带和滑床（见图3），相应的物理参数如表1所示。

3岩土体的屈服准则

在库水位变化和降雨共同作用下，滑坡前缘和中部主要发生剪切破坏，后缘出现张拉裂隙，发生张拉破坏。岩土体发生剪切破坏的屈服准则常采用广义米赛斯准则（D-P准则）和莫尔一库仑准则（M-C准则）。D-P准则在主应力空间上的屈服面为一圆锥面，在π平面上为圆形，不存在尖顶处的数值计算问题。M-C准则的屈服面为不规则的六角形截面的角椎体表面，在π平面上为不等角六边形，存在尖顶和菱角。运用现有技术，已能较好处理尖顶，同时M-C准则在边坡工程中得到广泛应用。因此，本文选用M-C准则来判定岩土体是否发生剪切破坏。岩土体发生张拉破坏常采用的屈服准则为最大拉应力准则。

复合准则在（σ1，σ3）平面上的描述如图4所示。图4中，A点到B点为M-C屈服准则f^s=0，其中f^s可以表示为：

4边界条件

采用饱和非恒定渗流与应力耦合理论进行三门洞滑坡流固耦合计算分析，渗流边界条件和应力边界条件的确定方法如下。

4.1渗流边界条件

（1）渗流初始边界条件。基于工程勘察资料确定坡体初始地下水位，然后采用Abaqus确定坡体内初始渗透场的孔隙水压力分布。如表1所示，采用不同岩土体渗透系数对应的渗流速度来模拟计算模型四周边界的透水条件。

（2）非恒定渗流边界条件。降雨引起坡面入渗及产流条件由降雨强度和坡面岩土体的人渗率决定。当降雨强度小于坡面岩土体的入渗率时，降雨产生的水流全部入渗至坡体;当降雨强度大于坡面岩土体的入渗率时，坡面上将产生径流。基于降雨在坡面人渗产流的主要特点，坡面入渗产流的边界条件可通过如下方法确定：

式中，v为雨水在坡面的渗流流速（由坡面岩土体的渗透系数决定）;q为降雨强度。

本文采用第一类边界条件来模拟库水位变化时坡体表面的水头变化情况。数值计算模型底边为不透水边界，四周及坡面不接触降雨或库水的边界为自由透水面。

4.2应力边界条件

数值计算模型的底部边界为法向约束，前后两侧及左右两侧均采用法向约束，坡体表面为自由边界。

最终建立的渗流边界和应力邊界如图5所示。

5滑坡变形现状

三门洞滑坡坡体上共布设ZG360-ZG365等6个GPS监测点，监测时间为2006年9月至2009年12月，监测结果如图6-7所示。

从图6可看出，三峡库区库水位每年均会经历不同程度的涨落，基本与各监测点的位移变形变化相互对应;库水位上涨时，滑坡变形变化不明显，而库水位下降时，滑坡变形变化大，且每次库水位下降均会导致滑坡累计位移曲线上扬，说明其对滑坡稳定性会产生很大影响。

从图7可看出，ZG360点累积水平位移为382.2mm，ZG361点累积水平位移为718.7mm，ZG362点累积水平位移为1205.3mm，ZG363点累积水平位移为433.7mm，ZG364点累积水平位移为413.9mm，ZG365点累积水平位移为1557.6mm。

基于监测数据资料，查明三峡库区水位调度过程（见图8）及当地50a内每月最大的连续3d降雨量（见图9）。

6计算工况及结果分析

6.1计算工况

为充分表现库水位变化和降雨的综合作用效应，将对比分析普通工况（工况1）和极端工况（工况2）条件下滑坡孔隙水压力分布、应力场、位移场和塑性区分布特征。其中，工况1为水库蓄水至175m水位;工况2为库水位从175.0m降至145.0m，并遇上50a一遇暴雨，荷载组合为自重和地表荷载。

6.2结果分析

（1）孔隙水压力计算结果及分析。水库蓄水至175m后，孔隙水压力如图10所示：滑坡前缘高程约140m处出现最大孔压，数值约为0.6MPa;滑坡后缘底部孔隙水压力主要受地下水位控制，最大孔压值为2.60MPa（见图10）。

随着库水位由175m降至145m，地下水位也相应下降但存在滞后，这是由于滑坡前缘孔隙水压力差消散需要一段时间;滑坡后缘地下水基本不受库水位变化影响，这是由于滑坡后缘离库水位变化带的距离相对较远，但受强降雨影响，滑坡后缘底部最大孔压值由2.60MPa增至2.83MPa（见图11）。

（2）应力计算结果及分析。水库蓄水至175m水位后大部分区域拉、压应力分布变化较小。但是，由于库水作用滑体中下部产生较大变化，最大主应力的最大值为-0.20MPa，最小主应力的最小值为-0.71MPa（图12-13）。

当库水位由175m降至145m，并遇上50a一遇暴雨时，滑坡整体应力增加，最大、最小主应力场整体变化较大。滑体中最大主应力的最大值为-1.51MPa，最小主应力的最小值为-2.60MPa;滑体中上部最大主应力的最大值为-1.10MPa，最小主应力的最小值为-3.14MPa。该工况应力场变化较大，主要原因是库水下降和降雨的共同作用引起渗流场的改变，而渗流场与应力场相互耦合作用，导致应力场发生改变（见图14～15）。

（3）位移计算结果及分析。水库蓄水至175m后，滑体下部发生较大水平位移，最大值为0.015m;滑体中部水平位移相对较小，最大值为0.004m;滑体上部水平位移较小，最大值为0.001m。这主要是由于库水对滑体产生浮托作用，使滑体抗滑力降低，导致滑体产生水平位移。在滑体上缘与滑床交界的部位出现一些突变位移，这是因为在此区域滑体与滑床岩土材料属性和地形均产生了不均匀特性（见图16）。此时，滑体上部发生较大向下的垂直位移，最大值为0.040m;滑体中部垂直位移相对较小，最大值为0.029m;滑体前缘产生向上的垂直位移，最大值为0.006m。这主要是由于库水对滑体前缘的浮托作用，使滑体产生向上的垂直位移（见图17）。

当库水位由175m降至145m，并遇50a一遇暴雨时，滑体下部发生较大水平位移，最大位移值为0.049m;滑体中部水平位移相对较小，最大位移为0.026m;滑体上部水平位移相对最小，最大位移为0.0013m。该工况水平位移较大，主要原因是库水下降使滑体内产生向外的渗透力，同时降雨作用使向外的渗透力增大，使滑体下滑力增大，导致滑体下部水平位移增大（见图18）。此时，滑体上部发生较大向下的垂直位移，最大值为0.102m;滑体中部垂直位移相对较小，最大值为0.063m;滑体下部垂直位移相对最小，最大值为0.052m。该工况的垂直位移较大，主要原因是库水下降引起滑体向下滑移，同时在自重和降雨作用下发生沉降，使滑体产生向下的垂直位移增大。滑体前缘产生向上的垂直位移，最大值为0.021m，主要原因是：库水对滑体前缘的浮托作用，使滑体向上的垂直位移增大（见图19）。

（4）塑性区计算结果及分析。工况1和2计算过程中滑坡未出现塑性破坏区（见图20～21），说明滑坡整体处于较稳定状态。

（5）边坡整体稳定性分析。对比分析两种工况的数值计算模拟结果和监测数据后发现，数值模拟结果与监测数据匹配度较高，证明了本文模拟结果的有效性。在强降雨和库水位缓慢下降综合作用下，三门洞滑坡在滑坡中后部和前缘分别出现局部的张拉和剪切破坏，滑体内渗流场发生很大变化，导致滑体内孔隙水压力、应力、位移和等效塑性应变整体增大，但并未出现塑性破坏区，边坡整体处于稳定状态。

7结语

本文以三峡库区三门洞滑坡为研究对象，建立了基于Abaqus的滑坡三维计算模型、库水位变化和降雨人渗的耦合边界条件、莫尔一库伦和最大拉应力准则相结合的复合准则，分析了滑坡在水库蓄水、库水位缓慢降落及强降雨过程中的流固耦合及三维稳定性问题。

三门洞滑坡的变形破坏是在三峡水库运行期水位波动条件下，耦合降雨过程长期发展的结果。由于篇幅有限，本文只考虑了水库蓄水、库水位缓慢回落和强降雨等工况。若要得到边坡变形破坏发展的真实趋势，须考虑库水位上升、下降、再上升、下降等反复循环的整个过程。