多联产鼓泡流化床内气固流动特性数值模拟研究

2022-08-31骆仲泱方梦祥王勤辉岑建孟

热力发电 2022年8期

王昊，骆仲泱，方梦祥，王勤辉，岑建孟

（浙江大学能源清洁利用国家重点实验室，浙江杭州 310027）

与常规煤粉炉或层燃炉等相比，流化床炉气化方式存在显著差别。其根本原因是基本流动、反应原理上的差别。流化床是一种气固两相流态化原理的技术，炉内存在大量的物料颗粒，这些颗粒不仅构成十分复杂且粒径分布十分广泛，从几微米到数十毫米不等，因此，形成了流化床独特的空气动力特性。正是由于其不同的空气动力特性，使得流化床在传热、燃烧、污染物排放、灰渣特性等许多方面，具有不同的特点和规律[1-3]。由于试验手段本身的限制，靠实验手段难以完全掌握流化床内的气固流动特性，尤其是流化床中颗粒的浓度以及运动分布规律（这些都是流化床的重要参数）。随着计算技术的快速进步以及人们对多相流颗粒模拟理论认知的进一步加深，使用计算流体力学（CFD）数值计算技术逐渐成为研究流化床内气固流动特性的一种重要的辅助手段[4-11]。通过数值计算，可以快速获得速度场、浓度场、温度场的数据；对数据分析处理后可以研究流化床内的气固流动规律。

图1是由浙江大学提出的循环流化床多联产系统[12-14]。该系统由蒸汽和再生煤气鼓泡流化床热解器及循环流化床气化炉组成。其工艺流程为：煤在热解炉内首先与燃烧炉返料器过来的高温循环热载体灰混合发生裂解反应，产生半焦和挥发分；析出的挥发分在旋风分离器中与细灰颗粒分离；经分离后的煤气进入煤气净化系统，依次经过洗涤剂塔、电捕焦油器；收集的焦油可以出售或提取高附加值化工产品，或者经催化加氢工艺制备高品位合成油；净化后的煤气部分通过再循环风机加压后送入热解器底部作为流化风，其余部分则进入装置进一步净化变成民用煤气或者作为合成气用于生产液体产品；热解炉的半焦和低温循环热载体灰进入燃烧炉；半焦燃烧的热能一部分用于把循环灰加热后再送回气化炉以提供气化炉所需热量，另一部分热能可以用于发电和供热。

图1 浙江大学1 MW循环流化床热电气多联产实验装置Fig.1 The 1 MW circulating fluidized bed heat-electricity-gas poly-generation test apparatus in Zhejiang University

MFiX软件是由美国联邦能源技术实验室（NETL）开发的气固两相流开源计算软件，采用Fortran语言编写，用于流固系统的传热、化学反应以及流体动力学计算。该软件包括双流体模型TFM和DEM模型，主要解决流化床中的颗粒运动和燃烧问题[15-23]。该软件收录了大量关于气固两相流研究的计算模型，同时又预留接口，方便研究人员对源程序进行二次开发。

Wang等人[24]研究了单固体颗粒相时二维、三维鼓泡流化床内的气固流动特性，得到了流型、气体速度、颗粒速度等数据，与实验结构吻合较好。Wu等人[25]研究了不同曳力模型下的流化床燃烧过程，考察了不同的曳力模型对反应速率的影响。Li等人[16]采用了开源软件MFiX对中试循环流化床系统进行了模拟，得到数据与实验结果相吻合，但也发现了采用多固体颗粒相时会使计算量大大增加的情况。多固体颗粒相的多相流多颗粒相模拟求解困难是MFiX软件存在的问题。本文基于MFiX软件，采用Fortran语言二次开发将其拓展至多固体颗粒相系统模拟，再采用了OpenMPI并行技术，最终实现了流化床三维多颗粒相的热解模拟。

1 数学模型

1.1 连续性方程

式中：α为气固相质量分数；ρ为相密度；u为速度；t为时间；下标i、j为坐标方向；下标g和s分别表示气相和颗粒相；下标gs表示从气相到颗粒相，下标sg表示从颗粒相到气相。

1.2 动量方程

纳维-斯托克斯（Navier-Stokes，N-S）方程是瞬时速度的动量守恒方程，在数值计算的双流体模型中，通过速度时均化，得到了N-S方程的时均动量方程：

式中：Igs=-Isg，为气体和固体颗粒相间的曳力；τ为应力张量，主要由垂直于面的弹性力和平行于面的剪切力构成。

1.3 能量方程

流动过程中通常包含热量的生成、传递和耗散，气体和颗粒相的能量方程可用同一形式表达：

式中：m、n为第m或n相；Hm为物质的焓，可以通过定压比热容计算化合物的焓值；λm为气相和颗粒相的导热系数，气相的导热系数可以通过不同组分导热系统和含量推算而出；Qmn为气体和不同颗粒相之间对流换的换热量热量；γ为辐射系数；Smn为化学反应产生的源相。

1.4 组分输运方程

本文主要模拟的是煤的热解和燃烧，它们都是多种化学反应的集合。气相内部，气相和颗粒相间都存在着组分的运输，均由多组分构成。气相包括CH4、H2、H2O、CO、H2O、CO2、N2和焦油气8种，固相为固定碳、挥发分、水分和灰分4种。各组分满足质量守恒方程：

式中：X为组分的质量分数；R为组分输运方程的源相，表示均相或异相反应的生成或消耗速度；J为组分扩散流量。

1.5 颗粒温度方程

式中：J为与碰撞有关的颗粒温度的消耗；κs为颗粒温度的扩散系数；ΠΘ为气固相之间的相互作用引起的颗粒相颗粒温度的耗散。

1.6 化学反应模型

化学反应模型采用单元平均法来计算。具体方法为：通过使用网格单元中流体相的平均性质来计算均相和非均相的反应，同时忽略气体扩散速率对系统空间内反应的影响，即认为反应速率完全由反应动力学控制。

为简化计算，煤热解归纳总结为3个异相反应和5个均相反应。异相反应方程由煤的热重曲线拟合得到，适用于500～1 200 ℃。气体均相反应采用文献归纳方程[13]。

1）挥发分析出反应：挥发分=0.235 8CO+0.109 5CO2+ 0.473 5CH4+ 0.806 9H2+1.229 3H2O+0.092 1Tar（Tar表示焦油）。

式中：R为颗粒相体积密度；Xs,V为颗粒相挥发分质量分数；VVM为中间变量。

2）半焦的氧化反应：Cs+0.5O2=CO。

式中：kfm为膜层阻力系数；p(O2)为氧分压；D(O2)为氧扩散系数；Shm为舍伍德数；Sc为施密特数；kam为灰层阻力系数；krm为表面反应的阻力系数。

3）水分析出反应：(H2O)s= H2O。

式中：Xs,m为颗粒相水分。

4）焦油裂解反应：Tar=0.18CO+2.88CH4+7.74C。

式中：C为对应物质的摩尔浓度；εg为气体体积分数。

5）焦油氧化反应：Tar+13.59O2=10.8CO2+5.76H2O。

7）H2燃烧反应：H2+0.5O2=H2O。

2 物理模型

在三维笛卡尔坐标系中对鼓泡流化床进行模拟，网格采用结构化网格，边界网格使用实际几何形状切割结构化网格，通过非正交项和内部插值法计算出边界处的速度节点位置。模拟网格如图2所示。采用20倍煤颗粒直径作为最大网格尺寸，逐渐减小进行网格无关性验证。网格尺寸为0.02 m×0.02 m×0.03 m时冷态轴向压降曲线趋于一致，此时总网格数为35×35×250。同时，近壁面网格采取细化处理。

图2 CFD模拟网格示意Fig.2 Schematic diagram of the CFD simulation mesh

恰当的边界条件和初始条件的设定是获得有效数值解的重要前提。本文模型中所有入口条件均设置为质量入口条件，速度垂直于入口界面；出口均设置为压力出口边界条件，压力为背景大气压力。气体壁面边界条件设定为无滑移边界条件，即气体壁面切向和法相的速度均为0；颗粒相壁面边界条件设定为自由滑移边界条件，即颗粒相壁面的法相速度为0，切向速度不受壁面力的影响，由动量方程求得。模拟采用的关键模型见表1。

表1 模拟关键模型Tab.1 Key models of the simulation

鼓泡流化床内初始床料高度为1.4 m，床料的体积分数为0.5，最大堆积体积分数为0.5，计算流体区域的气相和颗粒相速度均为0。计算时间步长为10-4s，每个步长的迭代直到收敛为止，但不超过20次，不收敛时自动减小时间步长再次迭代直到收敛，再过渡到下一步计算时间。在模拟时间5 s后，对动量和质量偏微分方程采用二阶差分方式计算，提高计算精度。模拟时间总共为30 s。

本文模拟了热解炉实验装置在500、600、700 ℃热解温度工况下的反应气固流动特性，共模拟了3组工况，每组工况参数见表2。

表2 鼓泡流化床煤热解模拟工况Tab.2 The simulation conditions for coal pyrolysis in bubbling fluidized bed

由于多联产实验中热解炉采用再循环热解气作为热解炉的流化风，本文中热解炉的模拟采用简化方案。首先采用惰性气体N2作为流化风，待模拟至炉膛出口气体组分分布稳定，即热解过程达到稳态时，采用此时炉膛出口CH4、H2、CO、CO2（除N2、H2O、焦油）4种热解气成分进行归一化处理后作为新的流化风成分。采用新的流化风模拟至流动稳定的状态来实现再循环热解气作为流化风的模拟过程。煤的工业分析和元素分析见表3。

表3 工业分析和元素分析 w/%Tab.3 Industrial and elemental analysis results of the coal

3 模拟结果及讨论

3.1 气体组分分布

图3—图5给出了热解温度为500、600、700 ℃时热解炉内煤气各组分分布。这是根据模拟计算所得到的气体组分分布，排除焦油和水蒸气组分后，对CO、H2、CO2、CH4采取体积归一化处理的结果，主要用于和实验得到的数据做进一步的对比。

图3 500 ℃时均截面气体组分分布Fig.3 The gas composition distribution on time averaged cross section at 500 ℃

图4 600 ℃时均截面气体组分分布Fig.4 The gas composition distribution on time averaged cross section at 600 ℃

图5 700 ℃时均截面气体组分分布Fig.5 The gas composition distribution on time averaged cross section at 700 ℃

从图3—图5可以看到，所有温度段的热解产物中的CO、H2、CO2在密相区的质量分数均高于稀相区，而CH4在密相区的质量分数低于稀相区。CH4不同于其他气体组分分布特性，主要是模拟中考虑了煤焦油二次裂解所致。同时，也可以看出反应器内煤焦油的二次裂解主要发生在稀相区内。煤焦油二次裂解的产物为CO、CH4、固定碳，其中CO的化学计量数较小，所以煤焦油二次裂解生成的CO较少，因而对稀相区CO的影响较少，所以CO和H2、CO2的组分分布规律具有一定的相似性。

从图3—图5还可以看到，CO、H2、CO2、CH4在密相区底部，尤其是靠近给料入口有一个明显的高质量分数区域。在该区域内，煤与床料混合均匀，反应剧烈；同时在该区域内也可以观察到焦油和热解水的生成，显示出煤热解过程的异相反应均发生在环-隙结构的底部。

随着热解温度的升高，反应速度加快，热解炉内气体组分达到稳态的速度随之加快。500 ℃热解时，气体组分的稳定区域在炉膛上部；600 ℃时，组分稳定面上升到了炉膛中部的渐扩段；而700 ℃时，所有的反应基本在炉膛下部的密相区内已经完成，此时热解产生的气体可以部分充当流化风的作用。这也解释了为何700 ℃时实验所需流化风量有明显的减少。

图6—图7为煤气中CH4、H2体积分数随热解温度的变化。由图6—图7可以看到，本文建立的多联产系统模型可以很好地预测煤气中CH4的体积分数且变化趋势也与试验一致。在热解温度为600、700 ℃时，对于煤气中H2体积分数的模拟数据与试验数据相关性好；500 ℃时H2体积分数的预测偏差较大，可能和挥发分反应方程在低温下的适用性有关。

图6 热解气中CH4的体积分数Fig.6 Volume fractions of CH4 in the pyrolysis gases

图7 热解气中H2的体积分数Fig.7 Volume fractions of H2 in the pyrolysis gases

将模拟获得的热解炉出口气体组分除去水蒸气后换算为标准状况下的数据，得到热解炉煤气中焦油成分的质量浓度。图8为多联产系统模拟热解炉出口煤气中焦油质量浓度与多联产试验的对比。由图8可以看出，建立的热解模型能较好地预测热解气中焦油的质量浓度。

图8 热解气中焦油的质量浓度Fig.8 Mass cocentrations of tar in the pyrolysis gases

3.2 流化床截面颗粒体积分数及速度分布

流化床不同截面（选定4个截面对应高度分别为0.2、0.6、1.0、1.5 m）颗粒体积分数及速度分布如图9—图10所示。

图9 颗粒时均体积分数Fig.9 The time averaged volume fraction of the particles

图10 颗粒时均速度Fig.10 The time averaged speed of the particles

由图9、图10可以观察到，热解流化床炉不同高度处径向截面上均出现了中心区域低颗粒体积分数和近壁面区域高颗粒体积分数的现象，即存在中心稀疏的环-核结构。此外还可以发现，环-核结构随流化床高度的增加也有变化。在密相区0.2 m和0.4 m高度处，高体积分数颗粒区域略有增大；0.6 m高度处，高体积分数颗粒区域的面积大幅减小。而在渐扩段上方1.5 m处，高体积分数颗粒区域面积又有较大增加，直至0.2 m处颗粒相已经极难观测到。这主要是由于渐扩段的存在，使一定质量的床料长时间地存在于渐扩段及其上部区域形成内部回流流动，同时阻止了下部密相区沿边界进入渐扩段的颗粒量。新的环-核结构的形成也是由于相同的原因。因此，由于渐扩段的存在流化床内形成了以渐扩段为界的相连的2个环-核结构。