颗粒物质在竖直振动U 形管中迁移的离散元方法模拟

郭 宇，凡凤仙，白鹏博，刘 举

（上海理工大学 能源与动力工程学院，上海 200093）

颗粒物质是由大量离散的固体颗粒相互作用而形成的具有内在联系的复杂体系，作为一种特殊的物质形态，其在受到周期性振动激励时呈现出复杂而奇特的行为[1]，如隆起和对流[2-3]、表面波和斑图[4-5]、尺寸分离[6-7]及毛细效应[8-10]等。近年来，颗粒在竖直振动U 形管内发生迁移的现象引起了研究者的关注[11-20]，其表现为：对填充有颗粒物质的U 形管施加竖直正弦振动，在适当的振动强度下，经过一段时间之后，U 形管内颗粒运动状态将达到稳定，此时U 形管两分支内颗粒柱形成一个稳定的高度差，该稳定高度差与初始高度差无关。

早在1976 年，Gutman[21]就发现了竖直振动U 形管中颗粒的迁移现象，并将其归因于U 形管底部水平方向上间隙气体压力梯度作用；1991 年，Rajchenbach[2]在真空条件下进行实验，发现颗粒迁移现象仍然发生，这表明间隙气体不是颗粒迁移的必要条件，进而提出了颗粒迁移的对流机理。即由于颗粒的随机速度梯度引起颗粒由疏松区向密实区运动，形成颗粒对流，导致了颗粒的迁移。然而根据该机理，U 形管一个分支内颗粒柱高度增加，另一个分支内颗粒柱高度下降，直至降为0，颗粒迁移才会停止，这与实际情况不符。1998 年，Ohtsuki 等[11]以插有竖直隔板而底部连通的颗粒床为研究对象，实验发现了隔板位置、颗粒粒径、振动频率、振动强度对颗粒迁移具有重要影响，同时，利用二维离散元方法（discrete element method，DEM）模拟获得了颗粒迁移受颗粒量、颗粒与管壁间摩擦系数的影响特性，展现了颗粒对流现象以及管壁对颗粒的平均剪切应力；此外，基于稳定状态下竖直方向颗粒微元受力方程，推导出隔板两侧颗粒表面高度差的理论表达式。2001 年，Akiyama 等[12]在不同气压下针对隔板两侧等宽的颗粒床进行实验，发现增加气压对颗粒迁移起促进作用，但气压不是颗粒迁移的必要条件；进而在忽略环境空气条件下开展DEM 模拟，给出了颗粒运动状态达到稳定时颗粒的平均速度分布，但缺少对颗粒迁移动态过程的展示和分析。近年来，随着振动颗粒行为研究兴起，围绕U 形管内颗粒迁移的报道增多。2007 年，King 等[13]给出了浸没在水中的颗粒物质在插有隔板的容器内的迁移特性，并基于压力梯度作用机理对颗粒迁移进行了解释；在此基础上，2009 年，Clement 等[14]对这一系统中颗粒物质的迁移进行了较为全面的研究。2009 年，Sánchez 等[15]通过实验研究给出了不同水平振动分量下U 形管两分支颗粒柱高度差随时间的演变情况；同时，基于循环流化思想建立了高度差增长模型，该模型适用于细颗粒、低频率和有空气存在的情况，粗大颗粒或无空气时的迁移机制和规律仍有待探讨。2011 年，Darias 等[16]通过实验发现，随着振动强度的增大，U 形管内颗粒物质呈现出三种不同的行为，即振动强度很小时颗粒物质表现出固体特性，不发生流动；在中等振动强度下，U 形管两分支颗粒柱高度差随时间呈指数增长，指数增长率随振动强度的变化曲线呈S 形；当振动强度很大时，颗粒表现出黏性流体行为，指数增长率突变为0。同年，Pérez 等[17]实验研究了颗粒与管壁间摩擦系数对U 形管两分支颗粒柱高度差指数增长率的影响，发现随着摩擦系数的增大，指数增长率减小，颗粒发生迁移所需的振动强度下限增大；将实验结果与基于循环流化思想的理论模型[15]预测结果相对比，发现该模型在摩擦系数较低时和实验符合较好，但在摩擦系数高时和实验的差异很大。2012 年，Sánchez 等[18]对U 形管中颗粒迁移进行了实验和二维DEM 模拟，指出颗粒迁移可归因于两类机理：压力梯度作用机理与壁面摩擦作用机理，两者在不同的条件下起主导作用。2013 年，Darias 等[19]利用实验和二维DEM 模拟展示了间隙气体作用可以忽略的情况下，竖直振动的准二维U 形管两分支内颗粒柱高度差随时间的演变过程，并通过二维DEM 模拟探讨了管宽度、颗粒密度、颗粒与壁面间恢复系数对稳定高度差的影响。2015 年，Sánchez 等[20]考虑到颗粒与U 形管之间延迟耦合，对基于循环流化思想的模型[15]进行了改进，得到了能够预测出大振动强度下高度差指数增长率突变为0 的理论解析式，并据此研究了管内颗粒总质量、黏性系数、颗粒密度、U 形管截面积、振动强度和频率对指数增长率的影响规律。

综上所述，竖直振动U 形管中颗粒迁移的相关研究主要集中在通过实验、理论和数值模拟探讨颗粒迁移的宏观效果（如高度差增长率和稳定状态下的高度差），虽然提出了颗粒迁移的间隙流体作用、对流、循环流化和壁面摩擦作用机理，但这些机理是不同的研究者在不同的研究条件下提出的，各机理的相互关联和作用范围尚不清楚。特别是对流和壁面摩擦作用机理均是为了解释无间隙流体存在的情况下颗粒的迁移现象而提出的，颗粒对流与壁面摩擦的关联有待探究。受测量手段、实验条件和颗粒物质对外部激励的非线性响应特性的限制，这些问题难以通过实验和理论研究得到完全解决，而DEM 模拟能够跟踪到每个颗粒的运动，且在参数设置上具有出色的灵活性，便于从颗粒尺度上探讨不同条件下颗粒物质的行为规律，是实验和理论研究的有力补充。虽然已有研究中也涉及到DEM 模拟，但这些模拟均为二维模拟，即在建模时将球形颗粒按片状处理，使得数值模拟结果难以准确反映实验现象，并且研究中只关注高度差演变的宏观效应，缺少对颗粒尺度动力学行为的探讨。基于此，本文对无间隙流体时竖直振动U 形管中颗粒的行为规律开展三维DEM 模拟，并结合文献中的实验结果，验证模拟的正确性，进而对有、无摩擦情况下U 形管两分支内颗粒动力学行为进行对比分析，为深入理解振动颗粒行为提供基础。

1 DEM 模型与数值计算方法

1.1 DEM 模型

根据牛顿第二定律，颗粒系中任意颗粒i 的运动方程可写为

式中：mi和Ii分别为颗粒i 的质量和转动惯量；和分别为颗粒i 的速度和角速度；t 为时间；g为重力加速度；N 为与颗粒i 相接触的颗粒的个数；和分别为颗粒i 和与它相接触的颗粒j 之间的法向和切向作用力；为颗粒j 对颗粒i 的切向作用力产生的力矩；为颗粒j 对颗粒i 的滚动摩擦力产生的力矩。

由于颗粒形状和物性的差别，描述接触颗粒的相互作用力的模型有多种[22-24]。本文利用黏弹性接触模型[25]描述法向作用力，采用修正的Cundall-Strack 模型[26]描述切向作用力，切向作用力由切向变形位移和切向动能损失计算，但切向作用力的极值受到颗粒摩擦系数与法向作用力乘积的限制，当切向作用力大于该极值时，两颗粒在接触面发生滑动。因此，有

此外，式（4）的积分路径path 为两颗粒接触期间在接触点的相对位移。

根据定向恒转矩模型[30]，可写为

式中：μr为滚动摩擦系数；为颗粒i 与j 的相对角速度，。

对于颗粒i 与管壁的相互作用，按照上述模型将管壁视为粒径无限大的颗粒处理。

1.2 数值计算方法

基于DEM 模型，借助开源颗粒系统数值模拟软件LIGGGHTS[31]对竖直振动U 形管中颗粒物质的动力学行为进行数值模拟。在数值计算中，考虑摩擦时，采用与Darias 等[19]的实验一致的U 形管结构（图1）、颗粒物性参数、振动条件，具体参数如表1 所示；不考虑摩擦时，则令摩擦系数为0。Darias 等[19]的实验系统中，U 形管内填充有聚甲醛树酯球形颗粒，信号发生器与音频放大器相连，驱动扬声器产生振幅为（3.9±0.1）mm，角频率为（73.5±0.1）rad/s 的竖直正弦振动；利用数码照相机记录实验照片。表1 中颗粒物性参数——杨氏模量、泊松比、恢复系数和摩擦系数，依据实验颗粒类型从文献[32-33]中获取。

图 1 U 形管示意图Fig. 1 Schematic diagram of the U-tube

表 1 数值模拟参数Tab.1 Parameters used in numerical simulations

数值模拟时，首先将一定数目的颗粒随机地加入计算区域，然后令颗粒在重力作用下沉降，由于重力作用、管壁的约束、颗粒与管壁的相互作用和颗粒之间的相互作用，U 形管内颗粒总动能先增大后减小并最终颗粒趋于总动能为0 的松弛状态，如图2 所示。在颗粒达到松弛状态后，在竖直方向对U 形管施加振幅为A，角频率为ω 的正弦振动，使得U 形管的竖直方向位移az随时间t 的变化关系满足az=Asin（ωt），并将施加振动的时刻作为计时0 点，即t=0。

图 2 施加振动前颗粒总动能随时间的变化关系Fig.2 Total kinetic energy of the particles as a function of time before applying vibration

对于给定的颗粒系统，DEM 模拟的准确性及其时间成本主要受时间步长的影响，时间步长Δt 可依据颗粒碰撞时间tc进行选择。对于无阻尼、非黏附碰撞，碰撞时间tc可由式（8）估算[29]。

式中：meff为有效质量，；νimp为碰撞速度。

在DEM 模拟中，时间步长Δt 常取为tc/50 ～ tc/10。本文取νimp=1 m/s 估算碰撞时间，得到tc≈ 3.2×10-5s，因此，时间步长取为Δt = tc/30 ≈ 10-6s。

2 结果与讨论

2.1 竖直振动U 形管中颗粒物质的行为模式

图3 给出了不同情况下竖直振动U 形管中颗粒物质的行为模式。图3（a）为文献[19]的实验结果，图3（b）～3（d）依次为采用与实验相同的参数（有摩擦），颗粒与管壁间摩擦系数为0，颗粒与颗粒间摩擦系数为0 时的数值模拟结果。图3 中，T 为U 形管振动周期，T= 2π/ω= 0.085 5 s。数值模拟结果中，颗粒颜色表征振动起始时刻颗粒位置。对比图3（a）和图3（b）可见，DEM 模拟结果和实验结果吻合良好，两者都表明，对初始时刻两分支颗粒柱等高的U 形管施加竖直振动后，一个分支颗粒柱高度逐渐上升，另一个分支颗粒柱高度逐渐下降，从而两分支颗粒柱出现高度差，且高度差逐渐增加，直至最终达到一个稳定值。与实验相比，数值模拟可以展现更为丰富的颗粒流动过程信息，具有实验无法比拟的优越性。由图3（b）可以看出：a. U 型管的振动带动颗粒运动，颗粒失稳而发生重排，竖直分支内颗粒进入水平段，由于两分支颗粒堆积的细微差别，从其中一个分支进入水平段的颗粒数目较另一个分支多，因此，两分支颗粒柱出现高度差；b. 在U 形管右分支颗粒向左分支迁移进而形成高度差的过程中，右分支内侧壁面附近的颗粒下降速度更为迅速，左分支中心区域颗粒的上升速度大于两侧壁面附近颗粒的；c. U 形管两分支呈现出截然不同的颗粒对流现象，相对右分支，左分支内颗粒对流较弱。图3（c）和3（d）则反映出无论是颗粒与管壁间摩擦系数为0 还是颗粒与颗粒间摩擦系数为0，竖直振动U 形管中颗粒定向迁移并最终稳定的现象消失，这表明摩擦是颗粒在U 形管内发生迁移的必要条件。对比图3（c）和3（d）发现，两种情况下颗粒的行为模式差别很大。图3（c）表明，在颗粒与管壁间无摩擦的情况下，两分支均是靠近内侧壁面的颗粒下降，靠近外侧壁面的颗粒上升，形成颗粒对流。可见，U 形管内颗粒对流现象的产生与器壁摩擦没有必然的联系。图3（d）显示，在颗粒与颗粒间摩擦系数为0 的情况下，无颗粒对流出现。此外，需要指出的是，在多次实验中由于颗粒初始堆积状态的随机性引起颗粒迁移方向的随机性，文献[19]中多次重复实验表明，U 形管左分支颗粒柱高度Hl大于右分支颗粒柱高度Hr的概率与U 形管右分支颗粒柱高度大于左分支颗粒柱高度的概率相等，但颗粒柱高度差的最终稳定值（|Hl-Hr|）不变；在数值模拟中，笔者将颗粒随机加入到U 形管两分支内不同高度范围的区域，并使颗粒沉降，以产生不同的初始颗粒堆积状态，得到了和实验一致的结论。

图 3 不同情况下竖直振动U 形管中颗粒物质的行为模式Fig.3 Behavior modes of granular matter in the vertically vibrating U-tube under different conditions

2.2 竖直振动U 形管中两分支颗粒柱高度差随时间的演变

为了定量描述施加竖直振动后U 形管左、右分支内颗粒柱高度Hl和Hr随时间的演变，将U 形管竖直分支沿其高度方向（z 向）划分为若干个连续的单元，各单元高度Δz = 2dp，利用下式计算各单元内颗粒的容积份额

式中：φ（z）为中心高度为z 的单元内颗粒的容积份额；Np（z）为中心高度为z 的单元内颗粒数目；Vp为颗粒体积，；Vc为单元体积，。

取满足φ（z）＞0.1 的单元的中心高度最大值z 为颗粒柱高度。据此，获得Hl和Hr，并采用ΔH=Hl-Hr表示U 形管两分支颗粒柱高度差。

图4 给出了数值模拟得到的U 形管左、右分支颗粒柱高度及高度差随时间的变化关系。由图4（a）可以看出，在有摩擦的情况下，U 形管两分支颗粒柱高度的演变过程可分为两个阶段：在0 ＜ t ＜ 300T 时，随着振动时间的延长，U 形管左分支颗粒柱高度逐渐增加，而右分支颗粒柱高度逐渐减小，从而两分支颗粒柱高度差增加，可称为增长阶段；t ＞ 300T 时，U 形管两分支颗粒柱高度达到稳定，两分支形成一个稳定的颗粒柱高度差，可称为稳定阶段。由图4（b）则可以看出，在颗粒与管壁间摩擦系数为0 时，U 形管两分支颗粒柱高度在初始高度附近上下波动，使得高度差时而为正，时而为负。图4（c）的结果表明，在颗粒与颗粒间摩擦系数为0 时，U 形管两分支颗粒柱高度几乎保持不变，高度差约为0。

2.3 颗粒竖直方向平均速度分布

为了探讨U 形管两分支内颗粒的运动特性，将两分支沿水平方向划分为等体积的5 个区域，分别计算施加竖直振动后特定时间内不同区域颗粒竖直方向平均速度。

式中：N'为计算时间范围内保存数据文件的个数，本模拟中每个周期内以相等时间间隔保存40 个数据文件；和分别为第i 个数据文件中计算区域内颗粒平均速度和颗粒数目。对于颗粒与管壁、颗粒与颗粒间均有摩擦的情况，时间选择为0～300T，即U 形管两分支高度差增长阶段；对于颗粒与管壁或颗粒与颗粒间无摩擦的情况，时间选择为0 ～ 500T，即整个计算时间。

图 4 不同摩擦系数条件下U 形管两分支颗粒柱高度及高度差随时间的变化关系Fig. 4 Granular column heights in two branches of the U-tube and the height difference as a function of time under different friction coefficients

图5 给出了不同摩擦系数条件下U 形管左、右分支颗粒竖直方向平均速度沿水平方向的分布情况。图中，横坐标为相对位置，即各区域中心距U 形管相应分支内侧壁面的距离与U 形管分支宽度的比值，其中，xin为内侧壁面横坐标，xout为外侧壁面横坐标，xout-xin=Wtube。由图5（a）可知，在颗粒与管壁、颗粒与颗粒间均有摩擦的情况下，左分支中紧靠内侧壁面的颗粒具有向下的平均速度，且速度值很小0.05 m/s），紧靠外侧壁面的颗粒平均速度几乎为0，中间的3 个区域颗粒平均速度均向上，且中心区域颗粒平均速度最大；右分支中除中心区域颗粒和紧靠右侧具有向上的平均速度外，其余区域颗粒的平均速度均向下，两分支内颗粒竖直方向平均速度特性与左分支内颗粒柱高度增加、右分支内颗粒柱高度减小的现象相符。由图5（b）可以看出，颗粒与壁面间摩擦系数为0 时，U 形管两分支内对应区域颗粒竖直方向平均速度较为接近，都表现出距离内侧壁面较近的区域内竖直方向颗粒平均速度向下，而其余区域竖直方向颗粒平均速度向上，这种速度分布特性使得U 形管两分支呈现类似的颗粒对流现象。由图5（c）可以得到，颗粒与颗粒间摩擦系数为0 时，U 形管两分支壁面邻近的颗粒竖直方向平均速度几乎为0；中心区域的颗粒竖直方向平均速度向上，与中心区域邻近的两区域颗粒具有竖直向下的平均速度，3 个区域内颗粒的净质量通量约为0，因此，U 形管两分支颗粒柱无明显高度差。此外，竖直振动U 形管中的颗粒受到3 种力的作用，分别是管壁对颗粒的作用力、颗粒之间的相互作用力以及颗粒自身的重力。图5（c）中，U 形管通过竖直壁面的摩擦作用和水平壁面的法向力作用对系统输入动能，引起U 形管内颗粒发生运动。但由于颗粒之间不存在摩擦，U 形管中心区域颗粒因摩擦而耗散的动能为0，与壁面附近颗粒以及颗粒之间有摩擦的情况相比动能耗散较小。因此，中心区域颗粒速度达到最大值，且该最大值较图5（a），5（b）中速度最大值大。

图 5 不同摩擦系数条件下颗粒竖直方向平均速度随相对位置的变化关系Fig. 5 Average vertical velocity of particles as a function of relative position under different friction coefficients

3 结 论

利用DEM 方法对竖直振动U 形管中颗粒的迁移运动进行三维数值模拟，将有摩擦、无摩擦情况下颗粒的动力学行为进行对比分析，得到以下结论：

a. 利用基于DEM 方法的三维数值模拟，再现了实验中得到的U 形管两分支颗粒柱高度差随时间演变的历程，展现了颗粒的对流现象等颗粒尺度动力学信息。

b. 当颗粒与管壁间摩擦系数为0 时，U 形管两分支颗粒柱高度交替增减，两分支内颗粒对流现象类似；而当颗粒与颗粒间摩擦系数为0 时，U 形管两分支颗粒柱高度差几乎为0，颗粒对流现象消失。

c. 摩擦系数均不为0 时，U 形管两分支颗粒的竖直方向平均速度分布差异很大；当颗粒与管壁间或颗粒与颗粒间摩擦系数为0 时，U 形管两分支内颗粒的竖直方向平均速度存在一定的一致性。