基于孔压-应力耦合下考虑渗透系数各向异性对库水位骤降边坡的影响

2021-12-06戚海棠郁舒阳任旭华张继勋

石河子大学学报（自然科学版） 2021年4期

戚海棠，郁舒阳，任旭华，张继勋

(河海大学水利水电学院，江苏南京 210098)

库水位的上升导致滑坡土体内部孔隙水压增大，土体的强度参数及有效应力降低[1-3]；此外，库水位的下降导致滑坡内部孔压变化滞后库水位，对周围土体骨架产生向临空面及向下的渗流动水压力(动态扩张力)，牵引带动坡面土体[4-6]，因此，对库水位变动下的滑坡失稳规律及影响因素的把握是正确认识库水位渗流作用机理及防治滑坡灾害的关键所在。

国内外学者针对库水位波动下的滑坡渗流力学响应展开了大量研究，其中，IQBAL J等[7]分析研究滑坡与周围环境因素的关系，LUO F Y等[8]研究发现随着坡度的增加、库水位下降下的边坡滑移面逐渐变深，HU X L等[9]利用非饱和渗流理论及有限元原理对黄土坡滨江滑坡稳定性进行了数值模拟分析，JIAN W X等[10]对不同降雨强度及库水位骤降速率下的千江坪滑坡的渗流稳定性进行了数值模拟，但是这些研究未考虑边坡土体的渗流各向异性。宋云奇等[11]研究表明黏土等土质因其絮状微观结构导致渗透系数存在各向异性，同时土质的渗透系数各向异性受其干密度及冻融循环等因素的影响巨大[12]；YEH H F[13]考虑土质边坡的渗透系数各向异性程度，对其渗流特性及局部安全系数进行数值模拟，但忽略了各向异性方向的影响。

由2013年6月至2015年12月三峡库区的水位波动及滑坡数量统计示意图[14-15]可知，三峡库区库水位变动是导致三峡库区滑坡的重要诱因。针对以往研究的不足，本文以三峡库区蔡坡堆积体为例[16]，对砂土与黏土2种滑坡土质，考虑到库水位下降下不同渗透系数各向异性程度kr及方向α情况的渗流、变形及稳定性进行有限元分析，从而为全面认识边坡渗流各向异性规律及对滑坡灾害的防治提供了一定的参考。

1 非饱和孔压-应力耦合计算模型

非饱和渗流各向异性孔压-应力耦合及稳定性计算利用Geo-studio中的SIGMA/W与SLOPE/W模块进行耦合分析，SIGMA/W对有限元网格的每一个节点创建3个方程，其中二个是平衡方程(位移)，另一个是渗流连续方程(孔隙水压力)，同时，求解3个方程能够认识位移跟孔隙水压力的变化。SLOPE/W可以根据SIGMA/W计算得出的孔压进行非饱和安全系数的计算。

1.1 非饱和土渗流有限元计算理论

根据饱和及非饱和达西定律[17]推导出渗流计算有限元公式中的渗流控制方程，其二维表达式为

(1)

式(1)中x、y分别为x、y方向的位置坐标，kx为x方向的渗透系数，ky为y方向的渗透系数，H为总水头，Q为施加的边界流量，t为时间，mw为储水曲线的斜率，γw为水的容重。

1.2 渗流-应力耦合计算理论

对于二维问题，非饱和土介质的增量应变-应力关系可以表达成为：

{Δσ}=[D]{Δε}-[D]{mH}(ua-uw)+{Δua},

(2)

式(2)中ε为法向应变，γ为工程剪应变，σ为法向应力，τ为剪切应力，ua为孔隙气压力，uw为孔隙水压力，[D]为排水本构矩阵，{mH}可以表达为：

(3)

式(3)中H′为与基质吸力(ua-uw)有关的函数。

1.3 非饱和安全系数

边坡安全计算系数采用基于极限平衡理论的Morgenstern-Price法，此法严格满足力平衡及力矩平衡，计算精度较高。

2 计算模型及边界条件

计算模型选择三峡库区蔡坡堆积体[16]为例，蔡坡堆积体位于鹤峰县燕子乡境内，距坝址上游约为1.32 km，其平面形态近似矩形，顺岸坡展布，地形坡角15°～40°，前缘高程约160 m，后缘高程220 m，总体东高西低。滑坡厚度为2～4 m，总面积约1.2万m2，体积约3.6万m3。为了减小边界条件对计算结果的影响，将坡脚和坡顶的范围进行延伸，整个模型共划分为1 084个节点、1 001个单元,计算模型如图1所示。

图1 计算模型及边界示意图

为了综合研究边坡土体各向异性对边坡不同位置的渗流、变形特性的影响，设置上部(x=100 m)、中部(x=200 m)、下部(x=300 m)3个监测面反映边坡不同位置的渗流、变形特性，监测面贯穿滑坡体。初始条件为AB 190 m、DF 175 m水头计算所得的稳定渗流场；边界条件设置如下：AB为190 m水头边界，DEF为175～145 m库水位变动边界，AF、CD为不透水边界。

3 材料参数及计算工况

3.1 非饱和计算参数

土-水特征曲线(SWCC)采用Fredlund & Xing模型[18]进行估算，该方法的控制方程见文献[18]。

边坡材料选取典型的砂土(渗透系数较大)及黏土(渗透系数较小)[19]，进行边坡非饱和渗流分析时的计算参数设置[20]如表1所示，同时，砂土与黏土的土-水特征曲线如图2所示。

表1 二种非饱和土体的计算参数[20]

图2 砂土与黏土的土-水特征曲线

3.2 物理力学参数

砂土与黏土的力学参数根据室内试验确定，采用Mohr-Coulomb弹塑性准则，其物理力学性质参数[20]见表2。

表2 砂土与黏土的物理力学参数[20]

3.3 渗流各向异性定义及计算工况

以往理论、试验及数值研究中大多忽略边坡土体的各向异性程度及各向异性方向，然而事实上无论是黏土，砂土等边坡土质均存在各向异性，对于渗透矩阵[C]而言，一般表达式可以表达为：

(4)

式(4)中C11=kxcos2α+kysin2α，C22=kxsin2α+kycos2α，C12=C21=kxsinαcosα+kysinαcosα。其中，kx、ky、α根据图2定义，kx为水平向渗透系数，ky为垂直向渗透系数，α为ky与y轴的方向夹角。当α=0°时，[C]退化为：

(5)

式(5)是各向异性渗透矩阵，定义各向异性程度kr=kx/ky，但是各向异性不仅仅存在水平竖直渗透系数不一致的问题，同时存在各向异性方向α的问题。

为综合讨论边坡不同土质的各向异性，设置以下工况：滑坡体土质采用黏土及砂土，各向异性程度kr=1、10、50、100,同时各向异性方位角α=0°、15°、30°、45°、60°、75°及90°。库水位从175 m高程变化到145 m高程，下降速率为1.2 m/d，同时考虑到库水位下降停止后的55 d，总共计算时间为80 d。计算工况如表3所示。

表3 计算工况

4 计算结果与分析

根据表3的计算工况对砂土与黏土进行56组工况的计算，依照砂土及黏土分类分析各向异性程度kr变化对边坡渗流稳定的影响。

4.1 各向异性对渗流特性的影响

4.1.1 砂土

砂土边坡不同各向异性程度kr及方向α下的浸润线变化如图3所示。由图3可见：库水位下降停止时刻，在接近库岸部位的浸润线存在上凸现象，即岸坡内部的浸润线变化要滞后于库区水位的变化。不同各向异性方向α对浸润线变化影响巨大，随着各向异性角度α的增大，库水位下降停止时刻(第25 d)岸坡内部的浸润线高程逐渐抬高，但是增加高程却逐渐衰减，对于α大于45°情况下浸润线几乎一致。

图3 砂土各向异性程度及方向对浸润线变化影响

不同各向异性程度kr对浸润线的影响也较大。随着各向异性程度kr的增大，浸润线抬高迅速，对于kr大于50的情况，浸润线高程几乎一致。

4.1.2 黏土

黏土边坡不同各向异性程度kr及方向α下的浸润线变化如图4所示。由图4可见：黏土边坡的整体浸润线比砂土边坡更高，同时，不同各向异性程度及方向下的浸润线差异较之砂土边坡更小，各向异性程度kr越大，各向异性方向角α越大，浸润线的高程越高。

图4 黏土各向异性程度及方向对浸润线变化影响

4.1.3 浸润线滞后高程分析

通过4.1.1及4.1.2节可以发现，无论是砂土还是黏土边坡，库水位骤降下岸坡内部浸润线变化均滞后于库岸的水位，这种滞后性会导致局部的渗透坡降陡增，渗流力指向坡外，加剧岸坡失稳。因此为定量化研究这种滞后效应，定义库水位骤降下的浸润线滞后高程为图5a中下部监测面与库水位骤降停止(第25 d)的浸润线的交点高程至死水位145 m高程的高度。图5b、图5c分别为不同渗流各向异性程度及方向下的砂土边坡与黏土边坡的浸润线滞后高程。

总体上看，随着各向异性程度kr及方向角α的增大，砂土边坡与黏土边坡的浸润线滞后高程逐渐增大至极大值。这是因为各向异性方向角α增大，对于kr>1的情况而言，与x轴平行方向的渗透系数越来越小，坡体内部的水很难排渗到库区，因此水位滞留现象越明显，因此浸润线滞后高程越大。同时，各向异性程度kr越大，土体的竖直向的渗透系数也越小，库水位较难下渗，因而浸润线滞后高程也越大。对比砂土与黏土边坡，可以发现渗流各向异性特性对砂土边坡的浸润线滞后高程影响巨大，仅考虑各向异性程度kr而忽略各向异性方向α情况下最大与最小浸润线滞后高程相差6.51%，而同时考虑各向异性程度kr及方向α情况下最大与最小浸润线滞后高程相差35.53%。对于黏土边坡，渗流各向异性特性影响较小，仅考虑各向异性程度kr而忽略各向异性方向α情况下最大与最小浸润线滞后高程相差4.7%，而同时考虑各向异性程度kr及方向α情况下最大与最小浸润线滞后高程相差9.94%。

图5 不同各向异性程度及方向砂土与黏土边坡的浸润线滞后高程

4.2 各向异性对滑坡不同位置变形的影响

4.2.1 砂土

基于孔压-应力耦合方法，得到了砂土不同各向异性下滑坡不同监测面(上部监测面、中部监测面及下部监测面)的水平位移分布规律，如图6所示。

由图6可见:库水位骤降停止时刻不同滑坡位置的位移变化各不相同，对于上部监测面而言，水平向位移随高程呈现先增大后减小的变化规律，其最大位移发生在距地表约为5 m处。而中部与下部监测面的水平位移则是随着高程逐渐减小。距离库岸越近，监测面的整体位移越大，结合4.1节分析可知，库水位下降靠近库岸岸坡内部水位线滞后，对周围土体骨架产生向临空面及向下的渗流动水压力(动态扩张力)，牵引带动坡面土体[6-8]，因此可以判定库水位下降坡趾最容易发生破坏。不同各向异性特性对库水位停止时的水平位移影响较大。上部监测面位移数值差异较小，各向异性方向角度α越大，下部监测面的整体位移越大，但是中部监测面的位移却越小；各向异性程度kr越大，中部监测面水平位移越小，而下部监测面的水平位移越大。这可能是因为各向异性方位角α越大，各向异性程度kr越大，坡趾位置的渗透坡降逐渐增大，而坡中位置的渗透坡降逐渐减小所致。

图6 不同各向异性特性下砂土水平位移变化

4.2.2 黏土

黏土边坡不同各向异性特性对水平位移影响规律如图7所示。

图7 不同各向异性特性下黏土水平位移变化

由图7可见：不同各向异性特性对黏土边坡的水平位移影响规律与砂土边坡类似，但是值得注意的是，各向异性程度kr及各向异性方向α对黏土边坡的影响比砂土边坡影响要小，体现在不同各向异性程度kr及各向异性方向α下的水平位移差异较小。与砂土边坡对比可知，黏土边坡的上部监测面与中部监测面的水平位移值整体上均小于砂土边坡，而下部监测面的水平位移值却大于砂土边坡，这也是容易理解的，因为黏土边坡的坡降突变集中于坡趾，库水位下降对坡中及坡顶的影响不大，而库水位骤降不仅仅影响到砂土边坡的坡趾，对坡中还具有一定的影响。

4.2.3 最大水平位移分析

为定量化研究不同各向异性特性下砂土与黏土边坡不同位移的最大水平位移，对不同工况的最大水平位移进行了统计，其统计结果如图8所示。

图8 不同各向异性特性下砂土与黏土的最大水平位移。

由图8可见：不同各向异性特性对最大水平位移的影响随监测面位置的不同而不同。下部监测面的位移差异最大，中部监测面次之，而上部监测面则最小。仅考虑各向异性程度kr但忽略各向异性方向α下砂土上部、中部及下部监测面的差异分别为69.5%、20.7%、6.9%，黏土上部、中部及下部监测面的差异分别为75.05%、77.2%、21.4%；同时考虑各向异性程度kr及方向α情况下砂土上部、中部及下部监测面的差异在80.0%、94.6%、45.0%，黏土则为75.1%、87.1%、32.0%。

4.2.4 库水位下降对边坡的渗流稳定性的影响

图9是各向异性程度kr=1，各向异性方向α=0°情况下砂土边坡的库水位动态变化图。图10显示：在库水位下降过程中，土体的持水性使浸润线的下降滞后于库水位，因此对边坡的稳定性影响主要体现在以下三个方面，一是坡体内具有较高的孔隙水压(超孔隙水压力)且不能及时消散，形成非稳态渗流场，对周围土体骨架产生向临空面及向下的渗流动水压力(动态扩张力)，同时由水头差产生的静水压力相当于对前缘施加了一个推力，运动着的地下水“拖拽”土体或土粒向渗流方向前进，这是对边坡稳定的不利方面；二是库水位在较高水平情况下库区水压力较高，对库岸滑坡起到阻滑的作用，但是当库水位下降至低水位时，上游水压力卸载，使得阻滑力降低，此为边坡稳定的不利影响；三是库水位降低导致岸坡内部的孔隙水压力降低，土体的有效应力及强度提高，此为边坡稳定的有利影响。综上所述，库区岸坡稳定性受以上三个因素的综合影响，其中，边坡土质、渗流各向异性特性则是影响边坡稳定的重要因素。

图9 库水位骤降不同时刻浸润线变化

由本文数值模拟可知库水位骤降下坡趾的水平位移最大，这同样也是由于水位线“滞后”导致土体内部指向坡面外的拖拽力导致。本文基于孔压-应力耦合方法得到库水位骤降全过程的水平位移云图(图10a)显示：滑坡坡趾处首先产生位移滑动，随着库水位的进一步下降，坡中及坡顶位移逐渐增大，可见库水位下降滑坡破坏首先发生在坡趾。江强强等[21]对库水位下降下的滑坡模型的室内试验发现，库水位骤降下滑坡的坡趾处首先产生贯穿裂缝，随后发生局部流滑动最终发生滑坡前缘的整体破坏，这与本文的数值模拟的结论是一致的，从而也验证了本文数值模拟及孔压-应力耦合方法的正确与合理性。因此，在实际库区滑坡治理工程中，应该加强滑坡坡趾处的加固措施。

图10 本文数值模拟与文献[21]试验结果对比

4.3 各向异性对滑坡安全系数的影响

4.3.1 砂土

渗流各向异性对砂土边坡的安全系数影响规律如图11所示，由图11可见：安全系数随库水位骤降呈现先减小后增大的规律(除α=45°情况下的kr=10及kr=100)，值得注意的一点是渗流各向异性特性不仅仅对边坡的安全系数大小有影响，同时对最小安全系数的时刻具有影响。渗流各向异性角度α越大，各向异性程度kr越大，边坡的整体安全系数越小，这与4.1节讨论的库水位下降下边坡内部的浸润线滞后效应有关。同时，渗流各向异性角度α及各向异性程度kr的增大同时也使得边坡最小安全系数出现的时刻提前，这是由于库水位下降情况下存在三个作用过程，一个是库岸坡面的水压力卸载，使得岸坡阻滑力降低，此为边坡稳定的不利因素，一个是岸坡内部浸润线滞后，导致水力坡降增大，此亦为不利因素，而另一个则是库岸内部的水位线降低，土体的有效应力及土体强度升高，此为边坡稳定的有利因素，结合4.1.1节分析结果可知，当渗流各向异性角度α及各向异性程度kr较小时，边坡的排渗能力较强，此时库岸坡面的水压力卸载及水位线滞后的影响要小于库岸内部的水位线降低导致土体的有效应力及土体强度升高的影响，随着各向异性角度α及各向异性程度kr逐渐增大，边坡的排渗能力降低，所以最小安全系数出现的时间逐渐滞后。

图11 砂土渗流各向异性对安全系数的影响规律

4.3.2 黏土

渗流各向异性对黏土边坡的安全系数影响规律如图12所示，由图12可见：各向异性程度及方向对黏土边坡的影响规律与砂土边坡类似，但是与砂土边坡相比，渗流各向异性特性对黏土边坡的影响更小。值得注意的是，黏土边坡的安全系数整体上要比砂土边坡更低，这也是由于黏土边坡库水位下降时的浸润线滞后效应更加明显导致。

图12 黏土渗流各向异性对安全系数的影响规律

4.3.3 最小安全系数

渗流各向异性对砂土与黏土边坡的最小安全系数影响规律(图13)显示：随着各向异性程度kr及各向异性方向α的增大，最小安全系数逐渐走向最低。值得注意的是，黏土边坡在库水位骤降下的整体最小安全系数要小于砂土边坡。仅考虑各向异性程度kr但忽略各向异性方向α下砂土边坡的最小安全系数相差2.72%，黏土边坡则相差9.29%，而综合考虑各向异性程度kr及各向异性方向α砂土边坡的最小安全系数相差11.58%，而黏土边坡则为12.31%。