化学链分级气化燃料反应器气固流动数值模拟1)

郭 月 崔 颖 宋 涛 卢 平 殷上轶 ,2)

*(南京师范大学能源与机械工程学院，南京 210023)

†(东南大学能源与环境学院，南京 210096)

生物质具有资源丰富、挥发分含量高等优点，常被选作气化的燃料。常规的气化技术存在合成气热值低等问题，针对此不足，国内外研究者们将化学链与气化相结合，由此形成了化学链气化技术[1]。化学链气化技术通过载氧体（金属氧化物或非金属氧化物）中的晶格氧为燃料提供氧来源，控制晶格氧和燃料的比值使生物质等固体燃料部分氧化，进而获得合成气。因此，其具有能够自热运行[2]、获得高品位的合成气[3]等优点。

在生物质气化过程中，不可避免会产生副产物焦油，若其凝结在下游设备上，会造成堵塞[4]，影响设备的安全运行。在化学链技术中，燃料反应器的选型以鼓泡流化床为主[5]，且载氧体颗粒选型为Geldart B类，而该种颗粒的典型流化特性是极易产生气泡，因此随着气泡在床层内产生、生长直至在床层表面附近破裂，极易造成气体逃逸。以上提及的燃料特性、载氧体选型、反应器设置，造成焦油大分子逃逸，导致合成气产率降低。为解决该问题，两级/多级气化的方式被提出，实验证明分级气化方式能够有效降低焦油含量[6]。目前多级燃料反应器设置有两级鼓泡床[7]、五级鼓泡床[8]等。但鼓泡流化床的缺点是压降较大，造成风机能耗高，尤其体现在较大规模的装置中。相比于密相流化床，喷动床气固接触良好，压降低[9]。因此为降低压降，本文提出一种新型化学链分级气化装置，其中，燃料反应器由下部的气化反应器（设计为鼓泡流化床）和上部的重整反应器（设计为喷动床）串联而成，中间由矩形狭缝连接，如图1所示。矩形狭缝的尺寸和结构影响着两床之间的颗粒连通量，进一步影响到燃料的气化效率、焦油的重整反应以及整个装置的载氧体循环利用。因此需要探究颗粒连通量受反应器结构、操作工况等的影响规律。

图1 化学链分级气化装置Fig.1 Chemical looping graded gasification unit

反应器中发生的复杂的多相流动不能完全由实验手段测得，而数值模拟由于低成本、获得的信息更丰富，能够帮助正确设计反应器结构和理解反应器中的气固流体动力学，近几十年来得到了较多的应用。常采用的模拟多相流的方法有以离散单元模型（discrete element method，DEM）为代表的欧拉–拉格朗日法和以双流体模型（twofluid model，TFM）为代表的欧拉–欧拉方法，主要特点是：计算精度高但受颗粒数目的影响[10]，能够模拟大型设备内的气固流动但计算精度受网格数目的限制[11]。由于流化床内的颗粒数量庞大，跟踪每个颗粒并不能保证在合理的计算时间内得到结果，因此近些年来出现了多相粒子单元法（multiphase particle-in-cell，MP-PIC）。该方法将具有相同物理性质的颗粒打包形成包裹，同时引入固体应力模型模拟颗粒碰撞模型，有效减小了计算量，在流化床的数值模拟研究中得到了广泛的应用[12]。

对本文提出的存在矩形狭缝的串级反应器结构，常规的模型选择、参数设置是否适用于小尺寸装置有待探究，并且上下两级反应器的颗粒连通量受操作参数和结构的影响规律的获得，有利于后续对反应器结构的优化和改进。

1 计算模型

1.1 数学模型

在MP-PIC中，气体动力学在欧拉框架下求解，气体湍流用大涡模拟方法描述；固相采用PIC法建模，颗粒运动、颗粒碰撞应力和气–颗粒反应分别用颗粒分布函数、空间梯度法描述。

气相的控制方程

式中，εg为单位网格空隙率，ρg为气相密度，ug为气相速度矢量，p为气相压力，τg为气相应力张量，F为气相和颗粒相间单位体积动量交换频率，Dp为曳力系数，up为颗粒相速度矢量，ρp为颗粒相密度，fmp为概率函数。

湍流模型采用大涡模拟[13]。

固相的控制方程如下。

粒子相的动力学是通过求解f(xp,mp,up,t)的输运方程得到

式中，fD为局部平衡状态的颗粒分布函数，τD为颗粒碰撞的松弛时间。

颗粒相动量守恒方程为

式中，εp为颗粒体积分数，τp为颗粒法向应力，u¯p为单位网格内平均颗粒速度矢量。

常用的曳力模型有Wenyu–Ergun模型和EMMS模型等。

Wenyu–Ergun模型

其中，DWenyu和DErgun表达式为

其中，θp是颗粒体积分数，θCP是计算粒子的颗粒体积分数，Re是雷诺数，rp为颗粒半径，θg是气相体积分数，ug是气相速度。

EMMS模型[14]

式中，μg为气体动力黏度，系数fe公式为

式中，ε为颗粒空隙率。

对于颗粒–颗粒碰撞，颗粒接触法向应力τp由欧拉网格上的空间压力梯度计算。Harris 等[15]提出的法向应力模型为

式中，Ps是单位为Pa且大于0的常数，本文取1；εb是堆积态时的颗粒浓度；γ是范围2～5之间的常数，本文取3；θ是数量级为10–7的常量。

对于气相，采用无滑移边界条件；对于离散相颗粒与壁面，采用不完全弹性碰撞

式中，ewn为颗粒–壁面法向碰撞恢复系数；ewτ为颗粒–壁面切向碰撞恢复系数。

1.2 计算模型

本文以分级的燃料反应器作为数值模拟的对象，如图2(a)所示。采用笛卡尔网格进行网格划分，划分的网格如图2(b)和图2(c)所示。初始条件：颗粒均匀堆积在反应器上下部，喷动床中的颗粒静止堆积高度为H0，颗粒浓度εp= εcp，本文取0.6；空气均匀填充其他计算区域。边界条件：将气化反应器底部作为均匀进气速度入口；将反应器顶部作为压力出口，设为一个大气压。此外，在距离气化反应器底部0.4 m处建立实时颗粒流量监测面，以获得从重整反应器掉落至气化反应器的颗粒流量MⅡ-Ⅰ。时间步长为0.1 ms，总的计算时长为20 s。

图2 计算模型示意图Fig.2 Schematic diagram of the calculation model

2 结果与讨论

2.1 模型参数对计算结果的影响

首先，采用不同的网格数量（58 740，33 696和4 050）对模型进行网格无关性验证的计算，得到的装置压力分布基本一致，且与实验数据[16]的最大误差为1.29%。综合考虑计算速度和准确性，采用网格数量为33 696进行后续计算。

2 .1.1 曳力模型

曳力是影响气固流动最重要的力之一，关系着反应器中的气固流动行为预测的准确性。Wenyu–Ergun和EMMS曳力模型对鼓泡流化床中的气固流动均有较好的预测性[17-18]，因此在本研究中采用这两种模型进行对比计算。图3显示了在表观气速Ug= 0.5 m/s时的瞬时颗粒体积分数εvf以及瞬时颗粒浓度εpc云图。

图3 曳力模型计算结果的对比Fig.3 Comparison of drag force models

两模型均能较好地显示床内气泡的生成、长大、破裂，但区别在于：在同样工况下，Wenyu–Ergun使得床内颗粒处于鼓泡流态化，而EMMS使得床内颗粒处于湍动流态化；并且，在同样的模拟时刻（t= 15 s），选择同样的模型参数，采用EMMS模型计算使得颗粒更易掉落。与实验拍摄图片进行对比，Wenyu–Ergun模型计算的结果更准确，因此后续的模拟采用Wenyu–Ergun模型进行计算。

2 .1.2 壁面边界条件

不同的壁面边界条件会对重整反应器掉落至气化反应器的颗粒连通量MⅡ-Ⅰ的计算结果产生较大的影响。在MP-PIC中，颗粒与壁面的相互作用主要体现在两个参数：法向和切向颗粒–壁面碰撞恢复系数ewn和ewτ，这两个参数表征着颗粒与壁面碰撞后的能量保留程度。ewn和ewτ的取值范围均在0～1之间，取值越接近1，颗粒和壁面碰撞后保留的能量越多。其中，当ewτ=0时，表示无滑移边界条件，当ewτ=1时，表示全滑移边界条件。

图4展示了颗粒–壁面法向碰撞恢复系数ewn对颗粒连通量MⅡ-Ⅰ和压力分布的影响（负号表明颗粒掉落的方向朝下，下同）。从图中可以看出，ewn对MⅡ-Ⅰ的影响较小，因此，对系统的压力分布影响也不明显。Yin等[19]在模拟循环流化床提升管时也得到了类似结论，这意味着该参数对循环流化床整个装置中的气固流动影响不大。将其与实验数据进行对比，选择ewn=0.3进行模拟计算。

图5展示了颗粒–壁面切向碰撞恢复系数ewτ对颗粒连通量MⅡ-Ⅰ（a）和压力分布（b）的影响。ewτ对MⅡ-Ⅰ影响较大，尤其当ewτ越接近全滑移边界条件时，其取值每变化0.1，颗粒越易掉落，气化反应器中的颗粒量增加，对气相阻力增加，因此压降增大，压力降低，计算结果产生了较为明显的差异。该影响规律与模拟提升管获得的结果差异较大[19-20]，这表明在化学链分级气化装置的模拟中，ewτ的取值需要进行择中处理，以确保固体循环通量的模拟准确性。与实验数据进行对比，本研究采用ewτ= 0.4进行后续的模拟计算。

图5 切向动量保持量Fig.5 Tangent-to-wall momentum retention

2.2 颗粒连通量

2 .2.1 表观气速的影响

图6展示了颗粒连通量MⅡ-Ⅰ受表观气速Ug的影响。表观气速Ug对颗粒连通量MⅡ-Ⅰ的影响较明显。随着Ug的增大，颗粒连通量MⅡ-Ⅰ逐渐减小。这是因为随着Ug的增加，气相对颗粒相的曳力增加，因此颗粒连通量MⅡ-Ⅰ减少。

图6 颗粒连通量受表观气速的影响Fig.6 Effect of gas velocity on particle flux

2 .2.2 矩形狭缝宽度的影响

图7展示了颗粒连通量MⅡ-Ⅰ受狭缝宽度z的影响，从图中可以发现，随着矩形狭缝宽度z的增大，颗粒连通量MⅡ-Ⅰ增加，且掉落的速度从慢变快。尤其当z= 0.025 m时，重整反应器中的颗粒在17 s左右已完全掉落至气化反应器中。当z= 0.02 m时，掉落速度与在z较小时不完全一致，其原因可能是：当狭缝宽度在某一临界值（15 mm ＜z＜ 20 mm）时，初始颗粒掉落主要受重力影响，速度较快；紧接着在气流的作用下颗粒开始形成非标准的喷动，掉落速度减缓。当z进一步增大至25 mm时，颗粒主要受重力作用，掉落速度较快。此外，相较于调节表观气速，改变狭缝大小对调节颗粒连通量更有效。

图7 颗粒连通量受狭缝宽度的影响Fig.7 Effect of slit width on particle flux

3 结论

本文以由气化反应器和重整反应器串联而成的燃料反应器作为计算对象，构建了三维模型，验证了曳力模型、模型参数（法向和切向颗粒–壁面碰撞恢复系数ewn和ewτ）对计算结果的敏感性，并与实验值进行对比，确定了合适的模型及模型参数，探究了表观气速、狭缝尺寸对颗粒连通量MⅡ-Ⅰ的影响规律。为类似存在渐缩渐括结构的流化床模型参数的选取提供了参考，并且为化学链气化循环装置的固体循环通量的控制提供了调节依据。获得的主要结论如下。

（1）Wenyu–Ergun曳力模型能够较好地模拟鼓泡床中的气泡，以及喷动床中的三区：喷射区、喷泉区、环隙区。

（2）法向颗粒–壁面碰撞恢复系数ewn对燃料反应器中的气固流动的影响不大，但切向颗粒–壁面碰撞恢复系数ewτ对气固流动的影响非常明显。确定的合适的模型参数为ewn= 0.3和ewτ= 0.4。

（3）随着表观气速Ug的增加，颗粒连通量MⅡ-Ⅰ逐渐降低。

（4）随着矩形狭缝尺寸z的增加，颗粒连通量MⅡ-Ⅰ逐渐增加。