流化床生物质气化过程的CFD–DEM模拟研究1)

2022-08-19孔大力樊建人

力学与实践 2022年4期

孔大力王帅罗坤樊建人

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

全球变暖和气候变化促使人们减少温室气体排放，其潜在途径是减少化石燃料的消耗，发展替代性的可再生能源。生物质由于挥发分含量高、储存量广、环境友好等优点，适合进行气化处理[1]。生物质气化主要在鼓泡流化床(bubbling fluidized bed, BFB)反应器中进行，其性能高度依赖反应器内的流体动力学。因此，深入了解炉内现象（如多相流、传热传质、化学反应）对于鼓泡流化床反应器的设计和优化具有重要意义。然而，复杂的热化学过程（如颗粒混合和碰撞、气泡形成和破碎、相间耦合以及化学反应）使得通过实验手段研究鼓泡床反应器内的生物质气化过程变得极为困难。

计算流体动力学(computational fluid dynamics, CFD)提供了一种极富前景的替代方案以描述复杂的稠密气固反应流。现有的CFD方法一般分为两类：欧拉—欧拉方法和欧拉—拉格朗日方法。前者将气体和颗粒均视为连续介质，在欧拉框架下进行求解。该方法能以较低的计算成本准确预测大型设备中的气固特性，并在过去几十年中被广泛用于模拟稠密气固反应系统[2-3]。然而，欧拉—欧拉方法无法解析颗粒的离散特性（如旋转、碰撞、耗散、缩核等），这使得人们无法获得稠密气固流中的微观尺度信息[4]。欧拉—拉格朗日方法可以很好地解决这一难题。作为此框架下最典型的一种方法，计算流体力学—离散单元法(CFD—DEM)将气相和颗粒相分别视为连续相和离散相，追踪每个颗粒的轨迹和其离散属性（如粒径、温度、组分、密度等）[5]。因此，该方法可以提供颗粒尺度的详细信息，例如颗粒间/相间的相互作用、传热特性和反应特性。

近年来，CFD—DEM方法逐渐被采用探索流化床反应器中的生物质气化特性。Yang等[6]发现单侧燃料入口会导致产物在反应器内的非对称分布，提高温度、蒸气/生物质比和当量比会增强颗粒受力和颗粒耗散性能。Wang等[7]观察到生物质颗粒在气化过程中倾向于迁移并聚集在床层表面，异相反应主要发生在密相区，而均相反应通常发生在稀相区。Ku等[8]比较了生物质和烘焙生物质之间的气化性能，并表示烘焙生物质需要更长的转化过程。Wang等[5]测试了接触和阻力模型对生物质气化过程热化学结果的敏感性，证明了使用线性弹簧阻尼接触模型的计算效率是使用非线性赫兹接触模型的6.8倍。综合来看，上述研究只关注了流化床反应器内的宏观尺度信息（如流动模式、合成气产量、压降、碳转化率和气体组分），而鲜有报道微观尺度的颗粒行为（如颗粒混合、传热、传质），这极大地限制了对生物质气化过程中气固流动特性和热化学特性的理解。

因此，本文发展了一种考虑颗粒粒径多分散性和热化学效应的CFD—DEM模型，在颗粒尺度上研究BFB反应器中生物质气化过程中的气固流体动力学和传热传质行为。本文的结构如下：第1节详细介绍了数学模型；第2节给出了模型设置；第3.1节从传热和热化学性质方面进行了模型验证；第3.2节展示了气固流动特性；第3.3节分析了生物质颗粒中不同传热模式的贡献以及不同操作参数下对四种传热模式的影响；第3.4节给出了气体产物浓度和反应速率的分布，以及它们之间的对应关系。

1 数学模型

1.1 控制方程

气相的质量、动量、能量和组分守恒方程分别为

式中，εg，ρg，pg，ug，Tg分别为气相的空隙率、密度、压力、速度、温度；δm˙p表示由异相反应引起的气体种类的体积消耗或生成速率；τg为气相应力张量；Cpg和κg分别为气体比热容和热导率；Qgp为气体和颗粒间的对流传热速率； ΔHrg为反应热；γRg代表气相的辐射传热系数；TRg为气相辐射温度；t为时间；下标g代表气相，m代表固相；Xn和Dn分别表示第n种气体的质量分数和扩散系数；Rgn是由均相反应引起的质量源项。Igm为相间动量交换源项，表达式为

式中，Np为网格c中颗粒总数，fd,i是气相对单个粒子施加的力，Vc和Vi分别为网格c和颗粒i的体积，xi为颗粒位置。气相作用力包含了压力梯度力和曳力，∇Pg为压力梯度力，up为颗粒速度。β为曳力系数，采用文献[9]计算。

固相的平动、转动方程为

式中，mp,i，vp,i，ωi和Ii分别为颗粒i的质量、平动速度、转动速度和转动惯量。L为接触点到颗粒质心的距离，n为法向单位向量。颗粒受到重力mp,ig，流体作用力fd,i和碰撞接触力fc,i的作用。k为颗粒i周围检索范围内的颗粒总数，fct,ij表示颗粒i和颗粒j的切向碰撞力。其中，碰撞接触力采用基于软球碰撞框架下的LSD（linear spring-dashpot）模型进行计算，其能在保持计算精度的前提下大幅提升计算速度[5]。

1.2 曳力模型

曳力在稠密气固流动的模拟中起着重要作用。在生物质气化过程中，生物质颗粒的性质（如颗粒质量、粒径等）会随着化学反应的进行而发生变化。化学反应引起的颗粒质量变化导致粒径分布存在较大差异。在多分散系统中，基于均匀流动系统开发的传统曳力模型[10-11]往往低估了大颗粒的曳力，高估了小颗粒的曳力[12]。因此，本工作采用了考虑多分散性效应的曳力模型[10,12]

式中，c为系统中颗粒的总数，dp,i为颗粒i的直径，μg为气相黏度，Rep为颗粒雷诺数，Ni表示与颗粒i具有相同大小的粒子数。

1.3 传热模型

颗粒i的能量守恒方程为

其中，Cp,i是颗粒的比热系数，Tp,i为颗粒温度。颗粒通过以下几种机制进行传热：颗粒—气体对流传热qgp,i，颗粒—颗粒传导传热qpp,i，颗粒—流体—颗粒传导传热qpfp,i，辐射传热qrad,i和化学反应热ΔHrs。各种传热机制的计算公式如下。

其中，Pr为普朗特数，κp为颗粒的热导率，Rc,ij为颗粒i和颗粒j接触区域的半径，Rp为颗粒半径。Rin和Rout为对流传热区的上下边界，Tenv为环境温度，Tg,Ω为在Ω区域内的气相温度，Np,Ω为在Ω区域内的颗粒个数。详细推导过程参见文献[13]。

1.4 反应模型

颗粒由于热解和焦炭转化等化学反应而发生质量损失。颗粒i的质量变化表示为

其中，Rsn,i为生物质颗粒中各个组分的质量变化，Ns为生物质颗粒中的组分个数。

采用粒径转化模型描述颗粒质量损失过程中的粒径变化，在反应过程中，颗粒密度不变，而粒径随着颗粒质量的变化而减小，表示为[14-15]

其中，ρp,i为生物质颗粒密度。

生物质进入流化床后热解生成固定碳、可燃性气体和焦油，随后焦油迅速热解。由于水蒸气的催化重整及高温状态作用加强了焦油热解反应程度，使得在生物质热解释放的大部分焦油充分裂解，生成可燃性气体。生物质的热解过程表示为

式中，系数αi根据生物质的元素分析和工业分析进行确定。由于生物质中的N元素和S元素含量非常低，此处不予考虑。采用基于单步的阿雷尼乌斯定律计算热解反应速率，表示为[16]

其中，mvolatiles为生物质中挥发分的质量，R为理想气体常数，Tp表示颗粒温度。

热解后，生物质颗粒中残留的焦炭与H2O和CO2反应生成CO和H2。除了与煤焦气化有关的异相反应外，反应器中还有均相反应，包括水煤气变换反应和甲烷气化反应[17-20]。焦炭和H2O及CO2发生异相反应（R1和R2）

其中，R1和R2的反应速率可以表示为

其中，对于R1和R2，指前因子（An）分别为45.6 s/(m·K)和 8.3 s/(m·K)，活化能（En）均为4.37×107J/kmol。rkin,n和rdiff,n分别为动力学速率和扩散速率。Cn为常数5.0×10—12s/K0.75。Ap为颗粒的表面积，Pgn为气体分压，Tg为气体温度。

气体间会发生均相反应

其中，R3和R4的反应速率可以表示为

式中，[gas]表示气体浓度，R4为可逆反应。

1.5 数值格式

在每个流体时间步，根据颗粒初始信息（如位置、速度、加速度、温度以及组分）计算气固相互作用源项，从而求解流体守恒方程。得到每个网格内的流场信息，然后根据颗粒时间步长通过一阶显式积分获取更新后的颗粒信息。再将这些颗粒信息反馈给网格进行下一轮的迭代求解。流体速度和压力耦合采用SIMPLE算法进行处理。流体相时间步长为1×10—5s，颗粒相时间步长为最小碰撞时长的1/50[21]。

2 模拟工况

模拟对象为一个BFB，如图1所示，其宽度（X）、厚度（Z）、高度（Y）分别为1 500 mm，230 mm，4.5 mm。在X，Y，Z三个方向划分网格为46，200，1，平均网格尺寸为4.5 ～ 7.5 mm。对于速度，入口设置为均匀流量入口边界条件，出口设置为零速度梯度边界条件，壁面设置为无滑移边界条件；对于压强，入口和壁面设置为零压力梯度边界条件，出口设置为恒定压力边界条件（101 325 Pa）；对于温度，入口和出口设置为零温度梯度边界条件，壁面设置为等温边界条件。床料和生物质粒径均假设为1.5 mm，则网格和粒径比值符合CFD—DEM计算的要求[22]。初始时刻，床料颗粒堆积在BFB下部，堆积高度为350 mm，床料颗粒数为27 045个。生物质颗粒以恒定质量流速从底部进入反应器。模拟工况及气固参数列于表1。生物质的工业分析和元素分析列于表2中。

图1 BFB几何尺寸（单位：mm）Fig.1 Geometric dimensions of BFB (unit: mm)

表1 模拟参数Table 1 Operating parameters

表2 生物质的元素分析及工业分析Table 2 Proximate analysis and ultimate analysis of biomass

3 结果与讨论

3.1 模型验证

在流化床内稠密气固流动系统的数值模拟中，CFD—DEM方法已经获得了广泛的验证[23]。本文针对流化床内颗粒冷却过程和流化床内生物质气化过程进行了验证，检验了传热和化学反应模块的可靠性。

传热验证参考由Patil等[24]开展的颗粒冷却实验。如图2所示，图中，dp表示颗粒粒径，m表示颗粒总质量，Ubg表示进气速度，exp为实验值，sim为模拟值。总质量为125 g，粒径为1.0 mm的颗粒堆积于宽、厚、高为80 mm，15 mm，250 mm的床体下部。床体和颗粒的初始温度为90℃。将20℃的气体以1.2 m/s和1.54 m/s的速度从底部引入床体冷却颗粒。具体的气固物性参数参见文献[25]。模拟所得的平均颗粒温度和实验数据吻合良好，表明模型可以用来准确描述流化床内颗粒传热过程。

图2 平均颗粒温度对比Fig.2 Comparison of average particle temperature

在对流体力学和传热进行验证后，对化学反应模型进行了验证。图3显示了BFB出口处气体组分瞬时质量流率的时间演化曲线。图中Tbed表示操作温度，S/B表示蒸气—生物质比。当流态化气体从底部分配器进入BFB时，N2首先从出口逸出，然后是H2O。随着化学反应的进行，其他四种气体（即H2，CH4，CO，CO2）在约0.6 s时逸出出口。在约5 s时，所有气体组分的质量流量随时间趋于稳定，表明系统达到动态平衡状态。因此，采用5 s后的数据作为统计样本，计算时间平均结果，以避免启动过程的影响。

如图3所示，H2和CO的预测结果与实验数据[26]吻合良好，而CO2和CH4的预测结果略有差异。在生物质气化过程中，存在数百个均相和异相反应（例如干燥、热解和气化），在模型中实现此类复杂反应是不切实际的。因此，化学反应被简化为几个整体反应。此外，反应动力学在化学反应中起着关键作用，根据不同的文献[27-29]，化学反应速率具有不同的值。数值结果与实验数据之间的差异是不可避免的，但差异在可接受的范围内。因此，本模型能够可靠地预测BFB反应器中生物质气化过程的热化学性质。

图3 出口气体组分流量随时间变化及模型验证Fig.3 Time-evolution of outlet gas flow rate and model validation

3.2 气固流态

图4显示了反应器内气固流态达稳定后随时间的演变。如图所示，小气泡在反应器底部不断生成，在气泡上升过程中，气泡不断合并，颗粒被气泡裹挟，在反应器内做混沌运动。在床层表面，气泡破裂，大量颗粒被抛洒至稀相区，随后沿着壁面下落。气泡动力学（如形成、上升、合并、喷发）引起粒子的混沌运动，加剧了气体—颗粒/颗粒—颗粒的混合，显著增强了传热过程，决定了床中的流动状态和流化质量。

图4 不同时刻颗粒运动及气泡演变（EP_G：空隙率， Usy：颗粒垂直方向速度）Fig.4 Time-evolution of motions of particles and bubbles (EP_G: porosity, Usy: particle vertical velocity)

3.3 传热占比

良好的混合和传热性能是BFB反应器的显著优势。具体而言，传热影响化学反应速率，并影响最终气体产物的质量。因此，研究BFB反应器内的传热机理具有重要意义。有四种传热模式，包括传导、对流、辐射和反应热。目前，对生物质气化过程中BFB中传热机理的探索鲜有报道。

图5给出了生物质颗粒四种传热模式占比的时间演化曲线（conductive, convective, radia-tive, reactive 分别表示传导，对流，辐射，反应传热模式）。在5～10 s内，对反应器中所有生物质颗粒的传热速率进行平均。在初始阶段，传导、对流和辐射产生的热量出现峰值，随着时间的推进逐渐降低，直至达到动态平衡状态，而反应热的启动滞后于其他三种传热模式。这表明低温生物质颗粒进入反应器后经历快速加热过程，达到所需活化能后开始化学反应。在达到平衡状态后，传导、对流、辐射和反应传热模式下的颗粒平均传热速率分别为0.013 W，0.106 W，0.05 W和0.041 W，这四种传热模式分别占总传热速率的6.2%，50.48%，23.81%和19.52%。因此，在生物质颗粒的传热过程中，对流传热占主导地位，其次是辐射传热和反应传热。传导传热所占比例最小。在BFB反应器中，导热发生在颗粒碰撞过程中。颗粒间的碰撞是随机的、不连续的，这导致系统中的所有颗粒不能同时发生传导传热模式，从而传导传热模式产生的总热量是所有传热模式中最小的。

图5 不同传热模式下生物质颗粒传热速率随时间变化曲线Fig.5 Time-evolution of heat transfer rate of biomass with different heat transfer modes

图6显示了不同操作温度和生物质/水蒸气比对不同传热模式贡献占比的影响。操作温度的升高导致生物质颗粒和气体/沙子颗粒之间的温差增大，从而使得传热过程更为剧烈。较高的温度促进吸热反应并增加反应消耗的热量。增加S/B会增加四种传热模式产生的热量。较大的S/B会增加H2O的浓度并促进气化反应（即R1，R3和R4），因此在反应过程中消耗更多的热量。随着S/B的增加，生物质温度降低，生物质颗粒与操作温度之间的温差增大，从而导致传导、对流和辐射传热模式产生的热量增加。

图6 不同工况下生物质颗粒不同传热机制占比Fig.6 Heat transfer contribution of different heat transfer modes under different operating parameters

3.4 反应特性

图7和图8给出了反应达稳态后反应器内时间平均反应速率和气体产物的分布。图中*_rates和x_*分别表示各个反应的反应速率以及各个气体组分的浓度。生物质颗粒进入反应器后，在高温的密相区迅速发生热解反应，对应于高热解速率（见图7(a)）。此外，密相区中部的热解速率远高于近壁区。在稀相区，热解速率几乎为零，表明生物质热解过程主要在密相区域完成。异相反应速率的分布与均相反应速率的分布有很大不同。异相反应（如R1，R2）主要发生在反应器的密相区。近壁区较高的反应速率表明此处的焦炭浓度较高。均相反应速率的分布与气体产物组分（如CO，CH4）的分布一致。床面附近存在较高的均相反应速率，表明均相反应主要发生在该区域，这与前人的研究结果一致[13]。高气体浓度对应于高均相反应速率。CH4，CO和H2主要由热解过程产生，这三种气体的分布大致相同。CO2主要分布在反应器上部，其分布趋势与H2O的分布趋势相反，说明CO2主要来源于水煤气转换反应。

图7 反应速率的时间平均空间分布Fig.7 Time-averaged spatial distribution of reaction rate

图8 气体产物组分的时间平均空间分布Fig.8 Time-averaged spatial distribution of gas product components

图9给出了不同操作参数下BFB反应器出口处的气体产物占比。随着工作温度的升高，H2和CO浓度升高，CO2和CH4浓度降低。由于生物质气化过程中的吸热反应，提高操作温度可提高反应速率。具体而言，对于可逆反应（即R4），增加操作温度可促进正反应。在较高温度下，R1，R3和R4反应增加H2浓度而减少CH4浓度。高温促进焦炭转化会导致CO2消耗量和CO产量的增加。随着S/B的增加，BFB中H2O的浓度增加，促进了气化反应（R1和R3）以及R4向正向进行，导致CO和CH4的消耗增加以及CO2和H2的产量增加。