刘曙光，钟文琪，陈曦

（能源热转换及其过程测控教育部重点实验室，东南大学能源与环境学院，江苏南京 210096）

引言

流化床锅炉由于燃料适应性广、污染物排放低等优点被广泛应用于工农业生产中[1-6]。在流化过程中，物料的混合特性[7-11]影响着颗粒间热质传递，进而决定了流化床整体设备的运行效率。布风板上方射流现象主导着物料间的混合行为，因此对布风板射流特征（如射流长度、直径和体积）开展研究至关重要。

流化床内部“射流”存在多种不同的定义，本文遵循Rowe 等[12]的阐述，即“射流”是在布风板上方形成的稳定空隙。以往流化床内射流现象的相关研究取得了重要进展[13-17]，例如Merry[18]采用视觉观察法对二维流化床中射流的长度特征开展了大量实验研究，并成功拟合出平均射流长度关联式，但是视觉观察的方式误差较大，获得的关联式准确性较差。Vaccaro 等[19-20]和Wen 等[21]分别采用压力探针方法和光学探针方法对三维流化床中的射流几何特征进行测量，虽然成功获得了射流几何特征的具体参数，但这两种方法都是侵入式测量技术，无法观察射流的精确形状特征并且会对流动过程产生影响。Rowe 等[12]和Cleaver 等[22]采 用X 射线透视成 像方法对三维流化床中的射流开展透视成像实验，得到射流的2D图像，该方法虽然避免了侵入式方法对流动过程的影响，但局限于单个射流。由于实验手段和测量方法的不足，传统对射流的相关研究大部分集中于单个水平或垂直射流，缺乏对流化床多孔布风板射流现象的认识，另外传统射流研究多集中于二维床或三维床近壁面区域[23-25]，难以获得三维床内部射流的结构特征，制约了对流化床内部流动机理的深层认识，从而阻碍了流化床装置的结构优化设计与高效稳定运行。

本文通过构建X 光层析成像（XCT）测量系统解决了传统手段难以获得流化床射流三维介尺度流动结构的难题，并采用该测量系统对流化床内布风板射流开展实验研究。首先基于锥束CT 三维重建算法（FDK）开发了CT 三维重建软件，实现了流化床内气固流动结构的三维重建。在此基础上设计了射流识别和量化算法，对布风板上方的气固射流进行统计和分析，获得了射流的形态结构及一系列量化特征，包括射流长度L、射流最大直径D和射流体积V。最后通过大量实验研究获得了平均射流速度Uj、床料粒径dp、布风板孔口直径d0和布风板孔口均分面积A0对射流形态结构和几何特征的影响规律。

1 实验装置与3D重建算法

1.1 X光层析成像实验装置

图1 为X 光层析成像实验装置示意图，该实验装置由流化床、X 光源、X 光探测器、空气压缩机、质量流量计和数据处理计算机组成。流化床为实验装置的主体部分，床室内径和高度分别为100 mm和900 mm，风室高度为100 mm，其中风室内填充了直径为10 mm的玻璃珠以满足均匀布风的要求。空气压缩机用来给流化床供气，质量流量计与数据处理计算机连接，通过计算机控制质量流量计实现空气流量调控。X 光源和X 光探测器是该实验装置的核心组件，X 光源发射的X 射线穿透流化床后投射在X 光探测器上，由X 光探测器拍摄并保存图像到数据处理计算机。流化床可以绕轴心360°旋转以完成X 光层析扫描成像，并由数据处理计算机实现流化床射流的三维重建。

图1 X光层析成像实验装置示意图Fig.1 Schematic diagram of X-ray tomography experimental setup

1.2 床料特性和布风板结构参数

实验采用三种不同粒径的石英砂作为床料，根据Geldart 颗粒分类法，三种床料均为B 类颗粒，石英砂的具体物性参数如表1所示。

表1 床料颗粒物性Table 1 Properties of bed material particles

布风板直径均为100 mm，其中孔口按照圆形阵列的方式排列。如图2 所示，A 类型布风板有37 个孔口，其中包含1 个中心孔口和三圈数量分别为6、12 和18 的孔口，相邻两圈孔口之间距离14 mm。B 类型布风板有61 个孔口，其中包含1 个中心孔口和四圈数量分别为6、12、18 和24 的孔口，相邻两圈孔口之间距离11 mm。A 类型布风板包含1、1.5 和2 mm 三种孔口直径，B 类型布风板孔口直径为1 mm，孔口均分面积为布风板总面积均匀分配给每个孔口周围的面积（A0=），布风板具体参数如表2所示。

图2 A类型布风板（a）；B类型布风板（b）Fig.2 A-type aeration plate(a);B-type aeration plate(b)

表2 布风板结构参数Table 2 Structural parameters of aeration plate

1.3 3D重建算法及几何位置误差分析

本文基于锥形束CT滤波反投影算法（FDK）[26-30]采用自主开发的软件实现流化床的3D 重建。不同于平行束和扇形束CT 重建算法（FBP）基于物体的一维投影，锥形束CT重建算法是根据物体的二维投影重建图像，这种重建方式一次投影的数据采集量更大，效率更高，而且更加符合实验所用锥形束X光的条件，最终获得的重建图像质量更高。

在进行X 光层析成像（XCT）实验时，X 光源焦点、流化床旋转中心和X 光探测器中心理论上应该保持在一条直线上，但是在调整三者位置关系时，由于人为校准误差的限制，很难保证三者精确地位于一条直线上，最终导致3D重建图像存在几何位置误差。3D 重建算法的几何位置误差分析主要是衡量目标物体重建图像的几何尺寸相对于真实尺寸的偏离程度。如图3 所示，目标物体采用3D 打印的立方体塑料框架，塑料框架边长为60 mm，每条边长由三个长度为20 mm 的小立方体的边长组成，总共包含64个交点。误差计算方法为：设塑料框架的一个角点为参考坐标点，其坐标为（0，0，0），其他交点根据与该角点空间相对位置确定相应坐标值。对该立方体框架进行CT三维重建，根据像素点统计重建图像中立方体框架64 个交点的坐标值，设第i个交点的真实坐标为（xi，yi，z）i，重建坐标为，则误差计算公式为：

图3 立方体塑料框架Fig.3 Cube plastic frame

经计算，该立方体框架的空间几何误差值为e=1.04%，满足实验要求。

2 3D图像处理方法

2.1 空隙率和CT图像

在流化床流化过程中，由于气泡穿过床体，气泡的上升、聚并和破裂使床体不同空间位置处的空隙率处于动态变化中，而布风板射流处的空隙率则始终比较稳定，比其他位置空隙率高，因此在流化床3D 重建图像中，布风板上方射流较为容易识别。采用已开发的重建软件重建出的原始3D 图像的体素值代表CT强度，流化床内某点对X射线的衰减度越强，重建出的该点CT强度值越小，本文基于空床、静态堆积床和流化床重建图像的CT 强度值之间的理论关系将流化床重建图像的CT 强度值转换为局部时均空隙率值，根据局部时均空隙率大小识别射流。

局部时均空隙率可通过流化床流化状态下X射线衰减系数（μ）f、床料颗粒X射线衰减系数（μp）和气体X 射线衰减系数（μg）来确定，由于X 射线衰减系数与CT强度值成正比，因此可以根据流化床流化状态下的CT 强度值（I）f、床料颗粒CT 强度值（Ip）和气体CT强度值（Ig）来计算局部空隙率。

由于单个床料颗粒的粒径较小，且是非均匀的，因此无法确定单个床料颗粒的CT 强度，但是，根据静态堆积床空隙率计算公式，床料颗粒的CT强度值可以由静态堆积床的空隙率（εb）、静态堆积床的CT 强度值（Ib）和气体的CT 强度值（Ig）来表示。

最终，流化床在流化状态下的局部时均空隙率（εg）可由流化床流化状态下的CT 强度值（I）f、静态堆积床的CT 强度值（Ib）、气体CT 强度值（Ig）和静态堆积床的空隙率（εb）计算获得。

在进行X 光层析成像（XCT）实验时，设定流化床每次旋转1°，两次旋转间隔时间为1 s，旋转一周大约需要6 min，并产生360 个投影方向，每个投影方向X 光探测器连续采集5 张X 射线衰减图像，并取其平均值作为该投影方向的最终投影图，因此3D 重建图像所得到的流化床内部形态结构和流动特征均为时均值。重建区域为流化床布风板至一个床高（100 mm）的区域，该区域具体尺寸为100 mm×100 mm×100 mm，该区域3D 重建体素（空间三维像素，对应平面二维像素）数量为200×200×200，每个体素边长为0.5 mm，因此最终重建出的流化床区域内部空间分辨率为0.5 mm。在最终的图像数据分析方面，选择射流3D 重建图像和流化床3D 重建图像的x切片进行研究，在3D 重建图像中，可以清晰地看到射流在整个布风板上的分布情况和三维形态特征，3D 重建图像x切片穿过流化床中心线，可以展示射流垂直方向截面的二维形态特征。

2.2 射流识别和量化

在流化床流化过程中，由于布风板上方射流的空隙率较高，所以本文采用目标体素空隙率与该体素邻域内空隙率平均值比较法来识别射流体素。如图4 所示，以布风板中心射流体素为例，对于3D空隙率图像的每一个横截面，以目标体素为圆心，半径为10 mm 的圆形区域作为该体素的邻域，计算该邻域内空隙率的平均值，邻域半径取10 mm 是为了最小化周围射流体素的影响而选取的最大邻域，这种邻域选取方法更能体现出射流体素值与邻域体素平均值在数值上的差异，然后比较目标体素值与邻域体素平均值，如果目标体素值高于邻域体素平均值并且高出10%以上，该目标体素即定义为射流体素，以白色体素表示，否则定义为非射流体素，以黑色体素表示。遍历3D 空隙率图像的每一个体素并对其进行射流识别，获得3D 射流二值图像。在二值图像中，可以清楚地看到每个射流，但由于实验过程中X 光源产生的量子噪声以及图像采集过程中产生的热电子噪声等一系列不可控因素，二值图像中除了射流体素外，还存在许多其他噪点，本文采用图像形态学膨胀、腐蚀等操作，在不影响射流体素的情况下去除其他小区域非射流体素，图5（a）为流化床3D 重建图像x切片原始灰度图，图5（b）为经射流识别算法处理获得的射流二值图像，从二值图中可以清晰地看到射流的形态结构。

图4 射流横截面Fig.4 Jet cross section

图5 射流灰度图（a）；射流二值图（b）Fig.5 Jet gray image(a);Jet binary image(b)

射流量化之前首先要定义射流的形状，图6 为通过X 光层析成像（XCT）实验重建的3D 射流，根据该3D射流的形状特征，定义射流近似为圆锥形。本文针对射流长度L、最大直径D、体积V三个特征进行量化。射流长度定义为布风板至单个射流所能达到的最大高度之间的距离，射流最大直径定义为单个射流最大横截面面积的等效直径，射流体积定义为单个射流所有体素体积之和，最终对布风板上所有射流的量化特征取平均得到某流化状态下的平均射流长度、平均射流最大直径、平均射流体积。

图6 射流3D重建图像Fig.6 Jet 3D reconstruction image

3 结果与讨论

本文采用已构建的X 光层析成像（XCT）测量方法对流化床布风板射流特征开展实验研究，分析不同流化风速下床料粒径dp、布风板孔口直径d0和布风板孔口均分面积A0对射流长度L、最大直径D、体积V和膨胀角θ的影响规律。对于粒径0.20～0.35 mm 床料以及0.35～0.45 mm 床料，实验采用的流化特征数Uj/Umf=1、1.5、2、2.5、3、3.5、4、4.5、5、5.5、6、6.5、7，对于粒径为0.45～0.60 mm 的床料，实验采用的流化特征数Uj/Umf=1、1.5、2、2.5、3、3.5、4、4.5、5。床层初始高度与床径相同，为100 mm。所有实验流化风速均采用降速法，从最大流化特征数对应的工况点开始，待该流化状态稳定后，对其进行X 光层析成像（XCT）实验，然后流化风速降低到下一个工况点进行实验，直至流化风速降到实验设定的最小流化风速。

3.1 床料粒径对射流特征的影响

图7 和图8 分别为布风板类型A-1，d0=1 mm 时三种床料粒径的射流x切片灰度图/二值图和射流3D 图，图中横轴为三种床料粒径（0.20～0.35 mm、0.35～0.45 mm、0.45～0.60 mm），纵轴为三种流化特征数U/Um（f1、3、5）。从图中可以明显看到均匀分布于布风板上方的射流现象，观察图7 和图8 的每一列可以发现，随着流化特征数U/Umf的增加，射流长度L、最大直径D和体积V都不断提高。对比图7 和图8 的每一行射流形态可以发现，在相同的U/Umf下，不同床料粒径产生的射流有明显的差异，床料粒径越小，射流形态越细长，具体表现为D和V越小，这是因为床料粒径越小，最小流化速度越小，在相同U/Umf下平均射流速度越小。当U/Umf较小时，例如当U/Umf=1 时，L分布极不均匀，大部分L很短，只有少部分射流可以渗透到流化床内部较深处，这是因为流化风量较小时，布风板孔口风量分配不均，大部分风量集中在少部分孔口，其余孔口的风量较小。观察图7 中的射流灰度图可以发现，部分邻近射流在上升至一定高度时相互靠近，最终合并形成气泡上升，这表明布风板上方邻近射流具有相互靠近并融合的倾向，同时也解释了流化床内初始气泡的产生。

图7 不同床料粒径下射流x切片灰度图/二值图Fig.7 Jet x-slice gray image/binary image with different bed material size

图8 不同床料粒径下射流3D图Fig.8 Jet 3D image with different bed material size

图9 横轴为平均射流速度U（j布风板孔口风速），纵轴为平均射流长度L。从图中可以看出，在相同Uj下，相对于较大的床料粒径，较小的床料粒径所形成的L更长。随着Uj的增加，三种床料粒径的L先呈上升趋势，并且在Uj较小时L增加较快，当Uj增加到一定程度时L趋于稳定，这是因为射流在达到一定高度后，射流顶部开始形成气泡脱离射流，L不再增加。对于0.20～0.35 mm 的床料粒径，L并不是一直处于上升趋势，在流化特征数U/Umf=2 时，L增加到最大值，随后L减小，然后上升，最后趋于稳定状态。图10 表示平均射流最大直径D随平均射流速度Uj的增加而增大，并且与相同U/Umf下床料粒径越小，D越小的规律不同，在相同Uj下，不同床料粒径所对应的D差异较小，床料粒径减小时，D轻微增加。随着Uj的增加，射流体积V逐渐增加，在相同Uj下，床料粒径越小，V越大，在床料粒径为0.20～0.35 mm 时，由于存在L的局部最大值，同样存在V的局部最大值。图11显示了平均射流体积V随平均射流速度Uj的变化趋势。

图9 平均射流长度L随平均射流速度Uj的变化趋势Fig.9 Variation of average jet length L with jet velocity Uj

图10 平均射流最大直径D随平均射流速度Uj的变化趋势Fig.10 Variation of average jet diameter D with jet velocity Uj

图11 平均射流体积V随平均射流速度Uj的变化趋势Fig.11 Variation of average jet volume V with jet velocity Uj

3.2 布风板孔口直径对射流特征的影响

图12 表示的是布风板类型A-1，0.35 mm＜dp＜0.45 mm 时不同孔口直径下射流3D 图，观察每一行图像中射流的形态可以发现，在相同的流化特征数U/Umf下，孔口直径越大，射流长度L、最大直径D和体积V就越小，这是因为在相同的流化风量下，孔口直径越大，平均射流速度Uj和通过布风板孔口的气体动量越小，形成的L、D和V就越小。当U/Umf=1 时，相比于流化床壁面附近的射流，布风板中心区域的射流现象更加明显，而当U/Umf=3、5 时，可以观察到布风板上靠近壁面的两圈射流现象更加明显，这表明当流化风速较小时，大部分流化风趋向于从布风板中心区域的孔口通入流化床，随着流化风速的提高，布风板中心孔口的风速增大，空气流动阻力增大，导致部分风量向布风板外圈孔口区域扩散，并最终在风室壁面的限制下，从风室壁面附近的孔口通入流化床，导致布风板外圈孔口风量增大，形成更为明显的射流现象。

图12 不同孔口直径下射流3D图Fig.12 Jet 3D image with different orifice diameter

图13～图15表示的是不同孔口直径下平均射流长度L、最大直径D和体积V随平均射流速度Uj的变化趋势，从图中可以发现，L、D和V随平均射流速度Uj的增加不断提高。孔口直径越大，L、D和V就越大，这是因为在相同的Uj下，孔口直径越大，流化风量就越大，因而产生的射流几何特征就越大。

图13 平均射流长度L随平均射流速度Uj的变化趋势Fig.13 Variation of average jet length L with jet velocity Uj

图14 平均射流最大直径D随平均射流速度Uj的变化趋势Fig.14 Variation of average jet diameter D with jet velocity Uj

图15 平均射流体积V随平均射流速度Uj的变化趋势Fig.15 Variation of average jet volume V with jet velocity Uj

3.3 布风板孔口均分面积对射流特征的影响

图16 表示的是d0=1 mm，0.35 mm＜dp＜0.45 mm时不同孔口均分面积下射流3D图，Ⅰ列图像对应的是A-1 型布风板，孔口数量为37，孔口均分面积为212 mm2，Ⅱ列图像对应的是B-1 型布风板，孔口数量为61，孔口均分面积为129 mm2。对比Ⅰ、Ⅱ两列图像可以发现，在相同的流化特征数U/Umf下，当孔口均分面积降低时，射流几何特征（L、D、V）明显减小，这是因为流化风量相同，孔口均分面积较小的布风板孔口数量较多，因此每个孔口的Uj更小。观察Ⅱ列图像可以发现，对于B-1 型布风板，当U/Umf=1 时，布风板最外圈孔口并没有产生射流现象，当U/Umf=3、5时，布风板最外圈孔口虽然有射流现象产生，但射流现象并不明显，相比于其他射流，布风板最外圈射流的几何特征（L、D、V）较小，这表明B 型布风板最外圈射流区域可能会产生死区，导致颗粒间不能充分混合，在实际工业流化床中会影响颗粒间的热量传递。

图16 不同孔口均分面积下射流3D图Fig.16 Jet 3D image with different orifice average area

图17～图19 表示的是不同孔口均分面积下L、D和V随Uj的变化趋势，从图中可以发现，L、D和V随Uj的增加不断提高，在同一Uj下，孔口均分面积越大，L、D和V就越大。观察图17、图19可以发现，当Uj较小时，不同孔口均分面积下的L和V差距较小，当Uj较大时，L和V差距较大，这表明在孔口均分面积较大的情况下，当Uj增加时，L和V增加的速度较快。

图17 平均射流长度L随平均射流速度Uj的变化趋势Fig.17 Variation of average jet length L with jet velocity Uj

图18 平均射流最大直径D随平均射流速度Uj的变化趋势Fig.18 Variation of average jet diameter D with jet velocity Uj

图19 平均射流体积V随平均射流速度Uj的变化趋势Fig.19 Variation of average jet volume V with jet velocity Uj

3.4 射流长度关联式

根据已有的实验数据，采用最小二乘法拟合流化床布风板射流长度L关联式。根据实验结果可以发现，射流长度L与布风板孔口直径d0，布风板孔口均分面积A0，平均射流速度Uj和床料粒径dp密切相关，所以关联式以d0，A0，Uj，Umf为自变量，其中Umf为床料最小流化速度，用来代替床料颗粒直径dp、密度ρp、流化风密度、黏度等参数，并以无量纲形式给出具体函数形式：L=F（d0，A0，Uj，Um）f。基于13 种流化特征数（U/Umf=1、1.5、2、2.5、3、3.5、4、4.5、5、5.5、6、6.5、7），两种床料粒径（0.35 mm＜dp＜0.45 mm、0.45 mm＜dp＜0.60 mm），三种布风板孔口直径（d0=1、1.5、2 mm），两种布风板孔口均分面积（A0=212、129 mm2）开展了55 组不同工况下流化床X 光层析成像实验，获得不同工况下布风板上方射流的具体长度值，最终拟合出流化床布风板射流长度L关联式：

式中，g为重力加速度，m/s2。

定义变量的指数因子为该变量的关联系数，分析射流长度L关联式可以发现，床料最小流化速度Umf的关联系数为-0.08，这表明射流长度L与床料最小流化速度Umf成反比，而最小流化速度Umf的主要影响因素是床料粒径，即床料粒径越小，射流长度L越大，与实验结果相一致。在射流长度L关联式中，平均射流速度Uj、孔口直径d0和孔口均分面积A0的关联系数分别为0.32、0.48 和0.20，这表明射流长度L与平均射流速度Uj、孔口直径d0和孔口均分面积A0成正比，符合实验结果。图20～图22 为射流长度L散点图和拟合关联式曲线图，从图中可以明显观察到，射流长度L曲线较好地拟合了实验数据。图23为在0.35 mm＜dp＜0.45 mm，d0=1 mm，A0=212 mm2时，实验数据、拟合关联式、Blake 等[25]和Rees 等[13]给出的关联式对比图，从图中可以明显看到，Blake 等和Rees 等两种关联式预测的射流长度值小于本文中实验测量的数值，不能很好地预测射流长度值。本文借鉴Rees 等关联式的形式，对变量的参数进行修正，使得关联式能较好地拟合实验数据，进一步提升了射流长度关联式的准确性。

图20 不同床料粒径下射流长度散点图和拟合曲线Fig.20 Jet length scatter plot and fitting curve with different bed material size

图21 不同孔口直径下射流长度散点图和拟合曲线Fig.21 Jet length scatter plot and fitting curve with different orifice diameter

图22 不同孔口均分面积下射流长度散点图和拟合曲线Fig.22 Jet length scatter plot and fitting curve with different orifice average area

图23 不同射流长度关联式对比Fig.23 Comparison of different jet length correlations

4 结论

采用X 光层析成像测量方法（XCT）对不同平均射流速度Uj、床料粒径dp、孔口直径d0和孔口均分面积A0下流化床中布风板射流的形态结构和几何特征开展了实验研究，得出以下结论：

随着平均射流速度Uj的增加，平均射流长度L、最大直径D和体积V呈上升趋势。当平均射流速度Uj较小时，不同射流间的形态结构和几何特征差异较大，少部分射流喷射速度较大，动量较大，可以渗透至流化床较深处，大部分射流由于动量较小，只能喷射至布风板上方较短距离，随着平均射流速度Uj的提高，不同射流间形态结构和几何特征趋于一致。布风板上方邻近射流具有相互靠近并合并形成气泡上升的倾向。

在同一平均射流速度Uj下，平均射流长度L、最大直径D和体积V与床料粒径dp成反比，与孔口直径d0和孔口均分面积A0成正比。通过实验数据成功拟合了无量纲射流长度关联式，能较好地预测不同平均射流速度Uj、床料粒径dp、孔口直径d0和孔口均分面积A0下的射流长度L。