库水位波动与降雨作用下的滑坡流固耦合分析

卿 菁， 朱 蕾， 陈 昊， 谢 鹏

(湖北工业大学土木建筑与环境学院， 湖北 武汉 430068)

随着我国工程建设的不断发展，诸多地质灾害问题也日益显著。三峡库区沿岸滑坡与水的影响密不可分，以库水位升降、地下水作用和库区降雨为主要影响因素的滑坡灾害多年来一直都是岩土工程界的重点研究对象。多种因素影响下的滑坡变形破坏是一个较为复杂的演变过程，这些过程里既存在有岩土体物理参数影响且改变周期相对较长的内因，同时也存在水位变化、降雨、人工影响等综合作用且改变周期相对较短的外因[1-3]。根据调查发现，绝大部分滑坡变形破坏起因都是因为地下水位的变化[4-7]。

在信息高速发展的今天，众多学者广泛采用了实验分析、数值模拟等方法对滑坡水的作用进行细致调查和剖析。例如，Rubin[8]针对二维渗流的相关问题，以有限差分法的手段求解了Richards方程的数值理论解；王恭先[9]对三峡地区滑坡的各种类型及其对应的分布的特征行为进行了归纳总结，同时对于滑坡的形成机理与稳定性状态进行了分析；王思敬[10]根据水岩作用及水环境改变的类型将库区滑坡分为水岩作用滑坡与天然形成的滑坡两大类；而张桂荣[11]针对三峡库区对某滑坡分析了水作用下的十种工况，并采用SEEP/W软件模拟进行分析。除此之外，三维有限元方法也日渐成熟，学者王媛[12]提出在Biot理论上，提出了一种假设孔隙水压力与位移量来确定出渗流场并与应力场耦合的方法。后来，陈庆中等[13]探究了渗流场计算方法，建立出新的流固耦合数学模型；刘建军[14]针对地下水与地层结构的相互作用关系，建立地下水流固模型来求解相关渗流问题。

现阶段大多数的滑坡破坏机理研究只进行了二维数值模拟或者是三维流固耦合分析，本文在前人已完成的众多工作的基础上，以秭归县卢家沱滑坡为例，结合二维分析有限元分析软件Geo-Studio数值模拟分析和基于Abaqus建模的三维流固耦合分析，进一步开展水作用对滑坡相关的影响行为研究，对该滑坡进行破坏机理和致灾原因等方面的剖析，旨在为今后开展相关滑坡防治的研究工作者们提供进一步的参考。

1 滑坡地质特征

1.1 基本概况

卢家沱滑坡前缘高程135 m，后缘高程370 m，前缘175～195 m高程带为一平台，临江处为坡度45°的陡坎，公路从滑体320 m高程处穿过。卢家沱滑左右边界以基岩山脊为界，整体坡度约23°。滑体纵向长约540 m，平均宽约280 m，均厚约22 m，面积约1.34×105m2，体积约2.95×106m3。具体滑坡平面情况详见图1。

1.2 物质组成及结构特征

卢家沱滑坡下部岩层组层多为泥岩，在库水位升降及降雨等多方面因素影响下，易产生滑移变形，为动水压力型滑坡。滑体物质分为两部分：上部物质为碎石土，褐黄色或灰褐色，碎石成分为灰岩，粒径一般为30～90 mm，最大250 mm，块石呈中～强风化状态。土质为砂质粘土，土石分布不均匀，可塑，稍密-密实，土石比为6∶4～4∶6，表层土较松散，向下碎块石含量增大。下部物质为块石土，块石成分为泥质灰岩，块径大小不一，通常为30～120 mm，块径最大约为260 mm，土质为土石比5∶5～3∶7的可塑性粉质粘土层。

滑带为角砾土，灰褐色，岩芯呈土柱状，岩芯内含角砾，砾径为1～4 mm，土质为石比为8:2的可塑性粉质粘土层，地表未见滑带露头。

滑床大部分表现为中风化状态，主要由三叠系巴东组泥岩、砂岩、泥灰岩组成。产状均为140°∠30°，通过勘查知道：岩石内部有正在发育的裂隙且有钙质填充，岩石强度较高，内见泥岩条带，岩石层面清晰。详细剖面情况如图2所示。

图 2 卢家沱滑坡地质剖面图

2 滑坡二维数值模拟分析

2.1 滑坡二维计算方法

以Geo-Studio软件中的SEEP/W模板对卢家沱滑坡二维模型进行不同水位工况组合分析，进而确定出对应的渗流场分布情况；在得到渗流场的分布之后，再由SLOPE/W模板来施加滑体自重及库水压力大小，并结合Morgenstern-Price对滑体进行稳定性计算。

根据卢家沱滑坡的地质结构，选取I-I′主剖面为计算剖面并建立出网格划分模型，确定了单元数为1493，节点数为1544，如图3所示。

图 3 卢家沱滑坡网格图

依照滑坡地勘报告中试验得到物理参数及相关建议值，本文确定出卢家沱滑坡计算参数如表1所示。

表1 卢家沱滑坡计算参数

针对卢家沱滑坡的计算工况考虑水位下降和降雨组合两类，根据三峡库区相关技术要求通常175 m水位下降至159 m速率为0.13 m/d，而}水位降至145 m时下降速率可在0.6～1.2 m/d中进行调控，同时结合降雨条件，本文采用的计算工况如表2所示。

表2 计算工况

2.2 滑坡工况对比分析

根据SEEP/W分析出渗流分布情况后，用摩根斯坦-普瑞斯(Morgenstern-Price)法进行计算，得出滑坡的稳定性系数分布，并在此基础上得出基于两种工况组合的卢家沱滑坡稳定系数。

从工况1的组合中可以看到：当水位由175 m跌落至145 m的阶段里，滑坡的整体稳定性系数随库水位的跌落而持续降低，且在145 m水位线处达到最小值。同时，滑体的稳定性也受到不同水位下降速率的影响，降速越大线性变化幅度越大，但从总体数据上看，稳定系数值变化的范围影响较小，在0.005内波动。由此可看出，滑坡的稳定性受库水位降幅变化的影响不大。

从工况组合2中可以看出：在库水位从175 m跌落的过程中，对施加降雨强度21.7 mm/d的工况进行分析，当在水位处于155～151.4 m区域时，施加降雨强度为55.33 mm/d的50年一遇暴雨进行工况组合，直至145 m水位线处。受到不同水位降速的影响，卢家沱滑坡稳定系数变化程度也均不相同，与工况1一样，曲线的降幅与水位下降速率相关，降速越大降幅越大。当水位线处于155～151.4 m区间中，由于施加50年一遇暴雨强度组合条件，曲线陡然下跌，待降雨停止后，曲线减小幅度明显降低，由此可判断滑坡稳定性受到降雨因素的影响，且具有一定的滞后性。由整体曲线分析，稳定系数值变化的范围影响较小，在0.005内。可得出结论：滑坡的稳定性受库水位降幅变化的影响不大。

工况组合2中叠加50年一遇三日暴雨，其稳定性系数比不叠加降雨的稳定性系数有所降低，可以确定出在暴雨条件下，滑坡的稳定性系数会有一定程度的降低，但影响较小，且随着水位降速增大的同时施加降雨，会明显降低滑坡的整体稳定性，所以当水位降速为1.2 m/d时叠加强降雨为最不利组合情况，该情况下的滑坡稳定性系数为1.057，说明卢家沱滑坡仍可处于基本稳定状态。

3 滑坡渗流场与应力场流固耦合分析

3.1 滑坡物理参数确定

由相关勘察及监测资料可以判断，卢家沱滑坡的形成是由于地表特征、岩层结构、地形地貌等内因与水位升降、降雨、人工活动等外因共同作用下的结果，因此在数值模拟的过程中应充分考虑相关因素影响下的物理参数取值，综合分析，确定住有限元参数取值如表3所示。

表3 卢家沱滑坡有限元计算物理力学参数取值范围

3.2 岩土体本构模型与屈服准则

滑坡体在受到水位升降、人工活动及降雨影响下，会发生各种程度的弹塑性变形，而在变形过程中会伴随着张拉和压剪破坏的产生。本文针对滑坡体在水位变化和降雨条件下的破坏模式，在考虑岩土体的屈服破坏状态时，选择采用拉破坏准则及Mohr-Columb准则共同来判断。

3.3 边界条件

根据饱和非恒定渗流与应力耦合理论来针对卢家沱滑坡进行流固耦合分析，所建立的边坡边界上的应力条件及渗流条件如图4所示。

图 4 卢家沱滑坡边界条件示意图

3.4 计算工况及荷载组合的确定

为重点体现出水位波动与降雨条件下变化过程，本文仅采取普通工况和极端工况(表4中工况1、2)对滑坡进行流固耦合数值模拟分析。

表4 计算及物理模型试验工况及荷载组合

按照滑坡对应的地形平面及剖面(图1、图2)条件，确定了滑体、滑带及基岩的分布情况，建立出卢家沱滑坡三维模型尺寸如图5所示，划分出了4004个单元，以及12 418个相关节点数。

图 5 模型计算网格

3.5 滑坡工况计算

为更清楚的分析滑坡的应力应变等物理变化情况及规律模式，可以采用Abaqus软件来进行不同工况组合下的流固耦合数值模拟分析。现取两种工况进行详细的分析，主要从这两种相对特殊的工况来对比研究滑坡稳定性的状态。

1)应力计算结果及分析

从图6中工况1的应力计算结果中可以发现：当处于175 m水位时，拉应力与压应力的分布区域较为明显。对于滑体而言，最大主应力位于滑体前缘部分，分布区间为-0.986～0.010 MPa；最小主应力位于滑体后缘部分，分布区间为-2.286～0.014 MPa。对于滑带而言，最大主应力位于滑带中后方区域，分布区间为-1.078～0.839 MPa；最小主应力位于滑带后部区域，分布区间为-2.386～0.123 MPa。

从图7中工况2的应力计算结果中可以发现：水位由162 m跌落至145 m同时遭遇50年一遇强降雨条件下，滑坡应力场发生了明显的改变。其中，前缘滑体主要受水位升降作用的影响较大，应力大小发生了改变；中后缘滑体主要受降雨条件的影响较大，拉应力与压应力的分布均发了较大的变化。对于滑体而言，在降雨停止之后，最大主应力位于滑体中前缘部分，分布区间为-1.007～0.075 MPa；最小主应力位于滑体后缘部分，分布区间为-2.341～0.021 MPa。对于滑带而言，最大主应力位于滑带中部段区域，分布区间为-1.103～0.833 MPa；最小主应力也位于滑带中部段区域，分布区间为-2.537～0.077 MPa。由于应力场与渗流场存在相互耦合的关系，因此导致应力发生变化的重要影响因素就是水位升降联合降雨作用促使滑体内部渗流场发生了变化。

图 6 工况1条件下主应力场分布 kPa

图 7 工况2条件下主应力场分布 kPa

2)位移计算结果及分析

对于图8中工况1的位移组合情况，可以得到：在库水位保持175 m时，水平位移最大为2.568 cm，出现在滑体中后部范围内。其中，后缘滑体的水平位移区间为0.297～2.568 cm，中段滑体的水平位移区间为0.287～2.448 cm，前缘滑体的水平位移区间为0.153～0.183 cm。

对于图8中工况2的位移组合情况，可以得到：最大水平位移位于滑体中后缘区域，大小为3.253 cm；在降雨条件停止后，后缘滑体的水平位移区间为0.266～3.189 cm，中段滑体的水平位移区间为0.271～3.253 cm；前缘滑体的水平位移区间为0.153～1.834 cm。经上述结构可以得出：相对于降雨作用，库水位升降对滑坡的位移影响结果相对较小。随着水位下降，水平位移增量先增加后处于减小趋势，水平位移加速度逐渐减小，最后趋于稳定，垂直位移增量变化以及垂直位移变化速率与水平方向位移变化规律一致，而垂直位移加速度先增大后不断减小。

(a)工况1

(a)工况2图 8 位移场分布 m

3)塑性区计算结果及分析

从图9a中工况1塑性区分布情况可以看到滑体不存在明显的塑性破坏范围，基本可以确定当蓄水至175 m水位时，滑坡没有发现显著的破坏迹象,一直处于稳定状态。从图9b中工况2塑性区分布情况可以看到，当水位从162 m降至145 m且施加50年一遇强降雨工况条件下，滑体不存在明显的塑性破坏范围，此时滑体没有发现明显的破坏迹象，处于基本稳定状态。

图9 塑性区分布图

4)孔压计算结果及分析

由图10a可知：175 m水位以下范围内，水位变化直接影响着孔压的分布，前缘最大孔压值大小为0.353 MPa，分布在滑体140 m高程处；地下水位对滑坡后缘孔压影响较大，最大孔压值为2.286 MPa。由图10b可知：库水位升降及降雨对滑坡前缘地下水的分布影响较大，当水位区间在162 m至145 m跌落时，地下水位也相对降低且具有一定的滞后性特征；库水位的升降作用对后缘滑坡的影响程度不大，但是受到降雨条件的影响，后缘孔压值会有一定的增加。

图10 孔压值显示图

5)滑坡整体稳定性

当库水位升降和大气降雨联合作用条件下，滑坡前缘至后缘均出现了一定程度的变形迹象，不同条件下的应力分布及位移量均存在差异，局部变形较为明显，可能出现局部滑动的情况，但在滑体内部未出现大面积的破坏区域，说明卢家沱滑坡目前仍保持基本稳定状态。

4 结论

以卢家沱滑坡为研究对象，通过Geo-Studio软件计算在不同水位、下降速率及降雨等组合工况条件下的滑坡稳定性系数分布；同时根据Abaqus软件模拟了卢家沱滑坡在库水位变动联合降雨组合工况下的渗流场、应力场、位移场及塑性区等流固耦合及三维稳定性情况，得到具体结论如下。

1)从二维稳定性系数曲线图可以看出:滑坡稳定随库水位下降速率增大而不断降低,叠加降雨条件情况下滑坡稳定性出现最小值,滑坡前缘坍岸将会加剧,但由于滑坡前缘土体为中-强透水性,前缘内的地下水渗出较为容易,滑坡稳定性受水位下降影响不大,失稳的可能性较小。

2)从三维数值模拟获得的位移、应力、孔压及塑性区在各工况下的计算结果表明:在库水位和降雨耦合作用下，滑坡会产生一定的变形,但整体仍处于基本稳定状态,仅局部区域(滑体后部及滑体前缘)在库水作用下变形相对较大,且同时考虑降雨影响,滑坡整体位移量相对较大。

3)由二维模型计算结果从宏观方向上分析，可以得出在库水位升降作用下，卢家沱滑坡的稳定性逐渐降低，但水位下降速率对滑坡的影响不大，若联合降雨条件作用下，对滑坡整体稳定性影响相对较大；由三维流固耦合模型结果从微观上分析，可以看出库水位变化联合降雨条件将会改变滑坡内部的渗流场，导致滑坡发生较大的变形，若遭遇不利暴雨条件，卢家沱滑坡存在局部失稳的可能。