基于介尺度曳力模型双组分颗粒混合的CFD模拟

宋晓皎， 王 帅

(1. 山西能源学院能源与动力工程系， 山西晋中 030600； 2. 哈尔滨工业大学能源科学与工程学院， 黑龙江哈尔滨 150001)

流化床反应器被广泛应用于催化裂解、 CO2捕集和煤粉颗粒燃烧等工业化生产中[1-3]。颗粒直径和密度等物理特性的不同，会导致多组分气固流化床中颗粒的混合分离现象的产生[4]，对气固反应特性及燃料的转化率产生直接影响，因此，对气固流化床中双组分颗粒分离特性的深入研究，将有利于获取系统最佳操作条件，进而提高燃料利用率。

近年来，数值模拟成为研究多组分颗粒流化床系统的重要方法[5-7]。Feng等[8]利用离散颗粒模型对流化床内双组分颗粒的分离和混合过程进行数值模拟， 并指出气固相间的作用力引发了流态化和颗粒间分离现象的产生。 Owoyemi等[9]利用数值模拟的方法， 对三维流化床内金红石颗粒的混合和分离过程进行研究， 并通过实验证实模拟结果的正确性， 指出多组分颗粒混合系统中固-固相间作用力和固相剪切力需进行进一步研究。 Mazzie等[10]利用矩阵直接求解法， 分析了流化床内多组分颗粒混合现象， 并对数值扩散对模型预测的影响进行了讨论。 Sharma等[11]讨论了颗粒密度和直径等物理特性对生物质和生物炭颗粒混合特性的影响， 结果显示生物质颗粒密度对混合物的分离和混合特性影响最为显著。

颗粒分离现象主要受流化床内气泡移动和颗粒湍动的影响。基于麦克斯韦速度耗散假设，非均一颗粒拟温度颗粒动力学理论被提出，并成功应用于双组分颗粒的鼓泡流化床的数值研究中[12-13]。Iddir等[14]基于能量非均一性考虑，进一步发展、 完善了混合物颗粒动力学理论，并指出该理论能较好地应用于剪切流和反应流系统中。Chao等[15]考虑到颗粒间的摩擦影响，得到双组分颗粒动力学理论，模拟分析显示双组分上升管中主要呈现出连续性流体的脉动现象。

双组分气固曳力模型的构建对于多组分气固流化床系统的数值模拟具有十分重要的作用。Beetstra等[16]基于格子玻尔兹曼模拟的方法获得应用于多组分系统的校正因子。Olaofe等[17]使用此校正因子对双组分流化床内颗粒的混合和分离特性进行数值研究，发现源于直接数值模拟(direct numerical simulation, DNS)方法的曳力关系式可以获得更好的模拟预测结果。Zhang等[18]应用离散颗粒方法，通过比较颗粒分离程度和气泡生成频率对多组分曳力模型的影响进行了研究；但此方法的应用受到颗粒直径和系统尺寸等多方面的限制，因此，非常有必要发展一种考虑到床层非均一性、适用于大尺寸系统粗网格模拟的曳力模型[19-20]。

Shi等[21]提出了以气泡为介尺度的能量最小介尺度曳力模型，并将其应用于鼓泡流化床的数值研究中。Lunge和Wang等[22-23]将曳力模型拓展应用于双组分气固流化床中，并且对双组分混合系统中的颗粒分离混合现象以及反应特性进行研究，但忽略了其他颗粒相的影响。Zhou等[24]针对2类颗粒聚团现象对曳力模型进行改进，并将改进后的模型应用于循环流化床提升管的数值研究中，预测结果与实验结果较吻合。然而，基于能量最小多尺度 (energy minimization multi-scale，EMMS)曳力模型、同时考虑乳化相中双组分颗粒的相互影响及气泡影响并应用于鼓泡流化床中的双组分气固曳力模型的相关研究较为少见。

本文中提出基于乳化相中颗粒相互影响并考虑以气泡为介尺度的曳力模型，结合颗粒动力学理论和欧拉(Euler)多流体模型，对具有双组分颗粒的三维鼓泡流化床中颗粒的混合和分离特性进行模拟研究，并通过实验数据验证双组分气固曳力模型的合理性，并进一步分析研究外界操作条件和颗粒物性参数对双组分颗粒的分离和混合过程的影响。

1 数学模型

假定2种固相分别为单一直径和密度的球形颗粒处于等温流场内，建立基于欧拉多流体模型和颗粒动力学理论、 以气泡为介尺度、 考虑乳化相中两相颗粒相互影响的双组分气固曳力模型。在验证双组分气固曳力模型准确性的基础上，可利用此曳力模型对双组分颗粒分离和混合过程进行数值模拟。

1.1 控制方程

在Euler多流体模型[23]中，固体颗粒和气体被作为拟流体和连续性介质进行处理。其质量守恒方程为

(1)

(2)

式中：φg为气相体积分数；ρg为气相密度， kg/m3；Ug为气相表观速度， m/s；xi为三维空间中某一方向的坐标，i=1，2，3；m、n为反应物质编号；Ng为气相反应物种类总数；Rgn为气相间转化的单位体积质量流率， kg/(m3·s)；φm为固相体积分数；ρm为固相密度， kg/m3；Umi为三维空间中某一方向的固相表观速度，i=1，2，3，m/s；Nm为固相反应物种类总数；Rmn为固相间转化的单位体积质量流率，kg/(m3·s)。

其动量守恒方程为

(3)

(4)

式中：Ugi和Ugj为三维空间中2个方向的气相表观速度，i=1，2，3，j=1，2，3，m/s；xj为三维空间中某一方向的坐标，j=1，2，3，m；Pg为气相压强，Pa；τgij为二维应力张量，Pa;gi为三维空间中某一方向重力加速度，i=1，2，3，m/s2；β为曳力系数；Umi和Umj为三维空间中2个方向的固相表观速度，i=1，2，3，j=1，2，3，m/s；Pm为固相压强，Pa；为向量微分算子符号。

其颗粒拟温度方程为

(5)

式中：Θm为固相颗粒拟温度，m2/s2；∏m为由于固相间滑移而产生的动能，J/(m3·s)；γ为因相间滑移而导致的能量碰撞耗散，J/(m3·s)；km为固相颗粒运动碰撞的热传导，J。

1.2 气泡介尺度曳力模型

本文中考虑气泡介尺度及乳化相中颗粒相的相互作用对鼓泡流化床内相间曳力的影响，在Shi等[21]提出的单组分气固曳力模型的基础上，构建适用于双组分混合物流化床系统的气固曳力模型。基于气泡内无颗粒相存在的假设，得到乳化相的平均物性[25]为

ρe=ρp1φep1+ρp2φep2+ρgφeg，

(6)

μe=μg[1+2.5(1-φep)+10.05(1-φep)2+0.002 73 exp(16.6-16.6φep)]，

(7)

Ue=(ρgUge+ρp1Uep1+ρp2Uep2)/ρe，

(8)

式中：ρe为乳化相的平均密度，kg/m3；ρp1和ρp2分别为乳化相中颗粒p1和p2的密度，kg/m3；φep1和φep2为乳化相中颗粒p1和p2的体积分数；μe和μg分别为乳化相和气体的动力黏度，N·s/m2；φep为乳化相颗粒总体积分数；Ue为乳化相表观气体速度，m/s；Uge、Uep1和Uep2为乳化相中气相、颗粒p1和p2的表观气体速度，m/s。

通过引入加速度项来表征非稳定特性，气泡相和乳化相中颗粒的力平衡公式为

(9)

(10)

式中：db为气泡直径，m；Cdb为传统曳力模型的气泡相有效曳力系数；Ueb为乳化相中气泡表观气体速度，m/s；ab为气泡相的加速度，m/s2；Cde1和Cde2分别为传统曳力模型乳化相中颗粒p1和颗粒p2的有效曳力系数；ae为乳化相的加速度，m/s2。

为了表征气泡相和乳化相的附加质量力和惯性的差异[26]，其差值为

(11)

式中：ρb为气泡相的平均密度，kg/m3；δb为气泡相体积分数；σ为局部颗粒体积分数浓度变化量。

基于气体和固体体积分数的定义以及气固质量守恒关系，可得到关系式：

(12)

Ug=Uge(1-δb)+δbUb，

(13)

(14)

式中：Usp1和Usp2分别为颗粒p1和颗粒p2的表观滑移速度，m/s。

为封闭模型求解，设定输送悬浮颗粒所需能量最小作为颗粒流动稳定准则[21]，即

(15)

式中：fb为气泡中气体体积与整个床层气体体积分数比值。基于以上方程的求解，可获得相关局部结构参数，进而得到双组分颗粒各固体相的曳力系数表达式

(16)

(17)

式中：Usp1和Usp2为2种颗粒表观滑移速度，m/s。基于此，为更好地表征考虑气泡介尺度后曳力模型的变化情况，引入非均质系数Hd， 气泡介尺度曳力系数计算式变为

(18)

Hd=βb/βw，

(19)

式中：Cd为单个颗粒的曳力系数；βb和βw分别为气泡介尺度曳力系数和传统WEN-YU曳力系数。

1.3 混合指数

为了更加深刻地分析多组分颗粒分离和混合特性，引用Owoyemi等[27]提出的混合系数定义式

M=wt/wo，

(20)

式中：wt为床层顶部(25%)下沉颗粒的平均质量分数；wo为整个床层中下沉颗粒的平均质量分数。

1.4 几何模型及计算参数

参考Sun等[28]模拟试验中的几何模型，建立的三维流化床反应器模型的结构示意图如图1所示。为了研究颗粒的混合和分离特性，以砂砾颗粒S1和稻壳颗粒S2作为床料，2种颗粒的平均密度分别为2 600、 950.6 kg/m3，平均直径分别为0.44、 1.54 mm；气相采用无滑移边界条件，固相采用Sinclair和Jackson边界条件[29]。

图1 三维流化床反应器模型的结构示意图Fig.1 Schematic diagram of three dimensional fluidized bed reactor model

三维模拟计算的初始条件以及物料的材料属性如表1和表2所示。

表1 三维模拟的初始条件

表2 材料属性

2 结果与分析

2.1 非均质系数

为更直观地呈现入口气体表观速度对两相颗粒非均质系数的影响，在Sun等[28]选取入口气体表观速度为0.79 m/s实验条件的基础上，添加入口气体表观速度为1.50、 2.00 m/s这2种工况进行分析对比。入口气体表观速度、床层空隙率和颗粒体积分数对非均质系数的影响如图2所示。由如图2可知，非均质系数随床层空隙率近似呈指数形式增长，不同颗粒相体积分数增长速率有所不同，这说明双组分颗粒的非均质系数不仅仅取决于床层空隙率，还会受到不同粒径的颗粒体积分数的影响，所以，利用单组分颗粒介尺度曳力模型对双组分颗粒系统进行数值模拟是不合理的；除此之外，从图2还可明显观察到，随着气体速度的增大，固相非均质系数呈逐渐减小趋势，说明较高气速时介尺度结构对曳力大小的影响会更显著。

a)砂砾b)稻壳图2 入口气体表观速度、床层空隙率和颗粒体积分数对非均质系数的影响Fig.2 Effects of inlet gas apparent velocity, bed voidage and particle volume fraction on heterogeneity coefficient

不同颗粒体积分数下非均质系数随床层空隙率的变化如图3所示。由图3可知，两相颗粒的非均质系数都随着颗粒直径的增大逐渐增大；随着床层空隙率的增大，砂砾颗粒的非均质系数增长幅度逐渐减小，而稻壳颗粒直径非均质系数增长情况则相反；较小直径颗粒的体积分数变化对较大直径颗粒的非均质系数变化趋势的影响较大，并且直径较小颗粒的非均质系数相对较小，说明较小直径或质量的颗粒对相间曳力的非均质性影响程度更大，相间曳力与非均质结构的关联性更强。

a)砂砾b)稻壳图3 不同颗粒体积分数下非均质系数随床层空隙率的变化Fig.3 Variation of heterogeneity coefficient with bed voidage under different particle volume fraction

部分流化床反应器并不总是处于常温工作状态，温度不同必然会影响到气相物理特性的变化，间接影响床层非均质特性及相间曳力。不同床层温度条件下非均质系数随床层空隙率的变化如图4所示。由图4可知，随着床层温度的升高，固相非均质系数逐渐减小，稻壳颗粒的非均质系数减小程度较小，这主要是由于床层温度的增大直接导致气体黏度增大，密度减小，使得相间曳力发生变化，导致床层温度变化对较轻颗粒的非均质系数影响较大，而对较重颗粒的非均质系数影响较小。

a)砂砾b)稻壳图4 不同床层温度条件下非均质系数随床层空隙率的变化Fig.4 Variation of heterogeneous coefficient with bed voidage at different bed temperatures

2.2 鼓泡流化床双组分颗粒混合模拟结果

参照文献 [28]中在入口气体表观速度分别为0.58、 0.79 m/s这2种工况下进行的双组分颗粒混合实验，利用气泡为介尺度的双组分气固曳力模型和传统Gidaspow曳力模型进行数值模拟，稻壳质量分数沿床层分布情况的模拟结果与实验结果对比如图5所示。由图5可知，在2种入口气体表观速度条件下，床层顶部的模拟结果均略微小于实验结果，但是整体上模拟结果与实验结果吻合度较高；在床层上部气泡介尺度曳力模型模拟结果比传统Gidaspow曳力模型模拟结果更接近实验结果，这主要是由于床层上部气泡尺寸更大，对相间曳力的影响更加明显，因此，以气泡介尺度的双组分气固曳力模型能够获得更准确的预测结果，后续研究均采用该模型。

考虑到2种颗粒动力学特性差异，为更直观地呈现入口气体表观速度对两相颗粒混合和分离特性的影响，选取入口气体表观速度分别为0.31、 0.79 m/s， 稻壳颗粒直径分别设为0.15、 0.44、 1.54 mm和床层温度分别为298、 398、 498 K为后续数值模拟的初始条件。

a)气体入口表观速度为0.58 m/sb)入口气体表观速度为0.79 m/s图5 稻壳质量分数沿床层分布情况模拟结果与实验结果的对比Fig.5 Comparison of simulation results and experimental results of rice husk mass fraction distribution along bed

在t=7 s(t为进气时间，以下同)时刻， 入口气体表观速度对瞬时颗粒浓度分布云图的影响如图6所示。 由图6可知， 入口气体表观速度为0.31 m/s时， 床层内部的气泡数量较多并且体积较小； 当入口气体表观速度增加到0.79 m/s时， 稻壳颗粒在床层上部的颗粒浓度相对高于床层底部的浓度， 但是两相颗粒分离程度不明显， 说明随着入口气体表观速度的增大， 气泡数量减少， 体积增大， 固相颗粒混合程度逐渐加大， 床层膨胀高度增大。 这是由于， 入口气体表观速度的增大使得气泡快速聚集成大直径气泡， 到达稀密两相分界面处时发生破裂， 对颗粒携带作用增强， 进而使得床内固相混合程度加大， 床层膨胀高度增加。

a)入口气体表观速度0.31 m/sb)入口气体表观速度0.79 m/s图6入口气体表观速度对瞬时颗粒浓度分布云图的影响Fig.6 Influence of inlet gas apparent velocity on instantaneous particle concentration distribution neogram

当入口气体表观速度为0.31 m/s时，稻壳颗粒直径分别设为0.15、 0.44、 1.54 mm，在t=7 s时刻，稻壳颗粒直径对瞬时颗粒浓度分布云图的影响如图7所示。由图7可知，随着稻壳颗粒直径的减小，颗粒最小流化速度逐渐较小，使得两相颗粒分离程度逐渐加大，双组分颗粒的分离和混合还取决于颗粒的物性。

在不同入口气体表观速度条件下， 稻壳颗粒直径对颗粒混合指数随的影响如图8所示。 由图8可知， 随着稻壳颗粒直径的减小， 混合指数逐渐减小， 分离程度逐渐加大； 当入口气体表观速度减小时， 混合指数也随之减小， 因此， 颗粒物性和外界操作条件的合理匹配将是实现双组分颗粒高效分离的关键之举。

a)dp2=0.15 mmb)dp2=0.44 mmc)dp2=1.54 mm图7 稻壳颗粒直径对瞬时颗粒浓度分布云图的影响Fig.7 Effect of rice husk particle diameter on instantaneous particle concentration distribution neogram

图8 稻壳颗粒直径对颗粒混合指数的影响Fig.8 Effect of grain diameter of rice husk on grain mixing index

流化床反应器工作环境有时处于较高的温度。 在入口气体表观速度为0.31 m/s时， 设定床层温度分别为298、 398、 498 K， 在t=7 s时刻的床层温度对瞬时颗粒浓度分布云图的影响如图9所示。 由图9可知， 随着床层温度的升高， 气体黏度增大， 颗粒相容易黏结， 进而导致两相颗粒分离程度逐渐减弱。

床层温度对颗粒混合指数的影响如图10所示。由图10可知，随着床层温度的升高，混合指数逐渐增大，分离程度逐渐减弱。 这主要是由于随着床层温度的升高， 气体黏度增大， 密度减小， 使得气固曳力系数减小，最终导致非均质系数减小； 两相颗粒的非均质系数下降幅度不同，使得差值逐渐减小； 同时，气体黏度增大，也使得颗粒较易黏结，不易分离，因此，除了通过改变颗粒物性和外界操作条件以外，还可以适当降低床层温度来提高颗粒分离程度。

a)298 Kb)398 Kc)498 K图9 床层温度对瞬时颗粒浓度分布云图的影响Fig.9 Effect of bed temperature on instantaneous particle concentration distribution nephogram

图10 床层温度对颗粒混合指数的影响Fig.10 Effect of bed temperature on particle mixing index

3 结论

以气泡为介尺度建立了双组分气固曳力模型，在三维鼓泡流化床反应器中模拟了双组分颗粒的分离和混合过程，得到以下结论：

1)非均质系数随床层空隙率近似呈指数增长；固相颗粒直径的减小、气体入口表观速度的增大以及床层温度的升高都将会加强两相间曳力的非均质性；相比于传统的Gidaspow气固曳力模型，双组分气固曳力模型的预测精度更高。

2)随着入口气体表观速度的增大，固相颗粒混合程度逐渐加大，床层膨胀高度增大。随着稻壳颗粒直径和入口气体表观速度的减小，两相颗粒分离程度加大，混合指数随之减小；颗粒物性和外界操作条件的合理匹配有利于双组分颗粒的分离；同时，随着床层温度的升高，两相颗粒分离程度逐渐减弱，混合指数逐渐增大。

在今后的研究中，将进一步优化双组分气固曳力模型，并将其拓展应用到反应性流化床系统中，模拟分析流化床内的气固反应和传质特性。