曳力模型对流态化两相流动数值计算影响的研究

，， ，

(1.上海电力建设启动调整试验所，上海 200031;2.动力机械与工程教育部重点实验室(上海交通大学),上海 200240)

0 引言

近年来，微颗粒流态化在煤粉流态化燃烧、物料干燥等诸多领域应用越发广泛,其所涉及的多相流动和传递特性成为国内外倍受重视的研究领域[1-2]。

诸多学者对流态化数值模型及不同类别颗粒的流态化进行了详细的研究与分析[2-8]：祁海鹰[3]等对流态化模拟中的曳力模型进行了论述，分析了EMMS理论在非均匀流动的颗粒团尺寸、内部固含率等关键参数上的缺陷、修正方法及实践验证结果。郑晓野[4]等采用改进的曳力模型对2D鼓泡流化床的流化特性进行了分析和验证，得出了更符合实际的模拟结果。郭雪岩[2]等采用欧拉-欧拉模型对Geldart D类颗粒气固流化床非定常传热流动进行了模拟，验证了6种气固传热模型的壁面传热系数实验关联式。晁东海[5]等模拟了大颗粒流化床在不同密度、布风装置及曳力模型下两相流动，得到了大颗粒流化特点、颗粒体积分率分布、床层压力瞬时变化、床层碰撞比以及颗粒径向速度和空隙率轴向分布规律。王会宁[6]等采用双流体模型就脱硫塔内气固两相流进行数值模拟，得到了塔内颗粒速度及浓度场的分布。汪新智[7]等和周明哲[8]采用数值方法对鼓泡流化床内颗粒分离行为进行模拟研究并结合实验数据对所选用的数值曳力模型进行对比和验证。刘富爽[9]等采用数值模拟方法对栅格型煤粉分布器气固两相流进行分析，为设计提供建议。

本文在先前研究的基础上，采用ANSYS FLUENT对2D流化床计算模型采用不同曳力模型和颗粒粒径的细颗粒进行数值模拟，旨在为后续研究提供参考。

1 数值模拟方法

1.1 数值计算理论

本次数值模拟研究采用欧拉-欧拉方法[10]，气固间动量交换采用经验曳力模型[3]进行描述,相应的控制方程组为

εg+εp=1

(1)

式中εg和εp——气相和固相体积分率。

(2)

(3)

(4)

(5)

(6)

(7)

式中CD0——单颗粒曳力系数。

气相湍流输运方程采用k-ε双方程模型[9,12]

(8)

(9)

其中(8)为k方程，(9)为ε方程，湍流粘性系数为

(10)

湍流产生率为

(11)

式中σk——湍流动能普朗特数;

σε——湍流动能耗散率普朗特数，相应模型系数取值：C1ε=1.44,C2ε=1.92,σk=1.00,σkε=1.30,Cμ=0.99。

1.2 数值网格和计算设置

数值计算采用宽度为0.186 m，高度为1.6 m的2D流化床简化模型，其计算网格和边界设置如图1所示。

图1 2D流化床模型流态化模拟试验网格及边界设置

该套网格有88 911个运算节点，计算域下部为速度进口边界，通过自定义函数设置其为抛物线分布，其中心最大速度为1 m/s，近壁面处速度为0.5 m/s。上部为压力出口边界，压力为标准大气压。其他壁面均为标准壁面边界。流态化过程不考虑传热、化学反应等过程作用。静止床层高度为0.382 m，固相体积分率考虑颗粒粒径在0.56～0.85范围内选取。数值模拟选取k-e计算模型，曳力模型分别采用Gidaspow模型、Syamlal-obrien模型和wen-yu曳力模型，颗粒碰撞还原系数为0.9。计算初定时间步长为0.000 5 s。

颗粒粒径参考照Geldart理论[14]进行选取[1]，如表1所示。

表1数值模拟的模型参数设置

颗粒种类D/μm密度/kg·m-3粒径/μm固含率5003 20050055503 200506553 200575

2 模拟结果与分析

2.1 不同曳力模型的分率随时间变化

在曳力模型的对比优选中，采用密度3 200 kg/m3，直径5 μm的颗粒为试验对象，分别选取syamlal-obrien、wen-yu、gidaspow三种曳力模型进行数值模拟计算。模拟设定初始固相体积分率均为0.75。

图2 Syamlal-obrien曳力模型下固相体积分率随时间变化

图2～图4为不同曳力模型的数值模拟结果图示。

图3 wen-yu曳力模型下固相体积分率随时间变化

图4 Gidaspow曳力模型下固相体积分率随时间变化

对比数值模拟结果，可以看出不同作用力模型下模拟出来的流态化过程和现象有着较大的差异：

Syamlal-obrien模型从气泡产生到流态化的整个过程中形态比较规整，模拟前期物相分布均匀对称，但其气泡的生成方式和形态同现实中有较大差异。此外在气泡破裂后的流态化过程中，模拟的流态化情况比较生动与真实。Wen-yu模型结果比较符合实验中气泡生长的形态[15]，但其流态化过程的模拟比较缓慢，在首次产生的气泡破裂后，才会产生平缓的流态化。Gidaspow曳力模型结果中气泡的成长与流态化过程良好，与真实实验结果较为接近。由于初始条件设置的颗粒较轻，固相区域上部的漂浮翻腾运动模拟也有较好效果。在细颗粒环境下气泡的产生与成长过程中，其产生方式、形态等特征都与现实中比较接近。

图5 典型的二维气泡及尾涡图示[15]

图5为实验中观察到的类似工况典型二维气泡及尾涡图示，可以看到其与Gidaspow曳力模型下获得的模拟结果较为类似。此外结合领域相关科研人员的研究推荐[2-5,16-19]，本文选取Gidaspow曳力模型进行后续D=500 μm和D=50 μm粒径的流态化差异的模拟计算。

2.2 固相体积分率随时间变化

图6、图7为模拟得到的颗粒粒径为500 μm和50 μm时流化床内固相体积分率随时间变化图示。

图6 粒径D=500 μm时固相体积分率的分布随时间的变化

图7 粒径D=50 μm时固相体积分率的分布随时间的变化

由图6可以看出，0 s时没有气流进入，床层处于静止状态。在0.1 s时有速度呈类抛物线分布的气流通过进气口进入床层，颗粒在气流的作用下上移小段距离后，在底部开始生成气泡；在0.3 s处，气泡逐渐长大，并在底部再次产生了两个近似对称的小气泡；在0.5 s时，最初产生的气泡不断长大并产生不规则运动；在0.6 s、0.7 s时，大气泡运动到静止床层的顶部，开始破裂。在1 s时，大气泡破裂，第二层两个气泡继续向上，床的底部不断产生新的气泡群。1.5s处，底部开始产生不规则的翻腾流化状态。

图7为粒径D=50 μm的流态化情况随时间变化的图示。可以看出其气泡生成的初期状态基本与颗粒粒径D=500 μm时的流态化情况基本一致，仅仅在气泡尺寸和形状方面有着一定的差异，但后期两者的流态化状态有着显著的不同，这与由于颗粒粒径以及固含率的设置有一定关系。

3 结论

本文以FLUENT为工具，对D=5 μm微颗粒在流化床中的流态化情况进行了模拟，获得了不同曳力模型下的固相体积分率随时间的图示，通过对比分析，采用优选的曳力模型进一步就D=500 μm和D=50 μm的颗粒流态化情况进行了数值模拟，开展了定性分析。研究表明syamlal-obrien、wen-yu、Gidaspow三个曳力模型在颗粒直径密度3 200 kg/m3，直径5 μm，固相体积分率均为0.75的条件下，Gidaspow曳力模型表现更佳。