基于动态双重网格下喷动床滞止区流动特性CFD-DEM数值模拟

2021-11-26王洪远纪律孟繁旭李斌杨建蒙陈海生

化工学报 2021年11期

王洪远，纪律，孟繁旭，李斌，杨建蒙，陈海生

（1 华北电力大学能源动力与机械工程学院，河北保定071003； 2 中国科学院工程热物理研究所，北京100190）

引言

喷动床作为流化床的一种，其应用领域十分广泛，如农业颗粒物料（稻谷、油菜籽、小麦等）的干燥和去皮、去壳，农作物种子以及制药工业中的包衣和化学肥料的造粒，以及固态物料的混合等[1]。喷动床在工作时床层受到通过小孔或喷嘴射入气流的影响，中央射流内的固体颗粒被气流抬升到一定高度，此时颗粒的重力大于射流作用在颗粒上的力，颗粒开始下降，之后又被气流重新朝上带起而呈喷泉状，从而导致颗粒在喷动床内往复运动。喷动床的工作效果主要受床体内的气固两相流运动影响，而床体内的气固两相流内部的流动规律又是一个极为复杂的过程组合，关于内在规律的研究还在进行。常规的实验观察对流化床内气固两相流动特性的研究难以获得颗粒详细的运动信息，而计算机模拟技术的发展，有助于获得微观层次上颗粒的运动信息[2-7]。

离散单元法(discrete element method,DEM)是目前较为先进的研究方法，凭借能追踪每个颗粒的位置、受力等运动信息的优点，其被广泛应用到气固两相流研究中[8-9]。国内外学者分析流化床内气固两相流动特性时侧重于CFD-DEM 和LBM-DEM 方法，但由于LBM-DEM 方法较于CFD-DEM 方法只能对床层尺度较小、颗粒直径较小的工况进行模拟，所以对于尺寸较大的喷动床常采用CFD-DEM方法开展。王敬哲[10]建立了气流-喷动床的三维数值模型，对床内的气固两相流动特性进行了模拟，较为完整地揭示了喷动的形成过程及气固两相在床内的速度分布情况。彭丽等[11]利用计算流体力学和离散元(CFD-DEM)方法模拟了拟二维矩形鼓泡床内的气固流动特性，分析了固相弹性、恢复系数对流场间歇性的影响。Olaofe 等[12]则模拟了颗粒在矩形喷动流化床内的混合过程。Hilton等[13]在CFDDEM 模型下，选取改进的压力梯度力模型与非球型颗粒（椭圆形、多边形、多球形）对气固喷动床的动力特性进行研究。杨建蒙等[14]基于MFIX 软件，利用MFIX 开源程序中的MFIX-DEM 模型，比较了Gidaspow 和Syamlal-O’Brien 2 种曳力模型下对压力源项处理的2种方式Model A和Model B的三维喷动流化床气固两相流动特性，并与同工况下的实验结果进行了对比。

诸多研究者在对气固两相流领域的研究中各有建树，他们的研究推动了气固两相流技术的进步[15-19]，但在现有的研究中均未对喷动床内滞止区现象进行深入研究，因此本文首先建立了基于动态双重网格技术的CFD-DEM 方法，然后进行了喷动床径向混合实验与模拟研究，同时与单网格模型进行对比，验证了动态双网格方法结果的准确性，然后利用动态双重网格方法建立了不同进口气速和不同初始堆积高度下的计算模型，同时对滞止区内颗粒进行示踪标记，分析示踪区域的颗粒流动情况，得出颗粒在不同模型的滞止区流动混合的规律。

1 数学模型

1.1 气固两相的运动模型

气固两相运动过程在流化床内非常复杂，而CFD-DEM 数值模拟技术已成为解决众多领域颗粒材料问题的有力工具。

颗粒之间不断碰撞并受气体夹带作用运动，同时近壁面颗粒还与床体碰撞，另外重力、电场力和磁场力等外界物理场产生的作用力也会对颗粒造成影响。因此本研究颗粒相的碰撞模型选择基于时间驱动原理的软球模型，其运动方程见文献[20]。气相运动模型则采用气固两相耦合的Navier-Stocks方程，湍流模型采用k-ε模型，详细的控制方程参见文献[20]。

1.2 气固两相间耦合

1.2.1 曳力模型气固相间曳力表达式通常表示为曳力系数乘以速度差的形式：

根据Gidaspow 的曳力模型，气相对颗粒的曳力系数[21]如下：

1.2.2 气固两相间动量耦合气固两相间的作用力应该遵循牛顿第三定律，颗粒受到了流体对它的曳力，相应的流体也会受到颗粒对它的反作用力，这个力应该与曳力大小相等，方向相反。所以气固两相间的动量耦合公式为：

1.2.3 动态双重网格模型图1为三种网格划分示意图，为了便于分清气相网格和颗粒相网格，用红色的网格代表颗粒相粗网格，用黑色的网格代表气相细网格。由图可见，在单网格模型下气相和固相的计算都在同一网格内进行，即气相和颗粒相的计算共用一套网格[22-27]。在固定双重网格模型下气相和颗粒相分别有一套网格，气相的细网格被包含在颗粒相粗网格中，所以如颗粒的速度、温度、床体空隙率等在颗粒相的粗网格中计算；而气相的流场和温度场等在细网格内计算，在计算过程中通过气固两相间的耦合使两套网格之间即时映射，在同一时刻下通过网格间数据的传递来交换相间能量和动量等信息。

图1 网格模型划分示意Fig.1 Schematic diagram of grid model division

动态双重网格模型与前两种网格划分情况不同之处在于，动态双重网格方法主要在颗粒相的网格的设置有极大改变，同样气相采用与固定双重网格相同的细网格，用来在计算过程中映射气相的计算信息，同样根据计算对象的不同可以改变气相网格的疏密程度；而颗粒相网格划分则不以传统网格划分方法为依托，其以颗粒的球形中心为网格中心，以图中颗粒i为例，在某时刻下以i颗粒的中心作为i颗粒网格的中心，同样其网格尺寸可以根据i颗粒的尺寸来改变，其网格尺寸设置值的范围为颗粒直径的3～4 倍[28]，同时判断统计位于i颗粒网格内的其他颗粒个数，以及气相网格编号和个数，在该时刻其他颗粒的网格也可以通过此方法划分，所以称随颗粒运动的颗粒相网格为动态粗网格，而气相网格称为动态细网格。在计算过程中通过气固两相耦合进行数据映射和相间的能量交换[29]。此外如何判定颗粒是否在所计算网格之内，是根据该颗粒的质心是否在计算的网格之内。只要该颗粒的质心在所计算单元之内，则可以认为该颗粒整体都属于此计算单元。

1.2.4 网格间数据传递在固定细网格中计算的气体速度等物理量可通过求解气相方程来得到[30]。而在动态粗网格中计算的颗粒相的速度等物理量则是由其包含的所有细网格中对应物理量的加权平均得到。

动态粗网格中气体速度计算公式如下：

粗网格映射给其包含的所有细网格空隙率为：

2 模型验证

为了探究喷动床混合规律和验证动态双重网格的准确性，实验过程选取直径4 mm的红黄两种颜色的亚克力球形颗粒作为实验物料，在喷动床流化过程中来观察颗粒混合特性。

图2 为进口表观气速1.5 m/s、初始堆积高度0.16 m的流化床内颗粒实验和模拟的径向混合序列图，图3 为数值计算的颗粒径向混合过程垂直速度分布云图。从图中可知，在t=0 s 时，两种颜色颗粒沿床高方向彼此分离，随后进口气流带动底部中间区域颗粒向上运动，而靠近壁面两侧区域的颗粒则向下运动，在颗粒往复运动过程中气流冲击床体内部而形成气泡，随后气泡不断扩张带动床体颗粒向上膨胀，红色颗粒开始向上混合形成一个“树”状图形，而黄色颗粒则向下运动来填补缺失的空间。当颗粒在扩散进程中自身重力大于气流的曳力时，颗粒便开始下行，同时与其他颗粒不断碰撞混合，缓慢下移至床层底部后，又被卷吸至喷口附近重新进入下一流化循环。颗粒在此循环运动中和其周围颗粒不断进行着混合与扩散，其混合程度也不断加深，最后随着时间的推移，会达到颗粒间混合的动态平衡。

图2 颗粒径向混合过程实验与模拟对比Fig.2 Comparison between experiment and simulation of radial mixing process of particles

图3 颗粒径向混合过程垂直速度云图Fig.3 Vertical velocity nephogram of particle radial mixing process

此外观察0～2 s内实验与模拟图片可以发现，在床层底部壁面与左右侧壁面由于空间位置的限制，存在一定的滞止区域，该区域内颗粒由于远离喷口位置，气流不能及时扩散，造成流动缓慢，在流化混合过程与壁面形成“三角”区域。该区域在流化时间内的面积随时间推移不断缩小，同时观察利用动态双重网格计算的模拟区域也可发现此规律，因此通过径向混合实验和模拟对比也论证了动态双重网格的计算结果与实验混合流化过程结果高度一致。在此模型基础上，本次研究还利用数值模拟中颗粒追踪技术开展滞止区颗粒流动特性的研究。

3 数值模拟及结果分析

模拟对象为单孔射流准三维矩形喷动床，其长宽高分别为150 mm、4 mm、900 mm。气流入口设在床层底部中心，宽度为0.01 m，喷动床模型如图4 所示。固相采用了1904 个球形颗粒，颗粒直径4 mm、密度2700 kg/m3。气相选择空气，密度和黏度分别取值1.205 kg/m3和1.8×10-5kg/(m·s)。具体模拟计算参数见表1。将流体在壁面处设置为无滑移边界条件，颗粒在壁面处设置为滑移边界条件。

图4 喷动床几何模型Fig.4 Geometric model of spouted bed

表1 数值模拟计算参数Table 1 Calculation parameter in numerical simulation

3.1 动态双网格模型与单网格模型对比分析

3.1.1 滞止区混合特性对比为了体现动态双网格技术的计算准确性，在计算过程中将床体内颗粒分为四个区域，利用黑色区域作为滞止区参考，其距离左侧和右侧壁面为0.008 m，红色区域作为喷动区参考，亮青色区域用来代表床层上半部分区域，绿色区域用来代表与滞止区相邻的区域。

图5和图6给出了初始堆积高度为0.18 m，在喷口气流速度为27 m/s时的单网格模型和动态双网格模型在2 s 内床体颗粒的混合流动情况。从图中可以看出，气流由底部进入床体，导致床层出现松动，在初始时刻由于空间狭小，导致喷口附近的局部压力瞬时增大，由此产生一条裂缝，推动颗粒往上运动。而整个床体中心处的颗粒速度较大，并且由于颗粒之间的动能传递，导致床体上部的部分颗粒跃起，并随着时间的推移，气流往上继续流动。在0.2～0.3 s 时刻附近，床层内部产生一个较大气泡，随后喷动区底部的颗粒在颗粒之间产生气泡的带动下一起朝上运动；到达床层表面附近即喷泉区后颗粒向床体壁面移动分散；由于位于环隙区上部的两侧壁面附近的气体瞬时速度较小，导致颗粒受到的曳力减小，在重力作用下沿着床体壁面向下部运动，并在这个过程中参与了环隙区内颗粒的混合，颗粒受气流卷吸影响，下降到床层底部后会重新进入喷动区，从而实现了颗粒在流化床内循环。

图5 单网格模型2 s内混合序列图Fig.5 Mixed sequence diagram within 2 s of single grid model

图6 动态双网格模型2 s内混合序列图Fig.6 Mixed sequence diagram within 2 s of dynamic double grid model

在这个循环过程中，由于壁面附近的颗粒远离喷口，其气流速度较小，致使床底与壁面附近的颗粒被压在床的底部角落处，形成流动停滞区，也称“死区”。由图5 可知，利用单网格方法在初始时刻，滞止区内黑色颗粒受到床层的膨胀后回落冲击，其高度随着时间推移逐渐变低，t=0.2 s 时,滞止区高度变为原来的一半；t=0.4 s 时,滞止区内颗粒完全随着喷口气流进入喷动区；t=0.6 s 时，观察可以发现在床层底部黑色颗粒所占比例大大降低；在随后的时间内整个床层内颗粒分布均匀。观察图6，则发现在0～2 s 内床层底部的滞止区示踪区一直存在，黑色颗粒在滞止区的含量随着时间的推移逐渐减少，t=0.5 s 时其滞止区黑色颗粒的高度变为原来的一半，t=2 s 时滞止区内仍然存在黑色颗粒，此时刻的滞止区内被其他颗粒所占据，重复此循环。

3.1.2 滞止区示踪颗粒流动高度对比为了对滞止区流动特性有更为清晰的描述，在研究中对黑色标记的颗粒进行统计追踪，取每个时刻的黑色示踪颗粒的平均高度作为参考，得到0～2 s内示踪颗粒的流动高度如图7所示。

图7 两种网格方法下滞止区示踪颗粒平均高度Fig.7 Average height of tracer particles in stagnation zone under two grid methods

图7 中初始床层高度为0.18 m,进口速度为27 m/s。可以看出，随着时间的推移，两种方法下的滞止区示踪颗粒平均高度不断下降，利用单网格方法比利用动态双网格方法滞止区示踪颗粒平均高度要低，并且利用单网格方法滞止区示踪颗粒平均床高变化较快，在t=0.5 s 滞止区内示踪颗粒平均床高趋于0，这说明在此后滞止区示踪颗粒含量很少，而利用动态双网格模型滞止区示踪颗粒平均床高变化较为缓慢，滞止区现象更为明显。此外，对动态双网格滞止区示踪颗粒平均高度在t=0.5 s 和0.9 s 出现波动情况进行分析，可能是由于床层波动造成挤压，致使示踪颗粒又往滞止区运动，造成滞止区示踪颗粒数量增多，进而使平均床高增加。

3.1.3 滞止区示踪颗粒流动宽度对比为了探究滞止区黑色示踪颗粒受喷口扰动的影响，在研究中取每个时刻所有黑色示踪颗粒与壁面的距离作为参考，得到0～2 s 内示踪颗粒的平均宽度如图8所示。

图8 两种网格方法下示踪颗粒平均宽度Fig.8 Average width of tracer particles under two grid methods

图8 中初始床层高度为0.18 m,进口速度为27 m/s。从图中可以看出，随着示踪颗粒的下滑，造成示踪颗粒在床层宽度方向逐渐延伸，从而逐渐接近喷口附近最后随同喷口气流往上运动。对比动态双网格方法与单网格方法可以发现，单网格方法示踪颗粒往喷口方向延伸速度更快，在0.8 s就延伸至喷口附近，而此时刻利用动态双网格方法其示踪颗粒才延伸至距左侧壁面0.04 m 的位置。而在2 s 时，两种方法的示踪颗粒平均宽度又接近一致，可能是由于随着时间的推移，两种网格方法下的示踪颗粒向喷口扩散的数量接近一致造成的。

通过对比单网格模型和动态双网格模型内关于示踪颗粒的运动可以明显发现，利用动态双网格模型的计算结果中滞止区区域更为明显，而在采用单网格模型模拟大颗粒稠密气固两相流时，由于气相网格不够精确，使得气流速度分布均匀，阻碍了该现象的形成。

3.2 不同速度下滞止区示踪颗粒流动高度分析

图9给出了基于动态双重网格方法下初始床层高度为0.18 m,进口速度为27、30、32 和35 m/s 时的床层滞止区示踪颗粒在2 s 时刻内的平均高度变化情况。从图中可以看出，随着时间的推移，滞止区示踪颗粒平均高度不断下降，四种速度下的床层滞止区示踪颗粒平均高度下降趋势基本保持一致,在0～0.5 s 内，可以看出随着进口速度的增加，其滞止区示踪颗粒平均高度下降速度也随之增快。此外，动态双网格滞止区示踪颗粒平均高度在t=0.5 s 和0.9 s 出现波动，可能是由于床层波动造成颗粒间挤压，致使示踪颗粒又往滞止区运动，造成滞止区示踪颗粒数量增多，进而使平均床高增加。整体来看，在计算时刻内，四种速度下床层内部都有良好的滞止区现象表现，速度的增加对滞止区颗粒高度下降影响不明显。

图9 不同速度下滞止区示踪颗粒平均高度Fig.9 Average height of tracer particles in stagnation zone at different speeds

3.3 不同速度下滞止区示踪颗粒流动宽度分析

图10 给出了基于动态双重网格方法下初始床层高度为0.18 m,进口速度为27、30、32 和35 m/s 时的床层滞止区示踪颗粒在2 s 时刻内的示踪颗粒平均宽度变化情况。从图中可以看出，随着示踪颗粒的下滑，示踪颗粒在床层宽度方向逐渐延伸，从而逐渐接近喷口附近最后随同喷口气流往上运动。从图中可以明显看出，在0～1 s 内，32 m/s 和35 m/s工况下，其延伸速度一致，较27 m/s和30 m/s的工况更快。32 m/s和35 m/s工况下在0.8 s就延伸至喷口附近，而此时27 m/s 和30 m/s 的工况才延伸至距左侧壁面0.04 m的位置。

图10 不同速度下示踪颗粒平均宽度Fig.10 Average width of tracer particles at different speeds

3.4 不同初始堆积高度下滞止区示踪颗粒流动高度分析

图11 给出了基于动态双重网格方法下初始床层高度为0.13、0.18、0.23 和0.28 m,进口速度为27 m/s时的床层滞止区示踪颗粒在2 s时刻内的平均高度变化情况。从图中可明显看出，随着床层初始堆积高度的增加，滞止区示踪颗粒平均高度下降速度也随之增加，但是初始堆积高度为0.28 m 的情况下在0～0.3 s下降速度最快，在0.3～2 s下降速度趋于平缓，这可能是由于初始堆积床层较大，而气体速度更接近最小流化速度造成的。同时对四种初始堆积床高进行对比，发现在0～0.3 s 快速下降后，在之后的下降速度都较为平缓，床层高度较大的床层，在0.3 s后存在较为明显的波动现象。

图11 不同初始高度下示踪颗粒平均高度Fig.11 Average height of tracer particles at different heights

3.5 不同初始堆积高度下滞止区示踪颗粒流动宽度分析

图12 给出了基于动态双重网格方法下初始床层高度为0.13、0.18、0.23 和0.28 m,进口速度为27 m/s 时的床层滞止区示踪颗粒在2 s 时刻内的示踪颗粒平均宽度变化情况。可以看出0～0.4 s，床层高度为0.18 m 的示踪颗粒宽度增长速度最大，床层高度为0.28 m 的变化则最慢。0.4～1.0 s，初始堆积高度为0.13 m 的宽度增长速度最快，床层初始堆积高度为0.28 m 的变化最慢。从四种初始堆积高度下其滞止区宽度变化情况来看，当初始堆积越高其宽度变化越慢，可能是由于床层高度增加导致滞止区颗粒下降速度增快，使其更接近喷口速度导致的。