基于离散元法的雪崩效应的仿真研究

冯靖禹 韩韧 张宇峰 辛昱鋆

摘  要： 伴隨着全球气候变暖，雪崩灾害的风险逐渐增加。因此，调查雪崩的发展规律和内在机理十分迫切和必要。而使用实验的方法来研究往往效率低下，并且只能获得宏观数据。为了更深入地了解雪崩效应，提出结合计算机仿真技术和转鼓雪崩实验，通过建立基于离散元法（DEM）的仿真模型来模拟发生雪崩时的颗粒流态。通过运用Hertz-Mindlin接触理论和牛顿第二定律塑造的仿真模型以及在相同配置下的物理实验的对比，证明了离散元法研究雪崩效应的有效性。基于该模型，我们跟踪每个颗粒的运动状态，并研究了自由表面上参与雪崩的颗粒数量的变化。参与雪崩的颗粒数量的增长速度经历了先增加，然后恒定，最后减少的过程。我们还发现不同区域内的颗粒平均速度随时间都呈现先增大后减小的规律。离散元仿真方法为雪崩颗粒流的研究提供了参考。

关键词： 水平转鼓;雪崩;仿真;离散元法

中图分类号： TP391. 9    文献标识码： A    DOI：10.3969/j.issn.1003-6970.2020.07.055

本文著录格式：冯靖禹，韩韧，张宇峰，等. 基于离散元法的雪崩效应的仿真研究[J]. 软件，2020，41（07）：264-268

Simulation of Avalanche Effect Based on Discrete Element Method

FENG Jing-yu， HAN Ren^*， ZHANG Yu-feng， XIN Yu-yun

（School of Optical-Electrical and Computer Engineering， University of Shanghai for Science and Technology， Shanghai 200093， China）

【Abstract】： With global warming， the risk of avalanche disasters is increasing. Therefore， it is urgent and necessary  to investigate the development rule and internal mechanism of avalanches. Studing by experimental methods is often inefficient and can only obtain macro-data. In order to learn more about the avalanche effect， a combination of computer simulation technology and drum avalanche experiments is proposed. The particle flow pattern during an avalanche is simulated by establishing a simulation model based on the discrete element method （DEM）. By comparing the simulation model using Hertz-Mindlin contact theory and Newton's second law with the physical experiments under the same configuration， the validity of studing avalanche effect by discrete element method is proved. Based on this model， we track the motion of each particle and investigate the quantitative change of particles involved in avalanches on the free surface. The growth rate of avalanching particles number experiences a process of increasing， then stabilizing， and finally reducing. We also find that the average velocity of particles in different regions increase at first and then decrease. The discrete element method provides a reference for the research of avalanche particle flow.

【Key words】： Horizontal rotating drum; Avalanche; Simulation; Discrete element method

0  引言

颗粒物质广泛存在于自然界，并与人类的日常生活以及生产活动密切相关。一般认为，单个颗粒的典型尺度在10^–6～10 m的范围内，其运动规律服从牛顿运动定律^[1]。颗粒系统在受到外力或内部应力状况变化时发生运动，表现出流体的属性，从而形成颗粒流。了解和研究颗粒流不仅对预防自然灾害，例如雪崩、山体滑坡和泥石流^[2-3]，而且对提高工业生产效率，包括制粒、研磨、干燥和混合过程^[4]都具有重要意义。

水平转鼓由于其简单性和可控性，是研究颗粒流动特性的理想工具^[5-6]。转鼓在低速旋转的情况下，转鼓中的颗粒流会表现出一系列离散的周期性雪崩^[7]。不少学者针对转鼓内雪崩问题做了大量实验研究，甚至借助正电子发射跟踪（PEPT）和磁共振成像（MRI）。虽然这些研究大大提高了我们对颗粒流运动行为的理解，但由于条件有限，此类实验往往是基于宏观观察，研究内容不够全面和深入。因此，在微观层面上缺乏实验数据不但意味着对实验结果的分析往往要凭经验确定，例如每个颗粒的位置和速度等，而且实验效率低，耗费时间长。

近年来，由于计算机的高速发展，把计算机仿真技术^[8-12]运用在分析颗粒流动行为成为一种可行且高效的研究方法。由于颗粒流的行为取决于粒子间相互作用的结果，单个颗粒的信息将有助于揭示整个颗粒系统的基本原理。因此，英国皇家工程院院士Peter Cundall博士提出的直接对固体颗粒进行跟踪的离散元法（DEM）得到了关注。这种方法特别适用于模拟离散颗粒组合体在准静态或动态条件下的运动及变形过程。运用上述这种建模思想，借助計算机强大的计算能力^[13，14]，能极大地提高研究速度。把离散元法仿真和转鼓雪崩相结合来研究具体有如下几点优势：

（1）使我们能够准确而快速地模拟出真实世界的颗粒系统。类似于雪崩、泥石流等大型自然灾害如果通过实验的方法复现，费用极其高昂，且存在诸多危险;对于生产环节，制备不同的颗粒材料和生产机器需要耗费大量的人力物力，且受到各种复杂环境因素的影响，实验难度很大。通过仿真，可以高效地搭建实验环境，迅速开展研究工作。

（2）离散元仿真能直接获得离散物质大量复杂行为信息以及不易测得的颗粒尺度行为的信息，为分析颗粒流的运动、受力、热量和能量传递等信息提供高级的解决途径。甚至对于连续介质理论无法解释的物质力学行为，也可以进行准确的预测和分析。

（3）对于实际的工作生产，通过计算机仿真的直观感受，能拓展工程师设计思路，加速在设计阶段找出设计的潜在问题，使用数据精确突出设计方案的优势，加速产品创新。

综上所示，运用基于离散元法的仿真模型来研究雪崩相较于单纯地实验观测有着巨大的优势，为研究转鼓内的雪崩现象开辟了一条全新的研究方向。

因此，本文基于离散元法进行雪崩效应的仿真研究，以进一步剖析雪崩颗粒流的运动特性。我们基于离散元法基本原理搭建迭代式计算流程框架。考虑到转鼓以及内部颗粒的物理特性，运用Hertz-Mindlin接触模型和软球离散元模型进行仿真建模。并根据理论基础对模型的材料属性进行参数优化。通过仿真结果和实验数据在表面角和颗粒速度方面的对比，验证仿真模型的可靠性和有效性。我们还对参与雪崩的颗粒数量以及不同区域的颗粒速度进行统计。为高效预测雪崩发展过程以及科学设计转鼓内的颗粒系统提供了一定的参考价值。

1  离散元法基本原理

离散元法自提出以来，在岩土工程和计算机仿真应用领域发挥了其他数值算法不可替代的作用。其基本原理是把不连续体分离为刚性元素的集合，使得各个刚性元素满足运动方程，用时步迭代的方法求解各刚性元素的运动方程，继而求得不连续体的整体运动形态。该方法允许单元间的相对运动，不一定要满足位移连续和变形协调条件，计算速度块，所需存储空间小，尤其适合求解大位移和非线性的问题。

在建模过程中，把物料中的每个颗粒单独作为一个粒子单元建立数学模型，并给定粒子单元的尺寸和物理性质，如质量、刚度和阻尼等，各个粒子之间存在接触和分离这两种关系。

在进行数值计算时，把总的仿真计算时间划分为一个个时间步长。选取的时间步长足够小，使得在一个单独的时间步长内，除了与选定单元直接接触的单元外，来自其他任何单元的扰动都不能传播开来，即在任意时刻单元所受到的作用力只取决于该单元本身及与之直接接触的其他单元。

当接触发生时，接触点处就会产生接触力和力矩，其大小可以根据接触力学模型求出。利用牛顿第二定律建立每个粒子的运动方程，并用中心差分法求解，整个介质的变形和演化由各个单元的运动和相互位置来描述。通过相邻颗粒之间的碰撞产生的接触力、力矩和不平衡力，计算每个颗粒在特定时刻的运动特征（速度、加速度等）。每一个时间步长进行一次迭代，并根据前一次计算所得到的单元的位置，作为下一次迭代的出发点，用以求出单元的新位置。通过遍历跟踪计算每个单元的微观运动，即能获得整个颗粒系统的宏观运动规律。其计算流程图如图1所示^[15]。

2  三维模型搭建

2.1  接触模型

接触模型是离散单元法模拟散体物料运动的重要基础，其本质是在准静态下分析散体颗粒接触时产生的接触力和弹塑性变形。接触模型中的分析计算对颗粒所受的力和力矩的大小起决定性的作用。在本研究中，为了最准确地仿真出现实世界的转鼓内雪崩现象，使用基于Hertz-Mindlin接触模型的软球离散元模型来进行研究。

2.2  模型搭建

仿真过程分为三个阶段：转鼓内填充颗粒、颗粒自然沉降和转鼓运动。在第一阶段，生成转鼓和颗粒。设置颗粒和几何模型的物理参数是利用离散单元法模拟颗粒运动的重要环节之一，参数设置的正确与否对仿真结果产生直接影响。本研究假定转鼓和颗粒为同种玻璃材料。鼓沿x轴方向施加周期性边界条件。如图2所示。这种处理可以在生成足够的数据进行分析的同时，使用较少的颗粒以缩短仿真时间。之后，在生成的鼓内填充52000个半径为1 mm的颗粒。填充方法为在鼓中随机位置判断有无填充空间，如果有则生成颗粒。在第二阶段，悬浮在转鼓中的颗粒在重力g= 9.8 m/s²的作用下完成自由落体，填充转鼓下半部分50%的空间。在第三阶段，在颗粒沉降稳定后，转鼓开始以0.19 rpm的速度旋转，带动颗粒床运动。模拟的时间步长设置为3×10^–7秒。总运行时间为300秒，模拟计算时间约为100小时。

3  仿真结果

3.1  仿真准确性验证

为了验证离散元法仿真的准确性，采用空间滤波测速法（SFV）获得实验数据^[16]。实验数据和仿真结果在相同配置的装置下进行。对于此转速，实验和仿真中都显示了离散的雪崩。

首先，我们在仿真和实验中研究表面角θ的变化。图3显示了θ随时间t的变化（取50个连续雪崩的平均值），其中红线表示仿真结果，蓝点表示实验数据。两条曲线都显示出相同的变化趋势，两者之间的最大差异也在1°以内。因此，仿真结果与实验数据吻合良好。

颗粒速度是雪崩的另一个重要特征。图4显示了在雪崩期间颗粒床表面的颗粒平均速度v沿y轴的分布（取50個连续雪崩的平均值）。仿真曲线与实验曲线一致，它们之间的平均差为0.78 mm/s，并且都在表面径向中心出现最大颗粒速度。此外，高速区的标准差值大于低速区的标准差。Lin等人也在在较高转速的滚动状态下观察到了相似的结果^[17]。

基于上述两种不同的全局属性的对比，我们验证了所提出的仿真系统可以准确地模拟旋转鼓中的雪崩现象。因此，离散元法仿真模型可以提供一种有效的方法来研究雪崩过程中的颗粒流行为。

3.2  参与雪崩的颗粒数量研究

参与雪崩的颗粒数量是表明雪崩规模的关键指标。我们观察到部分颗粒在滚动一定距离后会停止运动。因此，我们使用通过跟踪每个颗粒的运动轨迹来判断颗粒是否被卷入雪崩。雪崩开始时颗粒所处位置被标记为起点。在雪崩期间，当其位置距离起点大于一个粒径时被认为参与了雪崩。图5显示了一个粒径厚度的表面层在雪崩过程中参与雪崩的颗粒数量N随时间t的变化（取50个连续雪崩前5秒的平均值）。颗粒的数量在持续增长，这意味着雪崩的规模在持续扩大。基于不同时段的不同增长速率，我们将它分为3个时段，分别用T1-T3标记。在T1（0-0.2 s）期间，颗粒数呈现出加速增长的趋势。在T2（0.2-0.3 s）期间，颗粒数以均速增加（在图5中用实心黑线标记）。在T3（0.3-0.5 s）期间，增长速度减慢并最终变成0，颗粒数量增加到最大。

3.3  参与雪崩的颗粒速度研究

床表面不同位置的参与雪崩的颗粒速度是另一项重要参数。我们把颗粒床一个粒径厚度的表面层按y轴方向等分成5个区域，按所处y轴由小到大的位置标号为1-5号区域。图5显示了这5个区域内的颗粒平均速度v随时间t的变化规律。5个区域的颗粒平均速度都呈现先增大后减小的规律。把同一时刻不同区域的速度进行对比可以发现，相对处于中心位置的2-4号区域的颗粒平均速度更大，这意味着此区域内的颗粒运动最剧烈。而把同一区域在不同时刻的速度进行对比可以发现，相对处于转鼓更高位置的4-5号区域的最大速度出现的时刻更晚。

4  结语

雪崩严重影响着人们的生命安危，且传统的实验的宏观分析十分低效，因此通过计算机仿真的方法来研究雪崩势在必行。本文提出使用基于离散元法的仿真来研究转鼓中雪崩状态下的颗粒流动行为。详细介绍了所使用的离散元计算流程，颗粒接触时采用的Hertz-Mindlin模型以及模型的生成过程。通过将模拟结果与实验数据（如表面角和表面颗粒速度）进行比较，证明了离散元法仿真能实现转鼓内雪崩的有效仿真。此外，运用离散元法能跟踪每一个颗粒，对每个颗粒运动状态进行分析。结果表明雪崩规模的增长速度在每个时间段中是不同的：在T1期间加速，在T2期间以均速增长，最后在T3期间减速直至变成0。最后，我们还发现不同区域内的颗粒平均速度随时间都呈现先增大后减小的规律。这充分展示了离散元仿真技术在分析雪崩效应的可靠性和优越性，为雪崩的预测喝分析提供了有效的途径。

参考文献

1. 胡国明. 颗粒系统的离散元素法分析仿真： 离散元素法的工业应用与EDEM软件简介[M]. 武汉： 武汉理工大学出版社， 2010： 3-7.

2. H. Yang， B. F. Zhang， R. Li， G. Zheng， V. Zivkovic. Particle dynamics in avalanche flow of irregular sand particles in the

slumping regime of a rotating drum[J]， Powder Technology， 2017， 311（1）： 439-448.

1. E. Marteau， J. E. Andrade. A model for decoding the life cycle of granular avalanches in a rotating drum[J]. Acta Geotechnica， 2017， 13： 549-555.

2. J. Rajchenbach. Development of grain avalanches[J]. Phys Rev Lett， 2002， 89： 074301.

3. D. A. Santos， C. R. Duarte， M. A. S. Barrozo. Segregation phenomenon in a rotary drum： Experimental study and CFD simulation[J]. Powder Technology， 2016， 294： 1-10.

4. N. Sepúlveda， G. Krstulovic， S. Rica. Scaling laws in granular continuous avalanches in a rotating drum[J]. Physica A： Statistical Mechanics and its Applications， 2005， 356： 178-183.

5. R. Li， H. Yang， G. Zheng， Q. C. Sun. Granular avalanches in slumping regime in a 2D rotating drum[J]. Powder Technology， 2018， 326： 322-326.

6. 刘波， 郭平， 沈岳， 等. 农业文本特征重组向量矩阵算法优化与仿真[J]. 软件， 2015， 36（4）： 6-13

7. 曹妍， 陈伟， 徐森. 图像去噪方法研究与仿真[J]. 软件， 2015， 36（4）： 33-36

8. 聂敬云， 李春青， 李威威， 等. 关于遗传算法优化的最小二乘支持向量机在MBR 仿真预测中的研究[J]. 软件， 2015， 36（5）： 40-44.

9. 刘果， 陈凡， 李剑锋， 等. 构建SDN 仿真实验平台的探讨与实践[J]. 软件， 2015， 36（6）： 103-108.

10. 王燕妮， 吴文辉. Cloudsim 3. 0 仿真流程分析[J]. 软件， 2015， 36（6）： 109-113.

11. S. H. Lin， H. Yang， R. Li， G. Zheng， V. Zivkovic， Velocities of irregular particles in a continuously avalanching surface flow within a rotating drum[J]. Powder Technology， 2018， 338： 376-382.

12. 何牧君， 郭力， 严历. 大规模并行粒子模拟系统代码级优化研究与实现[J]. 计算机与应用化学， 2008， 09： 65-68.

13. 王国强， 郝万军， 王继新. 离散单元法及其在EDEM上的实践[M]. 西北工业大学出版社， 2010： 1-2.

14. 付宏， 乌兰， 黄万风. 基于图元的三维离散元法边界建模方法[J]. 计算机集成制造系统， 2008， 014（12）： 2328-2333.

15. 徐伟， 李玉忱， 王丽. 分子动力学并行算法研究[J]. 计算机工程与应用， 2002， 14： 218-221.