基于EDEM的碎屑流运动规律及冲击性能研究*

潘 青 张清照 李艺灵

(①同济大学土木工程学院地下建筑与工程系，上海 200092，中国)(②同济大学岩土及地下工程教育部重点实验室，上海 200092，中国)

0 引 言

随着国家大力推进公路建设发展进程，我国公路建设的重点也由东南沿海转向了西南欠发展地区，该地区内公路往往需要布设于沿河傍山地段，由于复杂的地质环境条件及人类工程活动，泥石流灾害十分发育，据不完全统计，我国西部几大交通干线公路受泥石流灾害影响多达3600余处(马东涛等，2003)。特别是每年的雨季，由于泥石流造成的公路阻断、堵车等问题层出不穷，以2019年8月20日四川阿坝州特大泥石流为例，泥石流的发生既对公路工程本身造成了损害，同时造成的公路中断也对受灾地区人民生命财产的及时抢救带来了影响，一定程度上加重了灾区灾情。目前对于泥石流的工程防治措施可分为跨越工程、穿过工程、排导工程、拦挡工程等，为了合理进行此类工程的设计与施工，需要对泥石流的运动规律及其冲击性能开展系统研究。

近年来，国内外专家学者分别通过试验研究、数值模拟以及理论分析等手段对泥石流的动力性能展开了大量研究，成果颇丰。何晓英等(2016)分别进行了清水流、水石流与泥石流工况下的泥石流冲击特性模型试验，模拟了沟谷的真实形态，考虑了浆体黏度与级配颗粒组合，采用小波理论对冲击荷载时程曲线及冲击荷载均值进行了分析；Scheidl et al.(2013)进行了小规模泥石流冲击特性模型试验，基于16组试验数据得出了泥石流最大冲击力的计算方法；刘兴荣等(2018)通过开展室内水槽试验，研究了不同粒径级配条件下工程弃渣泥石流启动机理。相较于试验研究方法，数值模拟由于其分析过程稳定可靠、成本低廉，已成为泥石流灾害研究的一种重要方式。李兆华等(2018)基于FEMLIP对泥石流固体、流体及固流转化3个阶段进行了模拟，并开展了泥石流对挡土墙的冲击效应的启发式研究；韩征等(2019)建立了整合HBP本构关系的稀性泥石流动力过程三维数值模型；马秋娟等(2015)，周健等(2015)均利用PFC软件对泥石流启动过程进行数值模拟研究；柳春等(2019)基于SPH-FEM耦合数值方法，考虑了大块石-泥石流浆体-拦挡坝的相互作用，分析了泥石流冲击被拦挡的过程、冲击力时程等力学行为；汤碧辉等(2018)提出了黏性泥石流沟内大颗粒的滚动启动模型，并用FLUENT对多种工况下的颗粒启动进行了数值模拟，验证了模型的适用性；王秀丽等(2015)采用流固耦合分析方法，将CFX与ANSYS相结合，将常规重力坝与新型带支撑拦挡坝进行了对比分析，得出了结构极限抗冲击承载特性；胡卸文等(2019)采用有限元与无限元耦合方法，系统分析了落石冲击荷载下桩板拦石墙的抗冲击特性。

从上述研究成果可以发现，泥石流的运动机理十分复杂，而其中固体颗粒间表现出的摩擦特性是泥石流-防治结构相互作用过程中能量耗散的关键，因此，从固体颗粒出发开展泥石流基本运动机理及其冲击性能的研究可为后续研究多相介质运动提供坚实的基础。本文基于离散单元法，通过EDEM 2018软件进行碎屑流冲击流槽试验模拟，综合考虑流槽坡度、防治结构角度、颗粒级配三方面因素，探明各因素对碎屑流运动过程及冲击性能的影响规律，为相关防治结构设计与施工提供理论依据。

1 数值模拟方法

1.1 EDEM软件介绍

EDEM作为一种DEM仿真商业软件，最初被用于工业生产中的颗粒处理及其制造设备的生产过程的仿真和分析。用户可利用EDEM轻松快速地创建颗粒实体的参数化模型，同时可以将力、材料特性和其他物理特性添加到EDEM中，形成完整的颗粒模型，为用户建立个性化的模型处理环境。

因此，借助EDEM的颗粒工厂技术(Particle Factory)，可很好地模拟碎屑流的完整运动过程，并且该软件具备高效的并行计算能力，可进行百万级颗粒运算，能够基本满足碎屑流冲击拦挡结构的计算要求。

1.2 模型设计与参数选取

本文旨在通过数值模拟方法进行碎屑流冲击流槽试验模拟分析，进而开展碎屑流运动规律及其冲击性能研究，因此，对比分析国内外学者已经开展的流槽试验模型，综合考虑其尺寸及参数完整性，最终选择以Jiang et al.(2013)所完成的流槽试验为参照进行数值模拟模型建立。其流槽试验装置如图1所示，其中流槽整体长2.93im，宽0.3im，高0.35im，流槽坡度角α为40°，底部拦挡结构角θ为50°，初始颗粒沉积区域长0.44im，宽0.3im，高0.150im。

图1 流槽试验装置示意图(Jiang et al.，2013)

图2即本文参照上述参数建立的流槽试验模型，主要包括槽底面、侧墙、拦挡结构、颗粒储存盒、触发门以及垂直与水平两个方向的加压面。设置加压面是为了保证各组试验的初始条件保持一致，即每组试验开始前，均需通过垂直与水平加压面对颗粒进行压密。上述流槽试验模型部件可由CAD导入，导入后需要对部件特性及部件运动进行指定，流槽试验模型材料参数如表1所示。需要注意的是，该流槽试验模型各部件除触发门外均为静态部件，即在模拟过程中不发生移动，因此只需要对触发门进行运动指定，本文通过触发门的线性旋转模拟触发门打开过程。

图2 数值模拟模型组成

表1 数值模拟中材料参数取值

在颗粒级配选取方面，Jiang et al.(2013)所选取的物质组成中粒径为10～20imm的颗粒占比在95%左右(图3)，已有研究表明对于粒径10～20imm的颗粒占比在95%以上的情况，由于粒径小于10imm的颗粒在初始颗粒样品中的质量百分比相对较低，忽略尺寸小于10imm的颗粒对计算结果的影响(Shen et al.，2018)，故本文直接选取10～20imm间的线性颗粒级配作为数值模拟中的物质组成如图3所示。

图3 颗粒级配曲线

此外，固体颗粒参数及接触参数如表1和表2所示，其中颗粒杨氏模量、泊松比、密度等参数主要根据颗粒介质数值模拟中的常用值进行设置(Salciarini et al.，2010)。对于接触参数，由于本文模拟的固体颗粒均为球形，而实际多为带棱角且表面不光滑的粗糙颗粒，因此在Jiang et al.(2013)的试验基础上，本文通过增大静摩擦系数，减小滚动摩擦系数来尽可能还原真实颗粒情况，具体取值如表2所示。

表2 数值模拟中有关接触的系数选取

1.3 工况设计

本文运用控制变量法，选取对碎屑流运动过程及冲击性能影响较大的3个因素作为自变量，即流槽坡度角α、底部拦挡结构角θ以及颗粒级配，具体工况设计如表3所示。

表3 工况设计

1.4 模型验证

为了验证本文所建立模型的可靠度，本节将数值模拟结果与Jiang et al.(2013)的试验数据进行了比较，图4、图5分别为Jiang et al.(2013)的试验与数值模拟得到的泥石流运动过程图像，可以看出两组颗粒运动形态基本吻合，表明该模型可以很好地模拟颗粒运动过程，并且通过总结发现颗粒运动过程主要可分为如下4个阶段：

图4 流槽试验中碎屑流运动过程照片(Jiang et al.，2013)

图5 数值模拟中碎屑流运动过程图像

(1)第1阶段，颗粒储存箱的触发门被移除，颗粒开始崩塌(图4a2和图5b2)；

(2)第2阶段，随着颗粒的持续滑动，滑体后缘逐渐离开其原始位置，前部滑体转变为薄层，而中段滑体扰动较小(图4a3和图5b3)；

(3)第3阶段，滑体前端撞击刚性拦挡结构并发生反弹，同时，部分颗粒开始沉积在流槽末端的拦挡结构处(图4a4和图5b4)；

(4)第4阶段，越来越多的颗粒到达已沉积颗粒上，直到所有的颗粒运动停止(图4a5和图5b5)，最终颗粒沉积形态如图4a6和图5b6所示。

此外，图6对比了试验与数值模拟中的冲击力时程曲线，基于离散单元法进行数值模拟时，当与拦挡结构接触的颗粒被困在相对长的力链中，该颗粒所产生的能量会显著增加，进而产生力的峰值，同时导致冲击力时程曲线发生明显波动。因此，本节采用中值滤波方法对原始时程曲线进行处理，经处理后的时程曲线与试验结果吻合，即该模型可以很好地模拟颗粒流的冲击性能。

图6 流槽试验与数值模拟中的冲击力时程曲线

2 碎屑流运动规律分析

2.1 碎屑流基本运动规律

各工况下的数值模拟结果所表现出的碎屑流运动过程基本相似，以S1-2(α=40°，θ=50°)为例，运动过程可分为如下3个阶段：

(1)第1阶段(冲击)：在颗粒撞击拦挡结构后，前端滑体形成薄层沉积体(图7a～图7b)，此时颗粒运动速度达到最大值；

图7 碎屑流固体颗粒运动过程(S1-2)

(2)第2阶段(爬升)：后端颗粒逐步向下运动，在前端既有沉积物上发生碰撞并逐步爬升，到达拦挡结构的更高处(图7c～图7e)。颗粒的不断爬升堆积使前端沉积高度逐步上升，而此时颗粒运动速度逐渐减小，靠近拦挡结构的颗粒运动速度均趋于0im·s-1。此外，该阶段的沉积形态具有一定的变化规律，即由起初的凹型(图7c)变为直线型(图7d)，再到最后的凸型(图7e)；

(3)第3阶段(堆积)：待所有颗粒完全静止且稳定后，其最终沉积形态如图7f所示，由于上层流动颗粒的加入，以及颗粒向上的爬升运动，沉积形状最终呈现为凸型。

图8则进一步展示了碎屑流固体颗粒速度场的变化过程，其中红色代表速度为3.08～3.96im·s-1，黄色代表速度为2.22～3.08im·s-1，绿色代表速度为1.35～2.22im·s-1，蓝色代表速度为0.49～1.35im·s-1，灰色代表速度为0～0.49im·s-1。从图9中可以发现静态沉积区在不断阻止颗粒的运动，导致颗粒减速向上爬升，同时，随着沉积体积的增大，静态沉积区逐渐向上游延伸。当既有沉积物坡度到达一个阈值后，随后的颗粒无法越过沉积物(图8a5)，只能冲击并堆积在现有的沉积物上。

图8 碎屑流固体颗粒速度场变化过程(S1-2)

图9 不同流槽坡度角下碎屑流运动速度时程曲线(S1-1，S1-2，S1-3)

2.2 流槽坡度角α对碎屑流运动规律的影响

如2.1节中所述，各工况下运动规律大致相似，但从图9中可以发现，随着流槽坡度角的增大，固体颗粒越快达到速度峰值，且峰值越大，当流槽坡度角α由35°分别增加至40°和45°时，对应的速度峰值涨幅分别为8.37%和22.62%。此外，颗粒运动速度从峰值到停止运动所需时间随着坡度的增大而减短，35°、40°和45°所对应的时间分别为1.50is、1.25is和1.10is。需要说明的是，为了保证各工况下初始颗粒运动状态一致，t=0is至t=0.5is期间在料箱内对颗粒进行预先压密处理，因此该区间内存在明显的速度变化。

有效残余力由最终沉积形状的几何形状决定(Jiang et al.，2013)，在流槽下端形成的速度为0的三角形区域被广大学者称为dead zone，该区域的长度在很大程度上取决于坡度，同时也影响着下端结构的平均受力。因此进一步观察各工况下的沉积形态(图10)可以发现，随着流槽坡度的增大，沉积体高度越低而长度越宽。

图10 碎屑流固体颗粒速度场变化过程(S1-1，S1-2，S1-3)

2.3 底部拦挡结构角θ对泥石流运动规律的影响

本节对比分析了底部拦挡结构角θ分别为30°、50°、70°和90°的工况，即S2-1，S3-1，S2-2和S2-3。数值模拟结果表明，在撞击拦挡结构前，各工况下固体颗粒的运动速度及运动规律大致相似，而撞击后的颗粒堆积情况如图11所示，其中a为侧视图，b为颗粒在拦挡结构面上的分布。

图11 碎屑流固体颗粒堆积形态(S2-1，S1-2，S2-2，S2-3)

图11利用不同颜色表示不同粒径的颗粒，其中红色粒径最大(18～20imm)，其次是黄色(16～18imm)、绿色(14～16imm)、蓝色(12～14imm)、灰色(10～12imm)。可以发现随着底部拦挡结构角θ的增大，颗粒在拦挡结构上的爬升距离逐渐降低，并且，当θ=90°时的固体颗粒的反序现象并不十分显著，而随着角度变缓反序现象愈发明显，即细颗粒主要分布于沉积物的底部，粗颗粒分布于顶部。这是由于随着固体颗粒爬升距离的增加，颗粒继续前进，此时，细颗粒仍可沿着粗颗粒向下移动，而粗颗粒则只能继续向前滑动，最终呈现为反序现象。

2.4 颗粒级配对碎屑流运动规律的影响

图12为S3-2下不同粒径颗粒的运动过程，其中2～5imm显示为黄色，5～10imm显示为蓝色，10～15imm显示为绿色，15～20imm显示为红色。当固体颗粒释放后直至沉积停止运动后，可观察到明显的反序分离现象，即粗颗粒主要分布于顶部与前部，而细颗粒主要分布于底部与后部。这是由于大多数粗颗粒在运动过程相互碰撞，飞跃后分离至顶部与前部，而细颗粒由于比表面积较大，相应的能量损失较大，运动相对较慢，因此在粗颗粒缝隙间穿过滑落至底部与后部。

图12 不同粒径颗粒的运动过程(S3-2)

图13进一步对比了S3-1，S3-2，S3-3工况下固体颗粒运动速度变化曲线，分别对应了粗颗粒含量由低变高的情况，发现随着粗颗粒含量的提升，固体颗粒峰值运动速度逐渐减小，表明细颗粒可一定程度提高碎屑流固体颗粒的运动速度。

图13 不同颗粒级配下碎屑流运动速度时程曲线(S3-1，S3-2，S3-3)

3 碎屑流冲击性能分析

3.1 流槽坡度角α对碎屑流冲击性能的影响

图14为不同流槽坡度角α下碎屑流固体颗粒对拦挡结构的冲击力时程曲线，观察到流槽坡度角α越大，泥石流冲击力则越快达到峰值，且峰值越大，同时后续也越快减小至残余冲击力，即沉积颗粒的重力分力稳定值，这是由于重力加速度在沿流槽方向上的分量发生了改变，进而影响了沿流槽方向的冲击力。

图14 不同流槽坡度下碎屑流冲击力时程曲线(S1-1，S1-2，S1-3)

3.2 底部拦挡结构角θ对碎屑流冲击性能的影响

图15为不同底部拦挡结构角θ下碎屑流固体颗粒对拦挡结构的冲击力时程曲线，发现各工况下，当泥石流到达拦挡结构(t=1.4is)后，冲击力迅速到达峰值，之后逐渐降低到残余冲击力，并且最大冲击力随拦挡结构角θ的增大而增大，这是由于当拦挡结构角θ增大后颗粒间的相互作用随之降低，最终导致了颗粒对拦挡结构冲击力的提升。

图15 不同拦挡结构角下碎屑流冲击力时程曲线(S2-1，S1-2，S2-2，S2-3)

3.3 颗粒级配对碎屑冲击性能的影响

图16为不同颗粒级配下碎屑流固体颗粒对拦挡结构的冲击力时程曲线，发现碎屑流最大冲击力随着粗颗粒占比的增大而显著增大。其中：S3-3的粗颗粒占比最大，其冲击力时程曲线的波动也更为剧烈，主要源于前端粗颗粒不断撞击拦挡结构以及后续颗粒碰撞前缘堆积体而间接冲击拦挡结构的过程，相应地，该组泥石流冲击力峰值最大；而S3-1中粗颗粒占比最小，主要以细颗粒为主，因此，颗粒间的摩擦耗能较为严重，此时泥石流对拦挡结构的冲击作用主要源于前缘堆积体的重力作用，同时细颗粒间的碰撞作用相对较小，导致相较于S3-3，其冲击力时程曲线波动小且冲击力峰值小。

图16 不同颗粒级配下碎屑流冲击力时程曲线(S3-1，S3-2，S3-3)

4 结 论

本文基于离散单元法，利用EDEM 2018软件，通过数值模拟手段分析了碎屑流固体颗粒的运动过程及其与拦挡结构间的相互作用，特别是碎屑流的冲击性能。主要考虑了流槽坡度角、底部拦挡结构角以及颗粒级配的影响，对碎屑流的运动规律及其冲击性能进行了系统研究，得到以下主要结论：

(1)流槽坡度越大，碎屑流固体颗粒越快达到其速度与冲击力峰值，且对应的速度与冲击力的峰值越大。

(2)拦挡结构角越大，与其相互作用的固体颗粒数量越多，因此碎屑流越快达到速度和冲击力的峰值，同时也越快停止滑动，且对应的速度与冲击力的峰值也越大。

(3)不同颗粒级配的碎屑流均呈现反序特性，即细颗粒主要分布于下部及后部，并且细颗粒的存在可显著提高泥石流的运动速度，而粗颗粒对运动速度影响不大。

(4)粗颗粒占比越大，冲击力主要源于前端粗颗粒对拦挡结构的撞击作用以及后续颗粒碰撞前缘堆积体而产生的能量，冲击力时程曲线波动大且冲击力峰值大，相反地，细颗粒占比越大，颗粒间的摩擦耗能越严重，冲击力主要源于前缘堆积体的重力作用，其冲击力时程曲线波动小且冲击力峰值小。