椭球颗粒材料在水平转筒内混合特性的超二次曲面离散元分析*

王嗣强 季顺迎

(大连理工大学,工业装备结构分析国家重点实验室,大连 116023)

转筒中的颗粒流广泛存在于工业生产中,颗粒形状是影响颗粒流动的重要因素.本文基于超二次曲面方程描述球体和椭球颗粒的几何形态,采用离散元方法对水平转筒中颗粒物质的流动和混合特性进行数值分析,并与椭球混合过程的实验结果进行对比验证.在此基础之上,进一步研究了圆筒转速、颗粒填充分数和颗粒长宽比对混合率的影响规律.计算结果表明：颗粒材料的混合率随着转速的增加或填充分数的减小而增加.在相同转速和填充分数下,椭球颗粒的混合率高于球形颗粒.同时,长宽比为0.75和1.5时椭球颗粒具有最高的混合率.当长宽比小于0.75时,混合率随着长宽比的增加而增加；当长宽比大于1.5时,混合率随着长宽比的增加而减小.此外,椭球颗粒出现更明显的速度分层现象.颗粒长宽比改变颗粒间的接触模式和系统的密集度,增加了颗粒系统的平动而限制了转动,在一定程度上提高了外部能量向颗粒系统转化的效率.

1 引 言

颗粒流是由大量离散固体颗粒组成,区别于传统流体介质的一种复杂流动[1,2].在外界激励下,颗粒物质表现出类似固态或液态的物理力学特性,并在一定条件下相互转化,是颗粒介质类固-液相变的宏观表现[3,4].颗粒间摩擦和非弹性碰撞可有效地衰减系统能量,从而引起非平衡相变的特殊转化规律[5,6].颗粒在流动过程中重新排列,毗邻颗粒相互接触形成强度迥异的力链网络,进而改变系统的流动状态[7,8].转筒内颗粒流的数值模拟是研究颗粒流动的重要手段之一,其混合特性的机理研究对颗粒材料在工业生产和设计中的应用具有参考价值[9,10].

离散单元法(discrete element method,DEM)是由Cundall 和 Strack[11]提出,并已成为模拟颗粒流物理特性并解决相关工程问题的重要工具[12].目前,水平转筒中颗粒材料的流动和混合特性研究主要通过对混合指数、持续时间和颗粒温度等讨论不同旋转速度、方向、尺寸以及颗粒系统的填充分数、弹性模量、摩擦系数等因素的影响[13-16].实验结果表明：颗粒摩擦和圆筒转速显著影响转筒内颗粒流的动力特性,并且随着摩擦系数和转速的增加,颗粒系统的动态休止角随之增加[17].进一步考虑容器形状对混合率的影响,结果表明在不同填充分数下4叶混合器具有更好的混合效果[18].然而,大量实验和数值模拟主要集中于球形颗粒,而对旋转圆筒内非球形颗粒流的研究则相对较少[19].此外,自然界或工业中普遍是由非球形颗粒组成的复杂体系,其流动过程中的动力学行为与球形颗粒具有显著差异[20].同时,球形颗粒系统得到的宏观和细观物理力学特性难以简单地推广到非球形颗粒系统[21].

为合理地描述非球形颗粒形状,基于粘结或镶嵌的组合颗粒单元[22]、基于几何拓扑的多面体单元[23]、基于Minknowski Sum理论的扩展多面体单元[24]、基于连续函数包络的超二次曲面单元[25]、基于球谐函数拟合的星形颗粒单元[26]等不同的构造方法不断发展和完善.通过SIPHPM方法模拟多面体颗粒在二维旋转圆筒内的混合行为[27],结果表明：与球形颗粒相比,三角形、四边形和六边形颗粒具有更高的混合度.采用超二次曲面方程构造椭球和杆状颗粒,研究不同长宽比对非球形颗粒系统混合过程的影响[28,29].结果表明：高长宽比颗粒的主轴方向基本平行于颗粒的流动方向.此外,采用球形镶嵌单元和平面-柱面组合单元模拟圆柱体颗粒的流动过程并与实验和理论结果进行对比,发现不同颗粒形状的构造方法显著影响颗粒流的动力学行为[30].因此,本文采用超二次曲面单元构造不同长宽比的椭球颗粒,并通过非线性迭代方法计算单元间的接触力,进而模拟椭球颗粒在转筒中的流动和混合行为.

本文基于连续函数包络的超二次曲面方程构造不同长宽比的椭球颗粒,采用DEM建立不同颗粒形状填充的颗粒床在水平转筒作用下的数值模型,并通过椭球颗粒混合过程的实验结果[31]对数值模型进行验证.在不同转动速度、填充分数和颗粒长宽比下,通过对Lacey混合指数以及混合率的计算,探讨颗粒形状对颗粒材料混合特性的影响规律,并通过颗粒系统的平动和转动动能揭示非球形颗粒混合特性的内在机理.

2 椭球颗粒的离散元方法

针对椭球颗粒单元的几何特点,这里采用基于连续函数包络的超二次曲面方程描述颗粒的几何形态,并采用相应的非线性接触模型计算单元间的接触作用,对旋转圆筒内椭球颗粒流的混合过程进行数值模拟.

2.1 基于超二次曲面方程的椭球单元

超二次曲面方程是数学意义上描述非球形颗粒的普遍方法,并由二次曲面方程经扩展得到的函数形式[32]：

式中,a ,b 和c 分别表示颗粒沿主轴方向的半轴长,n1和n2表示颗粒表面的尖锐度参数.当参数满足n1= n2= 2时,该方程可以构造球体或椭球颗粒.图1显示在超二次曲面方程中改变不同的半轴长得到不同的球体或椭球单元模型.

图1 球体和椭球的超二次曲面离散元模型 (a) a = b =2c ；(b) a = b = c ；(c) a = b = 0.5cFig.1.Super-quadric discrete element model of spherical and ellipsoidal elements：(a) a = b = 2c ；(b) a = b = c ；(c) a = b = 0.5c.

2.2 椭球单元间的接触搜索算法

由于椭球单元是通过超二次曲面方程进行构造,因此单元间的接触判断方法主要采用超二次曲面单元的求解方法[28].考虑将单元间最短距离的优化问题转化为求解四元非线性方程组进行求解,并通过牛顿迭代算法计算单元间的重叠量[33].因此,四元非线性方程组可表示为[34]

式中 X=(x,y,z)T,∇F(X) 表示超二次曲面方程的一阶导函数,μ2表示两个颗粒单元间的外法向平行且方向相反,下标i 和j 分别表示单元i 和j.这里,牛顿迭代公式可表示为[35]

式中X(k+1)=X(k)+dX,μ(k+1)=μ(k)+dμ,∇2F(X)表示超二次曲面方程的二阶导函数.如果中间点X0满足 Fi(X0)＜ 0 且 Fj(X0)＜ 0 ,则表明两个颗粒单元发生接触,接触法向定义为n =∇Fi(X0)/‖∇Fi(X0)‖或n =-∇Fj(X0)/‖∇Fj(X0)‖,如图2所示.颗粒表面点Xi和Xj满足：Xi=X0+αn 和 Xj=X0+βn ,这里未知参数α和β可通过非线性牛顿迭代数值求解和因此,法向重叠量可表述为

图2 超二次曲面单元间的接触检测Fig.2.Contact detection between super-quadric particles.

2.3 椭球单元间的非线性接触模型

椭球颗粒间的接触力F主要是由法向力Fn和切向力Fs组成,并考虑法向力和切向力不通过单元形心引起的力矩Mn和Ms以及由单元间发生相对转动时滚动摩擦引起的力矩Mr[36,37].颗粒间的法向力主要由弹性力Fcn和黏滞力Fdn组成,可分别表示为：

单元间的切向接触力主要有弹性力Fcs和黏滞力Fds组成,同时考虑基于Mohr-Coulomb摩擦定律.则切向弹性力Fcs为

式中μs为滑动摩擦系数,δt为切向重叠量,δt,max=μs(2-ν)/2(1-ν)·δn.

单元间的切向黏滞力Fds可表示为

式中Ct为切向黏滞系数,vt,ij为单元间切向相对速度.

当单元间发生相对转动时,由滚动摩擦引起的力矩Mr可表示为

式中,μr为滚动摩擦系数,为单位相对转速,可表示为

3 椭球颗粒材料混合过程的实验验证

为验证椭球离散元计算模型的可靠性,这里分别对转速ωr为20 r/min和40 r/min的扁长形颗粒混合过程进行数值模拟,并与文献[31]的实验结果进行对比.离散元模拟中水平转筒的内径D0为200 mm,长L0为20 mm.椭球颗粒的函数参数满足：2a = 2b = 5 mm,2c = 10 mm,n1= n2= 2,主要的离散元计算参数列于表1中.同时,Lacey混合指数被用于计算转筒中颗粒材料的混合度,可表示为[39]

总计1000个椭球颗粒在圆筒容器中考虑单元的随机位置和角度并在重力作用下实现堆积.当容器内颗粒无相对运动时,圆筒开始旋转并将离散元模拟的数值结果与实验结果[31]进行对比,如图3(a)所示.同时,将不同转速ωr下的Lacey混合指数进行对比,如图3(b)和图3(c)所示.可以发现,随着转筒圈数的增加,红色和蓝色颗粒呈螺旋形进行混合,且混合指数逐渐增加并最终达到稳定.尽管离散元计算的数值结果与实验结果存在一定的偏差,但在一定程度上可以很好地反映椭球颗粒的流动过程和混合特性,同时表明基于超二次曲面DEM模型的有效性.

表1 椭球颗粒离散元模拟的主要计算参数Table 1.DEM simulation parameters of ellipsoids.

图3 不同转动圈数下颗粒混合过程的实验结果[31]和离散元数值结果的对比 (a) ωr= 20 r/min的流动图案；(b) ωr=20 r/min的Lacey混合指数；(c) ωr= 40 r/min的Lacey混合指数Fig.3.Comparison of mixing process between experiment results[31]and DEM simulation results at different rotating speeds：(a) Mixing pattern at 20 r/min；(b) Lacey mixing index at 20 r/min；(c) Lacey mixing index at 40 r/min.

4 转动速度与长宽比对混合率的影响

颗粒形状、填充分数和圆筒转速是影响颗粒材料混合特性的重要因素.这里,水平转筒的内径D0为200 mm,长L0为50 mm,并且在z 方向施加周期性边界条件,如图4(a)所示.采用超二次曲面方程构造不同长宽比的椭球颗粒,函数参数满足：σ= c /a (=b)且σ∈[0.25,3.0],如图4(b)所示.不同形态的颗粒具有相同的质量,且等体积球体半径为2 mm,颗粒数量为10000个.颗粒间的摩擦系数μs= 0.65,颗粒与筒壁间的摩擦系数μws=0.85,阻尼系数Cn= Ct= 0.1,其余计算参数列于表1中.

4.1 颗粒材料混合过程的离散元模拟

图5(a)和图5(b)显示长宽比σ= 3.0时椭球颗粒在ωr= 20 r/min和40 r/min下的混合过程.可以发现,红色和蓝色颗粒的流动呈螺旋状并出现明显的分层现象.随着转动圈数的增加,这个现象逐渐消失并且红色和蓝色颗粒基本均匀地分布在圆筒内.同时,靠近圆筒的颗粒被连续提升并流动至自由表面,并在重力作用下呈现S形流动图案.随着转动速度的增加,颗粒混合速度变快且靠近圆筒的颗粒被提升的高度增加,并最终引起更显著的S形图案.图5(a),图5(c)和图5(d)显示当ωr=20 r/min时长宽比σ= 3.0,1.0和0.5的椭球和球体颗粒在不同时刻下的混合过程.可以发现,在相同转速下,椭球颗粒相比球形颗粒更快地达到均匀混合状态.下面将进一步通过改变转动速度和颗粒长宽比来定量分析颗粒流动性能的差异.

图4 三维水平圆筒的离散元数值模型(a)和不同长宽比的椭球模型(b)Fig.4.Schematic diagram of three-dimensional horizontal drum simulated by DEM model (a) and examples of ellipsoids with different aspect ratios (b).

图5 旋转速度和颗粒长宽比对混合过程的影响 (a) ωr= 20 r/min的椭球颗粒(σ= 3.0)；(b) ωr= 40 r/min的椭球颗粒(σ=3.0)；(c) ωr= 20 r/min的球体颗粒(σ= 1.0)；(d) ωr= 20 r/min的椭球颗粒(σ= 0.5)Fig.5.The influence of rotating speed and aspect ratio on the mixing process：(a) Ellipsoids (σ= 3.0) at 20 r/min；(b) ellipsoids(σ= 3.0) at 40 r/min；(c) sphere (σ= 1.0) at 20 r/min；(d) ellipsoids (σ= 0.5) at 20 r/min.

4.2 转动速度对混合率的影响

初始时刻下红色和蓝色颗粒完全分离,即初始混合指数I0为0.当圆筒转动充足时间后,颗粒系统最终达到完全混合状态,即最终混合指数If为1.对于球体或椭球颗粒,Lacey混合指数随时间的变化过程可以通过积分函数进行拟合,该函数可表示为[41]

式中,G(Rmt) 是积分函数,可表示为G(x)=这里,Rm是混合率,可通过拟合DEM模拟数据得到.图6(a)和图6(b)所示为长宽比σ= 0.5和3.0的椭球颗粒的Lacey混合指数随时间的变化过程.其中,实线表示采用 (11) 式对模拟数据进行拟合的结果.可以发现,对于不同长宽比的椭球颗粒,Lacey混合指数最终保持稳定并趋近于1.同时,(11) 式可以较好地反映颗粒系统混合指数的变化.因此,这里对球体和长宽比σ=0.5和3.0的椭球颗粒在不同转速下的混合率进行比较,如图6(c)所示.可以发现,随着转速的增加,颗粒混合速率随之增加.在相同转速下,椭球颗粒的混合率高于球形颗粒.同时,椭球颗粒的长宽比显著影响颗粒系统的混合率.

图6 旋转速度和颗粒形状对Lacey混合指数和混合率的影响 (a) 椭球颗粒(σ= 0.5)；(b) 椭球颗粒(σ= 3.0)；(c) 混合率Fig.6.The influence of rotating speed and particle shape on the Lacey mixing index and mixing rate：(a) Ellipsoids (σ= 0.5)；(b) ellipsoids (σ= 3.0)；(c) mixing rate.

4.3 颗粒填充分数对混合率的影响

颗粒填充分数Ψ是指颗粒材料总体积与转筒体积的百分比,这里分别取圆筒内径D0= 200,220,240,280和320 mm,进而对球体和椭球颗粒的混合率随填充分数在8%—22%范围内的变化关系进行离散元分析,如图7所示.可以发现,随着填充分数的增加,球体和椭球颗粒的混合率都逐渐降低.这主要是由于靠近转筒的颗粒数目随着转筒内径的增加而增加,从而提高了外部能量通过颗粒与转筒间的摩擦向颗粒系统转化的效率,这在一定程度上加快了颗粒的运动和混合行为.此外,在相同填充分数下,椭球颗粒的混合率高于球形颗粒.颗粒长宽比对混合率的影响为主要因素,而填充分数对混合率的影响为次要因素.下面将进一步通过改变颗粒长宽比来定量分析颗粒材料混合性能的差异.

图7 ωr= 30 r/min时椭球和球形颗粒的混合率随填充分数的变化Fig.7.Mixing rate of the ellipsoid and spherical particles varies with the fill level at ωr= 30 r/min.

4.4 颗粒长宽比对混合率的影响

图8(a)显示ωr= 30 r/min和50 r/min时椭球颗粒的混合率随长宽比的变化.可以发现,球体具有最低的混合率,而椭球在长宽比σ= 0.75和1.5时具有最高的混合率.当长宽比小于0.75时,混合率随着长宽比的增加而增加；当长宽比大于1.5时,混合率随着长宽比的增加而减小.这主要是由于颗粒系统的混合过程紧密联系于颗粒间堆积和接触状态,紧密的面-面接触和有序的堆积结构具有更好的混合效果[27].然而,颗粒长宽比显著影响颗粒间的接触模式和紧密堆积程度,如图8(b)所示.其中,C0为椭球颗粒系统的初始密集度.可以发现,球形颗粒的体积分数最低,而椭球在长宽比σ= 0.75和1.5时具有最高的体积分数,并且体积分数随长宽比的变化规律与文献[42,43]中基本一致.颗粒长宽比减小了椭球间的空隙,并形成排列紧密的颗粒系统,这增加了非球形颗粒间的接触时间,并在一定程度上提高了外部能量通过颗粒间的摩擦和非弹性碰撞向颗粒系统转化的效率.

图8 ωr= 30 r/min时椭球颗粒的混合率和初始体积分数的变化 (a)混合率；(b)初始密集度Fig.8.Mixing rate and initial packing fraction under various aspect ratios at 30 r/min：(a) Mixing rate；(b) initial packing fraction.

5 颗粒材料在混合过程中的动能分布特性

5.1 颗粒系统的速度分布

水平转筒通过颗粒间的摩擦和互锁效应将外部能量转化为颗粒系统自身的动能,进而驱动或加速颗粒的运动和混合.图9显示球体和长宽比σ=0.5和3.0的椭球在不同转速下颗粒系统的速度分布.可以发现,颗粒系统被分为三部分：靠近圆筒的流动层、颗粒床中间的静态层和颗粒床表面的滑动层.同时,随着转速的增加或球形向非球形的转变,这种颗粒速度的分层现象变得更加显著.一般而言,与球形颗粒相比,非球形颗粒更容易被圆筒提升,这使得颗粒床呈现更明显的S形状.同时,非球形颗粒床表面的颗粒具有更高的滑动速度,其流动层和滑动层的相对面积增加,而处于中间的静态层面积减小,这表明整个颗粒系统具有更高的动能.

5.2 颗粒系统的平动和转动动能

图10显示球体和椭球颗粒在不同转速下的平动动能Et和转动动能Er随时间的变化关系.可以发现,当圆筒开始旋转时,颗粒系统从初始堆积状态到混合状态需要吸收更多的动能,并且椭球颗粒系统吸收的动能高于球形颗粒系统.这主要是由于椭球颗粒间具有更紧密的接触模式和互锁效应,从而形成更稳定的堆积结构.因此,促使椭球颗粒系统开始流动需要更多的动能.当转动时间大于3 s时,椭球颗粒系统基本进入稳定流动和混合状态.然而,球形颗粒系统的平动动能随时间的变化具有较大的波动性.这主要是由于球形颗粒间更容易发生相对滑动或滚动,这不利于外部能量向颗粒系统的连续转化.此外,随着转速的增加,颗粒系统的平动和转动动能都随之增加.

图9 不同转速下球体和椭球颗粒的速度分布 (a) 椭球(σ= 0.5)；(b) 球体(σ= 1.0)；(c) 椭球(σ= 3.0)Fig.9.The velocity profiles of granular bed for differently shaped particles：(a) Ellipsoids (σ= 0.5)；(b) spheres (σ= 1.0)；(c) ellipsoids (σ= 3.0).

为进一步比较颗粒形状对系统动能的影响,这里将3—10 s内颗粒系统的平动和转动动能取平均值,并分别统计转速为30 r/min和50 r/min下球体和椭球颗粒系统的平均动能Et*和Er*,如图11所示.可以发现,在相同转速下,椭球颗粒系统的平动动能高于球形颗粒系统,而转动动能低于球形颗粒系统.同时,随着长宽比的增加或减小,系统的平动动能随之增加而转动动能随之减小.这主要是由于椭球颗粒在混合过程中其主轴方向基本平行于流动方向[28],有利于外部能量通过颗粒间的摩擦转化为系统自身的平动动能.然而,颗粒长宽比引起的互锁效应不利于颗粒间的相对转动,这在一定程度上降低了外部能量转化为系统的转动动能.

图10 在不同转速下球体和椭球颗粒的平动动能和转动动能随时间的变化 (a),(b) 椭球(σ= 0.5)；(c),(d) 球体(σ= 1.0)；(e),(f) 椭球(σ= 3.0)Fig.10.Translational and rotational kinetic energy at different rotating speeds for differently shaped particles：(a),(b) Ellipsoids (σ=0.5)；(c),(d) spheres (σ= 1.0)；(e),(f) ellipsoids (σ= 3.0).

图11 ωr= 30和50 r/min时球体和椭球颗粒的平动和转动动能 (a)平均平动动能；(b)平均转动动能Fig.11.Translational and rotational kinetic energy at 30 and 50 r/min for spheres and ellipsoids：(a) Average translational kinetic energy；(b) average rotational kinetic energy.

6 结 论

基于超二次曲面方程构造球体和椭球颗粒,采用离散单元法对颗粒材料在水平转筒内的流动和混合特性进行了数值分析.通过与椭球颗粒混合过程的实验结果进行对比验证,表明基于超二次曲面模型的离散元方法适用于分析椭球颗粒流的混合特性.在此基础之上,进一步研究了转动速度、填充分数和颗粒长宽比对混合率的影响.研究结果表明：颗粒材料的混合率随着圆筒转速的增加或填充分数的减小而增加.同时,相比表面光滑的球形颗粒,椭球颗粒具有更高的混合率.长宽比σ= 0.75和1.5时椭球颗粒具有最高的混合率.当长宽比小于0.75时,混合率随着长宽比的增加而增加；当长宽比大于1.5时,混合率随着长宽比的增加而减小.颗粒长宽比显著影响系统的密集度,紧密的接触模式和有序的堆积结构导致更高的混合率.此外,外部能量通过颗粒间的摩擦和互锁效应转变为颗粒系统自身的动能,同时椭球颗粒出现更明显的速度分层现象.颗粒长宽比提高了外部能量向颗粒系统平动动能的转化效率,而其引起的互锁效应限制了外部能量向系统转动动能的转化.