在燃煤电站运行过程中，煤粉颗粒的粒径和形状对电厂锅炉的运行和燃烧有重要影响[1]。不恰当的煤粉粒径和形状将直接增大不完全燃烧热损失或磨煤机能耗，降低锅炉燃烧效率，情况严重时还可能引起煤粉管路堵塞而影响电厂运行，因此，为了及时调整锅炉运行中的相关参数， 必须实时检测煤粉粒度。
在实际检测过程中， 煤粉颗粒粒径的测量精度受颗粒团聚等因素影响， 沉降法、 吸附法等传统的形态分析和表征方法难以表征煤粉颗粒的真实形状， 煤粉颗粒的形状通常按球形或双球形进行简化处理[2]。
对于煤粉颗粒形状的识别，最常用的方法是借助于图像识别对单个煤粉颗粒的形态特征进行分析。Xu[3]提出一种基于二维主成分分析算法的煤粉图像识别方法，有效地减少了计算的复杂程度，提高了煤粉颗粒识别精度。Pan[4]提出煤灰图像自动识别算法，采用阈值迭代和最大类间方差法(Otsu算法)相结合的单阈值法进行图像分割，进一步提取了煤灰图像的几何特征。
煤粉颗粒进入在线检测系统后，大部分在主流方向流动，少部分会在微弱的湍流卷吸作用下偏离主流方向而向其他方向扩散形成颗粒沉积。早在1992年，Tsuji等[5]就利用CFD-DEM法模拟鼓泡流化床的流动行为，得到了与实际相符的流化速度等床层信息。Jin等[6]针对煤粉浓相气力输送立管内填充床流、 塞流、 搅拌流和低浓相流，研究了4种流型的控制机理和多尺度统计规律，为密相气力输送流型转换的研究提供了理论支持。Lu等[7]、 Jin等[8]通过实验研究了水分对于煤粉流动性的影响，并提出了可预测模型来优化相关单元操作。张春燕等[9]对颗粒在不同弯径比的弯管输送中的运动特性进行耦合模拟，发现弯管弯径比增大时颗粒受到的离心力作用时间增多，颗粒不容易发生分散。
本文中将以在线检测设备中的煤粉颗粒作为研究对象，借助图像识别技术分析颗粒形态，采用EDEM软件建立煤粉颗粒模型，采用Fluent软件对颗粒相的运动进行数值模拟，通过对比不同类型煤粉颗粒的流动特性，为煤粉在线检测设备的长期运行提供数据基础。
1 煤粉颗粒特征的提取
使用EV-E1600C型号相机对煤粉颗粒进行静态拍摄，获得煤粉颗粒的照片，然后对照片进行灰度和二值化处理[10]，得到可分析煤粉颗粒形状的图像。煤粉颗粒的二值化图像如图1所示。由图可见，单张成像中存在大量颗粒，部分颗粒由于颗粒堆积使得成像颜色较深，不易进行颗粒形状分析，因此，在提取颗粒进行颗粒建模时，选择图像中相对空旷区域内且不受周围其他颗粒影响的、 与周围颗粒无黏连的单个颗粒进行分析。
图1 煤粉颗粒的二值化图像
Fig.1 Binary image of pulverized coal particles
为了对颗粒粒度进行处理，对图像中的不同颗粒进行标号和测量，统计每个颗粒所占像素点的总和，计算出颗粒面积(可视为颗粒最大投影面积)；对颗粒的边缘像素点作出标记，统计所标记的像素个数，计算出颗粒周长；对每个颗粒进行横向与纵向像素点统计，记录每个颗粒的统计数据，将统计数据的平均值作为颗粒的平均粒径；沿着颗粒边界像素点建立颗粒的最小外接多边形，以多边形的最小外接圆直径作为颗粒外接圆直径。
对于颗粒形状，采用长宽比、圆形度以及球形度来进行表示。长宽比a的计算公式[11]为
a=L/W，
(1)
式中： L为颗粒的最大长度， μm； W为颗粒的最小宽度， μm。
圆形度e的计算公式[12]为
(2)
式中： S为颗粒面积， μm2； C为颗粒周长， μm。
球形度φ的计算公式[13]为
φ=d/D，
(3)
式中： d为每个颗粒的最大投影面积对应的直径， μm； D为颗粒外接圆的直径， μm。
对所选二值化图像进行特征提取，由于煤粉颗粒数据量巨大，因此选取10个煤粉颗粒进行统计。10个煤粉颗粒的形状参数如表1所示。由表可知，随着颗粒直径的增大，颗粒的周长和面积也随之增加，但长宽比、圆形度和球形度的变化无规律，其中长宽比波动较大，球形度变化较小。
表1 10个煤粉颗粒的形状参数
Tab.1 Shape parameters of 10 pulverized coal particles
颗粒序号颗粒直径/μm颗粒周长/μm颗粒面积/μm2长宽比圆形度球形度115.29177.324423.5171.6030.9470.603220.204138.852786.5342.1740.7200.559335.812253.1732 420.5242.2490.6890.539441.323275.3653 297.2522.1180.7400.571551.615341.4584 840.0282.1030.7230.527660.416363.7316 503.5071.9130.7880.551772.931556.1869 831.2542.4270.6320.483887.231720.75214 429.276 2.6270.5910.4579101.367 534.24817 968.516 1.7620.8470.57910 111.518 710.38722 506.468 2.0270.7490.534
2 模型创建
2.1 颗粒的二维模型
模拟颗粒的流动性与颗粒间的碰撞时，需要保证颗粒模型的接触面为球面，因此，借助于重叠离散元簇理念，利用圆球来堆砌颗粒模型。
颗粒二维模型构建过程如图2所示。由图可知，首先，由某颗粒的二值化图像放大后得到颗粒轮廓；其次，根据粒子形状和边缘不均匀程度构建内部骨架线，骨架线中心及端点作为内部圆圈定位点，即圆心位置；最后，在每个定位点选择与表面形状内切的最大圆圈进行布置，使重叠后的圆圈外表达到与颗粒表面形状较高的吻合度。
图2 颗粒二维模型构建过程
Fig.2 Construction process of two-dimensional particle model
2.2 颗粒的三维几何模型
在二维图像中完成内部圆圈绘制后，采用EDEM中独特的用户自定义球面填充技术创建物料颗粒。用若干球面的组合来表示颗粒表面，既能体现颗粒的非球形特征，又可以使颗粒的接触满足球面接触的物理模型[14]。
将圆圈直径视为三维颗粒内部内切圆球的直径，进行三维颗粒模型创建。选取几种典型的颗粒形状进行堆砌，获得单球形、双球形以及2种异形颗粒，分别记为P1、 P2、 P3、 P4。EDEM软件中堆砌出的颗粒模型如图3所示。
图3 EDEM软件中堆砌出的颗粒模型
Fig.3 Particle model stacked by EDEM software
为了验证EDEM软件中堆砌出的颗粒模型的精度，将颗粒模型与采集的实物颗粒的球形度进行对比，颗粒模型的精度验证如表2所示。
表2 颗粒模型的精度验证
Tab.2 Accuracy verification of particle model
颗粒球形度实物颗粒颗粒模型误差绝对值/%颗粒球形度实物颗粒颗粒模型误差绝对值/%P10.831.0017.50P30.640.7918.90P20.610.7013.43P40.440.5519.45
由表2可知，颗粒模型与实物颗粒的球形度误差绝对值范围为10%～20%，因此，基于计算能力的考虑，可在模拟过程中将颗粒视为单球形或双球形颗粒。由于计算机算力不足，因此本文中未能对超细煤粉以及其他形状的煤粉颗粒进行大量模拟研究。
2.3 气固两相的数学模型
在CFD-DEM方法中，采用牛顿第二定律描述颗粒相的运动，采用Navier-Stokes方程描述气相的运动[14]。考虑到煤粉颗粒间碰撞过程中发生的形变和滑移，利用软球模型[15]处理颗粒间的相互作用。在计算相间耦合时，采用Beestra关联式[16]来计算气固相间曳力。
假定气力输送管道内的空气属于不可压缩的非稳态牛顿流体，根据流体力学[17]可知，质量守恒方程为
+
·(εg ρgug)=0，
(4)
动量守恒方程为
(5)
式中： εg为气相体积分数； ρg为气相密度， kg/m3； t为时间， s； ug为气相速度， m/s； 为向量微分算子； p为压强， Pa； Fsg表示气固相间作用力， N； τg为气相应力张量， Pa； g为重力加速度，m/s2。
对于离散相，颗粒运动受力方程[18]为
(6)
式中： ma为颗粒质量， kg； ra为颗粒半径， m； Va为颗粒体积， m3； β 为相间曳力系数； vs为单个颗粒的速度， m/s； Fca为颗粒碰撞时受到的接触力， N。
颗粒-颗粒间的碰撞受力为
Fabn=-knδnnab-ηnvabn，
(7)
δn=(Ra+Rb)-|ra-rb|，
(8)
式中： Fabn为颗粒碰撞时的法向接触力， N； kn为法向刚性系数， N/μm； δn为颗粒碰撞时发生的形变量， μm； nab为法向上的单位向量； ηn为法向阻尼系数，(N·s)/μm； vabn为法向相对速度， μm/s； Ra、 Rb分别表示颗粒a、 b的粒径， μm； ra、 rb分别表示颗粒a、 b的半径，μm。
2.4 检测装置的三维模型
针对煤粉在线检测装置建立简化模型，检测装置的结构参数和网格模型如图4所示。由图4 (a)可见，煤粉颗粒在主流通道中由上部入口进入，下部出口流出，检测元件置于检测通道的左右两侧，观察煤粉在检测通道中的扩散情况。由图4 (b)可见，采用Workbench软件对模型进行网格划分获得的模型结构较为简单；由于所需计算精度较高，网格的数量及质量直接影响计算的精度和收敛速度，因此采用正四面体来划分网格，并在进、出口边界及贴近壁面处设置网格局部加密。
(a)结构参数(b)网格模型图4 检测装置的结构参数和网格模型Fig.4 Structure parameters and grid model of detection device
2.5 模拟参数及边界条件
为了对气体速度分布、颗粒运动轨迹、颗粒速度等数据进行分析，研究不同形状颗粒在检测设备中的不同流动特性，气相运动采用ANSYS 19.0软件中的Fluent进行模拟，颗粒相运动采用离散元软件EDEM进行运动轨迹计算，二者通过Fluent的用户自定义函数UDF实现耦合。
在气相运动中，气体密度为1.225 kg/m3，动力黏度为1.789 4×10-5 Pa·s，假定气体在室温条件不可压缩；湍流模型选择RNG k-ε模型；入口边界选择速度进口条件，假设进气口速度均匀，另外考虑到相机曝光时间和移动距离，入口速度设定为0.02 m/s；出口边界选择压力出口条件。
采用随机轨道模型[19]作为颗粒相的湍流扩散模型。Hertz-Mindlin模型[20]适用于常规颗粒的接触作用，限制条件更少。该模型将颗粒间的接触按照静态弹性接触进行处理，从而得到颗粒间球形接触面积与弹性变形的关系，因此可以作为颗粒的基本接触模型。
由于煤粉颗粒质量分数较大，在测量前需要通过旋风分离器以及分级采样器对煤粉进行分级稀释，以便获取质量分数适中的煤粉进行测量。煤粉颗粒的密度为1 440 kg/m3，剪切模量为1×109 Pa，泊松比为0.3，依据煤粉图像中提取的数据将煤粉颗粒速率设置为5 000 s-1。检测装置的密度为7 800 kg/m3，剪切模量为7×1010 Pa，泊松比为0.3。
4种颗粒相的几何参数和时间步长如表3所示，碰撞系数如表4所示。
表3 颗粒相的几何参数和时间步长
Tab.3 Geometric parameters and time step of particle phase
颗粒序号形状粒径/μm时间步长/s颗粒序号形状粒径/μm时间步长/sS1单球形505×10-8S3双球形1005×10-8S2单球形1006×10-8S4异形1003×10-8
表4 碰撞系数
Tab.4 Collision coefficient
颗粒-颗粒恢复系数静摩擦系数滚动摩擦系数颗粒-壁面恢复系数静摩擦系数滚动摩擦系数0.5000.5450.0100.5000.5000.010
3 结果与讨论
3.1 速度场
为保证检测设备内模拟计算达到更好的效果，在加入颗粒前单独进行气相的流场模拟。流场充分流动后在Y=0处截面上的速度场如图5所示。
(a)速度矢量图(b)速度流线图图5 在Y=0处截面上的速度场Fig.5 Velocity field on cross section at Y=0
由图5 可以看出， 气相在煤粉颗粒的主流通道内流动比较均匀， 流速在靠近管道中心位置最大， 沿管径增大方向逐渐减小， 能够保证大部分煤粉颗粒在气力输送作用下经过检测区， 但在检测通道内有比较明显的流体旋涡生成。 主流通道中气相的最大流速为0.035 m/s， 通道交汇处平均流速约为0.025 m/s， 检测通道内平均流速为0.002 m/s，二者比值为12.5∶1。尽管旋涡速度与主流速度相比较小，但对于微米级的煤粉颗粒，在这种流体旋涡的作用下会有部分颗粒进入检测通道内，在旋涡的作用下扩散至边缘检测元件处，进而影响清晰成像。
3.2 颗粒相的流动性
在模拟出充分流动的流场后加入煤粉颗粒，发现形状不规则的颗粒S3、 S4在模型中的流动特性相似，因此选择颗粒S1、 S2、 S3进行对比。分别观察3种颗粒达到稳定状态后在水平检测通道内的扩散情况，取30 s为一个模拟时间段，分别在距检测元件10、 20 mm内的区域为计量区域1、 2，检测区域位置示意图如图6所示，由图可见，左、 右两侧的红框内共为计量区域1，蓝框内(包含红框区域)共为计量区域2。
图6 检测区域位置示意图
Fig.6 Location diagram of observation area
为显示清晰，将3种颗粒放大10倍后进行检测，检测通道内不同类型煤粉颗粒在不同时刻的位置分布如图7所示。由图7(a)、(b)可见，S1进入检测通道后少部分与两侧检测元件接触，大部分聚集于检测通道两端的底部形成颗粒堆积；S2沿检测通道的轴向方向均匀性分布，在检测通道下部形成沉积。这是由于，颗粒在检测通道内的运动轨迹由颗粒在壁面的湍流卷吸效应和颗粒的质量共同决定的，而且卷吸效应与距中心区域的距离呈负相关的关系；S1、 S2粒径不同，S2的质量是S1的8倍；S1经壁面反射时，入射角和反射角较大，并且不会接触上壁面，在反射2～3次后运动至检测通道两侧边缘位置；颗粒S2在反射时接近于垂直入射，多数颗粒在经第1次反射后会与上壁面接触，多次反射后形成在通道底部的均匀式堆积。由图7(c)可见，对于形状不规则的颗粒S3，受湍流卷吸的影响较小，进入检测通道两侧边缘位置的颗粒较少，运动轨迹与颗粒 S2类似，在检测通道底部沉积，不会形成在检测通道两侧边缘位置堆积的现象。
(a)颗粒S1(b)颗粒S2(c)颗粒S3图7 不同煤粉颗粒在不同时刻的位置分布Fig.7 Position distribution of different pulverized coal particles at different time
对扩散在计量区域1、 2内以及检测通道内的3种颗粒进行统计。计量区域内的颗粒数如表5所示。由表可知，在计量区域1内，颗粒S1扩散的个数最多为250，S2的个数为135，S3的个数为33；在计量区域2内，S1、 S2的扩散数量相近，分别为280、 316，S3仅有108；在水平检测通道内，S2的扩散个数最多为468，S1的个数为308，仍是S3的个数最少，仅有178，说明煤粉颗粒30 s内到达通道交汇处后，在气体二次流与湍流的共同作用下，相较于不规则的双球形颗粒S3，形状规则的单球形颗粒S1、 S2扩散到检测通道两侧检测元件附近的情况比较严重。
表5 计量区域内的颗粒数
Tab.5 Numbers of particles in measuring area
颗粒计量区域1计量区域2检测通道S1250280308S2135316468S333108178
沿煤粉主流通道方向(Z轴方向)选取不同截面，不同截面的颗粒速度分布如图8所示。由图可见，3种颗粒速度分布的总体趋势相同，在Z轴正向方向距离为0.11 m处，也就是通道交汇处，由于在检测通道内存在扩散的煤粉颗粒，该截面上的平均速度突然减小，3种颗粒速度最小均趋近于0；进入通道交汇处的S1、 S2、 S3的速度分别为0.58、 0.71、 0.83 m/s，在流动出口处，S1的速度最大为1.42 m/s，S2的速度为1.06 m/s ，S3的速度最小为0.91m/s，说明在相同的气力输运作用下，颗粒获得的速度与其质量相关，S1因质量最小而速度变化最为明显，在流出检测装置时的速度也最大，颗粒S3因绕流阻力较大而速度较小。
图8 主流通道上不同截面的颗粒速度分布
Fig.8 Distribution of particle velocity of different sections in mainstream channel
4 结论
首先，基于图像识别技术在煤粉二值化图像中进行煤粉颗粒基本特征的提取；其次，基于重叠离散元簇理念建立煤粉颗粒的二维模型，通过球面填充技术建立颗粒的三维几何模型；接着，在Workbench中建立检测装置的三维模型；然后，根据气固两相的数学模型设定模拟参数及边界条件；最后，采用CFD-DEM耦合方法模拟煤粉的气固两相流动，进行了气相速度场的分析以及颗粒相的流动分析，通过颗粒相不同时刻的位置分布研究了颗粒的扩散机理，获得了颗粒相不同位置处的速度分布和体积分数，归纳出颗粒的流动特性。
1)气相在主流通道内流动比较均匀，在靠近管道中心位置最大流速为0.035 m/s，沿管径增大方向逐渐减小；在检测通道内有比较明显的旋涡生成，湍流卷吸效应影响了煤粉颗粒影像的清晰度。
2)形状规则的球形颗粒S1、 S2更易偏离主流流动方向形成颗粒沉积；形状不规则的颗粒S3受气力输运作用影响较大，且不易分散，大部分在煤粉气流主流通道内流动；在计量区域1内，颗粒S1的扩散数量最多，S3的扩散数量最少，说明在颗粒受湍流卷吸效应偏离主流通道后，粒径较小的颗粒更易运动至检测通道两侧检测元件附近，长时间运行影响检测效果。
[1]葛凌峰. 燃煤电站煤粉粒径监测研究进展[J]. 能源与节能, 2021, 195(12): 189-190.
[2]赵霖, 高正阳. 基于OpenCV图像识别对一次风煤粉细度与形状的测量[J]. 能源与节能, 2018, 159(12): 156-157.
[3]XU Y.Research on coal and rock image recognition based on 2DPCA algorithm[J]. IOP Conference Series:Earth and Environmental Science, 2020, 514(2): 22004-22009.
[4]PAN J.Research on the image auto-recognition algorithm for the fusibility of coal ash[J]. Foreign Electronic Measurement Technology, 2008, 27(10): 1-3.
[5]TSUJI Y, TANAKA T, ISHIDA T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[J]. Powder Technology, 1992, 71(3): 239-250.
[6]JIN Y, LU H F, GUO X L, et al .Multiscale analysis of flow patterns in the dense-phase pneumatic conveying of pulverized coal[J]. AIChE Journal, 2019, 65(9): 16674.
[7]LU H F, GUO X L, JIN Y, et al. Effect of moisture on flowability of pulverized coal[J]. Chemical Engineering Research and Design, 2018, 133: 326-334.
[8]JIN Y, LU H F, GUO X L, et al. The effect of water addition on the surface energy, bulk and flow properties of lignite[J]. Fuel Processing Technology, 2018, 176: 91-100.
[9]张春燕, 马超, 晏飞. 基于CFD-DEM方法的不同弯径比弯管中气固两相流动特性[J]. 江苏科技大学学报(自然科学版), 2018, 32(6): 803-808.
[10]李颖. 基于深度学习的煤粉颗粒分割方法研究[D]. 徐州: 中国矿业大学, 2021.
[11]陈海洋, 汪稔, 李建国, 等. 钙质砂颗粒的形状分析[J]. 岩土力学, 2005, 26(9): 1389-1392.
[12]陈云波, 徐培涛. 水泥颗粒形貌的表征及其研究方法[J]. 水泥, 2003(2): 17-19.
[13]刘清秉, 项伟, BUDHU M, 等. 砂土颗粒形状量化及其对力学指标的影响分析[J]. 岩土力学, 2011, 32(增刊1): 190-197.
[14]ZHONG W Q, YU A B, LIU X J, et al. DEM/CFD-DEM modelling of non-spherical particulate systems:theoretical developments and applications[J]. Powder Technology, 2016, 302: 108-152.
[15]KOBAYASHI T, TANAKA T, SHIMADA N. DEM-CFD analysis of fluidization behavior of Geldart Group A particles using a dynamic adhesion force model[J]. Powder Technology, 2013, 248: 143-152.
[16]BEETSTRA R, HOEF M, KUIPERS J. Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres[J]. AIChE Journal, 2007, 53(2): 489-501.
[17]王松岭. 高等工程流体力学[M]. 北京: 中国电力出版社, 2011: 74-80.
[18]王婧, 王军武. 气固两相流中颗粒速度分布函数统计分析[J]. 中国粉体技术, 2018, 24(5): 1-5.
[19]郭佳才. 电旋风除尘器的气固两相特性研究与数值模拟[D]. 上海: 上海工程技术大学, 2020.
[20]王国强, 郝万军, 王继新. 离散单元法及其在EDEM上的实践[M]. 西安: 西北工业大学出版社, 2010: 17-18.