基于鼓泡流化床的新型曳力模型的验证分析

李慧君， 王业库， 张久意

（华北电力大学能源动力与机械工程学院，河北保定 071003）

鼓泡流化床由于具有气固两相接触与混合充分、温度分布均匀和传热传质效率高等优点，已被广泛应用于工业生产中［1-4］。流化床内的气固流动本质上是动态的，通过实验研究流化床费用较高，条件苛刻，很难获取有效的实验数据。随着计算机模拟技术的快速发展，计算流体力学（CFD）被广泛用于研究鼓泡流化床的气固流动特性。

在气固流动模型中，欧拉-欧拉双流体模型（TFM）因计算量较小且准确性较高而被广泛使用［5］。气固两相间的作用力非常复杂，一般情况下，相较于曳力，其他作用力可忽略不计，因此仅考虑曳力作用，进而曳力系数的选择对于流化床的模拟精度影响很大［6］。TFM包含常用的曳力模型，即Wen Yu模型、Ergun模型和Gidaspow模型等均匀模型，其中前两者均有其合适的颗粒体积分数范围；而Gidaspow模型虽然没有具体颗粒体积分数范围限制，但没有考虑气泡和颗粒团簇等非均匀结构。气泡和颗粒团簇的特性很大程度上会影响流化床的流体动力学和运行效率［7-10］，因此引出介尺度结构的亚网格曳力模型［11］。李静海等［12］在考虑颗粒团簇的基础上提出了能量最小多尺度（EMMS）模型，该模型以及以此为基础的其他模型均获得了很高的精度［13-15］。Shi等［16］基于鼓泡流化床提出EMMSBubbling模型，发现该模型能精准预测出床内的颗粒分布状态以及颗粒轴向和径向速度的变化。

EMMS模型虽然模拟精度较高，但受到颗粒本身物性参数的影响，且没有普适性，很难应用于实际工程。笔者根据不同曳力模型的特点，在不同颗粒体积分数区间选择不同的曳力模型。为克服新型曳力模型不连续性的问题，在颗粒体积分数分界点引入光滑函数使其连续，在此基础上耦合TFM，对使用典型B组颗粒的二维鼓泡流化床进行模拟，并将其模拟结果与Gidaspow模型的模拟结果以及实验值进行对比，选取最适合B组颗粒的修正因子。

1 数学模型

颗粒整体的运动表现出类似流体的运动规律，可以把颗粒相视为具有相同颗粒性质的一种拟流体［17］，采用TFM［18］结合颗粒动力学（KTGF）［19］对鼓泡流化床进行模拟，通过求解气固两相的各项连续性方程、动量守恒方程以及用来封闭方程的本构方程，获得模拟结果。

1.1 控制方程

气相质量方程［17］为：

式中：φ为体积分数；ρ为密度，kg／m3；u为平均速度，m／s；mgs为气固间的质量源项，kg／（m3·s）；t为时间，s；下标g、s分别表示气相和固相（即颗粒）。

气相动量方程为：

式中：τ为剪切应力，Pa；p为压力，Pa；β为曳力，kg／（m3·s）；g为重力加速度，m／s2。

气相能量方程为：

式中：λ为导热系数，W／（m·K）；H为比焓，J／kg；T为温度，K；hgs为气相与固相间的传热系数，W／（m2·K）。

固相质量方程为：

固相动量方程为：

固相能量方程为：

式中：α为空隙率。

式中：I为单位张量；kΘs为能量的扩散系数；γΘs为拟颗粒温度因非弹性碰撞引起的能量耗散系数；φgs为颗粒随机波动动能。

1.2 本构方程

各相剪切应力τi

［18］为：

式中：V为瞬时速度，m／s；ξ为体积黏度，Pa／s，对于气相而言，ξ为0 Pa／s；μ为剪切黏度，Pa／s；下标i可表示g或s。

固相压力ps为：

式中：go为径向分布函数；ess为颗粒弹性碰撞后的恢复系数。

固相体积黏度ξs为：

式中：d为直径，m。

径向分布函数为：

式中：Φ为内摩擦角，（°）；I2D为偏应力张量第二部变量，N／m2。

能量的扩散系数为：

拟颗粒温度因非弹性碰撞引起的能量耗散系数为：

颗粒随机波动动能为：

考虑到颗粒的聚团作用，将鼓泡流化床分为密相（φs＞0.4）、过渡相（0.1≤φs≤0.4）和稀相（φs＜0.1）。密相采用Ergun模型，稀相采用Wen Yu模型［19］，过渡相采用经过介尺度修正的Mckeen模型［11］。3种曳力模型计算式分别为：

式中：βErgun、βWen Yu和βMckeen分别为Ergun模型、Wen Yu模型和Mckeen模型的曳力，kg／（m3·s）；Re为雷诺数；C为修正因子，对于大颗粒常取0.5～0.99。

上述3种曳力模型在固相体积分数分界点处并不连续，因此引入光滑函数ψ1和ψ［19］2，使其成为连续的曳力模型。

利用式（23）和式（24）连接3种曳力模型，重新构建新型曳力βgs的数学模型，其计算式为：

鼓泡流化床常用B组颗粒作为流化颗粒，因此新型曳力的修饰因子分别选取0.7、0.8和0.9。同组颗粒具有相似的物理性质。通过将TFM与新型曳力模型和Gidaspow曳力模型耦合，对比不同曳力模型模拟结果与实验值的误差，得到B组颗粒的修正因子。

2 几何模型

根据文献［20］中的实验进行参数设定，以验证新型曳力模型的正确性。该实验将3个差压传感器置于距床层底部0.03 m、0.3 m和0.6 m的位置上，横跨气体分布器。在表观气速为0～0.8 m／s的条件下对总压降和床层膨胀率进行监测，当流化床流动达到稳态后，记录总压降和床层膨胀率的测量值。使用反射式光纤探头PC-4粉末空隙测量仪对床层空隙率进行测量，该探测器包含一束光纤，由于颗粒直径比光纤束直径小得多，光线被测量床层中的颗粒反射，从输出电压可获得瞬时固相体积分数。该实验在高为1.0 m、宽为0.28 m、厚度为0.025 m的二维有机玻璃柱中进行。球形玻璃微珠直径为250～300μm，密度为2 500 kg／m3，属于典型的B组颗粒，在环境条件下采用空气进行流化。静态床层高度为0.4 m，固相体积分数为0.6。采用二维简化鼓泡流化床进行模拟，模拟时间为8 s。二维鼓泡流化床模型示意图见图1。入口边界条件为速度入口，出口边界条件为自由出口，时间步长为0.001 s，最大迭代次数为20，收敛准则为0.001。具体给定的模拟参数见表1。

图1 二维床模型的示意图及网格Fig.1 Schematic diagram of the two-dimensional bubbling fluidized bed model with its grid

表1 主要模拟参数和条件Tab.1 Main simulation parameters and conditions

3 结果分析

3.1 网格敏感性分析

选取不同网格数进行模拟，其中静态床层区域的网格划分与其他区域相同，并采用Gidaspow曳力模型进行评估。不同网格数下轴向颗粒体积分数分布见图2。由图2可知，3种网格数下轴向颗粒体积分数分布趋势一致。网格数为2 800时，床层底部的轴向颗粒体积分数相比其他2种网格数较小，随着网格数的增加，床层底部轴向颗粒体积分数逐渐增大。与网格数为2 800时相比，网格数为11 200时轴向颗粒体积分数的相对偏差约为23.86%；网格数为28 000时轴向颗粒体积分数的相对偏差约为29.17%。考虑到模拟精度和计算时间，最终选用网格数为11 200。

图2 不同网格数下轴向颗粒体积分数分布Fig.2 Distribution of axial particle volume fraction under different grid numbers

3.2 鼓泡流化床稳定时域分析

在流态化初期，鼓泡流化床内流动状态不稳定，不利于流态化分析，因此需排除初始流态的影响。图3为采用Gidaspow模型和新型曳力模型获得的瞬时压差Δp1曲线。其中，瞬时压差Δp1为距床层底部0.03 m与距床层底部0.3 m压力监测线的压力差。由图3可知，在3 s后不同模型的瞬时压差基本稳定波动，表明流化床内颗粒已达到稳定流态化。其中，压差波动是由气泡穿过床层时发生分裂、碰撞和合并造成的。

图3 采用新型曳力模型和Gidaspow模型得到的瞬时压差曲线Fig.3 Instantaneous pressure difference obtained by new drag model and Gidaspow model

3.3 床层时均压差

采用Gidaspow模型和带有不同修正因子的新型曳力模型模拟床层时均压差Δp，并将其与实验值进行比较，结果见图4。其中，床层时均压差Δp为距床层底部0.03 m与0.6 m处压力监测线的压差。当表观气速较小（Ug＜1.5Umf）时，模拟值与实验值吻合度较差。其原因是在模拟过程中表观气速过小，颗粒没有被充分流化，气固两相间作用力以摩擦力为主，导致沿气体运行方向气体压差过大。在实验过程中，颗粒并非均匀分布，如果有颗粒堵塞在进气喷口附近，则进气口压力降低，导致床层时均压差实验值较小。当表观气速较大时，床层时均压差模拟值与实验值的吻合程度很高。表观气速为0.38 m／s时，修正因子为0.9的新型曳力模型的床层时均压差模拟值与实验值最接近，其误差约为1.8%；当表观速度达到0.46 m／s时，修正因子为0.8的新型曳力模型床层时均压差的精度较Gidaspow模型提高了约3.1%；新型曳力模型模拟得到的床层时均压差始终优于Gidaspow模型，其中修正因子为0.8和0.9的新型曳力模型预测所得床层时均压差精度相似，且优于修正因子为0.7的新型曳力模型。

图4 不同表观气速下床层时均压差模拟值与实验值的对比Fig.4 Comparison between simulated and experimental values of average pressure difference at different apparent gas velocities

3.4 床层膨胀率

采用床层膨胀率E［19］作为衡量新型曳力模型准确性的依据。

式中：H′为床层膨胀高度，m；H′0为静止床层高度，m。

不同表观气速下床层膨胀率模拟值与实验值的对比见图5。由图5可知，随着表观气速的增大，不同曳力模型下床层膨胀率模拟值均逐渐提高。表观气速为0.1 m／s时，Gidaspow模型和带有不同修正因子的新型曳力模型的床层膨胀率模拟值相差较小；表观气速为0.38 m／s时，修正因子为0.9的新型曳力模型的床层膨胀率模拟值相比Gidaspow模型提高了4.7%，在其他2种修正因子下新型曳力模型的床层膨胀率模拟值略大于Gidaspow模型；当表观气速为0.51 m／s时，修正因子为0.9的新型曳力模型的床层膨胀率模拟值相比Gidaspow模型提高了7.1%，修正因子为0.7和0.8时则分别提高1.8%和1.1%。随着表观气速的不断增大，新型曳力模型的床层膨胀率模拟精度不断提升，其中修正因子为0.9的新型曳力模型预测所得床层膨胀率精度最高，修正因子为0.7和0.8的新型曳力模型的预测值接近。当φs为0.1～0.4时颗粒聚团严重，经过介尺度修正的Mckeen模型相比Gidaspow模型模拟精度更高。

图5 不同表观气速下床层膨胀率模拟值与实验值的对比Fig.5 Comparison between simulated and experimental values of bed expansion rate at different apparent gas velocities

3.5 颗粒体积分数

图6给出了Ug=0.38 m／s时瞬时颗粒体积分数云图。由图6可知，不同曳力模型模拟所得瞬时颗粒体积分数与实验具有相似性。从图6（a）可以看出，实验中床层底部附近有小气泡，当气泡上升到顶部表面发生合并时，气泡变大。气泡的生长是由壁面效应以及与其他气泡的相互作用导致的。曳力模型模拟与实验观测到的气泡大小基本一致，但仍存在误差。这是因为在实验中气体分布器的设计对其附近的射流穿透和流体流动特性有一定影响，而模拟过程中设置为均匀进气，未考虑气体分布器对气泡的影响。

图6 瞬时颗粒体积分数云图Fig.6 Instantaneous particle volume fraction nephogram

从图6还可以看出，不同曳力模型下静态床层高度模拟值低于实验值。其原因如下：（1）实验时颗粒为非球形，导致实际曳力偏大；（2）TFM耗散相对更强，模拟的气固结构不易维持；（3）实验时颗粒直径并非定值，而模拟时取为定值。

图7给出了床层高度为0.2 m截面处不同表观气速下的时均颗粒体积分数变化。与Ug=0.38 m／s相比，当Ug=0.46 m／s时不同曳力模型的时均颗粒体积分数分布更加对称，但其与实验值之间的误差较大。当Ug=0.38 m／s时，孔隙度曲线出现轻微的不对称。这是由于床层中形成了一定的流型以及颗粒在床层上均匀分布的时间太短。当Ug=0.46 m／s时，在靠近床层中心处模拟和实验得到的时均颗粒体积分数曲线均较为平坦，说明表观气速越大，流动越均匀。由图7可知，新型曳力模型与Gidaspow模型得到的时均颗粒体积分数分布相似，其中修正因子为0.8的新型曳力模型的时均颗粒体积分数模拟值与实验值误差率最小，Gidaspow模型的误差率最大。当Ug=0.38 m／s时，Gidaspow模型的最大误差率约为21%；修正因子分别为0.8、0.7和0.9的新型曳力模型的最大误差率分别约为11.1%、16%和17.8%；当Ug=0.46 m／s时，Gidaspow模型的最大误差率约为33%，修正因子分别为0.8、0.7和0.9的新型曳力模型的最大误差率分别约为18.1%、19.1%和19.3%。综上，当Ug=0.38 m／s时，修正因子为0.8的新型曳力模型精度最高，比Gidaspow模型提高了约10%；当Ug=0.46 m／s时，修正因子为0.8的新型曳力模型精度最高，比Gidaspow模型提高了约15%。综合床层时均压差、床层膨胀率和颗粒体积分数，修正因子为0.8的新型曳力模型最适用于B组颗粒。

图7 不同表观气速下的时均颗粒体积分数变化Fig.7 Variation of time average particle volume fraction at different apparent gas velocities

4 结 论

（1）随着表观气速的增大，新型曳力模型得到的床层时均压差模拟值逐渐接近实验值。当表观气速增大至0.46 m／s时，修正因子为0.8的新型曳力模型精度比Gidaspow模型提高了3.1%。

（2）随着表观气速的增大，新型曳力模型得到的床层膨胀率模拟精度优于Gidaspow模型。表观气速为0.51 m／s时，修正因子为0.9的新型曳力模型的床层膨胀率模拟值与实验值之间的误差最小，相较于Gidaspow模型，其精度提高了7.1%。

（3）表观气速为0.38 m／s时，Gidaspow模型的时均颗粒体积分数最大误差率约为21%，而修正因子为0.8的新型曳力模型的最大误差率约为11.1%；当表观气速为0.46 m／s时，修正因子为0.8的新型曳力模型的时均颗粒体积分数分布曲线最优，其精度较Gidaspow模型提高了约15%。

（4）新型曳力模型的预测精度均高于Gidaspow模型。修正因子为0.8和0.9的新型曳力模型精度接近，且优于修正因子为0.7的新型曳力模型。修正因子为0.8的新型曳力模型最适用于B组颗粒。