滑坡碎屑流颗粒分选效应的数值模拟

李天话 樊晓一 姜元俊

(1.西南科技大学土木工程与建筑学院 四川绵阳 621010； 2.中国科学院成都山地灾害与环境研究所 四川成都 610041)

由不同形状、大小、密度的颗粒构成的颗粒物质系统在运动过程中会出现颗粒分选效应，如受到制药、运输、冶金等工业广泛关注的“巴西果”分离、“反巴西果”分离和“三明治”分离现象[1]。已有研究表明[2-3]，颗粒分选效应同样出现在滑坡碎屑流的运动过程之中，具体表现为：粗大粒径块石易向滑体前缘及表层聚集，而细小碎屑主要分布于滑体中下部，该效应最终导致了滑坡碎屑流反粒序堆积结构[4]的形成。

滑坡碎屑流颗粒物质的分选效应作为滑坡碎屑流动力学的主要表现特征之一，对沿程沟道床面物质的“铲刮”效应、滑体运动距离、影响范围、动能及冲击作用都有极重要的影响[3,5-7]，研究滑坡碎屑流颗粒分选效应和机理，对了解滑坡碎屑流运移堆积特征及滑坡冲击灾害的规划防治工作具有重要的科学意义。

尽管早期学者利用模型试验手段，通过惯性颗粒流理论、随机振动筛机理及挤压驱除机制等理论分析方法[7-9]，阐明了颗粒分选的形成机理，但各理论对颗粒流的密实度、离散度、黏性、级配等颗粒参数依赖较大，难以形成普适性的理论分析结果。

近年来，关于颗粒分选的研究多以现场调查[4,10]和模型试验[2,6,11]为主。虽然通过现场调查对碎屑流流动堆积区的粒组组成进行定量分析能够揭示滑体反粒序堆积特征，通过滑槽模型试验手段能揭示大、小颗粒含量、粒径差等因素对分选程度的影响，然而两种方法均难以实现对颗粒分选形成过程的直观分析。本文以模型试验[12]的资料和数据为依据，运用三维颗粒离散元素法( Discrete Element Method, DEM) 分析软件EDEM建立滑坡碎屑流颗粒流数值模型，分析颗粒流不同粒径组分在运动过程中分布位置及动能的演化过程，初步揭示滑坡碎屑流颗粒分选效应的形成过程，进而深入分析颗粒分选对滑坡碎屑流动力学特征的影响。

1 研究方法

1.1 离散元素法

离散元素法是基于粒子间接触的位置关系以及相互作用的力学关系来求得粒子的运动状态的一种数值计算方法[13]，该方法由Peter Cundall于1971年首次提出并应用于岩土力学的研究，经过40多年的不断深入与发展，DEM已被广泛运用于滑坡碎屑流运移冲击作用机制及承灾体动力响应等问题的相关研究，对于模拟大变形或大位移条件下的滑坡碎屑流运动具有较好的适用性[14]。

离散元素法利用牛顿第二定律和力-位移法则确定颗粒的运动及受力情况，其核心是颗粒接触特性，即接触本构模型[15]。接触模型的分析计算直接决定了颗粒所受的力和力矩的大小，对不同的仿真对象，需建立不同的接触模型，以提高仿真结果的准确性[16]。由于室内模型试验材料为无粘连的类球状碎石，具有散粒物料特性，颗粒表面黏附力很小，因此本研究采用Hertz-Mindlin (no slip)接触模型。

1.2 模型尺寸

如图1所示，以滑槽模型试验的模型尺寸为依据，建立数值模型。滑槽底板长2.19 m，宽0.3 m，侧板高0.35 m，用于重现碎屑流的流动环境。滑槽的倾角为35°；模型料槽长0.44 m，高0.2 m，宽0.3 m，用于装填滑体物料，碎屑流的启动通过一个触发门来控制；斜滑槽底端为刚性挡墙，拦挡滑体颗粒，便于直观获取滑体在斜滑槽上的堆积体形态，分析不同粒组的颗粒分选及运移堆积过程。

图1 滑槽模型示意图Fig.1 Sketch map of the flume model

1.3 数值模型参数的确定

离散元数值模拟中，材料的泊松比、剪切模量和密度等本征参数表示各种材料的自身属性，与外界条件无关。基于模型试验实测及计算结果，确定颗粒与几何体的材料属性如表1所示。

表1 材料属性表Table 1 Material properties

在离散元素法中，材料宏观的运动行为取决于颗粒间接触的微观力学参数，然而，物理模型试验中滑体颗粒的运动是多因素、多参数共同作用的复杂力学过程，数值模拟不能完全复原实际模型试验中的所有影响因素，因此，将数值模拟软件中有限的微观力学参数进行实测标定会导致得到的模拟运动过程与试验原型产生较大差异。对数值模拟材料参数的选择尚未形成一个完善成熟的理论方法，较多的方式是通过物理模型和数值模型结果的对比验证来确定相关力学参数。

图 2给出了不同时刻模型试验与数值模拟所得滑体形态的对比结果。图2表明，离散元数值模拟结果能够较好地与物理模型试验相吻合，进而表明表2所示的材料间的接触属性参数可用于进行下一步数值模拟研究，以弥补物理模型试验在仪器设备以及技术手段等方面的局限性，实现对颗粒分选形成过程的直观监测及分析。

1.4 试验设计

混合粒径颗粒流在运移堆积过程中会产生掺混和分选效应，导致滑体流动性增强[3,6-7]，运动堆积体呈现“上粗下细”的反粒序分布特征[2,4,11]。为了分析滑坡碎屑流颗粒分选过程中不同粒径滑体颗粒的动能及沿垂直于滑动方向(垂向)的分布位置演化等运移堆积特性，建立了如图3所示的M1，M2，M33组初始滑体。

2.1.3 使用者评价。在交流过程中，李主任表示，截至目前，该产品每天都在运行，情况总体良好，煤种和水分适应性不错，缩分精度、样品重量、样品各级粒度等各方面指标均达标。他特别指出，三德科技的该自动制样系统有几项关键技术解决了制样过程中的难题。①全通破碎机，该设备在破碎过程中采用了伞旋自清扫技术和给料皮带联动技术，解决了6mm湿煤破碎沾堵的难题，达到湿煤全通的效果；②风透低温快速在线干燥技术，有效解决了水分对制样环节带来的影响，保证制样过程不粘不堵；③高性能制粉技术，该自动制样系统制粉环节采用德国技术，制粉收集率达到95%以上，过筛率达到100%，远高于同类产品。

图3 3组滑体模型初始堆积状态Fig.3 Initial deposition forms of 3 samples

如图3及表3所示，初始滑体M1，M2和M3均由KL1至KL44个粒组组成，其中，以红色标记的KL1表示粒径为2～5 mm的滑体颗粒，以绿色标记的KL2表示粒径为5～10 mm的颗粒，以黄色标记的KL3表示粒径为10～20 mm的颗粒，以蓝色标记的KL4表示粒径为20～40 mm的颗粒。3组初始滑体M1,M2和M3中，4个粒组的质量百分比均相同，而各粒组的初始堆积位置不同。M1初始堆积体呈现反粒序堆积特征，具体表现为：堆积体沿垂直于滑槽方向，自下向上粒径依次增大(图3(a))；M3初始堆积体呈正粒序分布[17]特征：堆积体颗粒从下向上粒径依次减小(图3(c))；M2初始滑体中，4个粒组呈混杂排列特征(图3(b)及图4)，然而粗细颗粒在料槽中沉淀静置的过程中，受重力作用影响，小粒径颗粒会透过大颗粒之间的空隙，渗透到滑体中下部，导致初始滑体难以保持均匀的随机混杂状态。尽管如此，作为M1和M3之间的过渡模型，M2能够体现各粒组随机掺杂分布的特性。

表3 初始滑体级配组成表Table 3 Particle size distributions of the samples

图4 滑体M2各粒组初始堆积状态Fig.4 Initial particle deposition of M2

2 滑坡碎屑流颗粒分选效应

2.1 滑体颗粒垂向分布位置演化特征

国内外高速滑坡堆积体中已发现了大量的反粒序地质证据[2,4,11]。为分析颗粒分选效应影响下滑体反粒序分布结构的形成过程，研究滑体颗粒沿垂直于滑槽方向的分布位置演化特征，引入特征指标

(1)

式中，垂向分布位置指标HKLi(i=1,2,3,4)表示KL1-KL44个粒组距滑槽底板的平均高度，hj表示第j个颗粒距离滑槽底板的垂向高度(如图5)，n表示该粒组滑体颗粒的总数量。

图5 颗粒垂向高度示意图Fig.5 Sketch map of the vertical distribution of particles

图6给出了M1,M2和M33组初始滑体各粒组垂向分布位置指标HKLi随滑体运动时间的演化曲线，KL1-KL44条曲线均呈现先下降后上升的变化趋势，表达了初始滑体从开始释放沿斜滑槽下滑，到滑体长度不断拉长而深度不断减小，最后受到坡脚挡墙阻挡而壅高停积的运动过程。图7给出了M1,M2和M33组初始滑体在不同时刻的运动快照。

图6(a)及图7(a)表明，初始滑体M1呈反粒序分布，随着运动触发，滑体沿斜滑槽向下运动，位于初始滑体最上部的KL4率先运动到滑体前缘部位，且粗大颗粒间碰撞作用显著，导致滑体前缘密实度较低。随着运动的持续发展，滑体长度沿运动方向延伸，滑体平均厚度逐渐减小，各粒组滑体垂向分布高度逐渐达到最小值，伴随着挡墙处颗粒体的持续堆积壅高，滑体厚度逐渐增加，各粒组垂向分布高度增大，最终趋于稳定，滑体运动停止。在整个运动过程中，图6(a)所示4条曲线未出现相交，表明滑体反粒序分布结构显著存在于整个运动过程之中。

图6(c)及图7(c)表明，初始滑体M3呈正粒序分布，随着运动触发，滑体沿斜滑槽向下运动，位于初始滑体最上部的KL1率先运动到滑体前缘部位，且细小颗粒比表面积大，颗粒间摩擦作用显著，受表层及前缘细小颗粒的阻挡覆盖，滑体整体密实度较高。随着运动的持续发展，滑体长度沿运动方向延伸，滑体平均厚度逐渐减小，粗颗粒穿越表层细颗粒滑体，开始向滑体中上部运移，此时各粒组间掺杂均匀，滑体正粒序分布特征被打破，并开始向反粒序分布过渡。位于滑体前缘的细粒组分受约束相对较小，率先在挡墙处堆积，惯性作用更为显著的大颗粒比表面积较小，摩擦耗能相对较少，在碰撞分离过程中保持了较大能量，翻越细粒堆积体后覆盖停积，最终堆积体呈显著的反序分布特征。图6(c)所示4条曲线在0.79～1.03 s时间段内相交，表达了在颗粒分选效应的作用下，滑体正粒序分布向反粒序分布的过渡。

图6(b)及图7(b)表明,作为M1与M3的过渡模型，M2滑体初始状态各粒组间掺杂相对较均匀，但粗颗粒之间孔隙大，接触相对较松散，堆积高度相对较高。此外，沉淀静置过程中，细颗粒难以避免地穿过粗颗粒间的孔隙向下运移，导致较粗粒组KL3和KL4的平均分布位置位于初始滑体中上部，较细粒组KL1和KL2平均分布位置位于滑体中下部。运动触发后，反粒序分布特征进一步发展，随着滑体在挡墙处的逐渐停积，堆积体高度持续增加，显著的反粒序堆积特征形成。图6(b)显示，初始滑体M2中KL1与KL2间的平均分布高度差为2.00 mm，KL3与KL4间的平均分布高度差为26.46 mm。最终堆积体中，KL1与KL2间的平均分布高度差为27.48 mm，KL3与KL4间的平均分布高度差为46.25 mm。表明对于粗细颗粒掺杂的初始滑体，随着滑坡碎屑流运动过程中的颗粒分选作用，将形成显著的反粒序分布特征。

图6 各粒组运动过程中的分布高度Fig.6 Location-time curves of vertical distribution of KL1-KL4

图7 不同时刻滑体运动状态快照Fig. 7 Motion state of different particle size components of a mobilized particle flow at different times

图8 最终堆积状态Fig.8 Forms of final depositions

2.2 滑体动能演化特征

滑坡发生时，滑体从高势能位置向低势能位置运动并伴随着巨大的能量释放。滑体颗粒间的碰撞及摩擦等动力特征可通过滑体动能的增减表达。受颗粒分选效应的影响，同样规模的滑体沿相同坡度的沟道运动，其动能差别很大，对受灾建筑结构的冲击破坏作用亦有很大的差异[7]。图9给出了M1-M33组初始滑体各粒组动能随运动时间的变化曲线。

图9 各粒径组分运动过程中的动能时程曲线Fig.9 Kinetic energy-time curve of KL1-KL4

如图9(a)所示，对于呈反粒序分布的初始滑体M1，各粒组最大动能与粒径及始滑位置距底板距离呈正相关，粒径越大，始滑位置距底板距离越大，动能峰值越大，且获得最大动能所需时间越短。

如图9(b)所示，对于粗细颗粒混杂的初始滑体M2，各粒组最大动能依然与粒径呈正相关，粒径越大，动能峰值越大。

如图9(c)所示，对于呈正粒序分布的初始滑体M3，各粒组最大动能不再与其始滑高度或粒径呈正相关。随着运动的触发，位于初始滑体最表层的KL1率先达到动能峰值，位于初始滑体最底层的KL4最后达到动能峰值，但底层初始滑体KL4获得的最大动能(19.44 J)略大于表层初始滑体KL1(19.15 J)。位于初始滑体中下部的KL3的动能峰值为18.16 J，略大于中上部初始滑体KL2(17.72 J)。

图9表明，滑体所得动能峰值受滑体垂向分布位置和粒径共同影响。位于滑体表层的颗粒所受约束相对较小，能量损失相对较少，便于动能发育；大粒径滑体颗粒惯性作用显著，在频繁的碰撞作用下，促进了颗粒间的能量传递[18]，且比表面积小，摩擦耗能小，动能得以充分发育。

2.3 滑坡碎屑流颗粒分选机制

(2)

式中，v表示颗粒的运动速度，h0表示该粒组初始堆积的质心距滑槽底板的距离，h表示颗粒沿垂直于底板方向的运动高度，g表示重力加速度，n表示颗粒的数量。

M1和M3两组初始滑体均由4个粒组组成，每个粒组沿垂直于滑槽底板方向的分布层次不同，对应着不同的初始堆积高度。图10给出了M1和M3两组初始滑体对应分布层次的粒组的时程曲线对比情况。

图10 各粒组时程曲线Fig.10 Comparisons of kinetic energy-time curves between different particle size components

自然条件下，受地质环境条件、坡体结构及滑体岩性等因素的影响[19]，不同滑坡碎屑流运动过程中，岩体破碎程度不同，粗大块石与细小碎屑混杂。颗粒流数值模拟结果表明，滑坡碎屑流运动过程中，大粒径滑体颗粒惯性作用显著，颗粒间接触以碰撞作用为主，碰撞过程中滑体密实度降低，小粒径颗粒易穿过滑体间的缝隙向下部运移，此外，小粒径颗粒比表面积大，摩擦耗能显著，在上部滑体的挤压作用下，动能发育受到进一步抑制，导致这些颗粒向滑体中后部聚集。而大颗粒比表面积小，碰撞作用显著，动能得以充分发育[18]，运动过程中越过下部细小颗粒，向滑体前缘及表层聚集，前缘及表层大粒径颗粒所受约束相对较小，在频繁的碰撞过程中进一步促进了大颗粒的飞跃分离(如图11)。因此，滑坡碎屑流颗粒分选效应具体表现为，粗大块石运动速度较快并向滑体表层及前缘运移，细小碎屑运动速度较慢，并向滑体底层及后部聚集。滑坡碎屑流颗粒分选效应是由振动筛分、小颗粒耗能较大率先停积、大颗粒碰撞分离等动力学机理共同作用的结果。

图11 滑坡碎屑流大、小粒径组分运动状态示意图Fig.11 Different trajectories of large and small particles

4 结论

本文以模型试验的资料和数据为依据，运用三维颗粒离散元素法分析软件EDEM建立滑坡碎屑流颗粒流数值模型，从滑体不同粒径组分垂向分布位置和动能两方面研究了滑坡碎屑流颗粒分选效应的形成过程，探讨了颗粒分选效应的形成机制，得出以下主要结论：(1) 滑坡碎屑流颗粒分选效应会导致滑体颗粒向反粒序分布结构发育。呈反序分布的初始滑体在整个运动过程中会保持反序分布结构，粗细混杂或者呈正序分布的初始滑体会在颗粒分选效应的作用下形成反粒序分布结构。(2)滑体所得动能峰值受始滑深度和粒径共同影响，始滑深度越浅，颗粒粒径越大，动能峰值越大。(3) 初始堆积高度一定的情况下，粗大块石在运动过程中易获得较大的运动速度并运移至滑体表层及前缘，而细小碎屑运动速度较慢，并向滑体底层及后部聚集。(4) 滑坡碎屑流颗粒分选效应是由振动筛分、小颗粒耗能较大率先停积、大颗粒碰撞分离等动力学机理共同作用的结果。