颗粒毛细效应影响因素的离散元分析1)

张华腾 凡凤仙 王志强

(上海理工大学能源与动力工程学院,上海 200093)

(上海理工大学上海市动力工程多相流动与传热重点实验室,上海 200093)

引言

在外加扰动作用下颗粒物质的敏感性和非线性响应等使得颗粒物质表现出复杂而奇特的动力学行为,日益吸弓着不同领域科研工作者的兴趣[1-6].近几年,文献报道了颗粒的毛细效应[4-6],即将一根细管插入填充有颗粒物质的容器中并对管施加竖直振动时,颗粒物质在管内迅速上升并最终达到一个稳定的高度的现象.该现象为实现颗粒物料的逆重力输送提供了一种潜在的技术途径,弓起了研究者的高度关注.

自Liu 等[7]对颗粒毛细效应进行实验报道以来,研究者迅速开展了针对颗粒毛细效应的实验、理论和数值模拟研究.首先,在实验研究方面,获得了管振动强度、振幅、频率、管插入深度、管口形状、管内颗粒填充高度、间隙气体、管内径、管截面形状、颗粒直径对颗粒上升高度的影响[7-11];对比了管顶开口与封闭两种情况下颗粒的上升过程,发现管顶封闭对颗粒上升起促进作用,分析了管径及管长对这种促进作用的影响[12].其次,在理论研究方面,将管内颗粒视为一个整体,建立了颗粒受力模型,对颗粒的上升高度和速度进行了理论分析,并提出了颗粒上升的“空穴填充”机理[7];建立了颗粒上升高度随时间变化的半经验关系式,进一步分析了颗粒上升的“空穴填充”机理,并且基于能量平衡假设,推导出颗粒最终毛细上升高度的理论表达式[13];基于一个振动周期内管内颗粒的流进和流出特性,建立了颗粒毛细上升高度的理论模型,该模型能够在一定程度上解释振幅、频率和管径对颗粒最终毛细上升高度的影响,但无法解释较大振动频率下颗粒最终毛细上升高度随频率下降的现象,以及较小管径下颗粒最终毛细上升高度随管径增加的现象[14].再者,在数值模拟方面,利用离散元方法,获得了颗粒毛细效应全过程[15];揭示了颗粒毛细效应的对流机理,并发现了颗粒最终毛细上升高度随管径的变化特性[4];展现了颗粒在竖直方向的运动规律[5];获得了颗粒填充率的变化特性,以及振幅、振动强度、恢复系数对颗粒上升过程的影响规律[16];展示了不同初始填充高度下颗粒毛细效应的过程,探讨了颗粒毛细上升过程中管内颗粒柱高度、颗粒速度场和填充率分布随时间的演变特性,获得了颗粒毛细效应过程中由容器传输到管内的颗粒的占比分布[6].

虽然已有研究在颗粒毛细效应的宏观效果、动力学过程和内在机理上取得了诸多进展.然而,由于颗粒系统参数的复杂性和多样性,目前颗粒毛细效应研究尚处于起步阶段,宏观的实验研究仍不系统,尚缺乏对不同操作参数下管径影响的探讨,而这方面研究对推动颗粒毛细效应在颗粒逆重力输运上的应用具有重要意义.此外,已有研究中提出的颗粒毛细上升高度的理论模型都是建立在把大量颗粒视为整体的基础之上,事实上管内不同高度和径向位置颗粒的运动速度往往不同,因此现有理论模型无法对颗粒的输运提供有效的定量指导.离散元方法通过跟踪每一个离散单元(颗粒)的运动以实现对整个颗粒系统的模拟,是从颗粒尺度探究颗粒动力学行为的行之有效的方法[17-19].本文将在课题组前期工作的基础上[4-6,15],利用离散元方法,开展颗粒毛细效应过程数值模拟,获得不同容器宽度、管振幅和频率下颗粒毛细上升高度随管径的变化规律,以探明颗粒毛细效应的影响因素,并为颗粒物料逆重力输运方案的优化提供科学依据.

1 离散元模型与计算方法

1.1 离散元模型

在离散元模型中,颗粒系中任意颗粒i的平动与转动方程由牛顿第二定律给出,即

式中,mi和Ii分别为颗粒i的质量和转动惯量;vi和ωi分别为颗粒i的速度和角速度;t为时间;n为与颗粒i相接触的颗粒的个数;Fij,n和Fi j,t分别为颗粒i和与它相接触的颗粒j之间的法向作用力和切向作用力;Mi j,t为颗粒j对颗粒i的切向作用力产生的力矩;Mij,r为颗粒j对颗粒i的滚动摩擦力产生的力矩;g为重力加速度.

采用黏弹性接触模型[20]描述法向作用力,采用修正的Cundall-Strack 模型[21]描述切向作用力,则

式中,ρ 为弹性参数,其是杨氏模量Y、泊松比υ、有效半径Reff=RiRj/(Ri+Rj)的函数;Ri和Rj分别为颗粒i和j的半径;δi j,n为颗粒i与j之间的法向重叠量,δi j,n=Ri+Rj-|ri-rj|;ri和rj分别为颗粒i和j的位置矢量;An为法向耗散系数,根据其与恢复系数ε 之间的关系[22-24],可采用Pad´e近似方法进行计算[25];en为法向单位矢量,en=(rj-ri)/|rj-ri|;µs为滑动摩擦系数;G为剪切模量,G=Y/[2(1+υ)];At为切向耗散系数,At=AnY/(1-υ2)[26];vij,t为接触点上颗粒i与j的切向相对速度;et为切向单位矢量.此外,式(4)中的积分路径为两颗粒接触期间在接触点的相对位移.

弹性参数ρ 的计算式为

根据定向恒转矩模型[27],Mi j,r可写为

式中,µr为滚动摩擦系数;ωi j为颗粒i和j的相对角速度,ωi j=ωi-ωj.

对于颗粒i与管壁面及容器壁面之间的相互作用,将壁面视为半径无限大的颗粒,利用上述模型进行处理[28-30].

1.2 计算方法

基于DEM 模型,借助开源颗粒系统数值模拟软件LIGGGHTS[31]对颗粒毛细效应动力学行为开展数值模拟.本文数值模拟采用的颗粒系统示意图如图1 所示.其中,容器的横截面为正方形,其边长为w;容器底面位于z=-18 mm 处,容器高度足够高;管竖直插入容器中心,管底面中心位于坐标原点,管内径为D,管壁厚为0.3 mm,管长度足够长;颗粒直径d=0.6 mm.数值模拟采用的颗粒物性参数见表1.数值模拟过程为:首先将一定数目的颗粒随机加入容器内(含管内区域),颗粒在重力作用下发生沉降,由于重力作用、颗粒与壁面以及颗粒与颗粒之间的相互作用,系统内颗粒总动能先增大后减小,趋于总动能为0 的松弛状态;待颗粒达到松弛状态后,对管施加竖直方向的周期性正弦振动,使得管底面位置zb随管振动时间t的变化规律为

式中,a为振幅,f为频率.

数值模拟时,时间步长的选择非常重要,时间步长越小,计算精度越高,但计算时间成本也越高.因而,在选择时间步长时,应兼顾计算精度和计算成本.在离散元模拟中,通常选择颗粒碰撞时间的1/50～1/10 作为时间步长[32-33].为确定时间步长,采用无阻尼、无黏性碰撞模型估算颗粒碰撞时间[34]

图1 数值模拟采用的颗粒系统示意图Fig.1 Schematic of the granular system used in numerical simulations

表1 数值模拟采用的颗粒物性参数Table 1 Physical properties of the particles used in numerical simulations

式中,meff为颗粒对的有效质量,meff=mimj/(mi+mj);vimp为碰撞速度.取vimp=1 m/s 作为碰撞速度参考值,则有tcol≈6.0 × 10-5s.因此,选择时间步长Δt=tcol/30 ≈2.0×10-6s.在数值模拟时,发现采用更小的时间步长,数值模拟结果不发生改变,这表明时间步长的选择是合理的.

2 结果与讨论

2.1 不同管径下颗粒毛细效应过程

在容器宽度与粒径比w/d=40、管振幅与粒径比a/d=14.33、管振动频率f=12 Hz 时,得到的典型管径与粒径比(D/d)下,颗粒毛细效应过程快照如图2 所示.图2(a)～图2(c)分别为D/d=3.33,8.33,15时的结果.图中,T为管振动周期,T=1/f;Hc为毛细上升高度;竖直方向速度接近0(|vz| ＜0.05 m·s-1)的颗粒用蓝色表示,其余颗粒中,向上运动的用红色、向下运动的用绿色表示.由图2(a)可见,在D/d很小时,由于堵塞严重,毛细上升速度较慢,且管内颗粒发生结团现象,使得颗粒柱出现中断.需要说明的是:图中给出的是振动周期整数倍时刻的结果,此时管内颗粒处于最密堆积状态[4,6],为结团和中断最弱的时刻,非整数倍周期时刻结团和中断更为严重.由图2(b)可以看出,D/d增大到8.33 时,起初毛细上升高度增大迅速,随后毛细上升高度的增大逐渐减缓,直至达到一个稳定的最终毛细上升高度;由于D/d仍较小,受堵塞影响,管径方向几乎不存在速度梯度;在12T时刻,颗粒柱顶部的少量颗粒因相互作用较弱,自由程较大,运动行为不同于管内其余颗粒而向下运动,管内其余颗粒间频繁的摩擦与非弹性碰撞不断消耗颗粒系统的动能,使得颗粒聚集在一起,颗粒自由程较小,颗粒呈现随管向上运动的状态;图中24T～108T各时刻,管内颗粒均向上运动,这与颗粒体系的能量耗散有关,经过更多振动周期,颗粒间的摩擦与非弹性碰撞使得颗粒柱顶部颗粒呈密堆积状态,管内颗粒柱顶部与下部的颗粒运动协同起来,因而颗粒柱整体向上运动.由图2(c)可知,在D/d=15 条件下,颗粒毛细上升高度先迅速增大,而后在t≥24T时出现起伏,这与图2(b)差异显著.由图2(c)还可以看出,t≥24T时,管内颗粒柱分离为速度截然不同的两层:管以最大速度向上运动时,管内上层颗粒向下运动,与管的运动方向相反;管内下层颗粒随管向上运动,且管壁附近颗粒的速度大于管中心区域颗粒的速度,在管径方向存在明显的速度梯度.

图2 数值模拟得到的颗粒毛细效应过程快照Fig.2 Snapshots of the granular capillarity process obtained by numerical simulations

图3 不同管径下颗粒毛细上升高度随时间的演变Fig.3 Evolution of capillary height of the granular materials with time at different tube diameters

为定量描述颗粒毛细上升高度(即管内颗粒表面高度和容器内颗粒表面高度之差)随时间的演变,将管内z＞0 区域沿z轴划分为高度为5d的连续的单元,计算各单元内颗粒填充率φ,并定义满足φ ＞0.1 的单元的最大高度为管内颗粒表面高度.类似地,将容器内z＞0,r＞D/2 区域沿z轴划分为高度为d的连续的单元,计算各单元内颗粒填充率,采用相同的方法确定容器内颗粒表面高度.在此基础上,即可获得颗粒毛细上升高度.不同D/d条件下颗粒毛细上升高度随时间的演变如图3 所示.其中,图3(a)给出了D/d=5,6.67,8.33 时的结果,图3(b)给出了D/d=10,15,20,30 时的结果.从图3 可以清楚地看到,颗粒毛细效应过程中伴随着颗粒毛细上升高度的波动,且颗粒毛细上升速度和最终毛细上升高度受到管径的影响.由图3(a)可知,管径较小时,随着管径的增加,颗粒毛细上升速度更快,所能达到的最终毛细上升高度也更大.总体看来,管径较大时,随着管径的增加,颗粒毛细上升速度减小,所能达到的最终毛细上升高度也更小,如图3(b)所示.同时,需要指出的是:D/d=15 时,颗粒毛细上升高度呈现出反复波动的非稳定特征,这一特征产生的原因可能是竖直振动过程中颗粒的填充率分布和堆积结构、颗粒与颗粒以及颗粒与壁面之间的碰撞等因素影响了颗粒的自组织行为.复杂体系中颗粒物质的行为模式对其内部结构及扰动具有高敏感性,其物理机制十分复杂,有待进一步深入研究.

2.2 临界管径的影响因素分析

Fan 等[4]在仅改变管径而保持其他参数不变时,发现在毛细效应能够发生的管径范围内,颗粒最终毛细上升高度随管径的增大,先增大后减小.为了叙述的方便,本文将颗粒最终毛细上升高度随管径增大而增大的管径区间的上限称为临界管径.基于此,临界管径下颗粒最终毛细上升高度达到最大,一旦管径大于临界管径,颗粒最终毛细上升高度下降.临界管径影响因素研究对于利用颗粒毛细效应输运颗粒物料的效果优化具有重要意义.

2.2.1 容器宽度的影响

为探究容器宽度对临界管径的影响,在a/d=14.33,f=12 Hz 条件下,对容器宽度与粒径比w/d=40,60 时颗粒毛细效应过程开展数值模拟,获得颗粒最终毛细上升高度随D/d的变化关系,如图4 所示.为体现颗粒最终毛细上升高度随管的振动而上下波动的特性,图中利用点表示颗粒最终毛细上升高度的平均值,利用误差条表示标准差.由图可知,w/d=40 情况下,临界管径对应的D/d=9.17,在D/d≥10 时,颗粒表现出液体的性质,由于管内颗粒的流化,阻碍了颗粒的上升,使得在临界管径附近颗粒最终毛细上升高度急剧下降.确切地说,w/d=40 情况下,在D/d=9.17 时,在D/d=10 时,两者相比颗粒最终毛细上升高度下降了67.6 mm.在w/d=60 情况下,颗粒最终毛细上升高度随着D/d的增大,同样呈现出先增大后减小的特征,并且在管径超出临界管径时颗粒最终毛细上升高度急剧下降.然而,两种容器宽度下,临界管径并不相同,在w/d=60 情况下,临界管径对应的D/d=10,相应的颗粒最终毛细上升高度为138.1 mm;在D/d=11.67 时,颗粒最终毛细上升高度为60.7 mm,与D/d=10 时相比下降了77.4 mm.究其原因是:有摩擦的容器侧壁边界的存在对于诱导容器内颗粒的对流起到了不可缺少的作用[4],随着容器宽度的增大,容器侧壁的边界效应降低,容器内颗粒的对流强度随之减小,这造成了由颗粒对流传输至管口的颗粒的质量通量的减小,从而削弱了管内颗粒的流化,因此,颗粒由堵塞到流化的转折点延后,表现为临界管径随容器宽度的增大而增大.

图4 容器宽度对临界管径的影响Fig.4 Effect of container width on critical tube diameter

2.2.2 振幅的影响

图5 振幅对临界管径的影响Fig.5 Effect of vibration amplitude on critical tube diameter

为探究振幅对临界管径的影响,在w/d=40 的条件下对振幅与粒径比a/d=14.33 与20 两种情况下颗粒毛细效应开展数值模拟,得到颗粒最终毛细上升高度H∞c随管径与粒径比D/d的变化情况,如图5 所示.在a/d=20 情况下,临界管径对应的D/d=11.67,对比显示,该情况下的临界管径明显大于a/d=14.33 时的临界管径;同时,颗粒最终毛细上升高度H∞c的最大值可达到203.2 mm,也明显大于a/d=14.33 时H∞c的最大值133.5 mm.振幅增大弓起临界管径增大的原因是:振幅增大,管的振动强度增大,管内壁附近的颗粒受到了更大的壁面剪切力作用,使得更多颗粒能够跟随管运动,减小了管内颗粒自由流动的区域,相当于减少了有效管径,因而增强了管内颗粒的堵塞效应,使得颗粒堵塞能够在更大的管径下影响颗粒毛细效应.在D/d=5 时,a/d=20情况下颗粒最终毛细上升高度比a/d=14.33 情况下略低,这同样能够通过颗粒堵塞进行解释.由于振幅较大时,更易出现堵塞,管径很小时,管内颗粒堵塞严重抑制了颗粒的上升,使得最终颗粒毛细上升高度下降,特别是当管径过小时,颗粒出现完全堵塞,颗粒毛细效应将无法发生[4].此外,D/d＞5 时,a/d=20 情况下的颗粒最终毛细上升高度均大于a/d=14.33 时的值,这与文献中实验报道的颗粒最终上升高度随振幅的增加而增加的现象一致[11],该现象可解释为振幅增加弓起管的振动强度增大,使得管外壁对容器内的颗粒产生了更大的剪切力作用,促进了容器内颗粒的对流,增强了容器内颗粒向管口的传输,有效地提高了颗粒毛细上升高度.

2.2.3 频率的影响

图6 频率对临界管径的影响Fig.6 Effect of vibration frequency on critical tube diameter

在w/d=40,a/d=14.33 条件下,考察f=12 Hz与f=15 Hz 时最终稳定毛细上升高度随D/d的变化情况,结果如图6 所示.对比发现,D/d≤8.33 时,在较高频率下,略低;D/d≥9.17 时,在较高频率下,也较高.特别是当10 ≤D/d≤15 时,两种频率下颗粒最终毛细上升高度差异显著:f=15 Hz 条件下,颗粒毛细效应动力学处于受堵塞影响的区域,保持很高的值;f=12 Hz 条件下,颗粒毛细效应动力学处于不受颗粒堵塞影响的流化区域,较低.与振幅相比,频率对颗粒毛细效应的影响更为复杂,频率升高弓起振动强度增大、振动周期减小.由于本文研究中采用的两种频率相差不大,振动强度的增大对颗粒毛细效应动力学的影响起主导作用.较高频率下,振动强度也较大,靠近管内壁的颗粒因受到更大的壁面剪切力作用更难以与管壁分离,导致管内颗粒的有效流动区域变小,使得管内颗粒更容易发生堵塞.频率较高时,更为严重的堵塞效应造成较小管径(D/d≤8.33)下颗粒更难以上升,从而产生更小的最终毛细上升高度.同时,频率较高时,堵塞效应能够在更大的管径(D/d≤15)下发挥作用,使得临界管径增大.此外,由于频率增大,增加了管的振动强度,因而增强了管外颗粒在对流机制作用下向管口的传输,导致D/d＞15 时,较高频率下,颗粒最终毛细上升高度较大.然而,需要说明的是:若振动频率在很大的范围内变化,则管振动周期的改变将可能对颗粒动力学行为起主导作用.来自文献的实验结果表明,保持振幅不变时,随着频率的增大,颗粒最终上升高度先增加后减小,存在一个使得最大的最佳频率值[11].更大的频率变化范围内,临界管径的变化特性,有待进一步研究.

3 结论

基于离散元方法,对颗粒毛细效应动力学行为开展数值模拟研究,给出了不同管径下颗粒毛细效应过程中颗粒竖直方向速度和毛细上升高度随时间的演变规律,分析了容器宽度、管的振幅和频率对颗粒最终毛细上升高度随管径变化的影响特性.通过本文研究,得到如下结论:

(1)在容器宽度与粒径比为40、管振幅与粒径比为14.33、管振动频率为12 Hz 情况下,管径与粒径比D/d=3.33 时,管内颗粒堵塞严重,使得颗粒上升缓慢,并造成颗粒柱中断;D/d=8.33 时,起初颗粒毛细上升高度迅速增加,随后颗粒毛细上升高度增速减缓,管内颗粒在管径方向几乎不存在速度梯度;D/d=15 时,随着颗粒毛细上升高度的增大,管内颗粒柱分离为速度截然不同的两层.

(2)在颗粒毛细效应能够发生的管径范围内,存在一个对应于颗粒最终毛细上升高度最大值的临界管径,当管径小于临界管径时,颗粒最终毛细上升高度随管径的增加而增加,当管径大于临界管径时,颗粒最终毛细上升高度随管径的增大而减小;临界管径受到容器宽度和振动参数的影响.

(3)临界管径随着容器宽度的增加而有所增大;增加振幅和适当增大频率能够有效促进临界管径的增大,从而将更多的颗粒输运到更高的高度,这对于利用颗粒毛细效应实现颗粒物料的逆重力输运具有重要意义.