三维单旋导向挡板鼓泡流化床内气固流动的CPFD 模拟

曲 悦, 潘 腾, 杨越林, 成有为, 王丽军, 李 希

(1. 浙江大学 化学工程与生物工程学院, 浙江 杭州310027;2. 浙江大学衢州研究院, 浙江 衢州 324000)

1 前 言

气固鼓泡流化床是化学工业中常见的流化床型式，例如催化裂化和甲醇制烯烃装置中的汽提器和部分立式输送管路均在鼓泡流态化下操作[1-2]。气固鼓泡流化床中存在乳相和气泡相，乳相中的气体由颗粒之间的空隙通过颗粒床层，气泡相中的气体则以气泡的形式通过颗粒床层，因此，气泡的生成、发展和运动是影响鼓泡流化床反应器特性的关键因素[3]。大尺寸气泡的存在可能导致部分气相无法有效地与固相接触、加剧颗粒的返混和夹带等问题，从而使反应器内的化学反应难以达到预期的转化率和选择性[4]。为了提高鼓泡流化床的流态化性能，可以通过设置内构件来抑制气泡生长、破碎大尺寸气泡、强化气固混合等[5]。内构件技术作为一个简单、经济的改善流态化性能的手段，目前已成功应用于丁烯氧化脱氢、苯氧化法制顺酐等工业过程中[3]。

为了研究内构件对鼓泡流化床内气固流动的影响，许多学者已成功借助计算流体力学(CFD)方法获得了实验难以测量的微观流动状况，例如气泡的尺寸、形状和频率及气固混合特性等[6-11]。石孝刚等[6]使用双流体模型 (TFM) 研究带有倾斜挡板的二维鼓泡流化床内的气固流动特性，结果表明，挡板影响区域内气泡的平均尺寸降低、气泡数量增加；且气速越高，挡板的作用越强。Gao 等[7]使用TFM 模型计算对比了带有V 型挡板和不带内构件的FCC 汽提器中的气固混合特性，结果表明，V 型挡板汽提器可以有效提高气固混合效率、降低催化剂返混，从而提高汽提器的汽提效率。Zhu 等[10]利用计算颗粒流体力学(CPFD)方法对化学链燃烧过程中的鼓泡流化床进行模拟，探究了挡板的数量、开孔率等参数对鼓泡流化床中气泡尺寸、压力梯度、颗粒分布等方面的影响，结果表明，挡板显著改善了流化床中的气固接触。Yang 等[11]使用CPFD 方法对二维单旋导向挡板鼓泡流化床进行模拟，探究单旋导向挡板对床内固含率、压差波动、气泡行为的影响，提出了二维单旋挡板在较高表观气速下破碎气泡的机制。

CPFD 是一种结合欧拉双流体模型与拉格朗日离散模型各自优点的计算流体力学方法，通过引入“计算颗粒”这一概念，将多个真实颗粒打包进行计算，在保证计算精度的同时极大地降低了计算资源和成本。目前已有许多学者采用CPFD 方法研究无内构件鼓泡流化床内的气固流动特性[12-15]，但针对三维单旋挡板鼓泡流化床的研究相对较少。本研究使用CPFD 方法耦合Igci 曳力模型，对三维单旋导向挡板鼓泡流化床内的气固流动进行模拟，首先将轴向和径向时均固含率分布模拟结果与实验数据进行对比来验证计算方法的可靠性，随后考察有、无内构件时鼓泡流化床内的气泡平均当量直径分布和气体停留时间分布。

2 计算方法

2.1 C PFD 数学模型

使用CPFD Barracuda 软件对单旋导向挡板鼓泡流化床进行模拟。在CPFD 中，气相被视为连续相，通过Naiver-Stokes 方程来描述，颗粒被视为离散相，通过拉格朗日框架下的MP-PIC 方法描述。引入“计算颗粒”这一概念，将多个具有同样物化性质和运动状态的真实颗粒打包，可大幅度降低计算成本。关于CPFD 方法的详细描述及相关表达式见文献[16]。

2.2 气固曳力模型

气固曳力是气固两相流中最重要的作用力之一，因此，选择合适的曳力模型是准确模拟和描述鼓泡流化床内两相流动的关键。通常有两种类型曳力模型被应用于模拟气固两相流：均匀模型 (homogeneous drag model) 和非均匀模型 (heterogeneous drag model)。前者以Wen&Yu 模型[17]和Gidaspow 模型[18]为代表，在应用中往往会高估曳力从而导致床层过度膨胀。为了克服均匀曳力模型的缺点，多种非均匀曳力模型被提出，例如经典的过滤曳力模型和EMMS 曳力模型。过滤曳力模型[19]是通过细网格模拟捕捉小尺度的非均匀结构，并在时空尺度上平均化，从而建立曳力系数与气固滑移速度、局部固含率等参数的本构关联式。EMMS 曳力模型[20]采用多尺度分析方法并运用能量最小原理，建立气固曳力与鼓泡流化床整体、局部参数之间的关系。

本研究采用Igci 曳力模型即Igci 等[21]基于颗粒动理学和双流体模型的细网格模拟推导出的过滤曳力模型对单旋导向挡板鼓泡流化床进行模拟计算，并通过实验数据对模拟结果进行验证。以下是Igci 曳力模型的数学表达式。

在CPFD 中，单颗粒的曳力通过式(1)和(2)描述：

式中：F为单颗粒曳力，N；mp为颗粒质量，kg；gu和pu分别为气体和颗粒的速度，m·s-1；D代表曳力方程(drag function)；ρg和ρp分别为气体和颗粒的密度，kg·m-3；dp为颗粒直径，m；Cd为曳力系数，在Wen&Yu 模型中Cd通过式(3)描述：

式中：φg为气相体积分数，Re为雷诺数。在CPFD Barracuda 软件中，Re通过式(4)描述：

Igci 曳力模型在Wen&Yu 模型基础上通过非均匀系数HIgci进行修正，如式(5)和(6)所示：

式中：DIgci和DWen&Yu分别表示经Igci 和Wen&Yu 曳力模型计算得到的式(2)中的曳力方程；、H2D和ξ均为模型参数，如式(7)～(10)所示：

式中：Δfilter为过滤尺寸，m；ut为颗粒终端速度，m·s-1；μg为气体黏度，Pa·s；φp为固相体积分数，亦称为固含率。

2.3 模拟设置

模拟对象为文献中的一套带有两层单旋导向挡板的三维冷模实验装置[9,22]，其结构如图1 所示。该装置中流态化段主体部分为φ140 mm×1 000 mm 的有机玻璃圆筒，底部采用金属烧结板作为气体分布器，内置两层单旋导向挡板，其中叶片宽20 mm、叶片倾角45°、叶片间距14 mm、挡板整体厚度14 mm，图1 中箭头表示倾斜方向。

图1 单旋导向挡板结构图[9]Fig.1 Configuration of the louver baffle[9]

实验所采用的颗粒为玻璃微珠，颗粒密度为2 450 kg·m-3，堆积密度为1 395 kg·m-3，平均粒径为53 μm。该玻璃微珠属于A 类颗粒，其粒径分布如图2 所示。在模拟计算中采用此真实粒径分布。模拟过程中所设置的参数见表1，流态化介质选择空气，计算中空气由底部均匀进入鼓泡流化床 (见图3)。流化床底部为速度边界条件，顶部为压力边界条件。在CPFD 中使用挡板组件 (Baffles)对两层单旋导向挡板内构件进行建模，如图4 所示。在三维模拟计算中，时间步长初始为5×10-4s，随后由程序自动调整以保证收敛条件判断数介于0.8 与1.5 之间。总计算时间为70 s，取最后20 s 进行时均处理和分析，从而得到计算结果。

图2 玻璃微珠累积粒度分布[9]Fig.2 Cumulative particle size distribution of glass beads[9]

表1 模拟参数Table 1 Simulation parameters

图3 单旋导向挡板三维鼓泡流化床模拟示意图Fig.3 Schematic diagram of the simulated 3D bubbling fluidized bed with louver baffles

图4 CPFD 挡板设置示意图Fig.4 Schematic diagram of baffle arrangement in CPFD

3 结果与讨论

3.1 网格数和计算颗粒数无关性分析

为了考察网格分辨率和计算颗粒数对模拟结果的影响，采用三种不同精度的网格和计算颗粒数对两层单旋导向挡板鼓泡流化床进行模拟，表2 列出了所使用的三种精度的网格数和计算颗粒数配置。图5 为轴向截面时均固含率模拟结果与实验数据对比图，图中H为分布器上方垂直高度，m。由图5 可见上述三种精度的参数配置下的模拟结果区别不大。综合考虑计算成本和精度，采用其中网格数为159 840、计算颗粒数为2 131 200 进行后续计算。

表2 网格和计算颗粒数Table 2 Parameters of grids and particle clouds

图5 不同分辨率下轴向时均固含率模拟值与实验值对比Fig.5 Comparison of simulated value and experimental data of time-averaged axial solid volume fraction at different resolutions

3.2 轴向和径向时均固含率分布

采用Igci 曳力模型分别对有、无两层单旋导向挡板鼓泡流化床进行模拟，图6 为轴向分布模拟结果与实验数据对比，图7 为H分别为0.15、0.20 和0.25 m 处水平截面上径向分布模拟结果与实验数据对比，图中r/R为归一化半径。由图可见，由Igci 曳力模型得到的径向分布模拟结果与实验值相对接近，与Yang 等[8-9]用双流体模型 (TFM) 耦合结构曳力模型所得到的轴、径向分布模拟结果趋势基本相同。图6 显示，CPFD 模拟成功地捕捉到了挡板附近区域降低现象；图7 显示，CPFD 模拟基本上可用于描述径向固含率分布情况。通过轴、径向模拟结果与实验数据的对比可知，使用CPFD 方法耦合Igci 曳力模型对两层单旋导向挡板鼓泡流化床进行模拟是可行的。

图6 有、无两层单旋导向挡板时鼓泡流化床轴向时均固含率模拟值与实验数据对比Fig.6 Comparison of simulated value and experimental data of time-averaged axial solid volume fraction of the bubbling fluidized beds with/without two louver baffles

图7 有、无两层单旋导向挡板鼓泡流化床径向时均固含率模拟值与实验数据对比Fig.7 Comparison of simulated value and experimental data of time-averaged radial solid volume fraction of the bubbling fluidized beds with/without two louver baffles

3.3 气泡平均当量直径分布

使用Bakshi 等[23]提出的多相流三维检测和追踪算法 (MS3DATA) 对有、无单旋导向挡板时的鼓泡流化床内的气泡进行统计和分析。在模拟计算过程中，每0.01 s 记录1 次床内的瞬时流场状态，取50～70 s 内的2 000 帧数据作为样本追踪和统计气泡，设置气泡相的阈值为0.7 (气相体积分数)。采用Darton[24]和Rowe 等[25]的关联式来预测鼓泡流化床内的气泡当量直径，其表达式分别为式(11)和 (12)，使用式(13)[26]来计算最大气泡直径。

式中：Db为气泡平均当量直径，m；ug为表观气速，m·s-1；umf为起始流化速度，m·s-1，实验中所用的玻璃微珠的起始流化速度为 0.005 1m·s-1；h0为常数，取h0= 0.03 m[7]；g为重力加速度，m·s-2；Dbm为最大气泡直径，m；At为气体入口截面积，m2。

图8 为有、无两层单旋导向挡板时鼓泡流化床内的Db随H统计分析结果。由图可见，有单旋导向挡板时的床内气泡特征不同于无导向挡板即空筒鼓泡流化床时的气泡特征。在空筒鼓泡流化床中，随H增大Db不断增大，且未出现气泡破碎现象；而在有单旋导向挡板鼓泡流化床中，第一层挡板下方的Db大小以及变化规律与空筒鼓泡流化床相似，但经过第一层挡板时，Db变小，这说明单旋导向挡板对气泡有一定破碎作用；经过第一层挡板后，气泡又重新开始聚并，Db逐渐增大；同理，经过第二层挡板时，Db再次变小。由图还可见，大部分床高处的气泡尺寸模拟结果与Darton 关联式预测结果趋势相吻合，但在流化床上界面附近，有、无两层单旋导向挡板时鼓泡流化床内的Db均出现明显变小的趋势，这是由于MS3DATA 算法基于体积计算气泡当量直径，而在此区域气泡已经开始穿越流化床界面，从而导致该算法统计的气泡体积偏小。图8 表明，有挡板的鼓泡流化床内Db小于空筒鼓泡流化床的，但减小程度相对有限，这可能是由于实验中所使用的单旋导向挡板叶片数量较少、开孔率大，从而起到的气固搅拌作用相对较弱，因此，在实际应用中需要进一步对此单旋导向挡板的开孔率、叶片数量和角度等进行优化。

图8 气泡平均当量直径沿床高分布图Fig.8 Axial distribution of average equivalent bubble diameters

3.4 气体停留时间分布

使用脉冲示踪法研究有、无两层单旋导向挡板时鼓泡流化床内气体返混现象和气体停留时间分布(RTD)。计算中用空气作为示踪剂，在第20 s 至20.01 s 内注入示踪剂，令其流速与主风流速保持一致。在此时间间隔内，设置主风流速为0，以保证示踪剂的注入不会对床内的气固流动状态产生影响。记录每个时间步时流化床出口处示踪剂浓度，每200 个时间步内对示踪剂浓度取一次平均，得到有、无两层单旋导向挡板时鼓泡流化床内的气体停留时间分布曲线，如图9 所示，图中θ为对比时间(无因次)，E(θ)为停留时间分布密度函数。

图9 气体停留时间分布图Fig.9 Profiles of gas residence time distribution

由图9 可见，有、无挡板时鼓泡流化床内的气体停留时间分布趋势基本相同，有单旋导向挡板鼓泡流化床内的停留时间分布曲线相比于空筒鼓泡流化床略微向右偏移，这说明单旋导向挡板对大气泡的破碎作用在一定程度上抑制了气体的短路问题。进一步对停留时间分布的数学特征进行计算和分析，得到气体在鼓泡流化床内的平均停留时间和方差σ2列于表3，可见，有、无挡板的鼓泡流化床的返混程度均相对较小，与平推流理想反应器接近，两者的平均停留时间和方差基本相同，表明设置两层单旋导向挡板后没有明显影响流化床内的气体返混。

表3 停留时间分布的数学特征Table 3 Mathematic characteristics of RTD

4 结 论

采用计算颗粒流体力学 (CPFD) 方法耦合Igci 曳力模型对有、无两层单旋导向挡板的鼓泡流化床进行模拟，模拟计算得到的轴向和径向时均固含率分布与实验数据较为吻合，表明此方法用来模拟单旋导向挡板鼓泡流化床的三维气固流动可行。使用MS3DATA 算法统计和分析了该鼓泡流化床内的气泡平均当量直径沿床高的分布，结果表明，单旋导向挡板对大气泡有一定的破碎作用，有挡板的鼓泡流化床内气泡平均尺寸小于空筒鼓泡流化床。使用空气作为示踪剂考察该鼓泡流化床内的气体返混和停留时间分布，结果表明，两层单旋导向挡板的存在没有明显影响鼓泡流化床内的气体返混，有、无挡板的鼓泡流化床内的气体平均停留时间和方差基本相同，返混程度均与平推流理想反应器接近。