文章信息
基于离散元方法的U形管中颗粒不稳定性
周英贵1，高信2，白鹏博2，凡凤仙2,3
1盐城工学院环境科学与工程学院，江苏 盐城 224051；2上海理工大学能源与动力工程学院，上海 200093；3上海市动力工程多相流动与传热重点实验室，上海 200093
引用本文
周英贵, 高信, 白鹏博, 等. 基于离散元方法的U形管中颗粒不稳定性[J]. 化工进展, 2024, 43(11): 6083-6090.
DOI：10.16085/j.issn.1000-6613.2024-0785
摘要
对填充颗粒物料的U形管施加竖直振动时，颗粒将呈现由一个分支向另一个分支定向迁移的不稳定性,该现象对颗粒物料的输运产生重要影响。为探究U形管中颗粒不稳定性的动力学行为及机理，基于离散元方法对颗粒动力学过程进行可视化模拟，分析了初始颗粒填充状况对颗粒不稳定性的影响，揭示了颗粒不稳定性的发展与饱和阶段速度场演化特性，考察了壁面摩擦对颗粒不稳定性的影响。结果表明，颗粒不稳定性与自放大效应有关，其发展阶段受到颗粒初始填充状况的影响；而最终饱和阶段U形管两分支颗粒表面高度差与颗粒初始填充状况无关，这与实验报道一致。在颗粒不稳定性的发展阶段，U形管向下运动过程中，水平段颗粒向两竖直分支膨胀，颗粒传输到颗粒量较多而底部颗粒速度较小的分支；在饱和阶段，U形管向下运动时，颗粒量较多分支中的惯性效应引起颗粒底部速度较大，抑制了水平段颗粒向其中的传输，导致颗粒表面高度差的增长停止。壁面滑动摩擦系数过小或过大时，颗粒定向迁移的现象均无法发生。
颗粒物料的输运在化工领域普遍存在，这种输运过程不可避免受到机械振动的影响。因而，有必要对振动激励下颗粒物料的行为规律加以研究。在外加振动作用下，颗粒物料将发生颗粒固体向颗粒流体的转变，并伴随着复杂而奇特的动力学行为。其中，U形管中颗粒不稳定性吸引了诸多研究者的兴趣，其具体为：向竖直放置的U形管填充一定量的颗粒后，施加竖直周期性振动，当振动强度高于某一阈值时，颗粒失去稳定性，进而发生由一个分支向另一个分支的定向迁移，使得两分支颗粒表面产生一定的高度差。
U形管中颗粒不稳定性研究已有很长的历史。然而，迄今对其机理的认识尚存争议，对其影响因素的研究尚不充分，描述颗粒不稳定性的理论模型仍不完善。例如，在机理方面，Gutman认为颗粒不稳定性的成因是间隙气体作用；Rajchenbach通过无间隙气体情况下U形管中颗粒动力学行为实验，提出颗粒膨胀引起的对流效应是颗粒不稳定性的原因；Chung与Liaw基于分子动力学模拟，提出颗粒物料由固体向流体的相变导致颗粒不稳定性。再如，在影响因素方面，通过实验与数值模拟考察了微小水平振动扰动、竖直振动强度、壁面粗糙度以及U形管分支宽度、颗粒密度、颗粒恢复系数对颗粒不稳定性的影响。同时，Sánchez等利用循环流化理论，建立了描述U形管两分支颗粒表面高度差增长速率的理论模型，但该模型未能考虑壁面摩擦效应。此外，一些学者探索了U形分割容器内颗粒不稳定性，可作为U形管中颗粒不稳定性研究的重要参考，相关研究涉及U形分割容器内浸没于水中颗粒的不稳定性以及U形分割容器存在间隙气体时颗粒的不稳定性。
已有相关工作在U形管/U形分割容器中颗粒不稳定性方面取得了一些有价值的成果。然而，由于U形管中颗粒动力学过程中涉及参数和因素繁多，如U形管几何参数、颗粒物性参数、振动条件、有无间隙流体及间隙流体种类和压力等，U形管中颗粒不稳定性关键的动力学过程和特性仍不清晰，亟需开展针对性研究。本文在前期工作对离散元模型进行验证的基础上，通过离散元方法的可视化模拟，分析颗粒不稳定性的动态过程与成因，揭示颗粒不稳定性发展阶段与饱和阶段速度场的演变特性，考察壁面摩擦系数对颗粒不稳定性的影响，为深入理解U形管中颗粒不稳定性行为规律及机理提供基础。
1
模型与方法
1.1
物理模型
Darias等采用如图1所示的准二维U形管开展了颗粒不稳定性实验，并使用数码相机记录实验图像（见图2）。图中U形管几何参数为：宽度Lx=(107±1)mm、厚度Ly=(3.10±0.05)mm、高度Lz=(362±1)mm，分支宽度Wtube=(33±1)mm。向U形管填充直径为(3.00±0.05)mm、密度为(1.42±0.02)g/cm3的颗粒，颗粒数量为(1695±5)个。U形管的振动方程为式(1)。
图1 U形管示意图
式中，az为z向位移；a为振幅，a=(3.9±0.1)mm；f为频率，f=(11.7±0.016)Hz。
1.2
离散元模型
本文采用开源颗粒系统离散元模拟软件Liggghts（版本号：3.5.0）开展U形管中颗粒动力学数值模拟。自从Cundall与Strack提出描述颗粒运动的离散元模型以来，离散元方法作为颗粒动力学研究的有效工具，在颗粒行为及机理研究得到了广泛应用。对图2所示的颗粒系统，考虑颗粒的平动与转动，建立颗粒运动的离散元模型，见式(2)、式(3)。
图2 实验得到的颗粒不稳定性动力学过程
式中，m为质量；I为转动惯量；v为速度；ω为角速度；t为时间；g为重力加速度；F为接触力；M为力矩；N为与颗粒i接触的颗粒的数目。下角标i、j以及ij分别表示颗粒i、j以及颗粒i和j之间的值；n和t分别表示切向和法向分量；r表示滚动摩擦有关的值。
接触力的表达式为式(4)、式(5)。
式中，δ为重叠量；ρ为弹性参数；A为耗散参数；e为单位矢量；μs为滑动摩擦系数；G为剪切模量；υ为泊松比；Reff为有效半径；s为积分变量；path为由接触点相对位移确定的积分路径。
耗散参数的法向分量An可采用其与恢复系数的函数关系求出，在此基础上，可计算耗散参数的切向分量At。At、ρ、G、Reff的表达式分别为式(6)~式(9)。
式中，Y为杨氏模量；R为颗粒半径。
此外，滚动摩擦引起的力矩Mr,ij的表达式为式(10)。
式中，μr为滚动摩擦系数；ωij=ωi-ωj为相对角速度。
1.3
计算方法
数值模拟采用Lx×Ly×Lz=107mm×3.1mm×362mm、Wtube=33mm的U形管，U形管填充有1700个直径为3mm、密度为1.42g/cm3的颗粒；杨氏模量设置为2.8MPa，泊松比和恢复系数分别设置为0.35和0.95；颗粒之间的滑动和滚动摩擦系数分别为0.4和0.06。沿z方向对U形管施加式(1)所示的正弦振动，振幅a=3.9mm、频率f=11.7Hz。时间步长的选择策略为：利用无阻尼碰撞时间的理论表达式估算碰撞时间，从而选择一个远小于碰撞时间的时间步长。基于此，可采用10-6s作为时间步长；模拟中发现，若采用更小的时间步长，计算结果几乎不发生改变，这表明该时间步长是适宜的。
2
结果与讨论
2.1
颗粒不稳定性动力学过程
图3给出了U形管中颗粒不稳定性动力学过程。图中，T=1/f为振动周期，按初始位置对颗粒进行着色。模拟采用的颗粒与壁面间滑动与滚动摩擦系数分别为0.4与0.06。尽管文献曾给出U形管中颗粒表面高度差演变过程，但本模拟采用的颗粒初始条件与其不同。文献为了验证离散元模型，将等量的颗粒分别随机添加到左右分支中并使颗粒自由沉降得到颗粒堆积状态，从而产生与实验尽可能一致的初始颗粒条件；而本文为了探索颗粒不稳定性是否与自放大机理有关，通过自定义文件方式对颗粒进行初始化，使颗粒空间分布关于yoz平面（见图1）对称。陈伟中等利用自放大机理解释了竖直振动床中颗粒的对流和隆起，认为随机堆积颗粒表面具有微小的隆起，使得微小隆起下部颗粒受到比邻近颗粒更大的压力，当颗粒床下降时隆起底部的颗粒将下凸，从而颗粒层表面颗粒向隆起运动；当颗粒床上升时，隆起底部首先接触颗粒床底面，从而呈现更大的隆起。由图3可以看出，颗粒初始空间分布对称情况下，随着振动的进行，颗粒出现重排，颗粒堆积变得更为致密，U形管两分支颗粒表面高度均略有下降，同时颗粒分布失去对称性，颗粒不稳定性产生，见t=140T时的结果；随后，颗粒空间分布的不对称性加剧，颗粒由U形管一个分支向另一个分支发生明显的迁移运动，同时伴随着竖直分支内侧颗粒向下、外侧颗粒向上的对流运动，且颗粒表面下降的分支中对流更为剧烈，如图中t=280T、t=420T时所示；t>420T时，宏观上颗粒不稳定性进入两分支颗粒表面高度差基本不变的饱和阶段，然而，微观上颗粒运动并未停止，颗粒对流仍在发生。在文献中颗粒随机分布情况下，U形管两分支颗粒表面高度差迅速增大，在t=140T时，已形成了显著的高度差；而本模拟条件下，t=280T时U形管两分支颗粒表面高度差仍很小。这表明颗粒初始空间分布对高度差增长速率产生重要影响，对称分布能够延迟高度差的发展。其原因与颗粒物料对外界微小作用的敏感性、非线性响应等特性有关，数值计算过程的微小误差使得颗粒分布偏离对称性，这种偏离逐渐放大，从而引起颗粒不稳定性。此外，本模拟结果也表明，间隙气体作用不是颗粒不稳定性发生的必要条件。
图3 U形管中颗粒不稳定性动力学过程
2.2
初始颗粒填充状况对颗粒不稳定性的影响
已有文献报道了U形管中颗粒不稳定性饱和阶段的颗粒表面高度差只与系统参数和振动条件有关，与初始颗粒填充状况无关；但缺乏对不同初始颗粒填充条件下颗粒表面高度差演化过程的探讨，而这一信息对于认识颗粒不稳定性非常重要。为对颗粒表面高度差进行定量，沿竖直方向将U形管两分支z≥0部分划分成高度为6mm的若干个网格，由下而上计算各网格内的颗粒体积分数，取颗粒体积分数不大于10%的首个网格下方相邻网格中心的z坐标作为颗粒表面高度。图4给出了U形管两分支颗粒表面高度以及高度差随时间的演变。模拟采用的壁面摩擦系数与2.1节相同。图4(a)~(c)对应的初始颗粒填充状况分别为两分支颗粒表面初始高度差为0（与图3一致）、初始高度差约为最终高度差、初始高度差显著大于最终高度差。此外，图中H表示颗粒表面高度；∆H=Hleft-Hright表示高度差；Hleft与Hright分别为左分支与右分支颗粒表面高度；∆H0与∆H∞分别为初始与最终高度差；d=2R为颗粒直径。由图4(a)可知，当初始高度差为0且初始颗粒空间分布对称时，在t/T<250时，两分支颗粒表面高度差接近0；之后，随振动时间的增加，一个分支中颗粒表面上升，另一分支中颗粒表面下降，导致两分支颗粒表面高度差迅速增大，此阶段为颗粒不稳定性的发展阶段；直至t/T>500时，U形管左右分支颗粒表面高度基本保持不变，两分支颗粒表面高度差几乎恒定，此阶段为颗粒不稳定性的饱和阶段。计算发现，t/T<250时无量纲高度差∆H/d的平均值为1.4，标准差为1.2；在饱和阶段∆H/d的平均值为18.2，标准差为1.7。图4(b)中结果表明，当初始高度差约为颗粒表面最终稳定的高度差，即∆H0/d≈18时，施加振动过程中，两分支颗粒表面高度以及高度差基本保持不变。分析显示，∆H/d的平均值为17.9，标准差为1.0。由图4(c)可以看出，当初始高度差大于饱和高度差时，颗粒由初始颗粒表面较高的分支中向初始颗粒表面较低的分支中迁移，使得前一分支颗粒表面高度下降，而后一分支颗粒表面高度上升，从而两分支颗粒表面高度差减小，当t/T>300时，颗粒不稳定性达到饱和阶段，在该阶段∆H/d的平均值为19.9，标准差为2.5。三种初始条件下，颗粒表面高度差的偏差在-2.5%~6.6%之间，考虑到统计误差，可以得到U形管中颗粒不稳定性饱和阶段颗粒表面高度差与颗粒初始填充状况无关，这与实验报道的现象一致，进一步验证了离散元模型的可靠性。
图4 不同初始颗粒填充状况下颗粒表面高度及高度差随时间的演变
2.3
颗粒不稳定性过程中速度场演变
为了揭示颗粒不稳定性的发展阶段与饱和阶段的颗粒动力学特性及机理，针对初始高度差为0且初始颗粒空间分布对称的颗粒填充状况[见图3与图4(a)]，分别给出了颗粒不稳定性的发展阶段（t=280T~281T）与饱和阶段（t=780T~781T）一个周期内颗粒瞬时速度场的演变，如图5与图6所示。
图5 颗粒不稳定性发展阶段一个周期内速度场的演变
图6 颗粒不稳定性饱和阶段一个周期内速度场的演变
由图5可以看出，t=280T时，U形管由平衡位置以最大速度向上运动，受U形管竖直壁面的摩擦力与水平壁面的支持力作用，颗粒向上运动，随着高度的增加，颗粒速度增大；t=280.25T时，U形管处于最高位置且运动速度为0，U形管竖直分支中除表面处颗粒具有较大速度外，其余颗粒速度很小；t=280.5T时，U形管向下运动至平衡位置处并达到速度最大值，受U形管竖直壁面的摩擦力与水平壁面的压力作用，颗粒整体向下运动，且竖直分支中位置越低的颗粒速度越大，同时由于颗粒之间的相互碰撞及壁面约束，在U形管底部少量颗粒向上运动；当t=280.75T时，U形管向下减速运动至最低位置处，U形管水平段中的颗粒受到两竖直分支向下运动的颗粒挤压作用而向上运动，与上侧壁面碰撞，颗粒动能发生耗散，并向两侧膨胀，由于右分支内侧颗粒速度更大，水平段颗粒与右分支内侧颗粒碰撞的动能耗散更大，抑制了颗粒向右分支的传输，从而更多颗粒进入左分支；而后，U形管转而向上运动，当t=281T时，U形管重新回到平衡位置，颗粒速度场与t=280T时类似，但两分支颗粒表面高度差变大。根据颗粒不稳定性的膨胀机理，竖直振动使得颗粒体积发生膨胀，从而周围的颗粒可以自由运动至膨胀区域，改变颗粒的堆积结构，引发颗粒对流的同时，颗粒不稳定性现象发生。图5(a)对应于颗粒体积增大的膨胀阶段，颗粒运动有序性良好，而在图5(d)所示的颗粒体积减小的压缩阶段，颗粒呈现明显的由一个分支到另一分支传输的特征，这表明膨胀机理无法解释U形管中颗粒不稳定性。
由图6可知，t=780T时，U形管以最大速度由平衡位置向上运动，颗粒受U形管带动向上运动，且随着高度的增加，颗粒速度增大；t=780.25T时，U形管向上减速运动达到最高位置且速度减为0，颗粒速度方向向上，但速度值很小；t=780.5T~780.75T时间内，U形管由平衡位置，向下做减速运动，壁面的摩擦和压力作用使得t=780.75T时颗粒随U形管向下运动；t=780.75T时，U形管处于最低位置处，其速度减小为0，U形管左分支颗粒数目更多而具有更大的惯性，颗粒继续向下运动挤压水平段颗粒，使得其与右分支颗粒发生碰撞，造成右分支颗粒的动能耗散，从而右分支中除颗粒柱上部之外的区域颗粒速度均很小；当t=781T时，U形管向上加速运动至平衡位置处，其速度达到最大值，受壁面摩擦和水平段支持力的影响，颗粒总体呈现出向上运动的特性，颗粒速度场与两分支颗粒表面高度差均与t=780T时接近。
2.4
壁面摩擦对颗粒不稳定性的影响
为着重探讨U形管壁面滑动摩擦的影响，将颗粒之间以及颗粒与壁面之间的滚动摩擦系数设置为0，颗粒之间滑动摩擦系数保持为0.4，改变颗粒与壁面之间的滑动摩擦系数，在初始时刻两分支颗粒表面高度接近相等的随机堆积条件下进行数值模拟，获得不同壁面摩擦下两分支颗粒表面高度差随时间的演变，如图7所示。总体而言，当壁面摩擦系数μs,w=0.1时，颗粒先快速向左分支迁移，使得∆H>0，而后迅速向右分支迁移，引起∆H减小，直至∆H<0，∆H保持基本不变一段时间后，颗粒又开始向左分支迁移；μs,w=0.2时，∆H始终小于0且其值波动较大；μs,w=0.3时，∆H在0附近上下波动；然而，μs,w=0.4时，∆H为正，且先增大，后迅速减小为0，接着再次增大；μs,w=0.5时，与μs,w=0.3时类似，∆H在0附近上下波动，但波动幅度略小；μs,w=0.6、0.7、0.8时，颗粒发生由右分支向左分支的定向迁移，且能够维持一定的高度差，特别是μs,w=0.7时，最终高度差约为20d，显著高于其他摩擦系数对应的结果；μs,w=0.9与1.0时，∆H又呈现出在0附近上下波动的特性。其原因可解释为：当壁面摩擦系数很小时，颗粒与壁面之间的摩擦力不足以支撑竖直分支中颗粒重力，因而难以维持稳定的高度差；当壁面摩擦系数很大时，颗粒与壁面之间的滑动摩擦做功耗散的能量很大，使得振动传递给颗粒的能量被颗粒与壁面的摩擦和碰撞所耗散，也难以带动颗粒产生定向迁移。因此，颗粒不稳定性由发展阶段向饱和阶段的平稳过渡需要适中的壁面摩擦。同时，注意到图7中μs,w=0.4时颗粒不稳定性行为与图3及图4(a)中截然不同，这表明滚动摩擦是影响颗粒不稳定性的重要因素。再者，文献采用峰值粒径为270μm、标准差为50μm的颗粒研究了U形管中颗粒不稳定性，发现壁面越粗糙，两分支颗粒表面高度差增长得越缓慢，与本文结果不同。其原因可解释为：文献中颗粒粒径较小，黏附力不可忽略，壁面越粗糙，颗粒与壁面的黏附力越小，导致颗粒与壁面的作用减弱，不利于高度差的发展；然而，本文模拟采用的颗粒粒径为3mm，远大于文献中的颗粒，此时颗粒与壁面间的黏附力可以忽略不计。这表明，颗粒不稳定性动力学行为及影响因素的复杂性。此外，文献曾针对横截面为圆形、下端为弯头的U形管中颗粒不稳定性进行研究，考虑到这种U形管在工业过程中更为常见，未来的工作应针对该结构U形管开展研究，以为实际应用提供更为有效的指导。
图7 不同壁面摩擦下颗粒表面高度差随时间的演变
3
结论
利用离散元方法对U形管中颗粒不稳定性进行可视化模拟，展示了初始颗粒空间分布对称情况下颗粒不稳定性过程，分析了初始颗粒填充状况对颗粒不稳定性的影响，揭示了颗粒不稳定性的发展阶段与饱和阶段的瞬时速度场，考察了壁面摩擦对颗粒不稳定性的影响。通过本文研究，得到如下结论。
（1）颗粒分布对称能够减缓颗粒不稳定性的发展，但不能阻止颗粒不稳定性的发生；颗粒不稳定性的发展阶段与自放大效应有关，并受到颗粒初始条件的重要影响，但其饱和阶段两分支颗粒表面高度差与颗粒初始条件无关。
（2）在颗粒不稳定性的发展阶段，U形管向下运动至最低位置时，水平段中的颗粒受到两竖直分支向下运动的颗粒挤压作用而向上运动，并与上侧壁面碰撞而向两侧膨胀，更多颗粒传输到速度较小的较高颗粒表面分支；在颗粒不稳定性的饱和阶段，U形管处于同一位置时，颗粒表面高度较高的分支中颗粒整体惯性更大，颗粒继续向下运动而挤压水平段中颗粒，抑制了颗粒向该分支的迁移。
（3）U形管中颗粒不稳定性过程对壁面摩擦系数非常敏感，壁面滑动摩擦系数过小或过大时，颗粒难以发生稳定的定向迁移运动，壁面滚动摩擦系数对颗粒不稳定性产生显著影响。
作者简介
第一作者：周英贵，博士，副教授，研究方向为大气污染控制、颗粒流与多相流。
通信作者：凡凤仙，博士，教授，研究方向为颗粒流与多相流。
（扫码关注我们）
邮发代号：82-311
订阅热线：010-64519502
网址：http://hgjz.cip.com.cn
欢迎您分享、点赞、收藏、在看
