王述红, 何坚, 杨天娇
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-07-10
基金项目：国家自然科学基金资助项目(51474050，U1602232)；辽宁省高等学校优秀人才支持计划项目(LN2014006)。
作者简介：王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。

摘要：运用有限元软件ABAQUS模拟降雨条件下边坡渗流场和应力场耦合, 并运用强度折减法, 以监控点位移突变和边坡形成连续贯通的塑性变形为边坡失稳判据, 采用数值方法计算出耦合后的安全系数, 结合孔压演变分析, 综合评价稳定性变化情况, 从而研究降雨影响边坡稳定性的机理, 并通过含黏土层边坡工程实例进行验证.分析发现:降雨强度越大, 浅层土体形成饱和区的速度越快, 极易发生浅层滑坡, 而黏土层则会加速上述过程, 危害边坡稳定性; 在降雨24 h内, 安全系数降幅最大.该研究结果为降雨条件下边坡事故防治提供了参考和分析依据.
关键词：降雨入渗边坡稳定性有限元法流固耦合强度折减法
Numerical Analysis on Stability of Slope Considering Rainfall Infiltration
WANG Shu-hong, HE Jian, YANG Tian-jiao
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Shu-hong, E-mail: shwangneu@126.com
Abstract: Finite element software ABAQUS was utilized to simulate the coupling effect of seepage field and stress field in slope under rainfall conditions. Regarding the saltation of the displacement at a characteristic nodal and the formation of continuously plastic deformation as instability criterion, The safety factor of slope was numerically computed by strength reduction method to assess the slope stability. By analyzing safety factor and the evolution of pore pressure, how the rainfall infiltration affects slope stability can be explained. And the result was verified by actual slope case including the layer of clay. It was found that the higher the rainfall intensity is, the faster transient saturated zone will form, and the layer of clay can accelerate this process, threatening the stability of slope. Therefore, the shallow landslides may occur. The maximum decline of the safety factor occurs during 24 h rainfall. This research provides good reference to slope monitoring and slope anti-slide measures.
Key words: rainfall infiltrationslope stabilityFEMfluid-structure interactionstrength reduction method
降雨入渗是诱发边坡失稳最显著和最频繁的环境影响因素^[1].受气候和环境因素影响, 大部分边坡土体是非饱和土, 在雨水渗入条件下, 土体内部产生渗流场, 导致土体应力场改变, 两场彼此影响.大部分对降雨条件下边坡稳定性的研究主要讨论渗流场作用, 关于两场耦合的研究较少^[2-4].文献[5]采用基于Biot固结理论的直接耦合计算方式实现了两场耦合运算, 但没有进一步验证结果的准确性.随着计算机和软件技术不断成熟, 通过现场采集的数据建立模型, 并在隧道、边坡、岩体破碎等工程中得到了较好的运用和验证^[6-8], 数值分析结合工程实际的方法获得广泛认可.文献[1]采集实际工程数据, 运用ABAQUS软件, 结合强度折减法, 获得边坡安全系数, 并将这种方法与广泛使用的极限平衡法对比, 验证了采用ABAQUS模拟边坡受降雨影响的方法具有较高的准确性和参考价值.
在此基础上, 本文通过降雨入渗条件下边坡内瞬态渗流场和应力场的直接耦合运算, 获得耦合作用后边坡的应力和应变状态, 后续使用强度折减法计算边坡安全系数, 并结合坡内孔隙水压力演变，综合研究边坡稳定性，最终通过常见的含有黏土层的边坡实例验证了分析的准确性.
1 流固耦合与强度折减法1.1 流固耦合原理1.1.1 流固耦合影响机理雨水渗入到土质边坡形成的渗流场会随着降雨不断变化, 渗流场产生的渗流力作用在土体上引起应力场变化, 导致土体变形, 土内部孔隙体积变小, 孔隙率降低, 雨水在孔隙中渗流作用减弱, 从而又反过来影响渗流场; 渗流场以改变渗透能力的方式使其与应力场不断相互作用.渗流体积力是渗流场产生的, 也是影响应力场的主要因素.根据水力学原理, 渗流体积力会随着水力梯度成正比例变化, 单元的渗流体积力F_V为^[5]

(1)

式中:F_x, F_y是F_V在x和y方向上的分量; J_x, J_y是在x, y方向上的渗透坡降; r_w为水的重度.渗流场改变土体孔隙, 孔隙率对渗透系数影响显著, 在没有试验测定相关性时, 常采用A.Riverra等人提出的渗透系数随孔隙率变化经验公式^[5]:

(2)

式中:n₀为土体初始孔隙率; n为变化后的孔隙率; k, k₀是与n, n₀对应的渗透系数.
1.1.2 ABAQUS中流固耦合方程根据虚功原理和质量守恒定理分别获得渗流场与应力场耦合的平衡方程和渗流连续性方程, 再通过有限元离散, 结合孔压边界条件与流量边界条件, 使用Galerkin解法, 联立各式获得流固耦合方程, 即为文献[9]中所推导的流固耦合方程:

(3)

式中:K, C, E, G, 均为矩阵, t为时间; u为单元节点位移; p_w单元节点孔隙水压力.
1.2 非饱和土强度理论非饱和土的孔隙中含有气体, 无法应用传统的摩尔-库伦强度理论.在工程中应用最广泛的非饱和土强度理论是由Fredlund等提出^[10]:

(4)

式中:φ′为有效内摩擦角; c′为有效黏聚力; p_w为孔隙水压力；p_a为孔隙气压力；φ^b为随基质吸力(p_a－p_w)变化的内摩擦角, 为方便分析, 此处设为常数; τ_f为非饱和土抗剪强度.由式(4)可知, 基质吸力减小会引起土体强度大幅下降, 因此基质吸力是研究土体强度和边坡失稳重要参考因素和指标.由于坡面直接接触空气, p_a=0, 故负的孔隙水压力可以表示基质吸力.
1.3 强度折减法强度折减法能在有限元分析中真实模拟边坡失稳和塑性变形过程, 因此本文用其求解瞬态渗流场和应力场耦合作用后边坡的安全系数, 分析边坡整体稳定性.在外荷载保持不变的情况下, 折减系数F_r对边坡整体土体的c′与φ′进行折减运算, 当边坡因强度折减至下滑力等于土体抗剪强度时, 即边坡处于失稳临界状态, 此时的F_r可视为边坡安全系数.经过折减的抗剪强度参数为^[11]

(5)

(6)

在土体强度折减过程中, 判断边坡失稳临界状态的判据有多种, 其中边坡塑性贯通区形成并且监控点位移发生突变是最准确的边坡失稳判据^[12-14], 算例均采用上述判据.虽然安全系数能较好评价边坡整体稳定性, 但无法反映边坡稳定性受降雨影响的机理；而分析孔隙水压力或饱和度等, 能解释稳定性演化过程和分析边坡局部稳定性, 如坡面的稳定性.因此本文将两者结合, 共同评价边坡稳定性.
2 边坡稳定性影响因素分析2.1 计算模型及计算方案运用ABAQUS进行有限元数值计算的二维模型算例如图 1所示.土体的材料属性:干密度1.3 kg/m³, 黏聚力15 kPa, 有效内摩擦角30°, 弹性模量100 MPa, 泊松比0.3, 渗透系数5×10^-6 m/s, 初始孔隙比1.0.本构模型使用摩尔-库伦准则.EF, AB为不透水边界, 边坡CD面和BC面的雨水入渗强度分别为20和15.62 mm/h.网格单元采用CPE4P单元, 土体涉及渗透系数的变化特征采用文献[15]研究数据.
图 1(Fig. 1)

图 1 模型尺寸及地下水位(单位:m)Fig.1 Model geometry and the water level(unit:m)

2.2 降雨强度和降雨持续时间的影响为研究边坡稳定性在不同降雨强度下的变化, 分别设置降雨强度为10，15，20 mm/h三种方案, 并采用强度折减法数值计算出边坡的安全系数, 得到图 2所示曲线.分析可知:边坡安全系数均随着降雨持续而降低, 降雨24 h后安全系数的降幅最大; 安全系数的降幅会随降雨强度增大而增大.
图 2(Fig. 2)

图 2 降雨强度对边坡安全系数的影响Fig.2 Effect of rainfall intensity on safety factors of slope

为说明稳定性削弱机理, 选取边坡模型右端顶点为监控点, 获得浅层土体监控点基质吸力变化曲线如图 3所示, 与图 2的安全系数变化趋势对应.降雨初期24 h内, 由于地表水力传导率较大, 雨水全部入渗, 降雨强度越大, 饱和区在浅层形成的速度越快, 基质吸力迅速降低, 土体抗剪强度削弱, 严重降低边坡整体稳定性, 安全系数降幅最大; 降雨持续24 h之后, 浅层土体达到饱和, 基质吸力减小的速度减缓, 安全系数降幅逐渐减小.因此, 降雨对边坡浅层首先产生危害, 逐渐形成的饱和区将降低边坡的稳定性.
图 3(Fig. 3)

图 3 浅层土体监控点基质吸力变化曲线Fig.3 Curve of the matric suction on controlled nodal in narrow area of slope

边坡浅层土体最易受降雨影响, 分析孔隙水压力演变可以较好评价其稳定性.设置10，15，20 mm/h三种降雨强度, 每种降雨强度分别持续24, 48, 72, 96 h, 共有12种降雨方案, 由此获得边坡土体孔压在降雨过程中的演化曲线, 如图 4所示，可进一步研究边坡在降雨强度和降雨时长两种因素共同作用下的稳定性.
图 4(Fig. 4)

图 4 孔隙水压力与坡体深度的关系Fig.4 Relationship between pore pressure and slope depth

分析发现边坡浅层土体孔隙水压力受影响程度:雨强20 mm/h, 持续24 h >雨强15 mm/h, 持续96 h; 雨强15 mm/h, 持续24 h>雨强10 mm/h, 持续96 h.即对比低强度长时间降雨，短时强降雨雨水的入渗量更多, 浅层土体含水率更高, 土体重度增大, 同时由于土体抗剪强度降低, 极易产生浅层滑坡.由此可见短时强降雨造成的危害更大, 并且更加迅速.综合上述分析可知, 边坡浅层土体既受降雨强度影响，也受降雨时长影响, 而降雨强度的影响更大.
3 工程实例3.1 实例计算方案选取某含不透水黏土层基坑边坡进行实例分析, 基坑边坡底距地面7 m, 坡度比1:1.5, 地下水位线在地下13 m处.根据当地气象局数据, 选取2016年梅雨季节最大日降雨强度10 mm/h进行分析.基坑有三层土, 底层为不透水的黏土层, 中间层为粉质黏土, 最上层为杂填土, 各土层参数见表 1.各层土体渗透系数采用经验公式(2)计算确定.
表 1(Table 1)

表 1 土体参数表Table 1 Parameters of soil 参数 杂填土 粉质黏土 黏土 干密度/(kg·m-3) 1.340 6 1.399 8 1.533 0 黏聚力/kPa 10 16 31 内摩擦角/(°) 8 9 9 弹性模量/MPa 3 5 9 渗透系数/(m·s-1) 2.31×10-6 1.16×10-8 5.79×10-9 孔隙率 0.97 0.953 0.793 层厚/m 1.7 4.3 6

表 1 土体参数表 Table 1 Parameters of soil

3.2 边坡雨水入渗降雨0.5 h在沿坡面方向距坡脚2.44 m处, 即粉质黏土与黏土在坡面上的交界处, 产生正孔隙水压力, 高达228 kPa, 如图 5中红色区域所示, 而红色区域所在土层的平均孔压为-3.57 kPa.这说明在降雨过程中雨水逐渐入渗到粉质黏土层底部, 由于黏土层不透水, 雨水在黏土层上部聚集, 孔隙水压力在粉质黏土层底部与坡面交界处迅速增大.对应地, 降雨0.5 h后边坡安全系数由降雨前的1.95变为1.88, 降低了3.4%, 降低速度较快, 符合2.2节所述降雨初期浅层孔压迅速增大, 基质吸力降低, 边坡安全系数下降幅度较大的特征.
图 5(Fig. 5)

图 5 降雨0.5 h后孔压分布图Fig.5 Pore pressure distribution after 0.5 h rainfall

为进一步说明此现象, 选取降雨1 h后的孔压云图, 见图 6.图中两处红色区域, 即沿坡面方向距坡脚2.27 m处以及3.72 m处, 孔隙水压分别高达261 kPa和181 kPa.说明滞留在黏土层上的雨水以及沿坡面入渗的雨水在粉质黏土层底部与坡面交接处逐渐聚集, 并沿坡面往坡顶方向发展, 逐步形成饱和区, 即图 7中粉质黏土层所在坡面位置形成的橙红色区域.随着降雨持续, 饱和区沿坡面扩展, 边坡稳定性将持续降低, 安全系数逐渐减小.
图 6(Fig. 6)

图 6 降雨1 h后孔压分布图Fig.6 Pore pressure distribution after 1 h rainfall

图 7(Fig. 7)

图 7 降雨1 h后饱和度分布图Fig.7 Saturation distribution after 1 h rainfall

3.3 边坡浅层土体破坏降雨1 h后, 孔隙水压力为261 kPa的土体部位已发生较大塑性变形, 见图 8中出现的塑性变形区域.该区域因雨水软化, 导致土体强度不足, 而孔隙水压力增大, 发生了较大塑性变形.
图 8(Fig. 8)

图 8 降雨1 h后边坡塑性应变Fig.8 Plastic strain of slope after 1 h rainfall

在实际状态下, 受多种因素影响, 雨水还会沿粉质黏土层所在坡面位置渗出, 内部土体颗粒被冲刷出坡体, 土体易发生冲蚀或滑流.因此, 该区域土体可能发生破坏和一定程度滑移.随着降雨持续, 经过4.5 h，在黏土层上部, 坡面浅层土体发生塑性破坏并沿坡面向上发展, 最终形成连续塑性破坏面, 见图 9.由于该区域土体破坏, 其上部土体缺乏支撑, 坡面会产生浅层滑坡, 验证了2.2节中短时强降雨易发生浅层土体滑坡的结论.
图 9(Fig. 9)

图 9 降雨4.5 h后坡面浅层形成连续塑性破坏面Fig.9 Continue plastic failure surface after 4.5 h rainfall

3.4 不透水黏土层的影响经过强度折减的边坡因黏土层上部土体产生连续贯通的滑动面而失稳破坏, 见图 10中形成的滑动带.该过程中黏土层的变形不大, 而黏土层上部透水土体受强度削减首先破坏.结合渗流场和边坡局部塑性破坏的演化过程分析, 边坡底部不透水的黏土层受降雨影响并不明显, 但由于黏土层阻挡入渗雨水补给地下水, 使得水不易在坡体内消散, 雨水逐渐聚集在坡面和边坡内部, 加速上部透水土层形成饱和区, 威胁边坡安全.
图 10(Fig. 10)

图 10 边坡形成连续贯通塑性变形Fig.10 Perforation of plastic zone in slope

综合上述分析, 该基坑边坡对强降雨较敏感, 短时间内易发生浅层滑坡破坏.因此该边坡应在降雨前进行加固处理, 在梅雨季节施工时, 应警惕降雨, 及时采取相应的基坑降水措施.
4 结论1) 对比降雨时长, 降雨强度对边坡稳定性影响更大; 短时强降雨可迅速在边坡浅层土体形成饱和区, 并逐步从坡脚往坡顶延伸, 极易造成边坡浅层滑坡事故.
2) 最大的安全系数降幅发生在降雨初期24 h内, 24 h后边坡的安全系数降幅会随着降雨持续逐渐减小.
3) 边坡内不透水的黏土层在降雨过程中所受影响不大, 但黏土层的存在会加速上层透水土体形成饱和区, 影响边坡稳定性.
参考文献

[1]	李炎隆, 陈波, 马成成, 等. 基于ABAQUS的降雨入渗条件下基坑边坡稳定性分析[J].应用力学学报, 2017, 34(1): 155–204. ( Li Yan-long, Chen Bo, Ma Cheng-cheng, et al. Study on the stability of the foundation pit slope under rainfall infiltration conditions based on ABAQUS[J].Chinese Journal of Applied Mechanics, 2017, 34(1): 155–204.)
[2]	Liu Z Z, Yan Z X, Duan J, et al. Infiltration regulation and stability analysis of soil slope under sustained and small intensity rainfall[J].Journal of Central South University, 2013, 20(9): 2519–2527.DOI:10.1007/s11771-013-1764-8
[3]	Lan H X, Zhou C H, Lee C F, et al. Rainfall-induced landslide stability analysis in response to transient pore pressure—a case study of natural terrain landslide in Hong Kong[J].Science in China(Series E:Technological Sciences), 2004, 47(1): 1–8.DOI:10.1360/03ye0074
[4]	Collins B D, Znidarcic D. Stability analyses of rainfall induced landslides[J].Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(4): 362–372.DOI:10.1061/(ASCE)1090-0241(2004)130:4(362)
[5]	陈丽刚.基于ABAQUS渗流与应力耦合作用的边坡稳定性分析[D].郑州: 郑州大学, 2010. ( Chen Li-gang.The slope stability analysis of taking into account the seepage and stress coupling effect based on ABAQUS [D].Zhengzhou: Zhengzhou University, 2010.http://cdmd.cnki.com.cn/Article/CDMD-10459-1011017802.htm)
[6]	Wang S H, Ni P P, Guo M D. Spatial characterization of joint planes and stability analysis of tunnel blocks[J].Tunnelling and Underground Space Technology, 2013, 38(1): 357–367.
[7]	Wang S H, Huang R Q, Ni P P, et al. Fracture behavior of intact rock using acoustic emission:experimental observation and realistic modeling[J].Geotechnical Testing Journal, 2013, 36(6): 903–914.
[8]	Wang S H, Ni P P. Application of block theory modeling on spatial block topological identification to rock slope stability analysis[J].International Journal of Computational Methods, 2014, 11(1): 1350044(1-24).
[9]	林子杨.流固耦合下的边坡稳定分析[D].郑州: 郑州大学, 2012. ( Ling Zi-yang.The slope stability analysis of taking into account the seepage and stress coupling effect[D].Zhengzhou: Zhengzhou University, 2012.)
[10]	Fredlund D G, Rahardjo H. Soil mechanics for unsaturated soil[M]. New York: John Wiley & Sons, 1993.
[11]	Zienkiewicz O C, Humpheson C, Lewis R W. Associated and non-associated visco-plasticity and plasticity in soil mechanics[J].Geotechique, 1975, 25(4): 671–689.DOI:10.1680/geot.1975.25.4.671
[12]	刘金龙, 栾茂田, 赵少飞, 等. 关于强度折减有限元方法中边坡失稳判据的讨论[J].岩土力学, 2005, 26(8): 1345–1348. ( Liu Jin-long, Luan Mao-tian, Zhao Shao-fei, et al. Discussion on criteria for evaluating stability of slope in elastoplastic FEM based on shear strength reduction technique[J].Rock and Soil Mechanics, 2005, 26(8): 1345–1348.DOI:10.3969/j.issn.1000-7598.2005.08.035)
[13]	Liu Y, Wang C, Yang Q. Stability analysis of soil slope based on deformation reinforcement theory[J].Finite Elements in Analysis and Design, 2012, 58: 10–19.DOI:10.1016/j.finel.2012.04.003
[14]	Huang M S, Jia C Q. Strength reduction FEM in stability analysis of soil slopes subjected to transient unsaturated seepage[J].Computers and Geotechnics, 2009, 36(1/2): 93–101.
[15]	费康, 彭劼. ABAQUS岩土工程实例详解[M]. 北京: 人民邮电出版社, 2017. ( Fei Kang, Peng Jie. Application of ABAQUS in geotechnical engineering[M]. Beijing: Posts and Telecom Press, 2017.)