坡面泥石流启动机理的细观数值模拟研究

王乐荣， 石 崇， 翟淑花， 张一平， 马金城

（1.河海大学岩土力学与堤坝工程教育部重点实验室，南京 210024；2.河海大学岩土工程科学研究所，南京 210024； 3.北京市地质研究所，北京 100120）

泥石流是指在山区或者其他沟谷深壑、地形险峻的地区，因为暴雨、暴雪或其他自然灾害引发的山体滑坡并携带有大量泥沙以及石块的特殊洪流. 作为一种突发性自然灾害，泥石流暴发时往往会伴随着巨大的人员伤亡和财产损失［1-3］. 对泥石流的启动机理进行研究可为泥石流灾害的预测和预防提供参考.

现有的研究泥石流启动机理的方法主要有三种：现场勘察研究、室内模型试验、数值模拟. Zhao等［4］根据九寨沟泥石流沟的地形特征分析了泥石流灾害的历史和成因，提出了由稳定结构、阻塞结构和沉积结构组成的泥石流综合防治体系. 陈鹏宇和彭祖武［5］应用现场勘察研究的方法分析了栾川县境内泥石流的物源特征及启动模式. 甘建军和储小东［6］通过野外调查及室内模型试验分析了江西省德安县泥石流的发育条件和形成机理，发现影响泥石流形成的主要因素是强降雨、丰富的松散物源和陡峻的地貌，同时利用极限平衡法计算了泥石流的稳定性，并提出了治理的方法和对策. 胡林等［7］通过资料收集、现场调查等手段，计算了泥石流的流速、总量、爬高等动力学参数并按照物源供给对泥石流进行分类. 王丽娟等［8］通过野外调查分析了滑坡型坡面泥石流和冲蚀型坡面泥石流的形态特征和颗粒组成差异，并结合地形陡峭指数和堆积区流动路径分析了两种坡面泥石流的流动特征. 周健等［9］通过室内模型试验研究了不同降雨强度条件下细粒含量对泥石流形态的影响. 胡言鹏和刘恩龙［10］通过室内模型试验初步探究了降雨特性对黏性泥石流启动过程的影响. 高波等［11］通过室内模型试验探究了不同规模的泥石流启动时所需的降雨量. 张登等［12］通过室内模型试验研究了泥石流启动时间与坡度、降雨强度之间的关系. 张永双等［13］通过室内模型试验分析了暴雨作用下宁陕城区高陡斜坡带坡面型泥石流的发育特征和形成过程，探讨了泥石流启动的临界状态. 随着计算机技术的发展，越来越多的学者［14-16］开始采用数值模拟方法对泥石流的启动机理进行研究. 姚志雄［17］通过室内模型试验和有限元数值模拟法分析了细砂型泥石流的诱发机制. 张奋翔等［18］采用FLO-2D软件模拟了不同降雨频率以及极端降雨条件下泥石流的暴发过程. 左正金等［19］采用PFC2D软件模拟了泥石流在降雨条件下的启动机理及运动过程.

本研究在典型坡面泥石流室内模型试验的基础上，利用三维颗粒流数值模拟软件PFC3D建立了坡面泥石流的数值模型，并提出等效渗透压力的简化方法，然后对不同坡角、不同坡表饱和区厚度及不同土体抗剪强度条件下坡面泥石流的启动过程和启动机理进行了分析，以期为泥石流地质灾害的预防提供理论参考.

1 室内模型试验

1.1 试验材料

构成泥石流物源体的颗粒级配会直接影响土体的物理特性和水理特性. 为保证试验结果具有真实性和代表性，从北京西山某坡面泥石流处采取原状土试样，采用振筛机对试样进行颗粒筛分与级配实验. 筛分后的颗粒级配曲线如图1 所示.通过计算得到试样的不均匀系数Cu=114.226，曲率系数Cc=0.541. 根据土的分类标准可知，该试样为级配不连续的土.

图1 颗粒级配曲线Fig.1 Particle grading curve

1.2 室内模型试验模拟坡面泥石流的启动过程

采用泥石流模拟槽体装置进行室内模型试验，槽体角度可调节，试验斜坡体长2 m，宽1.5 m，用人工降雨机模拟降雨，试验中降雨以高强度连续降雨为主. 进行室内模型试验时，分别设置三种坡角的斜坡体，如图2 所示，左侧斜坡体坡角为30°，中间斜坡体坡角为40°，右侧斜坡体坡角为50°.

由图2 可知，在室内模型试验模拟的整个过程中，右侧斜坡体先出现坡脚临空面土体的局部坍塌，而后在持续降雨作用下，右侧斜坡体的上部土体率先出现裂缝，裂缝的出现加速了土体与雨水的耦合作用，随着降雨入渗量的增加，裂缝逐渐侵蚀、扩展，土体的含水率逐渐加大、孔隙水压逐渐升高、抗剪强度逐渐降低，当土体的下滑力大于抗滑力时，右侧斜坡体发生较大规模的滑塌. 在室内模型试验模拟的整个过程中，中间斜坡体不仅出现了坡脚临空面土体的局部坍塌现象，还出现了细颗粒运移、坡体表面粗化以及坡体坡度变缓等现象. 左侧斜坡体的整体稳定性较好，在降雨雨力溅蚀下，左侧斜坡体出现了坡体表面粗化以及坡体坡度变缓等现象，但未出现坡脚临空面土体的局部坍塌现象.对比三种斜坡体的坡面泥石流启动过程发现，坡角与坡面泥石流形成的难易程度呈正相关，且坡角越大，坡面泥石流造成的危害也越大.

坡面泥石流启动过程的剖面堆积示意图如图3所示. 在降雨初期，因为土体初始含水率低，所以雨水全部下渗进入土体，从而导致土体中孔隙水压力升高，土体抗剪强度降低. 在降雨中期，坡体上缘出现微小的裂缝，裂缝的出现加速了土水的相互作用，使得表层土体逐渐饱和，产生坡面径流，坡体表面的细颗粒被表面径流冲刷至坡脚. 随着降雨的持续，即在降雨后期，坡体产生的裂缝逐渐增多、扩展，最终贯通形成滑动面，当土体破坏面上的剪应力高于抗剪强度时，坡面泥石流启动.

图3 坡面泥石流启动过程的剖面堆积示意图Fig.3 Schematic diagram of profile accumulation of slope debris flow startup process

2 数值模拟试验

2.1 数值模型的建立

基于北京西山某坡面泥石流的室内模型试验，利用三维颗粒流数值模拟软件PFC3D建立坡面泥石流试验槽数值模型. 先建立滑槽，滑槽坡角α为30°，滑槽水平长2 m，宽1 m，如图4（a）所示；然后在滑槽内生成0.4 m高的土体颗粒，如图4（b）所示，土体颗粒的最小粒径为5 mm，坡面最大粒径为15 mm.

图4 滑槽和坡面泥石流试验槽的数值模型Fig.4 Numerical model of chute and slope debris flow test chute

2.2 细观参数的选取

在颗粒流方法中，通过球（ball）来模拟土体颗粒，并通过接触（contact）来模拟土体颗粒之间的相互作用，通过赋予接触不同的细观参数来模拟土体颗粒之间的不同力学性质.为确定数值模型中采用的细观参数的合理性，需先建立标准试样的数值模型，并设定数值模型的初始细观参数，在此条件下进行标准试样的直剪试验，从而可获得初始的土体抗剪强度，然后通过不断调整细观参数，使数值模拟中标定出的土体抗剪强度与室内模型试验确定的土体抗剪强度相对应，最终即可获得数值模型所采用的细观参数，如表1所示.土体颗粒间的接触模型均选用接触黏结模型.

表1 坡面泥石流试验槽数值模型中采用的细观参数Tab.1 Mesoscopic parameters used in numerical model of slope debris flow test chute

3 降雨影响的简化分析

在降水作用下雨水入渗至坡体表面的堆积体内，降雨过程涉及饱和-非饱和渗流过程.在颗粒流数值模型中，因降雨引起的饱和-非饱和渗流涉及水-力相互耦合作用，通常需要进行复杂的模拟及计算，为了简化分析，本研究仅考虑因降雨引起的坡表饱和区的渗流作用.

在降雨过程中，雨水下渗进入土体，当坡面的雨水入渗速度大于土体的雨水入渗速度时，会率先在坡表形成饱和区.颗粒流数值模拟方法中，考虑降雨对坡体的影响时，会通过施加附加力来考虑土体颗粒与流体的作用.在饱和土体内，由于水头压力差的存在，流体对土体产生动水压力和拖曳力，这两部分力可简化为等效渗透压力作用在土体颗粒上［20］，每个土体颗粒受到流体施加的等效渗透压力F为：

式中：γw为水的重度，N·m-3；i为水力梯度，其值为土体中两点水头之差（H1-H2）与其渗透距离L的比值；V为单个土体颗粒的体积，m3.

在坡表一定厚度的饱和区内可能形成如图5（a）和5（b）所示的两种基本的渗透压力形式，即矩形渗透压力和三角形渗透压力.渗透的边界情况不同产生的渗透压力形式也会不同，同时还会存在两种基本渗透压力的组合形式，即梯形渗透压力.对于多数的边坡而言，水有向坡脚运动的趋势，均会产生如图5（c）所示的梯形渗透压力.

图5 坡表饱和区的渗透压力形式Fig.5 The forms of seepage pressure in the saturated area of the slope surface

土体颗粒在饱和状态下还会受到浮力作用，其大小通过公式（2）计算得出：

式中：ρw为流体密度，kg·m-3；g为重力加速度，m·s-2；V为土体颗粒的体积，m3.

由于降雨入渗会引起土体的重度增加，因此降雨影响下的土体重度按照公式（3）计算：

式中：ρd为土体骨架密度，kg·m-3；ρw为流体密度，kg·m-3；θ为体积含水率，%.

除此之外，土体在饱和状态下，水对土体具有软化作用，会使其物理力学参数存在一定的折减，根据工程经验近似折减15%.

4 坡面泥石流启动过程分析

根据以上理论，基于如图4所示的坡面泥石流试验槽数值模型，采用PFC3D软件对坡面泥石流的启动过程进行数值模拟.采用施加流体作用力方法考虑降雨对土体的作用，并采用公式（1）～（3）计算土体的等效渗透压力、浮力及土体重度，同时对各参数进行15%的折减来考虑降雨作用的影响，饱和区渗透压力采用如图5（c）所示的梯形渗透压力，坡表饱和区厚度为0.05～0.15 m，坡角为30°，摩擦系数为0.3，数值模型采用的细观参数如表1所示，模拟结果如图6所示.

图6 坡面泥石流的启动过程Fig.6 Startup process of slope debris flow

从图6（a）可以看出，受到降雨影响的土体在流体的渗透压力作用下首先会在坡体上方形成少量裂隙，此时坡面形态完整，如图6（b）所示.随着降雨的持续影响，坡体产生的裂隙逐渐增多，且裂隙主要集中在坡顶和坡中，如图6（c）所示，裂隙的存在加快了水土作用，使得坡顶的土体向坡脚流动，并在坡底堆积，坡顶变薄，如图6（d）所示.降雨后期，这些裂隙逐步贯通，如图6（e）所示，在坡体表层形成滑动面，处于滑动面以上的土体在流体作用下向下滑动，最终形成坡面泥石流，如图6（f）所示.

5 坡面泥石流启动过程的影响因素分析

5.1 坡角对坡面泥石流启动的影响

为了分析坡角对坡面泥石流启动的影响，利用PFC3D软件分别设置坡角为20°、30°、40°、50°的坡面泥石流试验槽数值模型进行模拟，数值模型采用的细观参数如表1所示，饱和区渗透压力采用梯形渗透压力，坡表饱和区厚度为0.05～0.15 m，摩擦系数为0.3，数值模拟结果如图7所示.

图7 不同坡角对坡面泥石流启动的影响Fig.7 Effect of different slope angles on the startup of slope debris flow

从图7可以看出：坡角为20°时，土体的滑动体积为0.245 m3，平均滑动距离为0.113 m，此时坡体表层土体有向下滑移的趋势，坡顶土体位移很小，坡脚有少量土体颗粒溢出，坡面形态完整；坡角为30°时，土体的滑动体积为0.257 m3，平均滑动距离为1.053 m，坡顶土体颗粒向前流动，坡脚出现较多的土体颗粒堆积，坡面变化明显；坡角为40°时，土体的滑动体积为0.267 m3，平均滑动距离为2.243 m，坡顶土体颗粒被全部带走，坡面形态几乎被破坏；坡角为50°时，土体的滑动体积为0.297 m3，平均滑动距离为3.084 m，坡顶土体颗粒被全部带走，坡顶土体颗粒滑移后呈圆弧形.以上结果表明，坡角越大，越容易发生坡面泥石流，且土体的滑动距离越远、滑动体积越大.

5.2 坡表饱和区厚度对坡面泥石流启动的影响

为了分析坡表饱和区厚度对坡面泥石流启动的影响，利用PFC3D软件分别设置坡表饱和区厚度为0.02～0.06 m、0.05～0.15 m、0.10～0.30 m的坡面泥石流试验槽数值模型进行模拟，数值模型采用的细观参数如表1所示，坡角为30°，摩擦系数为0.3，饱和区渗透压力采用梯形渗透压力，当设定的坡表饱和区厚度超出堆积体厚度时，将渗透压力施加到滑面上，数值模拟结果如图8所示.

从图8可以看出：坡表饱和区厚度为0.02～0.06 m时，土体的平均滑动距离为0.759 m，滑动体积为0.268 m3；坡表饱和区厚度为0.05～0.15 m时，土体的平均滑动距离为1.002 m，滑动体积为0.279 m3；坡表饱和区厚度为0.10～0.30 m时，土体的平均滑动距离为1.752 m，滑动体积为0.283 m3.以上结果表明，坡表饱和区越厚，土体的滑动体积越大、滑动距离越远，即坡面泥石流的规模随坡表饱和区厚度的增大而增大.

图8 不同坡表饱和区厚度对坡面泥石流启动的影响Fig.8 Effect of different thicknesses of slope surface saturation area on the startup of slope debris flow

5.3 土体的抗剪强度对坡面泥石流启动的影响

土体的抗剪切强度用摩擦系数来表征，为了分析土体的抗剪强度对坡面泥石流启动的影响，利用PFC3D软件分别设置土体摩擦系数为0.1、0.3、0.5的坡面泥石流试验槽数值模型进行模拟，数值模型采用的细观参数如表1所示，坡角为30°，饱和区渗透压力采用梯形渗透压力，坡表饱和区厚度为0.05～0.15 m，数值模拟结果如图9所示.

图9 不同土体抗剪强度对坡面泥石流启动的影响Fig.9 Effect of different shear strengths of soil mass on the startup of slope debris flow

从图9可以看出：土体摩擦系数为0.1时，土体的平均滑动距离为1.350 m，滑动体积为0.279 m3；土体摩擦系数为0.3时，土体的平均滑动距离为0.990 m，滑动体积为0.257 m3；土体摩擦系数为0.5时，土体的平均滑动距离为0.020 m，滑动体积为0.001 m3.以上结果表明，土体的抗剪强度越大（即土体的摩擦系数越大），发生坡面泥石流的可能性越小，这是因为在降雨作用下，摩擦系数较大的土体产生的位移较小，坡体更稳定；土体的抗剪强度越小（即土体的摩擦系数越小），土体的滑动距离越远、滑动体积越大，发生的坡面泥石流规模也越大.

6 结论

结合室内模型试验，采用三维颗粒流数值模拟软件PFC3D建立了坡面泥石流试验槽数值模型，并通过颗粒流数值模拟方法分析了不同条件下坡面泥石流的启动过程和启动机理，得出结论如下：

1）基于渗流的流体作用过程，在数值模拟中，对饱和区的土体采用等效渗透压力的简化方法可以近似地模拟降雨对坡面泥石流的作用，并可以提高计算效率.

2）坡体在渗透力作用下发生坡面泥石流时，坡体上首先会形成一系列微小的裂隙，这些裂隙逐步侵蚀、扩展、贯通，在坡表形成滑动面，处于滑动面以上的土体颗粒在流体作用下向下滑动，形成坡面泥石流.数值模拟结果与室内模型试验结果大致吻合.

3）坡角越大，则坡体的稳定性越差，在相同降雨条件下更容易发生坡面泥石流.坡面泥石流规模随坡表饱和区厚度的增大而增大，坡表饱和区越薄，则雨水越容易形成地表径流，造成坡表冲刷.土体的抗剪切强度用土体的摩擦系数来表征，摩擦系数越小的土体越容易滑动，则坡面泥石流越容易发生，且发生的坡面泥石流规模越大.