基于拟静力法分析库水变化条件下地震荷载对某边坡稳定性影响研究

曹千红 伍 岳 李 耀 邓永华

（1.湖北省土地规划勘测院，武汉 430071；2.三峡大学 土木与建筑学院，湖北 宜昌 443002；3.宜都北斗星科技发展有限公司，湖北 宜昌 443300）

水库蓄水诱发地震在世界水利史上已不鲜见［1］，国内学者对水库诱发地震的机理及其对边坡工程的危险性评价和预测方法都作了较为深入的研究工作［2－6］，而地震作用对库岸边坡的变形破坏有何影响及影响程度有多大在这些研究工作中却极少提及.

本文采用有限元方法分别对受库水变化和地震荷载综合作用和仅受库水变化作用两种工况条件下某土质滑坡的位移和等效塑性应变进行了计算，通过两种工况计算结果的比较得出地震荷载对库岸边坡稳定性的影响.

1 地震荷载模拟方法

目前地震荷载模拟方法主要有拟静力法和动力分析法.动力分析法虽然计算精度高，但操作复杂，计算结果的评价在现有的规范中亦没有明确的标准［7］；拟静力法计算简单，使用经验丰富，抗震设计中被广泛应用［8－10］.

拟静力法是一种用静力学方法近似解决动力学问题的简易方法，它发展较早，迄今仍然被广泛使用［11］.其基本思想是在静力计算的基础上，将地震荷载简化为一个惯性力系附加在研究对象上，其核心是设计地震加速度的确定问题［12］.

本文将地震荷载简化为水平加速度，来模拟地震作用对滑坡稳定性的影响.根据GB50011－2010国家标准建筑抗震规范，设计基本地震加速度值与抗震设防烈度对应关系见表1.

表1 设计基本地震加速度值与抗震设防烈度关系表

计算中水平加速度选取抗震设防烈度为Ⅶ时对应的加速度0.1g.

2 渗流场－应力场耦合原理及双场耦合数值模型

渗流场和应力场的相互作用可表述为［13－14］：地应力作用于岩土体的空隙结构，改变地下水的运移通道，进而改变渗透系数、渗透力；岩土体内存在的地下水，通过物理、化学和力学等作用亦改变岩土体的结构，给岩土体施加静水压力和动水压力，改变岩土体应力场的分布.

根据渗流场与应力场的耦合机理，将连续介质渗流场数值模型、应力场数值模型及渗流与应力关系的经验模型［15］组合在一起，得到渗流场与应力场耦合的数值模型如下：

式（1）是渗流方程，［k］为渗透系数矩阵，｛P｝是渗透压力矩阵，｛Q｝为汇源项列阵，［S］为给水度列阵；（2）式与（3）式是应力场方程，｛K｝为刚度矩阵，｛U｝为位移列阵，｛F｝为外部载荷列阵，｛σ｝为应力矩阵，［D］为弹性矩阵，［B］为几何矩阵；式（4）为Louis在1967年由实验得出的表征岩体应力与渗流的耦合关系为节点i沿j轴方向的渗透系数，其中t表示迭代步数，n表示耦合节点总数.

3 滑坡概况

滑坡地处大巴山与秭归盆地的过渡地带，属于构造侵蚀中浅切割中低山河谷地貌.长江自西向东从滑坡前缘流经，滑坡呈弧形向长江凸出.滑坡地形为斜坡，北高南低，倾向长江，滑坡前缘临江一带及滑坡后缘（位于滑坡后部公路上下一带）地形较陡，地形坡度20～40°，滑坡中部地形较缓，地形坡度一般为14～20°，局部较陡，滑坡前缘高程125m，后缘高程365m，相对高差约240m.

滑坡为土质滑坡，其成因为崩塌滑坡堆积形成.滑体中上部块碎石以紫红色泥岩为主，夹少量泥灰岩，下部黄色土层块碎石为泥灰岩、灰岩，在滑坡后部地表可见大量石英砂岩孤块石；滑带一般位土岩接触面，其物质成分主要为含砾粉质粘土，多呈灰黄色，局部呈紫红色，厚10～150cm，滑带埋深40～75m，滑面形态呈折线状，后缘较陡，前缘较缓，坡角10～36°；滑床由三叠系巴东组泥灰岩、泥岩组成，岩层倾向340～20°，倾角16～26°，节理裂隙较为发育.

4 数值计算模型及计算参数

根据所选滑坡剖面的地质结构特征，建立该滑坡剖面二维数值计算模型.模型垂直于长江水流方向为407.7m，前缘高程78m，后缘高程288m.

计算域包含滑体、滑带和基岩，整个计算域剖分了6 079个六面体单元，共计12 414个节点，二维计算模型与网格如图1所示.

图1 模型计算网格

根据该滑坡地质勘查报告中提供的有关岩土体的物理力学性质试验数据及参数建议值为基本依据，并类比三峡库区其它类似滑坡岩土体的物理力学性质参数值获得该滑坡二维空间有限元分析计算参数取值见表2.

表2 滑坡有限元计算物理力学参数取值表

5 计算工况

按照三峡水库水位调度曲线，每年从1月1日到5月31日，库水位由175m水位降低至145m水位，历时约150d，下降速率约为0.2m／d；从9月30日到10月31日，库水位由145m水位上升至175m水位，历时约30d，上升速率约为1.0m／d.

根据此水位变化情况，结合地震作用制定如下计算工况（见表3）.

表3 计算工况

6 边界条件

1）渗流边界条件.数值计算模型四周及坡面不接触库水的边界为自由透水边界，接触库水的坡面为动水头边界，其水头高度等于库水水位.

2）应力边界条件.计算域的底部边界为固定支座约束，前后两侧及左右两侧均采用法向约束，坡体表面为自由边界.

7 计算结果及分析

采用ABAQUS有限元计算软件对滑坡剖面二维数值计算模型在两种工况条件下的位移和等效塑性应变分别进行了计算.工况1和工况2两种工况条件下的位移计算结果见图2，等效塑性应变计算结果见图3.

从图2（a）和2（c）所示的水平位移计算结算结果可以看出：1）在工况1条件下，滑坡的水平位移主要集中在滑体下部表层部分，量值约为0.14～0.19m；2）在工况2条件下，滑坡的水平位移主要集中在滑体中部表层部分，量值约为0.17～0.25m；3）通过两种工况计算结果进行对比，可以看出受到地震作用后，滑体水平位移集中区域由滑体下部转移到滑体中部，水平位移集中区域的面积有较为明显增大，水平位移最大值增量约为0.06m，其增幅约为32%.

图2 不同工况条件下位移计算结果

从图2（b）和2（d）所示的水平位移计算结算结果可以看出：1）在工况1条件下，滑体的垂直位移主要集中在滑体上部表层部分，量值约为0.25～0.33m；2）在工况2条件下，滑体的垂直位移主要集中在滑体中部表层部分，量值约为0.26～0.40m；3）通过两种工况计算结果进行对比，可以看出受到地震作用后，滑体垂直位移集中区域由滑体上部转移到滑体中部，垂直位移集中区域的面积有较为明显增大，垂直位移最大值增量约为0.07m，其增幅约为21%.

图3 不同工况条件下塑性应变计算结果

从图3（a）和3（b）所示的塑性应变计算结算结果可以看出：1）在工况1条件下，滑坡的等效塑性应变分布范围较小，集中在滑带前缘，等效塑性应变最大值为0.017；2）在工况2条件下，滑坡的等效塑性应变在滑带处较大范围出现，主要集中在滑带上部和中部，滑带前缘较小范围出现，等效塑性应变最大值为0.075；3）通过两种工况计算结果进行对比，可以看出受到地震作用后，等效塑性应变分布范围在滑带处明显扩大，等效塑性应变值明显增大.

8 结 语

通过对滑坡在2种工况条件下的计算结果进行对比分析表明：

1）在地震作用下，滑体位移有较明显增大，地震作用对滑体水平位移的影响大于对垂直位移的影响.

2）在地震作用下，滑坡等效塑性应变在滑带处的分布范围明显扩大，等效塑性应变值有较为明显增大.

3）地震作用是影响滑坡稳定性的不利因素，受到地震作用后滑坡失稳的可能性增大.

［1］ 易立新，车用太，王广才.水库诱发地震研究的历史、现状与发展趋势［J］.华南地震，2003，23（1）：28－37.

［2］ 曹建玲，石耀霖.河道型水库蓄水诱发地震的数值模拟［J］.中国科学院研究生院学报，2011，28（1）：19－26.

［3］ 黄振伟，黄海蛟，丁淑平，等.丹江口水利枢纽大坝加高工程水库诱发地震初步分析［J］.资源环境与工程，2011，25（5）：391－395.

［4］ 罗启祥，王 昆，高海军，等.云鹏电站水库诱发地震分析［J］.云南水力发电，2010，27（1）：30－34.

［5］ 张丽芬，廖武林，曾夏生，等.水库诱发地震对混凝土坝的影响及抗震设防［J］.长江科学院院报，2010，27（5）：71－75.

［6］ 尹岳明.印尼Jatigede大坝工程水库诱发地震的判定及评估［J］.四川水利，2011（1）：51－53.

［7］ 宫经伟，严新军，夏新利，等.拟静力法在溢洪道闸室抗震稳定计算中的应用［J］.人民黄河，2010，32（1）：107－108.

［8］ 郭诚谦.混凝土面板堆石坝坝顶结构设计及其拟静力法的抗震分析［J］.水利水电技术，2008，39（5）：26－30.

［9］ 於昌荣，黄 真，黄醒春，等.用拟静力法对斜隅支撑框架的地震反应分析［J］.铁道建筑，2008（12）：110－113.

［10］王运永，高永涛，蔡嗣经，等.考虑地震影响的地下连续墙稳定性拟静力法研究［J］.北京科技大学学报，2010，32（5）：559－564.

［11］刘如山，胡少卿，石宏彬.地下结构抗震计算中拟静力法的地震荷载施加方法研究［J］.岩土工程学报，2007，29（2）：237－242.

［12］包世华.结构动力学［M］.武汉：武汉理工大学出版社，2004.

［13］王晓鸿，忤彦卿.渗流场－应力场耦合分析［J］.勘查科学技术，1998（4）：3－6.

［14］李宗坤，陈丽刚，韩立炜.基于ABAQUS渗流与应力耦合的边坡稳定性分析［J］.人民黄河，2011，33（2）：103－104.

［15］梁 冰，鲁秀生.裂隙岩体渗流场与应力场的耦合数值分析［J］.水资源与水工程学报，2009，20（4）：14－16.