刘浩然, 李 阳, 王昌建

(合肥工业大学 土木与水利工程学院,安徽 合肥 230009)

生物质热解是生物质利用的重要方法,是生物质燃烧和气化过程中的重要步骤[1]。作为在厌氧或缺氧环境下将生物质转化为二次燃料的常见热过程,生物质热解是一个极其复杂的过程,包括一系列均相和非均相反应,并且伴随着质量、动量和能量的传递。因此从数值模拟的角度研究生物质热解的过程十分必要。

在生物质热解数值模拟的发展过程中会遇到许多挑战,例如反应物的尺寸、生物质的物理性质和化学成分、时间尺度和简化的动力学参数。此外,在之前的生物质热解数值模拟中,研究者通常假设热解的可燃产物可以作为单一组分,例如CO[2]。文献[3]强调,当建模不准确时,可能会由于提供不足的输入而破坏原本成功的详细模型。为了应对这些挑战,基于OpenFOAM平台开发了一个求解器,命名为biopyrolysisFOAM,用于模拟生物质热解过程及其子过程。在该求解器中,实现并整合了多孔介质内部与周围的瞬态流动、多孔介质周围和内部气体与固体之间的传热传质、生物质及其热化学性质的灵活定义、反应热直接定义或由生成的焓导出均匀和非均匀反应、可修改的热解和气化动力学机制。

本文选取松针作为典型生物质,开展热重实验。通过传统的热解动力学理论计算方法得到相关的热解动力学参数。随后利用优化算法对其进行优化,最后将优化后的热解动力学参数输入求解器中,将数值模拟计算结果与相关的实验数据进行对比验证。

1 实验部分

1.1 实验原料与仪器

本文实验中所使用的松针样品均来自于四川省眉山市。热解实验升温范围为300～1 300 K,升温速率分别为10、30、50 K/min,实验气氛为氮气。所有样品都被磨碎至80目(0.20 mm),然后放置在烘干箱中,在100 ℃的环境中进行24 h的烘干,目的是尽可能地去除样品中的水分,以防对热解实验结果的干扰。热解实验使用德国耐驰STA449F3型热分析仪,如图1所示。

1.2 热解动力学计算方法

对于相关的热解动力学参数计算选择Flynn-Wall-Ozawa(FWO)和Kissinger-Akahira-Sunose(KAS)这2种model-free方法以及分布活化能模型(distributed activation energy model,DAEM)和shuffled complex evolution(SCE)这2种model-fitting方法。在热重分析中,通常采用model-free和model-fitting相结合的方法来求解热解动力学参数。本文中FWO、KAS、DAEM 3种方法不需要确定具体的反应机理就可以计算参数,从而避免了反应机理的不恰当选择而引起的错误。SCE方法用于优化反应方案的动力学参数,可以得到更加准确和普适的参数。

(1) FWO法。计算公式为:

(1)

其中:β为升温速率;A为指前因子;E为活化能;R为理想气体常数;G(α)为反应机理函数;α为转化率;T为温度。根据(1)式,ln[AE/(RG(α))]独立于β,通过lnβ和1/T做一条直线,根据斜率-1.052E/R可以计算得到相应的活化能值。

(2) KAS法。计算公式为:

(2)

通过ln(β/T2)和1/T做一条直线,根据斜率-E/R可以计算得到相应的活化能值。

(3) DAEM法。计算公式为:

(3)

由(3)式可知,ln(β/T2)与1/T的关系给出了直线方程,ln(AR/E)+0.607 5可以提供方程的截距用来计算指前因子的值。

(4) SCE法。当使用复杂的反应模型时,大量未知的动力学参数会使在数值模拟的建模过程变得困难甚至不切实际。为了解决这一问题,一些新的优化方法被开发出来用于优化估算动力学参数。目前,主要应用到热重分析中的全局优化算法有遗传算法(genetic algorithm,GA)、SCE、simplex method 3种。文献[4]通过比较发现SCE全局优化算法的优化结果更加准确,整个优化过程更加快速。因此使用SCE全局优化算法进行相关的研究。

众所周知,生物质中3种主要的组分分别为半纤维素、纤维素和木质素,在前人研究的基础上,生物质热解模型可以视作是多组分的单步平行反应。随后,生物质热解后会产生焦炭和挥发分。松针热解的反应模型为:

半纤维素→vh焦炭+(1-vh)挥发分

(4)

纤维素→vc焦炭+(1-vc)挥发分

(5)

木质素→vl焦炭+(1-vl)挥发分

(6)

其中,v为热解产生的焦炭比例,下标h、c、l分别为生物质中的半纤维素、纤维素、木质素。

1.3 biopyrolysisFoam求解器的控制方程

生物质热转化的三维数学模型是建立在气相和固相守恒定律上的。固体是各向同性或各向异性的多孔介质。气相由不同来源的混合气体组成,主要包括最初填充反应器的气体、供应过程的气体、在生物质干燥过程中蒸发的水蒸气以及在固体热解过程中产生的气体。潜在的液相(水和液态碳氢化合物)被认为是蒸汽。这个假设对于大多数过程都是可以满足的。其他假设如下:在当前的设置中,忽略了多孔介质的运动;水分嵌入在多孔结构中,并且从固体中释放出来的水分会立即蒸发;气体在多孔介质中根据达西定律流动;从计算域的边界到固体表面存在辐射转移;化学反应和热过程可能导致热不平衡;只有固相存在时,整个结构域的空隙率γ为0,没有固相存在时,孔隙率为1;多孔介质在均质和非均质过程中发生质量损失。

气相的控制方程包括动量守恒方程、连续性方程、质量守恒方程和能量守恒方程。具体表达式如下:

(7)

(8)

(9)

(10)

Γ+SGrad

(11)

α∑(TG-TS)+Hγ+SSrad

(12)

1.4 网格无关性分析

以升温速率为50 K/min的热重实验为例,采用3种不同的网格精度对其进行模拟计算,网格数量分别为16 000、8 000、1 000个,计算域的尺寸均相同,计算的结果如图2所示。由图2可知,网格数为16 000、8 000个的模拟结果比较接近,而网格数为1 000个的模拟结果相差较大。因此考虑到模拟的精确度和计算的时间,选择网格数为8 000个来进行后续的计算。

图2 不同网格精度下质量损失的模拟结果对比

2 结果与讨论

2.1 松针热解动力学参数计算

根据热重实验数据可以得到在3个不同升温速率下松针的质量损失速率曲线,如图3所示。由图3可知,松针从450 K左右,开始发生热分解,直至850 K。从质量损失速率曲线中可以观察到2个峰。在不同转化率下,第1个峰对应的温度分别为594、600、626 K;第2个峰对应的温度分别为627、640、655 K。这2个峰分别对应半纤维素和纤维素的热解。

图3 松针的质量损失速率曲线

随后,选择转化率在0.05～0.90之间对松针热解活化能和指前因子的值进行计算,具体数值见表1、表2所列。

根据FWO法计算的斜率-1.052E/R和KAS法计算的斜率-E/R,在递进的转化率[5]下得到了不同的活化能值。由表1可知,FWO法和KAS法计算得到的活化能范围分别为63.90～205.23 kJ/mol和59.56～203.85 kJ/mol。2种方法计算的结果十分接近,且这2种方法的活化能分布基本相同,误差小于2%,这是由于不同方法使用不同的近似值造成的[5]。

表1 FWO法和KAS法计算的松针活化能

表2 DAEM法计算的松针活化能和指前因子

随着转化率的增加,活化能在转化率为0.05～0.40范围内逐渐增大,然后在转化率为0.45～0.70之间保持相对稳定,最终在转化率为0.75处再次增大。值得注意的是,不同的活化能值与反应速率和反应起始的难度有关[6]。因为不同转化率下活化能的变动可以归因于不同的反应,即样品中3种主要成分的反应,所以可以将松针热解分为3个阶段。

2.2 松针热解动力学参数优化

对于松针中3种组分的比例已经有学者开展相关的实验研究[7-9]。根据数据,选取最大值和最小值之和的1/2作为初始值。忽略除3种组分之外的其他组分,这是由于它们在初始的总质量中占据非常小的比例。于是可以得到3种组分的初始值分别为Yc,0=0.34,Yh,0=0.19,Yl,0=0.47。对于3种组分活化能和指前因子的初始值可以设为2.1节中计算得到的3个阶段的平均值。对于焦炭的初始值,可以参考文献[10]的研究。设置半纤维素热解生成的焦炭初始值vh=0.23,纤维素热解生成的焦炭初始值vc=0.06,木质素热解生成的焦炭初始值vl=0.46。所有参数的优化区间设为上述初始值的50%～150%。

对于反应级数来说,多数研究者将其定义为1。文献[11]对生物质废弃物的热解展开研究。使用1级级数反应模型和n级级数反应模型来预测样品的热重曲线。在其研究中发现,当反应级数在0～3之间,模型得到的结果与实验数据较为吻合。超过这个范围,模拟得到的热重曲线会出现不切实际的形状。因此,生物质3种组分的反应级数优化区间为0～3,初始值为1.5。

对这3组实验数据同时进行优化,将2.1节中提到的所有需要优化参数的初始值和优化区间列出,最终的优化结果见表3所列。随后根据优化后得到的热解动力学参数,预测松针在不同升温速率下的热重曲线,并与实验数据进行对比。

表3 松针SCE全局优化算法设置

SCE优化参数预测结果与实验结果对比如图4所示。

图4 SCE优化参数预测结果与实验结果对比

由图4可知,将所有升温速率下优化得到的热解动力学参数预测得到的热重曲线与实验数据高度吻合。肩状峰和尖状峰区域都可以被很好地捕获到,这说明通过SCE全局优化算法优化得到的参数具有可靠性。通过SCE全局优化算法还可以得到半纤维素、纤维素、木质素和焦炭的热重曲线,同样绘制在图4中。由图4可知,半纤维素和纤维素的热解峰分别对应于整体热重曲线的肩状峰和尖状峰区域,正好对应前人对生物质热解中峰值对应反应的解释。

2.3 数值模拟与实验结果对照

将数值模拟的结果与实验数据对比分析,结果如图5所示,在3个升温速率下的模拟结果与实验结果吻合度较高,但是存在一定的差异。主要体现在热解反应后剩余物的量以及热解反应发生的主要温度范围。造成这些差异的原因可能是不够完善和过于简化的热解动力学以及材料本身成分的变化导致的,具体的原因不能得到确定。值得注意的是,在模拟结果中,所有反应在大约800 K时停止,此时应该只剩下焦炭。但是在实际的实验结果中,样品还继续发生反应,可能是在较高的温度下仍然有少量的焦炭发生挥发,这并不影响已经选择的热解动力学。

图5 质量损失的实验结果与模拟结果比较

综上所述,基于3组分平行反应热解模型建立的biopyrolysisFOAM求解器,对热重实验的模拟与实验结果取得了较好的一致性。这充分证明了通过SCE全局优化算法优化得到的热解动力学参数的可靠性以及biopyrolysisFOAM求解器应用在生物质热解上的可行性。

3 结 论

本文通过热重实验对典型生物质松针的热解动力学展开研究,主要结论如下:

(1) 松针的热解温度范围为450～850 K,从质量损失速率曲线中可以观察到2个特征峰,分别对应半纤维素和纤维素的热解。

(2) 通过FWO法和KAS法计算得到不同转化率下的活化能分别为63.90～205.23 kJ/mol和59.56～203.85 kJ/mol。DAEM法可以计算相应的指前因子。

(3) 利用SCE法优化得到的热解动力学参数预测得到的热重曲线与实验数据高度吻合。

(4) 将优化后的热解动力学参数作为初始参数,输入biopyrolysisFOAM求解器中,对热重曲线进行模拟计算。结果显示,模拟结果与实验数据吻合度较高。充分证明了该求解器应用在生物质热解上的可行性。