降雨条件下北川斜坡破坏的数值模拟

王淑云，叶天力，鲁晓兵，聂雪媛

(中国科学院力学研究所流固耦合系统力学重点实验室，北京 100190)

0 引言

降雨是诱发斜坡破坏的主要因素之一。降雨所致的泥石流和滑坡灾害，会严重威胁人民生命安全和带来经济损失。近几年来，我国每年都会有较大的泥石流发生，尤其在滇西北和东北、川西、甘南及白龙江流域等地区。典型的例子是在2008年9月25日，四川省北川县城，因持续降雨引发多处山体滑坡和泥石流，30多人死亡或失踪，老县城一半以上被泥石流掩埋，最深处的污泥达40m以上。2010年9月9日，北川再遭泥石流冲击，桥梁被冲毁，雷禹路全线瘫痪。2013年7月16日，北川的桂溪乡陈家坝因暴雨突发泥石流，50余万m3泥石流冲入青林沟，将河床整体抬高近4m，形成堰塞湖。

降雨诱发泥石流是雨水不断入渗至地表土层内部的结果。在降雨过程中，地表土层的饱和度增加，孔压上升，基质吸力减少，抗剪强度下降。当降雨达到一定时间，雨水入渗到土层的深度增加，土体强度大大降低，土体失稳，泥石流起动。目前，在降雨引发滑坡和斜坡破坏研究方面，已有不少文献报道［1-9］。但因为每个地区的降雨强度和历时、现场地形和坡度以及土质情况的复杂性，对其在降雨下的斜坡破坏的模拟和分析结果会有所不同，尤其反应在定量方面。因此对特定地区土的斜坡破坏的分析和预测必须建立在其所处的外界环境和土本身的特性基础上。

本文采用FLAC两相流模块，针对北川地区地形(主要是坡角)和土质(现场土)，对北川地区土在不同降雨强度和降雨历时情况下的斜坡破坏进行了数值模拟，并对降雨过程中土体位移、速度矢量等变化规律进行了分析。同时开展了室内降雨模型实验，并将数值模拟计算结果与室内实验结果进行比较，验证了数值模拟计算和分析的可信度和可行性。

1 数值模拟和模型试验

1.1 数值模拟介绍

土体模型采用修正的Mohr-Columb模型，二维的饱和-非饱和土中水的渗流控制方程采用Richards方程。有限差分网格如图1所示，其中X轴为沿土坡方向，Y轴为垂直于坡面方向。模型形状为梯形体型。坡顶350mm，坡底520mm，坡面土层厚度100mm。在数值模拟中，土坡坡角分别取27°、30°、33°，本文分析以坡角27°时的数值计算结果为主。

图1 有限差分网格Fig.1 Element mesh

北川土的主要土性参数为:干容重1.58g/cm3，粘聚力20kPa，体积模量800MPa，剪切模量600MPa，渗透系数1.0×10-5m/s。土体初始饱和度设为0.5，并假定残余饱和度为0。降雨参数的选取基于北川当地降雨情况，分别采用5种雨强:1×10-5m/s、2×10-5m/s、3 ×10-5m/s、4 ×10-5m/s和 5 ×10-5m/s。

1.2 室内模型实验介绍

模型实验简图如图2所示，土坡坡度为27°，均匀铺土，土层长 ×宽 ×厚为1000mm×1200mm×100mm，土的颗分结果参见图3。

图2 模型实验示意图Fig.2 Sketch of model test

图3 北川土的颗分曲线Fig.3 Grain size distribution

实验选取的雨强分别为:1.75×10-5m/s、2.25 ×10-5m/s、3.25 × 10-5m/s、3.5 × 10-5m/s、4 ×10-5m/s、5 ×10-5m/s，记录降雨历时和斜坡破坏时间。

2 数值模拟和模型试验结果分析

2.1 数值模拟结果分析

2.1.1 土体水平位移分析

图4为雨强1×10-5m/s时土体水平方向位移降雨历时的变化情况。从图4可以看出，小雨强情况下，降雨初期土体的水平位移相对较小，最大水平位移发生在坡脚处，这与室内试验现象是比较一致，在试验中发现，在降雨初期，坡脚处的土体首先发生坍塌滑动。数值模拟结果证实坡脚处的这种坍塌滑动主要是由土体的水平方向的位移引起的。土体水平方向的位移主要发生在模型左侧，如果取一沿坡面方向的截面，则自左至右土体水平位移呈逐渐减小的趋势，最右侧土体的水平位移很小甚至为零。如果在坡面左侧部分取一竖直截面进行分析，则土体的水平位移自上而下先增加后减小，在该截面中部的某一点土体水平位移取得最大值，且随着该截面的向右移动，最大值点的位置将不断的向上移动，这些最大值点构成了一条线，这条线所处的位置即为潜在的滑移带。

图4 土体水平方向(沿坡面方向)位移分布图Fig.4 Distribution of horizontal displacement

2.1.2 土体竖直位移分析

图5为雨强1×10-5m/s时土体竖直方向位移随降雨历时的变化情况。从图5可以看出，小雨强情况下，土体发生了不均匀沉降。土体的竖向位移随降雨历时的增加不断增大，且土体上层的竖向位移远超过下层土体的竖向位移。土体竖向位移最大值出现在土体表层，且随降雨历时的增加，土体竖向位移的最大值位置逐渐向左侧移动，最终土体竖向位移的最大值点大约位于距坡顶三分之二处。结合前文对土体水平位移的分析，可以得出，土体之所以出现不均匀的沉降主要是由于土体发生了水平位移引起的。随降雨历时的增加，左侧土体不断的向左下方滑动，土体总体积是一定的，这就引起土体的不均匀沉降，出现图中的现象。正是由于土体竖向位移的不均匀性，随竖向位移最值点的位移值不断增大，该点土体将与上部土体产生错位，这为滑移带的产生在表层形成了突破口。

图5 土体竖向方向(垂直坡面方向)位移分布图Fig.5 Distribution of vertical displacement

2.1.3 土体速度矢量分析

图6为雨强1×10-5m/s情况下土体速度矢量随降雨历时的变化情况。从图6可以看出，在小雨强情况下，随降雨历时的增加，土体的速度矢量不断增大，且在降雨初期，土体速度增加的速率比较小，但随着降雨的持续，增加的速率也越来越大。从图7可以看到，土体速度矢量在水平方向分量的最大值位于坡脚处，这与前文分析中得到的土体的水平最大位移发生与坡脚处是一致的。土体速度矢量在竖直方向分量的最大值大约在坡体表层距坡顶大约三分之二处，这与前文的分析结果也是一致的。从整个坡体来看，速度矢量在竖直方向上的分量明显小于其在水平方向上的分量，也就是说土体的水平运动处于主导地位。从土体产生较大运动的面积上来看，随降雨历时的增加，面积逐渐减小，且产生较大运动的区域不断向左侧偏移。

图6 土体速度矢量与降雨历时关系Fig.6 Distribution of velocity vector

图7所示为土体速度矢量最大值随降雨历时变化曲线。可以看到，土体速度矢量的最大值也随着降雨历时的增加而增大，降雨历时越长，最大速度矢量值增长的越快。

图7 土体速度矢量最大值随降雨历时变化曲线Fig.7 Maximum velocity vector versus rainfall duration

2.2 斜坡破坏及与室内试验结果的分析和比较

图8、9为四种雨强下土体水平向最大位移以及竖向最大位移随时间变化的关系曲线。可以看出:相同雨强条件下，整个土体位移场的最大位移值随时间不断增大，在降雨初期土体位移最大值比较小，当降雨达到一定时间土体位移的最大值将有一突变。随雨强增大这种突变所需要的时间越短。这一现象也在坡脚30°和33°情况中得到了证实。

图8 不同雨强下土体水平向最大位移变化曲线Fig.8 Maximum horizontal displacement under different rainfall density conditions

图9 不同雨强下土体竖向最大位移变化曲线Fig.9 Maximum vertical displacement under different rainfall density conditions

对于坡角的影响，选取雨强为4×10-5m/s情况，针对坡脚27°、30°和33°进行模拟计算。图10、11 为不同坡角下土体水平向最大位移以及竖向最大位移随时间的变化曲线。从图中可以看出，相同雨强下，随坡角的增大，土体最大水平位移和最大竖向位移发生突变所需时间大大缩短。

图10 不同坡角下土体水平最大位移曲线Fig.10 Maximum horizontal displacement under different slope angle conditions

本文模拟主要采用坡面特征点的位移突变作为泥石流是否起动的判据，因为竖向最大位移发生在距坡顶三分之二处，该处是土体滑移带与坡面的交汇点，所以我们取如图1所示的A点作为特征点。

从图12中可以看出，在坡角27°时，雨强越小，斜坡破坏所需时间越长;当雨强从较大开始变小时，斜坡破坏所需时间急剧减小，随着雨强的进一步减小，斜坡破坏所需时间减小的速度逐渐变小，最终斜坡破坏时间趋于一个很大值。相同雨强下，坡角越大，斜坡破坏所需时间越短。这些规律跟室内的人工降雨试验是一致的。从图13可以看出，在坡角为27°时，A点水平速度随降雨历时的变化趋势与图12相同，降雨初期速度比较小，当降雨达到一定时间速度也发生突变.从曲线还可以看出，速度突变时间超前于相应的位移突变时间。

图11 不同坡角下土体竖向最大位移曲线Fig.11 Maximum vertical displacement under different slope angle conditions

图12 坡角27°时A点水平位移曲线Fig.12 Horizontal displacement of A when slope angle is 27°

图13 坡角27°时A点水平速度曲线Fig.13 orizontal velocity of A when slope angle is 27°

以特征点位移是否发生突变作为泥石流起动与否的判断标准，从相同坡角不同雨强情况下的特征点的位移曲线，可以得到不同雨强下所对应的斜坡破坏所需的降雨历时，并将其与试验所得结果进行对比，结果如图14所示。可以看到，数值计算得到的斜坡破坏-降雨历时曲线与试验得到的结果规律相同。即在坡角一定时，当雨强较小，斜坡破坏所需时间较长，随雨强的不断增大，斜坡破坏所需时间越来越短，且所需时间减小的速率越来越大。

图14 斜坡破坏时间随雨强变化关系(室内试验和数值模拟结果对比)Fig.14 Slope failure time versus rainfall duration

3 结论

初步结论如下:

(1)土体水平方向最大位移发生在坡脚处，且土体水平方向位移随降雨历时的增加不断增大。在降雨初期土体位移较小，当降雨达到一定时间，土体的位移将有一突变。无论是坡脚处的坍塌，还是斜坡破坏过程中滑移带的形成，土体水平方向的位移都起到了至关重要的作用。

(2)土体竖向位移具有不均匀性，在大约距坡顶三分之二处取得最大值，这种不均匀性主要是由于土体发生水平方向的运动引起的，也正是因为这种不均匀性，为斜坡破坏在表层形成了突破口。

(3)可以采用特征点位移是否发生突变作为斜坡破坏与否的判断标准，其数值模拟结果与室内模型试验结果吻合较好。

［1］Iverson R M.The physics of debris-flows［J］.Reviews of Geophysics，1997，35(3):245-296.

［2］Klubertanz G.，Laloui，L.，Vulliet，L.，Identification of mechanisms for landslide type initiation of debris flows［J］.Workshop on Mechanics and Velocity of Large Landslides，Engineering Geology，2009，109(1-2):114-123.

［3］Tang，C.，Zhu，J.，Li，W.L.，Rainfall-triggered debris flows following the Wenchuan earthquake［J］.Bulletin of Engineering Geology and the Environment，2009，68:187-194.

［4］CUI Peng，HU Kaiheng，ZHANG Jianqiang et al..Prediction of debris flow area by combing hydrological and inundation simulation methods［J］.Journal of Mountain Science，2011，8(1):1-9.

［5］文海家，张永兴，柳源.滑坡预报国内外研究动态及发展趋势［J］.中国地质灾害与防治学报，2004，15(1):1-4.WEN Haijia，ZHANG Yongxing，LIU Yuan.A home and abroad investigative-tendency of landslide forecast［J］.The Chinese Journal of Geological Hazard and Control，2004，15(1):1-4.

［6］李海亮，黄润秋，吴礼舟，等.非均值土坡降雨入渗的耦合过程及稳定性分析［J］.水文地质与工程地质，2013，40(4):70-76.LI Hailiang， HUANG Runqiu， WU Lizhou， et al.Rainfall infiltration coupling process and stability analysis ofa heterogeneous soilslope［J］.Hydrogeology＆Engineering Geology，2013，40(4):70-76.

［7］Guzzetti，F.，Peruccacci，S.，Rossi，M..The rainfall intensity-duration control of shallow landslides and debris flows:an update［J］.Lanslides 2008，5(1):3-17.

［8］Berti，M..Simoni，A..Experimental evidences and numerical modelling of debris flow initiated by channel runoff［J］.Landslides，2005，2(3):171-182.

［9］Tohari，A.，Nishigaki，M..Komatsu，M..Laboratory rainfall-induced slope failure with moisture content measurement［J］.Journal of Geotechnical and Geoenvironmental Engineering， 2007， 133(5):575-587.

［10］Major J J，Iverson R M.Debris-flow deposition:effects of pore-fluid pressure and friction concentrated atflow margins［J］.Geological Society of America Bulletin，1999，111(10):1424-1434.

［11］戚国庆，黄润秋.泥石流成因机理的非饱和土力学理论研究［J］.中国地质灾害与防治学报，2003，14(3):4-9.QI Guoqing，HUANG Runqiu.Study on genetic and mechanical analysis of debris flow based on unsaturated soils mechanics［J］.The Chinese Journal of Geological Hazard and Control，2003，14(3):4-9.

［12］陈小念，孙海洋.降雨条件下饱和-非饱和土坡的渗流分析［J］.科学技术与工程，2007，7(24):115-118.CHEN Xiaonian，SUN Haiyang.Saturated-unsaturated seepage analysis of slope under rainfall［J］.Science Technology and Engineering，2007，7(24):115-118.

［13］陈守义.考虑入渗和蒸发影响的土坡稳定性分析方法［J］.岩土力学，1997，18(2):8-12.CHEN Shouyi.A method of stability analysis taken effects of infiltration and evaporation into consideration for soil slopes［J］.Rock and Soil Mechanics，1997，18(2):8-12.

［14］朱文彬，等.降雨条件下土体滑坡的有限元数值分析［J］.岩石力学与工程学报，2002，21(4):509-512.ZHU Wenbin，LIU Baochen.Forming and development process of soil landslide during rainfall［J］.Chinese Journal of Rock Mechanics and Engineering，2002，21(4):509-512.

［15］Sun，Y.，Sakajo S and Nishimura ，M..Appliaction research on a numerical model of two-phase flow in deformation porous medium［C］.Int.Conf.Computer Methods and Advances in Geomech，Wuhan，1997，Ⅱ:1171-1176.

［16］Sun，Y.，Nishigaki，M.and Kohno，I..A study on stability analysis of shallow layer slope due to raining permeation［C］.Proceedings of the 1stInt.Conf.on unsaturated soils.Rotterdam Balkema， 1995:1135-1141.