从质量流向漏斗流转变过程中的动力学分析

球床模块式高温气冷堆(简称高温气冷堆)采用氦气为冷却剂，以全陶瓷型包覆铀燃料制成球形颗粒，与石墨颗粒一起堆积形成反应堆堆芯，意在以循环方式实现不停堆连续运行。堆芯颗粒流属于准静态流动，其力学现象非常丰富。确定单颗粒的动力学特征，确定介尺度结构特征，是国际学术界的竞争热点，如图1所示。采用筒仓颗粒流开展物理实验和数值计算，是研究堆芯颗粒流的主要方法。

图1

图1   堆芯颗粒流中的多尺度结构

Fig.1   Multiple scales in the flow of fuel element spheres in the pebble-bed high temperature reactor (HTR)

在重力作用下，筒仓内的颗粒流型包括质量流、漏斗流和混合流。理想筒仓内颗粒物质是按照质量流形式运行的，并遵循“先进先出”的原则[1]。然而，在实际运行中，大部分筒仓都是以混合流的形式进行物料运输和加工[2]，这是由于颗粒动力学的不均匀性引起的[3]。混合流会增加局部剪切速率、加剧颗粒间的磨损、危害筒仓的安全运行[4]。比如，在漏斗流区域与边壁区域的颗粒流型变化较大，与颗粒的旋转直接相关，统计颗粒的动力学特征，探究其变化规律，受到人们的普遍关注。

近年来，人们通过不同的实验方法测量了筒仓卸料过程中固体颗粒的平移速度[5-11]。Gentzler等[12]利用核磁共振成像技术测量了三维锥形筒仓中粉末颗粒的平移速度。结果发现，在漏斗流区域，筒仓中心线上的颗粒平移速度满足指数函数分布。孔口附近的颗粒具有较大的平移速度和较低的体积分数，这表明该区域内颗粒间的碰撞磨损较小。Guillard 等[13]通过基于连续 X 射线的成像技术和基于傅里叶变换的图像处理方法，测量了矩形筒仓卸料过程中椭球形颗粒的平移速度场。观察到漏斗流区域内静止颗粒的排列方向具有一致性，呈径向分布。而在孔口附近，颗粒的排列方向呈轴向分布。这种排列方式不同会导致颗粒的受力存在差异，从而造成不同区域内颗粒间的磨损是不相同的。颗粒的速度包括平移速度和旋转速度，为了了解筒仓内颗粒物料的动力学特征，研究它们在筒仓内的分布是很重要的。Zhu 等[14]利用42 mm 直径的球形探测器，实时测量了三维锥形筒仓卸料过程中的颗粒偏转角。发现筒壁附近的颗粒累积偏转角较大，这说明了边壁效应会导致靠近筒壁的颗粒产生更多旋转磨损。赵颖等[15]利用直径30 mm的无线探测小球作为示踪颗粒，测量了圆柱形筒仓卸料过程中颗粒的旋转速度，来量化颗粒间旋转产生的磨损。结果发现，在筒壁附近旋转速度较大，约为中心位置颗粒的两倍。表明了靠近筒仓中心位置处旋转磨损量小，而靠近筒壁处旋转磨损量大。

随着计算机技术的发展，离散单元法 (distinct element method, DEM) 越来越多地应用于颗粒流研究[16-18]。DEM 可以实时跟踪每个颗粒的运动，获得颗粒的旋转运动信息。Feng等[19]通过基于线性接触模型的DEM，模拟了方形筒仓内椭球形颗粒的混合流动。观察到在卸料过程中颗粒的旋转速度分布是稳定的，数值大小与颗粒的空间位置有关。在靠近筒壁的剪切区，由于颗粒与筒壁之间的摩擦作用导致了颗粒的平移运动缓慢，而颗粒间相对切向的旋转运动十分剧烈，这可能会造成颗粒的磨损增加。Saxén等[20]利用DEM方法模拟了漏斗筒仓卸料过程中的球形颗粒流动模式，讨论了边壁静摩擦因数对漏斗流型和边壁颗粒运动的影响。结果发现颗粒在除孔口外的整个筒仓区域具有相同的平移速度，并且在缺料过程中颗粒运动由旋转速度控制。随着颗粒与边壁间静摩擦因数的减小，边壁颗粒的平移速度和旋转速度数值均增大。

本工作通过DEM 模拟获得的数据来研究颗粒在不同流型区域中的平移和旋转运动。模拟仿真的详细参数和模拟过程在实验方法部分给出。根据仿真结果，划分了卸料过程中筒仓中的不同流型区域，分析了不同流型区域下颗粒平移速度和旋转速度的空间分布及其空间相关性，最后计算了颗粒的滚动贡献率，并探究其分布规律。

1 实验方法

1.1 装置系统

实验装置包括一个总高度为750 mm的半圆柱体筒仓和一个倾斜角为30°的锥形漏斗，筒仓的横截面直径为300 mm，漏斗孔口的直径为80 mm。选择80 mm的孔口尺寸是为了减少由于堵塞而导致流动中断的可能性。整个筒仓由聚甲基丙烯酸甲酯（PMMA）制成，边壁厚度为5 mm。卸料的球形颗粒由聚甲醛（POM）制成，其直径为9 mm，总填充高度为600 mm。图2(a)为筒仓的正视图。为了统计颗粒的平均运动信息，将筒仓表面切片划分为30×60个网格。同时为了保证每个网格内至少存在一个颗粒，每个网格的尺寸为10 mm×20 mm×10 mm。此外，边壁区域和中心区域的位置如图2(a)中标注。图2(b)为筒仓的三维视图。O点位于筒仓底部直径的中心，以 O 点为坐标原点建立笛卡儿坐标系。图2(c)为单颗粒运动，包括平移和旋转，影响颗粒之间相互作用形成的介尺度结构，进而影响宏观流动行为。

图2

图2   筒仓及初始床层的几何结构以及单颗粒的运动

Fig.2   The geometry of the silo with an initial packed bed and movement of single particles

1.2 仿真过程及参数

首先将筒仓孔口进行密封处理，在筒仓顶部处随机产生颗粒，在重力作用下沉降至筒仓中。当颗粒床层的高度达到600 mm时，加载结束，颗粒总数约为2.7×105个。之后等待颗粒达到静止状态，以避免加载过程中颗粒的运动对排放过程造成影响。达到静止状态的条件是筒仓中所有颗粒的平移速度为零，此时所有颗粒到达稳定位置，床高和床荷达到平稳状态。最后打开筒仓底部的孔口，开始对筒仓颗粒进行卸料，直到筒仓内所有颗粒都流出。

仿真过程中，设置采集颗粒运动数据的时间间隔为0.01 s。为了保证仿真的准确性，最小时间步长设置为2.71×10-5 s，小于颗粒的Rayleigh波传播时间5.82×10-5 s[21]。DEM模拟的其他系统参数如表1所示。其中，密度、杨氏模量和泊松比为材料的固定参数。颗粒间的静摩擦因数μp-p通过测量颗粒静止堆积的最大休止角确定[22]。恢复系数e 通过测量颗粒自由落体过程中的下落高度和回弹高度确定[23]。材料的滚动摩擦因数是通过球体材料在相应材料的板上自由滚动，并测量球体的速度以及覆盖的距离得到的[24]。

表1   模拟仿真中的物理参数及其值

Table 1  Physical parameters and their values of the simulations

新窗口打开| 下载CSV

1.3 接触模型

DEM模拟的关键是选择合理的物理接触模型。Hertz-Mindlin with RVD Rolling Friction模型调整了滚动摩擦的计算方法，以确保三个维度均具有适当的功能而不影响计算时间，特别适用于对颗粒旋转特性有着严格要求的强旋转系统[25]。因此，本文采用基于Hertz-Mindlin with RVD Rolling Friction的物理接触模型来模拟三维锥形筒仓中球形颗粒的卸料过程。在该模型中，颗粒之间的法向力计算基于Hertzian接触理论，而颗粒之间的切向力计算基于 Mindlin-Deresiewicz 接触理论[26]。根据接触颗粒之间的相对旋转速度，计算颗粒的滚动摩擦[27]如式（1）所示。

τ=-μrFnR*ω̂relFn=43E*R*δn32(1)

式中，扭矩τ表示颗粒的滚动摩擦；μr为滚动摩擦因数；Fn为颗粒间的法向力；R* 为颗粒间的等效半径；ω̂rel为相对旋转速度的单位向量；E* 为等效杨氏模量；δn 为法向重叠量。

等效半径R* 和等效杨氏模量E* 的计算如式（2）所示。

1R*=1Ri+1Rj1E*=1-νi2Ei+1-νj2Ej(2)

式中， i和j表示接触的两个颗粒；Ri 和 Rj 为接触颗粒的半径；νi 和νj 为泊松比；Ei 和Ej 为杨氏模量。

1.4 质量流指数

质量流指数（MFI）是用来划分筒仓内颗粒流流型的通用标准[28]，MFI=0.3作为区分质量流与漏斗流的临界值。当MFI>0.3时，颗粒运动满足质量流模式；当 MFI<0.3 时，颗粒运动满足漏斗流模式。 MFI的定义如式（3）所示。

MFI=vwZvcZ(3)

式中，vw(Z)为边壁区域颗粒的 Z 轴平移速度； vc(Z)为中心区域颗粒的 Z 轴平移速度。

1.5 滚动贡献率

颗粒间的相对切向运动，具体表现为颗粒间的相对旋转运动[29]。利用滚动贡献率Rr[30]可以表达颗粒旋转运动对总运动的贡献，也可以量化颗粒之间的磨损，计算公式如式（4）所示。

Rr=vrotationvtranslation+vrotation=ωRωR+v(4)

式中，ω是颗粒的旋转速度；R是颗粒的半径；v是颗粒的平移速度。当Rr>0.3 时表明此时颗粒由旋转造成的颗粒间相对速度比例较大，颗粒表面磨损较大，不可忽略其对宏观流型的影响。

2 实验结果与讨论

2.1 流型转换

为了验证模拟的有效性，对比了物理实验和数值模拟的流动过程，发现模拟结果与实验结果一致。图3为卸料过程中不同时刻的筒仓颗粒流型，卸料过程总时长为32 s。可以看出，筒仓中心区域的颗粒以较快的平移速度向下运动，而筒仓底部靠近边壁的颗粒以非常缓慢的平移速度向孔口运动。从颗粒流型变化来看，筒仓内的颗粒在卸料过程中处于混合流状态[31]。筒仓上部床层的颗粒流型较为均匀，属于质量流。随着床层高度的降低，质量流逐渐过渡到漏斗流。

图3

图3   筒仓卸料过程中的流型

Fig.3   The silo flow pattern diagrams during the discharge process

图4(a)为t=1 s时边壁区域和中心区域颗粒的Z轴方向速度的变化曲线。可以看出，当t=1 s时，筒仓内的流型转化高度为63 mm。进一步，取采样时间间隔为0.1 s，统计卸料过程中不同时刻的流型转换高度，如图4(b)所示。结果表明，稳定的混合流状态下，筒仓内的流型转换高度基本保持稳定[32]。值得注意的是，当t=12 s后，筒仓内的MFI值均小于0.3，这表明此时筒仓内的颗粒处于漏斗流状态。由于在2~6 s时间内，筒仓内的床层具有足够高度，且流量稳定。因此，为了探究混合流状态下颗粒的动力学特征，接下来仅关注2~6 s时间内筒仓内颗粒的运动信息。当高度在80 mm以上时，颗粒处于质量流状态；当高度在55 mm以下时，颗粒处于漏斗流状态。

图4

图4   卸料过程中的MFI与流型转换高度

Fig.4   MFI and the flow pattern transformation height during the discharge process

2.2 颗粒平移速度

图5为2~6 s内筒仓颗粒的平均平移速度分布，其值为平移速度的模长。如图5(a)所示，颗粒的平移速度在孔口区域达到最大值，并向上扩散减小。在筒仓底部靠近边壁的“滞留区”，颗粒平移速度几乎为0。取Z=0、50、100和400 mm为代表，统计筒仓内颗粒平移速度的径向分布，采用最小二乘法对原始数据进行拟合，拟合误差小于5%，如图5(b)所示。结果表明，颗粒平移速度的径向分布关于筒仓的中心轴对称，中心区域的颗粒平移速度具有最大值，并沿 Y 方向减小。同样地，取Y=0、50、100和150 mm为代表，统计筒仓内颗粒平移速度的轴向分布，如图5(d)所示。结果表明，在靠近筒仓边壁区域，随着质量流向漏斗流转变，颗粒的平移速度明显增大；在靠近筒仓中心区域，随着质量流向漏斗流转变，颗粒的平移速度略微下降。

图5

图5   2~6 s内颗粒平均平移速度在空间上的分布

Fig.5   Spatial distribution of average translation velocity in 2—6 s

2.3 颗粒旋转速度

图6为筒仓内颗粒的平均旋转速度分布，其值为旋转速度的模长。与平移速度类似，颗粒的旋转速度在孔口区域达到最大值，并向上扩散减小。值得注意的是，“滞留区”的颗粒并没有因位置相对固定而静止，相反地，此处的颗粒速度以极低的平移速度和一定的旋转速度向孔口滚动。取Z=0、50、100和400 mm为代表，统计筒仓内颗粒旋转速度的径向分布，如图6(b)所示。结果表明，在质量流区域，靠近边壁的颗粒旋转速度略小于中心附近；在漏斗流区域，靠近中心的颗粒旋转速度远大于边壁附近。但由于颗粒动力学的不均匀性，分布并不完全对称。取Y=0、50、100和150 mm为代表，统计筒仓内颗粒旋转速度的轴向分布，如图6(d)所示。结果表明，除了边壁区域的颗粒外，颗粒的旋转速度随高度位置变化呈指数分布，且在漏斗流区域急剧增大。在质量流区域，边壁附近的颗粒具有较大的旋转速度。

图6

图6   2~6 s内颗粒平均旋转速度在空间上的分布

Fig.6   Spatial distribution of average rotation velocity in 2—6 s

2.4 颗粒平移切向速度与旋转切向速度的空间相关性

图7为不同流型区域内颗粒间由平移运动产生的相对切向速度和由旋转运动产生的相对切向速度的空间相关性分布，采用泊松相关系数进行分析。相关系数ρ的计算公式如式（5）所示。

图7

图7   2~6 s内不同流型区域内颗粒间由平移运动和由旋转运动产生的相对切向速度的空间相关性

Fig.7   Spatial correlation of relative tangential velocities caused by translational motion and rotational motion among particles in different flow pattern regions within 2—6 s

ρ=∑vt-v¯tvr-v¯t∑vt-v¯t2∑vr-v¯t2(5)

式中，vt代表颗粒间由平移运动产生的相对切向速度，即颗粒的平移速度；vr代表颗粒间由旋转运动产生的相对切向速度，vr的计算公式如式（6）所示。

vr=ωR(6)

式中，ω是颗粒的旋转速度;R是颗粒的半径。

可以看到，相关系数随床层高度的增加而减小，总体上遵循指数函数分布。在质量流区域，颗粒间由平移运动和由旋转运动产生的相对切向速度在径向上的分布呈负相关；在漏斗流区域，颗粒间由平移运动和由旋转运动产生的相对切向速度在径向上的分布呈正相关。

2.5 颗粒滚动贡献率

图8为筒仓内颗粒滚动贡献率Rr的分布。由图8(a)可知，在筒仓底部以及靠近边壁的区域，Rr呈现出较大的数值，这表明此处的颗粒由于旋转造成相对切向速度比例较大，从而受到更严重的磨损。图8(b) 显示了 Rr在 Y 方向的分布。取Z=0、50、100和400 mm为代表，采用最小二乘法拟合原始数据。可以看出，对于不同的流型区域，Rr的分布表现出相似性，满足开口向上方向的抛物线分布。图8(d) 显示了Y=0、50、100和150 mm时，Rr 在 Z方向的分布。结果表明，颗粒的Rr值随高度的降低呈指数分布，在质量流区域分布较为均匀，在漏斗流区域急剧增大。在边壁区域，颗粒的Rr值远大于同高度的其他区域。因此在生产生活中，需要对筒仓的漏斗流区域及边壁区域加以保护，以减小由于颗粒相对切向运动引起的颗粒间及颗粒与边壁间的磨损。比如增大筒仓的质量流率、减小颗粒的半径等，这些优化措施将在后续的研究中进行验证。

图8

图8   2~6 s内颗粒平均滚动贡献率在空间上的分布

Fig.8   Spatial distribution of average particle rolling contribution rate in 2—6 s

3 结 论

本文利用基于Hertz-Mindlin和 RVD 滚动摩擦接触模型的DEM方法，模拟了三维筒仓中颗粒的卸料过程。为了掌握筒仓内的流型转化以及不同流型区域中颗粒运动的情况，划分了不同筒仓流型的区域，并定量研究了颗粒在不同区域中的平移和旋转特性。根据所得结果，得出以下结论。

（1）在卸料过程中，筒仓自上而下呈现出质量流向漏斗流过渡的混合流状态，质量流和漏斗流的转化高度可由MFI指数确定。在稳定的卸料状态下，筒仓的流型转换高度稳定在70 mm左右。

（2）在质量流区，颗粒间由平移运动和由旋转运动产生的相对切向速度在径向分布上呈负相关。 而在漏斗流区，颗粒间由平移运动和由旋转运动产生的相对切向速度在径向分布上呈正相关。

（3）稳态流动中，旋转造成相对切向速度比例较大的区域，主要集中在漏斗流区域与边壁区域。需要对该区域加以保护，以减小由于颗粒相对切向运动引起的颗粒间及颗粒与边壁间的磨损。

（4）在漏斗流区域与边壁区域的颗粒流型变化较大，是由于旋转造成相对切向速度较大引起的。

本文实验结果揭示了筒仓内颗粒平移运动和旋转运动的规律，为改进筒仓结构、提高筒仓的有效使用面积和使用寿命提供了参考数据。