渗流-应力耦合效应对非饱和土边坡稳定性影响分析

刘 俊 ，孙树林,2 ，刘 博 ，陈怿旸 ，王恩喜

(1.河海大学 地球科学与工程学院，江苏 南京 211100；2.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)



渗流-应力耦合效应对非饱和土边坡稳定性影响分析

刘 俊1，孙树林1,2，刘 博1，陈怿旸1，王恩喜1

(1.河海大学 地球科学与工程学院，江苏 南京 211100；2.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)

基于VG模型得到土-水特征曲线与渗透系数函数曲线，通过建立有限元数值分析模型，采用有限元软件SIGAM/W与SEEP/W分别进行渗流-应力耦合分析与只考虑渗流的非耦合分析，将得到的渗流场与应力场导入SLOPE/W计算安全系数并对土坡内应力应变进行分析，研究耦合效应对边坡稳定性影响的相关规律。结果表明，渗流-应力耦合效应使得土体中雨水下渗速度更快，土中基质吸力消散更快；降雨对土坡稳定有不利影响，耦合分析土坡安全系数比非耦合分析小，耦合效应对非饱和土边坡稳定性影响不可忽略；降雨入渗条件下，土坡浅层出现应力集中现象，湿润锋深度处出现应力集中带，土坡易发生浅层破坏；坡脚水平位移比坡顶大，土坡失稳多从坡脚开始。

非饱和土；渗流；渗流-应力耦合；边坡稳定

由于非饱和土中存在基质吸力，与饱和土相比，增加了土的抗剪强度，提高了土边坡的稳定性。降雨入渗过程中，雨水通过土中孔隙不断下渗，土体从坡面开始逐渐趋于饱和，在地下水位以上形成暂态饱和区，土体基质吸力消失，边坡稳定性降低。在渗流的作用下，土体中的位移场和应力场会发生改变。应力场的改变会使土的孔隙率、孔隙比、渗透系数等参数发生变化，从而影响到水在土中的渗流速度，即渗流场改变[1]。因此，在进行非饱和土边坡稳定分析时，渗流-应力耦合作用是必须要考虑的因素，它对土坡稳定性的影响不可忽略。随着渗流-应力耦合理论的逐渐成熟，土体中渗流与应力的耦合成为了研究热点。高小育[2]建立了渗流场和应力场耦合作用下的微分方程，为今后的研究工作提供了基础；徐晗[3]对非饱和土渗流场和应力场进行有限元计算，得到了非饱和土边坡变形和应力的若干规律；张社荣[4]进行了瞬态渗流条件下饱和-非饱和渗流场和应力场耦合分析，得到渗流-应力耦合作用下破坏接近度评价方法；田东方[5]提出了非饱和渗流与应力耦合的计算方法，并验证了方法的可行性；H. Ghasemzadeh[6]提出了各向同性荷载条件下水力耦合模型，该模型能充分描述加载时土体中吸力与应力之间的相互影响。耦合效应作为影响边坡稳定的重要因素，得到了越来越多学者的重视。如O. Davies[7]通过建立水力耦合模型对边坡稳定性进行分析，发现孔隙水压力是导致边坡失稳的重要因素；王柳江[8]研究了降雨过程中土坡的渗流变形耦合特征，阐明了降入过程中基质吸力与边坡变形之间的相互作用；WANG H.B[9]发现边界条件与渗流-变形耦合效应对非饱和土体的渗流有重要影响。虽然许多学者对土中渗流-应力耦合效应进行了广泛研究，但是关于耦合效应对土坡稳定性影响的研究还不够深入。本文用有限元软件SIGAM/W与SLOPE/W对遭受降雨的非饱和土坡进行渗流-应力耦合分析，用有限元软件SEEP/W与SLOPE/W进行只考虑渗流作用的非耦合分析，将两者计算结果进行对比，通过分析坡内孔隙水压力、安全系数、应力和位移变化，对耦合效应的影响机理进行详细解释。

1　土-水特征曲线与渗透系数函数

在非饱和土力学中，基质吸力是其重点研究对象之一，也是非饱和土性质区别于饱和土的最主要原因。基质吸力与土的含水量密切相关，其大小是含水量的函数，随含水量的变化而变化。基质吸力与土体含水量之间的关系可用土-水特征曲线表示。土-水特征曲线对研究非饱和土的物理力学性质十分重要，利用土-水特征曲线可以确定非饱和土的强度与渗透系数等有关参数[10]。

渗透系数是表征土体渗透性的重要指标。非饱和土体中，部分孔隙被气体占据，导致土的渗透系数减小，其减小的程度与土的含水量有关。非饱和土的含水量是基质吸力的函数，因此其渗透系数也可以表示成基质吸力的函数，即渗透系数函数。

1980年，Van GENUCHTEN提出经典的VG模型[11]，该模型能很好地描述土水特征曲线与渗透系数函数曲线，形式如下：

(1)

K=Ksψ1/2[1-(1-ψ1/m)m]2

(2)

式中：θ为体积含水量；θs为饱和体积含水量；θr为残余体积含水量；h为负孔隙水压力；a,n ,m为曲线拟合参数，n=1-1/(1-m)；K为非饱和土渗透系数；Ks为饱和土渗透系数；ψ=(θ-θr)/(θs-θr)。VanGENUCHTEN模型具有3个拟合参数，拟合效果好，因此被广泛使用。

本文计算中的土-水特征曲线与渗透系数函数曲线见图1。

2　渗流-应力耦合分析基本理论

2.1 非饱和土渗流控制方程

忽略气相的影响，假设非饱和土孔隙中气体与大气相连，即孔隙气压力ua=0。水的流速与水头之间的关系可以用达西定律表示如下：

(3)

式中，vi为i方向的流速；ki为i方向的渗透系数，与饱和度s有关；uw为孔隙水压力；γw为水的重度； y为位置水头。

非饱和土的水流连续方程中假设水不可压缩，不计土颗粒的变形。取土单元体作为研究对象，根据水流连续方程和达西定律，得出二维渗流连续方程为：

(4)

式中t为时间。

用加权余量的Galerkin方法对渗流连续方程进行离散，得到简化的二维有限元渗流方程[12]：

[K]{H}+[M]{H}={Q}

(5)

式中[K] 为单元特征矩阵；{H}为总水头；[M]为单元质量矩阵；{Q}为单元上施加的流量矢量。

2.2 渗流-应力耦合数值分析

在非饱和土有效应力原理的基础上，考虑土单元应力平衡，建立应力平衡方程，利用虚功原理，将渗流方程表示为孔隙水压力与体积应变的有限元形式，采用向后差分的数值分析方法，得到饱和-非饱和固结的有限元方程[13]：

[Kd]{Δδ}=[Ld]{Δuw}={ΔF}

(6)

式中：[Kd]为结构整体刚度矩阵；{Δδ}为位移增量矩阵；[Ld]为结构整体的孔隙水压力刚度矩阵；{Δuw}为孔隙水压力增量矩阵；{ΔF}为结构整体的节点外荷载增量。

式(5)与式(6)联立就构成了非饱和土的渗流-应力耦合模型，加上渗流边界条件和位移边界条件后即可进行有限元求解。

3　算例分析

3.1 计算模型

某边坡高12m，其几何形态见图2。材料类型为弹塑性，其物理力学参数见表1。在降雨入渗条件下，边坡表面会在很短时间内由非饱和变为饱和状态。为了在数值分析上处理这种急剧的边界变化，确保计算结果的准确，SEEP/W会在急剧变化边界创建表层，通过划分精细网格解决问题[14]。表层可在其它模块通用，厚度可以任意确定。本次研究的边坡，假设表层厚度2m，土体性质与下部土体一致。表层被划分为单独一层，其网格为沿垂向节点的四边形网格。因为在降雨入渗过程中，初始的水力梯度在垂直地表的方向上更陡，在表层内使用三角形网格会导致计算结果波动，而且四边形网格能大大减少单元的数目，因此使用四边形网格效果会比三角形网格好。表层使用优化的网格使得计算结果更准确，而下部的网格则可由软件自动划分，以便于计算。

进行渗流-应力耦合分析时，既要确定水力边界条件，也要确定力学边界条件。对于非耦合分析，只需定义水力边界条件即可。假定初始地下水位水平，距离坡底7m，坡内初始孔隙水压力沿高程线性分布，地下水位以上为非饱和区，孔隙水压力为负值。侧向边界上地下水面以下为定水头边界，地下水面以上的为不透水边界，因为地下水面以上的侧向流量对边坡稳定影响不大。边坡底部亦为不透水边界。边坡底部的垂向与水平向位移约束，两个侧向边界允许发生垂向位移。

表1　土的物理力学参数

降雨视为附加于边坡顶部与坡面的流量边界条件。假设坡面的降雨入渗系数为1。降雨的入渗量取决于降雨强度q与饱和土渗透系数ks和土壤允许入渗的容量f之间的大小关系[3]。当qf时，大于f部分的雨水在边坡表面形成径流，此时流量边界将转变成坡面单元节点孔隙水压力为零的压力边界，SEEP/W会自动考虑这种变化。本次研究降雨历时15h，降雨强度恒定为1×10-2m/h，时间步长为3h。饱和土渗透系数为2.0×10-6m/s。

3.2结果分析

(1)孔隙水压力分析

在本次研究时间范围内，由于入渗雨水主要集中在边坡浅层，坡下部土层受降雨入渗影响较小，因此重点对表层内的孔隙水压力变化进行分析。图3表示水平坐标为7 m断面上表层内耦合与非耦合的孔隙水压力变化曲线。

从图3可以看出，耦合分析与非耦合分析的初始孔隙水压力一致，沿高程线性变化，表层内孔隙水压力为负值。降雨开始后，坡面孔隙水压力小于零，这是因为降雨强度与饱和土渗透系数十分接近且小于土壤允许入渗的容量。入渗的雨水使得表层上部土体趋于饱和，孔隙水压力迅速增大，基质吸力逐渐消失。在一定深度时，耦合与非耦合分析中的孔隙水压力曲线都有一段较平缓，然后出现急剧变化，这个深度就是降雨入渗时湿润锋的深度。随着降雨的进行，湿润锋深度不断增加，同一深度位置的孔隙水压力随时间增加而增大。与非耦合分析相比，耦合分析的湿润锋下降速度明显比非耦合分析快，相同时间同一位置耦合分析的孔隙水压力比非耦合分析大，这种差异在湿润锋深度处最为明显，表明耦合分析基质吸力消散速度比非耦合分析快，入渗雨水对基质吸力消散有重大影响。湿润锋位置以下土层内，耦合分析与非耦合分析的孔隙水压力与初始值一致，说明下部土层几乎不受降雨的影响。

(2)安全系数分析

图4为基于VANAPALLI强度理论[15]计算得到的边坡整体安全系数与时间的关系。一种用M-P法计算得到的非耦合分析边坡安全系数，第二种用有限单元应力计算得到的耦合分析边坡安全系数，第三种用有限单元法计算的安全系数，其中应力来自于耦合分析，孔隙水压力来自非耦合分析。

从图4可以看出，在耦合与非耦合分析中，边坡的安全系数随降雨的持续不断减小，非耦合分析得到的安全系数比耦合分析的安全系数大，耦合分析得到的边坡稳定性更低。因为降雨时，雨水从坡面开始往下渗，坡内基质吸力不断消失，土体抗剪强度变低。耦合分析时雨水下渗速度更快，相同时间内基质吸力消散范围更广，使得边坡稳定性更低。非耦合分析时，基质吸力较耦合分析大，土体抗剪强度大，因此用耦合方法计算非耦合条件下的边坡安全系数比正常耦合计算得到的安全系数大。相较于非耦合分析，耦合分析时湿润锋位置更深，上覆土层更厚，滑动力更大且基质吸力更小，因此湿润锋位置处的安全系数也更小。

(3) 应力分析

图5为降雨15h后坡内最大剪应力云图。从图中可以看出，边坡表层应力等值线密集，尤其在表层下部最密，表层内出现应力集中带。结合孔隙水压力变化曲线分析得出的湿润锋位置可知，应力密集带出现在湿润锋深度附近，表明边坡易发生沿湿润锋的浅层破坏。

(4)位移分析

图6为降雨15h后边坡位移云图。从图可以看出，水平位移随着深度的增加而增大，在坡内形成环形等位移圈。表层水平位移从坡顶到坡脚逐渐增大，在坡脚处达到最大值。坡顶表层位移以竖向位移为主，几乎不发生水平位移。坡顶竖向位移为负值，表现为沉降；坡脚竖向位移为正值，表现为隆起。降雨后边坡顶部开始沉降，边坡整体发生水平向的挤压，其中坡脚受挤压最强烈，并伴随有局部隆起，表明坡脚最易发生破坏，这也解释了为什么降雨条件下土坡的破坏多从坡脚开始。

4　结论

1)渗流-应力耦合效应使土坡内雨水下渗速度变快，孔隙水压力消散速度更快，土体强度降低，土坡稳定性变差。

2)耦合分析与非耦合分析得到的安全系数都随时间降低。耦合分析得到的安全系数比非耦合分析小，渗流-应力耦合效应对边坡稳定有不利影响，在进行边坡稳定分析时应该考虑耦合效应的影响。

3)降雨使边坡表层的孔隙水压力迅速变化，在表层形成应力集中区，湿润锋位置处形成应力等值线密集带。在降雨影响下，土坡易发生浅层破坏。

4)降雨条件下坡顶发生沉降，坡脚受挤压发生大的水平位移和隆起，坡脚易发生破坏。在进行土坡安全治理时，尤其要注重坡脚的加固。

[1]柴军瑞,仵彦卿.均质土坝渗流场与应力场耦合分析的数学模型[J]．陕西水力发电．1997，13(3)：4-7.

[2]高小育,廖红建,丁春华.渗流对土质边坡稳定性影响[J]. 岩土力学，2004，25(1):69-72.

[3]徐晗，朱以文，蔡元奇，等.降雨入渗条件下非饱和土边坡稳定分析[J].岩土力学，2005，26(12)：1957-1962.

[4]张社荣，谭尧升，王 超，等.强降雨特性对饱和一非饱和边坡失稳破坏的影响[J].岩石力学与工程学报，2014，33(3)：4102-4122.

[5]田东方，刘德富，王世梅，等.土质边坡非饱和渗流场与应力场耦合数值分析[J].岩土力学，2009，30(3)：810-814.

[6]H.GHASEMZADEH,S.A.GHOREISHIAN AMIRI.A hydro-mechanical elastoplastic model for unsaturated soils under isotropic loading conditions[J].Computers and Geotechnics, 2013,51:91-100.

[7]O.DAVIES,M.ROUAINIA,S.GLENDINNING,et al.Investigation of a pore pressure driven slope failure using a coupled hydro-mechanical model[J].Engineering Geoglogy,2014,178:70-81.

[8]WANG LIUJIANG, LIU SIHONG, FU ZHONGZHI, et al. Coupled hydro-mechanical analysis of slope under rainfall using modified elasto-plastic model for unsaturated soils[J].Journal of Central South University，2015，22：1892-1900.[9]WANG HUABIN, LI JIWEI. Analytical Solutions to the One-Dimensional Coupled Seepage and Deformation of Unsaturated Soils with Arbitrary Nonhomogeneous Boundary Conditions[J].TRANSPORT IN POROUS MEDIA, 2015, 108(2):481-496.

[10]栾茂田，李顺群，杨 庆.非饱和土的理论土-水特征曲线[J].岩土工程学报，2005，27(6)：611-615.

[11]VAN GENUCHTEN M TH. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980,44(5): 892-898.

[12]付宏渊，曾 铃，王桂尧，等．降雨入渗条件下软岩边坡稳定性分析[J].岩土力学，2012，33(8)：2359-2365.

[13]张芳枝，陈晓平.非饱和堤岸的渗流特征及其稳定性研究[J].岩土力学，2011，32(5)：1561-1567.

[14]中仿科技公司.非饱和土体渗流分析软件SEEP/W用户指南[M].北京：冶金工业出版社，2011.

[15]VANAPALLI S K，FREDLUND D G，PUFAHL D E，et al. Model for the prediction of shear strength with respect to soil suction[J].Canadian Geotechnical Journal,1996,33(3):379-392.

(责任编辑李军)

Analysis of Seepage and stress coupling effect on stability of unsaturated soil slope

LIU Jun2, SUN Shulin1,2, LIU Bo2, CHEN Yiyang2, WANG Enxi2

( 1College of Earth Science and Engineering, Hohai University, Nanjing 211100,China ;2State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098,China)

Soil-water characteristic curve and hydraulic conductivity function are determined based on VG model, and then by way of constructing finite element numerically analysis model, both seepage-stress coupling analysis and uncoupling analysis which only considering seepage are conducted by using finite element software SIGAM/W and SEEP/W respectively, factors of safety is calculated by software SLOPE/W using seepage field and stress field from coupling analysis and uncoupling analysis and analyze stress and strain to research the principle of coupling effect on slope stability. The results show that the rate of rain infiltration advance is faster due to seepage-stress coupling effect which results in faster loss of suction compared with uncoupling analysis. Rainfall has negative influence on slope stability ,the factors of safety from coupling analysis is lower than that from uncoupling analysis, thus the influence of coupling effect on slope stability cannot be neglected. Stress concentration appears in surface layer and stress concentration belt appears at the depth of wetting front on the condition of rainfall infiltration, which means that slope is prone to shallow failure. Slope toe is the first region to failure because horizontal displacement in slope toe is higher than that in slope top.

unsaturated soil; seepage; seepage and stress coupling; slope stability

2016-04-25

河海大学水文水资源与水利工程科学国家重点实验室开放研究基金项目(2005408911)；留学回国人员科研启动基金项目(20071108)。

刘 俊(1993-)，男，湖南邵阳人，硕士研究生，主要从事地质灾害研究工作。

1673-9469(2016)03-0038-05

10.3969/j.issn.1673-9469.2016.03.008

TU431

A