降雨条件下膨胀岩边坡失稳数值模拟研究

徐 晗，饶锡保，汪明元

（长江科学院水利部岩土力学与工程重点实验室，武汉 430010）

1 概 述

膨胀岩边坡失稳是工程中经常遇到并难于处理的问题，正在设计的南水北调中线工程的膨胀岩土渠段长达340 km，如何正确进行边坡的稳定性分析是合理设计的前提。在实际工程中，膨胀岩边坡大多是在降雨时发生破坏的，因此有必要研究降雨条件下膨胀岩边坡的失稳过程。目前数值模拟降雨入渗引起边坡失稳较为常见的方法是先用渗流理论求解瞬态的渗流场，再用技术已经成熟的极限平衡法进行土坡暂态稳定性分析［1－4］；但降雨条件下边坡的失稳是一个应力场和渗流场不断变化、相互作用的渐进过程，因此，在考虑降雨入渗边坡稳定性分析中，引入土体的本构关系，用非线性有限元流固耦合的分析方法模拟边坡的稳定性，具有重要的意义。

经过长时间的研究与探讨，作者认为对膨胀岩边坡的稳定性数值分析应正确模拟膨胀岩的膨胀性、裂隙性，以及非饱和膨胀岩的吸湿软化特性。目前多数学者在膨胀岩边坡稳定分析技术中，不考虑其非饱和强度基本特性及膨胀性、裂隙性［5］，使得计算结果偏离实际，因此，有必要对膨胀岩边坡稳定技术进行改进，使之能考虑膨胀岩与普通土的不同之处，揭示膨胀岩边坡破坏的浅层性与渐进性特点。为此，采用流固耦合的非线性有限元分析方法，对膨胀岩边坡降雨失稳现场试验进行了数值模拟研究，改进了膨胀岩边坡稳定分析技术，使之能考虑膨胀岩的膨胀性、裂隙性，以及非饱和膨胀岩的吸湿软化特性。

2 流固耦合的分析方法

采用Galerkin有限元格式，将节点位移和孔隙水压力作为节点自由度进行空间离散，则力学平衡方程及渗流连续方程可写成如下矩阵形式［6］：

平衡方程为

式中：［K］为通常的刚度矩阵；｛Δ珋δ｝为位移增量；［L］为节点孔隙水压力所对应的节点力；｛Δ珋p｝为孔隙水压力增量；｛F｝为节点外荷载；｛I｝为增量迭代过程中上一增量步中的不平衡力。

渗流连续方程为

直接解耦合方程，引入时间积分，差分公式为

式中：ξ为一系数，0≤ξ≤1；Δt为时间增量。

为了确保数值稳定，可以选择ξ＝1（后向差分法），此时方程（3）可写为

因此，在t＋Δt时刻的渗流连续方程为

其中｛R｝＝Δt［－｛Q｝t＋Δt＋］T｛珋v｝t＋Δt＋［H∧］T｛p｝t＋Δt］，为流体体积变化修正量；［B］为节点变形对应的流体体积改变；［H］为孔压变化对应的流体体积改变。

3 膨胀岩的特殊性质

3．1 吸湿强度软化特性

近年来，非饱和土力学的发展为定量计算因水分入渗而引起的土体软化的强度变化提供了理论基础。Fredlund（1978年）提出一个抗剪强度公式为

式中：τf为非饱和土抗剪强度＇，φ＇分别为有效凝聚力和有效内摩擦角；σn为法向应力；φb为随吸力变化的内摩擦角。根据非饱和膨胀岩三轴试验，在400 kPa的基质吸力范围内，非饱和原状黏土岩的抗剪强度随着基质吸力的增大而线性增加，该线段的斜率为基质吸力引起的强度增量tanφb，其中＝20°。

因此，可在400 kPa的基质吸力范围内取φb＝20°，但非饱和土的强度不可能随基质吸力无限制提高，设定当基质吸力大于400 kPa时，抗剪强度按400 kPa选取。

3．2 土水特征曲线

泥灰岩的土水特征曲线采用室内试验成果。水分运动数值分析中的有关非饱和土参数将利用实测的土水特征曲线确定，并采用VG模型进行计算。Van Genuchten模型能够在较大的水头范围内表征水分特征数据，应用较为广泛

式中：Ks为饱和渗透系数；Se＝θ－θr／（θs－θr）为有效饱和度；α，m和n为模型的形状参数；参数m和n的关系可表达为m＝1－1／n。

将泥灰岩原状样非饱和试验得到的土水特征曲线进行拟合，所得公式的具体形式为

式中：θr＝0．12，θs＝0．266 5；α＝0．061；m＝0．14；n＝1．163；h的单位为 kPa。

将黏土岩原状样非饱和试验得到的土水特征曲线进行拟合，所得VG模型公式的具体形式为

式中：θr＝0．12，θs＝0．312 016；α＝0．002 179；m＝0．178 1；n＝1．216 7；h的单位为 kPa。

3．3 吸湿密度增加

对任意时刻计算出来的单元积分点的饱和度，根据土的三相比例指标换算关系可以求出该情况下的含水量，然后根据干密度可以求出该时刻非饱和土的湿密度ρ。

式中：w为含水量；e为孔隙比；Gs为相对密度；Sr为饱和度；ρd为干密度。

3．4 体积应力、含水量变化相关的膨胀性

假设膨胀应变是由于应力第一不变量的改变所引起的，膨胀应变仅与体积应力有关，而与应力偏斜张量无关，采用GDS公司的动三轴仪器，对中膨胀土试样开展了三轴膨胀实验研究，得到了以体积应力、含水量变化为主要变量的膨胀应变公式，公式如下。

式中：ε为膨胀体积应变，按百分数表示；Δw为含水量变化百分数；p为体积应力；Pa为大气压力。

3．5 裂隙性

膨胀岩的裂隙性是影响膨胀岩边坡稳定失稳的最直接最主要的原因，且受含水量影响显著，在膨胀岩边坡稳定性分析中，可采用用残余强度作为风化层的最主要的强度参数，这为膨胀岩边坡稳定评价提供了一定的强度储备［7］。但在边坡一定深度范围内，如超过大气影响深度，膨胀岩可视为保持原状土的性质，在选择强度参数时，可采用峰值强度。

4 膨胀岩现场试验数值模拟

4．1 数值分析的基本思路

（1）根据现场试验报告描述建立自重作用下的初始应力场，设定地下水位，建立初始含水量场。

（2）进行降雨数值模拟，期间考虑如下过程（整个过程为连续计算）：入渗引起的边坡含水量变化场；随着边坡含水量的变化引起的土体重度增加（由湿重度变为饱和重度）；土体非饱和强度指标随含水量增加而软化降低；孔压与应力场的耦合作用。

（3）对于膨胀力的考虑，根据三轴膨胀试验提出的本构模型，分阶段或者仅在计算的最后阶段考虑膨胀力的影响（当计算达到某一阶段后，根据当前计算点的应力与含水量场的变化，假定在该阶段膨胀过程中应力场不变，采用三轴膨胀试验提出的本构模型计算膨胀引起的体积应变，转化为温度荷载施加到模型上）。

（4）边坡稳定计算，对任一降雨时刻，提取该时刻的单元的强度指标、重度、孔压、含水量变化引起的膨胀应变，采用强度折减法来计算边坡稳定。

4．2 现场试验简介

长江科学院岩土力学与工程重点实验室在新乡膨胀岩进行了大型的膨胀岩降雨失稳现场试验研究，该边坡从坡底到坡顶17 m高，一级马道以下高9 m。降雨试验主要在一级马道以下进行，坡比为1∶1．5，降雨范围具体如图1斜方向箭头所示。一级马道以下大致可分为5层，位于顶部为壤土层，厚约0．2 m；下为厚约6．8 m的泥灰岩、其间有0．8 m厚黏土岩夹层，最下层为黏土岩。除表层壤土外，各地层自由膨胀率分别为21%～29%，63%，40%～54%，48%～70%。

图1 边坡有限元网格及降雨示意图Fig．1 Slopemesh and rainfall schematic diagram

该试验从2008年9月24日开始，在3 d的正式监测期间，共进行了5次降雨过程，现场试验降雨强度为16 mm／h，共历时16 h，在数值计算中假定降雨为连续不间断过程。

当现场试验累计降雨约6 h左右时，一级马道以下开始发生大面积滑坡，发生整体的塑性变形破坏，导致表层水平位移急剧增加至12 mm，在滑坡后的第2 d继续降雨1 d后，水平位移持续以较大的幅度增加。

4．3 初始状态考虑

在计算中将地下水位线设在渠坡底上1 m左右。由于该边坡开挖后表层裸露在大气作用下半年后才开始降雨，因此，根据泥灰岩、黏土岩入渗深度设定表层2 m范围内左右为裂隙发育区，强度参数按照有裂隙情况下的残余强度选取。表层2 m以下为未风化层，强度参数按照峰值强度选取。

4．4 计算参数取值

计算参数均根据泥灰岩与黏土岩原状样室内试验统计选取。如表1与表2所示。

表2 壤土与砂岩的材料参数Tabale 2 Material parameters of loam and sand rock

5 计算结果分析

5．1 初始状态

由于该边坡高度仅17 m，因此可将饱和度初始状态设定为沿地下水位的自重分布，如图2所示。位移场采用自重作用下的位移分布，水平位移与竖向位移如图3所示。可知一级马道坡面的水平位移约在2 mm左右，竖向位移在1～3 cm之间，处于稳定状态。

图2 降雨试验开始前饱和度Fig．2 Saturation degree before rainfall

表1 材料参数取值（由室内试验统计得到）Table 1 Material parameters（acquired by indoor test）

图3 降雨试验开始前水平位移、竖向位移Fig．3 Horizontal displacement and vertical displacement before rainfall

5．2 变形分析

为了监视边坡在降雨过程中发生的变化，在坡面设置了3个监控点，分别在坡脚A、坡面B、坡顶处C。具体示意图如图1所示。图4为降雨后边坡饱和度分布图，可知在边坡2 m范围内接近于饱和状态。

图4 降雨后饱和度Fig．4 Saturation degree after rainfall

图5 （a）为降雨引起的等效塑性应变，图5（b）为考虑膨胀力后的等效塑性应变，可知膨胀力会加剧塑性区的扩散作用。图6（a）至图6（d）分别为降雨后未考虑膨胀力与考虑膨胀力位移分布图，基本规律为坡顶水平位移最小，坡底水平位移最大，坡面中间渐变。可知未考虑膨胀力坡底水平位移为9．3 cm，在考虑膨胀力后，位移场基本上分布在边坡坡面附近，坡底水平位移达38 cm，表明边坡已发生较大的滑动变形。

图5 等效塑性应变Fig．5 Equivalent plastic strain

图7 （a）为降雨过程中监控点饱和度过程曲线，可知8 h左右（30 000 s）边坡坡面会形成一层饱和状态，此时降雨入渗率会逐步下降并维持稳定。图7（b）为降雨过程中监控点（未考虑膨胀）的合位移过程曲线，可知在降雨过程中坡脚A点位移变化迅速，从刚开始的1 cm左右迅速变化到12 cm以上；位于坡面中心B点位移变化比坡脚小，但也达到了8 cm；位于坡顶的C点位移变化很小，说明边坡的破坏过程是由坡脚逐步向坡顶延伸的渐进过程，具有明显的渐进破坏特点。

根据变形分析可知：

（1）降雨对膨胀岩边坡产生较大的影响，降雨过程中坡脚范围会先饱和并软化，而且坡脚是应力集中的部位，因此在该处首先出现较大的塑性应变。

（2）随着降雨的进行，饱和度逐渐向边坡内部延伸，重度进一步增加，下滑力逐渐增大，同时，其抗剪强度又因含水率增加而急剧下降，位移也因此非线性地增大；塑性应变逐步沿着坡脚向上部扩展开来，这是降雨引起边坡渐进性破坏的主要原因。

（3）当考虑膨胀特性后，等效塑性应变范围由坡顶到坡底形成贯穿的区域，表明膨胀力会加剧边坡下滑趋势，是滑坡的重要影响因素。

图6 降雨后位移Fig．6 Displacements after rainfall

5．3 稳定性分析

图8 （a）、（b）分别为初始状态、降雨后的滑弧示意图，可知初始状态安全系数为1．45，处于稳定状态，而降雨后处于滑动状态。

由稳定性分析可知：

（1）只有考虑膨胀岩的非饱和强度基本特性，才有可能保证边坡存在浅层裂隙而未降雨情况下仍处于稳定状态，否则边坡在初始状态下就会失去稳定，不可能呈现为深层滑弧。

图7 监控点饱和度和合位移的过程曲线Fig．7 Saturation degree and displacement process curves atmonitoring points（without considering expansion force）

图8 滑弧示意图Fig．8 The sketch of the sliding curves

（2）考虑了膨胀岩的膨胀特性后，边坡的安全系数会进一步降低，加剧边坡稳定破坏。

（3）在降雨后，由于浅层存在裂隙，入渗较快，强度指标急剧下降，同时浅层重度也在不断增加，致使下滑力逐步增大的同时抗滑力逐步减小，从而引起安全系数急剧下降，稳定由深层滑弧逐渐变化到浅层滑弧。如果没有浅层裂隙的存在，雨水难于入渗，也不会因含水率增加而引起抗剪强度下降。因此大气影响带膨胀岩的风化、干湿循环影响的浅层裂隙，是膨胀岩边坡浅层滑动的根源。

6 结 语

（1）膨胀岩边坡在雨水入渗的情况下，重度增加，下滑力增大，同时，其抗剪强度又因含水率增加而急剧下降，这是降雨引起边坡渐进性破坏的主要原因，而膨胀特性进一步加剧了这种影响。

（2）大气影响带膨胀岩的风化和干湿循环影响导致的浅层裂隙，是膨胀岩边坡浅层滑动的根源。

（3）本文的方法全面考虑了膨胀岩的基本特性与特殊性质，可以完整地描述膨胀岩边坡降雨引起滑坡发生的渐进破坏过程。通过对膨胀岩现场试验的数值模拟，完整再现了膨胀岩边坡降雨失稳的整个过程，且材料参数均为室内试验得到的成果，故分析成果具有相当的实用价值。

［1］ 林鲁生，蒋 刚．考虑降雨入渗影响的边坡稳定分析方法探讨［J］．武汉大学学报（工学版），2001，34（1）：42－44．

［2］ 陈善雄，陈守义．考虑降雨的非饱和土边坡稳定性分析方法［J］．岩土力学，2001，22（4）：447－450．

［3］ 魏保义．用非饱和渗流理论研究降雨人渗的特性及对边坡稳定的影响［D］．南京：河海大学，2001．

［4］ 刘艳华．膨胀土边坡的稳定性研究［D］．武汉：长江科学院，1999．

［5］ 韦立德，杨春和，徐卫亚，等．考虑饱和一非饱和渗流场和应力场耦合的三维强度折减有限元程序研制［J］．水文地质工程地质，2006，33（3）：16－20．

［6］ Hibbitt，Karlson＆Sorensen，Inc．ABAQUSManualsversion 6．3［M］．USA：Hibbitt，Karlson＆Sorensen，Inc．，2003．

［7］ 陈善雄，陈守义．膨胀土判别与分类方法探讨［J］．岩土力学，2005，26（12）：1895－1900．