太阳能热化学反应器内流动和反应的数值模拟

翟占虎，李 巍，刘国栋

(1.哈尔滨工业大学能源科学与工程学院，黑龙江 哈尔滨 150001；2.海洋石油工程股份有限公司，天津 300452)

太阳能因其储量丰富、清洁可再生等特性，成为当前研究的热点问题，但是由于其热流密度低、具有周期性且不易储存，目前还只是作为辅助能源.近年来太阳能聚焦技术取得长足发展，将光热转换和高温热化学反应相耦合的光化学利用技术引起重视，该技术以氧载体材料为依托，在高温下将 CO2、H2O 等热分解产生H2和CO等碳氢燃料.

围绕着太阳能反应器，国内外学者开展大量的研究工作.瑞士PSI太阳能技术实验室Steinfeld[1]教授领导的团队设计并试验的太阳能热化学反应器，能够将CH4-ZnO在高温下分解为Zn和合成气.Chueh等[2]设计的太阳能反应器，参考植物将太阳能转化为燃料的过程，研究了CeO2在热化学循环中生产燃料的热力学和动力学性质以及反应的稳定性，在不进行废热回收的前提下，该反应的热化学学循环效率理论上可达到16%～19%.Meier等[3]设计并实验的太阳能间接加热回转窑，可以在1 200 K～1 400 K温度区间内实现高温煅烧反应，并获得高纯度石灰.吴玺[4]针对ZnO与CH4-Ar的重整制氢反应，数值模拟研究了吸收腔内温度场分布情况，并分析了腔体温度、反应物料进口速度和CH4-Ar体积比对重整制氢反应效率的影响.王富强等[5]采用数值模拟对比研究了辐射模型多孔介质太阳能接收器传热性能的影响.黄兴等[6]研究了催化剂NiFe2O4的热解机理并应用于太阳能热化学反应器，采用Ansys Fluent研究了反应器内速度场和温度场，分析了颗粒直径和气体流速对反应器内部的温度分布影响.Lougou等[7]数值模拟研究了太阳能热化学反应器作用下的甲烷重整反应，分析了NiFe2O4催化作用下，混合气体流速、操作压力、催化剂装载量和辐射温度分布等参数对反应转化率的影响.张浩等[8]采用数值模拟结合反应动力学，研究了太阳能反应器中的甲烷重整反应和积碳特性，得到多孔介质的最佳空隙率和颗粒直径.Bader等[9]设计了太阳能高温热化学反应器，以二氧化铈(CeO2)为固定床催化剂，利用该金属氧化还原作用的分解CO2，该反应器中催化材料与吸收腔不直接接触，CeO2相对于Fe3O4和ZnO等催化材料，具有储放氧性能高、比表面积大、传热传质速率高和结构稳定等优点[10]，避免了催化剂的污染和分离回收，同时独立的反应器提升了系统稳定性和操作性，适用于存在一定温度波动的高温反应器.

太阳能热化学反应器反应效率受热流密度分布和金属氧化物催化剂的氧化还原过程的影响，实验直接获得相关的参数比较困难，因此本文数值模拟研究Bader所设计的太阳能热化学反应器，反应过程中氧空位浓度δ的变化通过求解编写的UDS(用户自定义标量方程)实现，质量和能量变化通过UDF(用户自定义函数)加载到多孔介质区域，最后分析氧化还原反应阶段的温度、氧空位以及产物浓度分布情况.

1 反应器结构与数学模型

1.1 反应器结构及物理模型

反应器整体结构如图1所示，由一石英窗、光圈、保温层、圆筒形吸收腔及腔内环形排列着的七个反应器构成，吸收腔内部长为347 mm，半径为152.5 mm，每个反应器由同轴圆管构成，反应器外管的外/内半径分别为34.9 mm/31.8 mm，内管的外/内半径分别为22.2 mm/19.1 mm，圆筒壁面厚为3.1 mm，两圆筒夹层厚为9.5 mm，夹层中放置718 g固定的金属氧化物催化剂二氧化铈(CeO2)，该催化剂为颗粒直径为5 mm，催化剂颗粒的孔隙率为70%，固定床层的空隙率为40%.

图1 反应器几何结构图

反应器工作流程如下，经过聚焦后的高通量太阳辐射能(3 kW)，透过光圈射在吸收腔壁面和七个圆管反应器的表面，吸收了太阳辐射能的氧化铝管壁将热量传入固定床层的二氧化铈，催化材料二氧化铈(CeO2)其表面及体相的晶格氧原子能够直接参与反应并被消耗[11]，同时形成氧空位(δ)，产出化学燃料.在反应器工作过程中，周期性循环地向反应管内通入N2和CO2，通入N2时长为100 s，质量流量为1.1 g·s-1，该阶段发生还原反应，出口产物有N2和O2；通入CO2时长为100 s，质量流量为1.7 g·s-1，该阶段发生氧化反应，出口产物有CO2和CO.反应方程式如下所示：

还原反应(Reduction)

(1)

氧化反应(Oxidation)

CeO2-δrd+(δrd-δox)CO2=CeO2-δox+(δrd-δox)CO

(2)

二氧化铈颗粒及其所构成的固定床，在模拟中被简化为多孔介质，根据相关公式计算出的输运性质如表1所示.

表1 二氧化铈和固定床层输运性质

为简化计算，取模型的七分之一进行模拟，首先模拟得到吸收腔内热流密度分布，然后将热流密度加载到单个反应器，再进行反应器内流动、传热和化学反应模拟.

1.2 控制方程

反应过程中满足连续性方程、动量方程、组分输运方程、能量守恒方程.

连续性方程

(1)

公式中：ρ、u、Ri分别为密度、速度和质量源项；下标f、s、i分别为流体、固体和组分i.

动量守恒方程

(2)

公式中：φ、p、μ、K、CE分别为空隙率、压力、动力粘度、渗透率和惯性阻力系数.

组分输运方程

(3)

公式中：φeff、Yi分别为有效空隙率和组分i的质量分数.

能量守恒方程

(4)

(5)

(6)

DO辐射方程：

(7)

由于反应器吸收腔内温度较高，需要考虑其内部的辐射特性，在模拟吸收腔内部热流密度分布时采用DO辐射模型.

1.3 化学反应动力学方程

目前的文献中，没有准确描述二氧化铈催化分解CO2的反应动力学表达式，本研究中采用的动力学表达式主要参考Bulfie和Rohini Bala的文章[14-15].在二氧化铈参与氧化还原反应的过程中，其非化学计量数δ决定反应进行速率的大小，以下为δ随时间变化的微分方程，以及各组分质量源表达式.

δ与t的关系：

(8)

(9)

在氧化与还原过程中，CeO2既能释放O2，也能吸收O2，因此氧气质量源项的表达式如下：

(10)

研究表明[17]，高温下CO2的直接分解作用也很明显，且为可逆反应，各组分质量源项的表达式如下：

(11)

公式中：i=O2，CO，CO2，相应地vi=0.5，1，-1.

O2的净产量为CO2分解所释放的量与CeO2氧化消耗的量之和，如下所示：

RO2，net=RO2，diss+RO2，ceria

(12)

预计热分解产物CO会在短时间内和O2重新结合成CO2，故本研究认为氧化阶段产生的O2，一部分被CeO2吸收来填补氧空位，其余的全部与CO重新结合.

本研究以上述的非化学计量数δ的微分方程、各组分的质量源项表达式为理论基础，来实现上述化学反应.

2 模拟结果分析

Bader[9]文中给出的太阳能热化学反应器为概念设计，并未进行实验测量，反应器在设计运行参数下，工作温度约为1 800 K，CO产率在预估为3.9×10-4mol·s-1，模拟得到的温度分布如图4和图11所示，以及出口CO摩尔产率如图17所示，较为接近，可认为模拟所采用模型是合理的.

2.1 吸收腔内热流密度分布

吸收腔内部热流密度分布情况，如图2所示.可以看出圆管反应器壁面的热流量明显高于吸收腔壁面，经计算圆管反应器壁面的平均热流密度为4 600 W·m-2，总热流量为318.5 W，吸收腔壁面的热流密度为181.5 W，热量损失占36.3%.沿着圆管反应器轴的热流密度分布曲线，如图3所示.沿着轴向先增高后减少，并且在x=0.12 m达到最大值6 900 W·m-2.

图2 反应器热流密度分布云图图3 沿X轴向热流密度

2.2 还原阶段分析

(1)反应器温度分析

t=1 s和100 s时反应器温度分布云图，如图4所示.从左侧图可看出反应器壁面温度最高为1 880 K，t=100 s时超过1 860 K温度的范围明显少于t=1 s时的，这是因为CeO2的还原反应是吸热反应，随着反应的经行，反应器温度下降.

图4 t=1 s和100 s时反应器温度分布

反应器壁面轴向温度在不同时刻的温度分布和多孔介质的体积平均温度随时间的变化，如图5、图6所示.壁面温度和体积平均温度都随着反应的进行，不断减小，进一步验证了CeO2的还原反应是吸热反应.

图5 反应器壁面轴向温度在不同时刻的温度分布图6 多孔介质体积平均温度随时间的变化

(2)氧空位浓度δ与O2产率

t=1 s时和t=100 s多孔介质区域的δ分布云图，如图7所示.在t=100 s时δ最高到0.041，且最高位置在出现在刚进入多孔介质区域的地方.这是因为气体刚进入多孔介质时，由于不含氧气，故此处氧气的分压为零，远低于平衡压力，CeO2中的氧原子在高温条件下能够以氧气的形式释放出来，故此处的氧空位最高.从图7中还可看出，靠近出口处多孔介质的δ最小，这是因为氧气的不断释放，氧气在下游的质量分数达到最高，即氧气分压也达到最高，故此处的CeO2中的氧原子受氧气分压的抑制，释放速率较慢.

图7 不同时刻δ在多孔介质区域的分布云图

沿轴向δ是不断减小的，且随着时间增加，δ的增长变缓，多孔介质上游区域的δ增长率明显高于下游，如图8所示.由图9可知，氧气的质量分数沿轴向不断变大，随之时间增加，氧气的质量分数不断降低.这是因为当CeO2刚开始反应时，由于多孔介质区域各处的氧气分压均为零，O2的产率很高，随着反应的进行，氧气分压升高，特别是多孔介质下游区域，其氧气分压明显高于上游，更接近平衡压力，当氧气分压越接近平衡压力，δ增长越慢，氧气的产率越小，相反地上游区域氧气分压始终较低，δ增长始终较快，所以x=0 m处的δ值明显高于x=0.374 m处的.

图8 轴向氧空位浓度δ分布(不同时刻)图9 轴向氧气质量分数分布(不同时刻)

出口氧气的摩尔流率，如图10所示.在刚通入N2的时候，出口氧气的摩尔流率可达到9.42×10-5mol·s-1，且随着反应时间的增加，氧气流率逐渐下降，这是因为总体上氧气的分压越来越接近平衡压力(大约为240 Pa)，δ增长减缓，从CeO2中释放出的氧气变少.

图10 出口氧气的摩尔流率

2.3 氧化阶段分析

(1)反应器内温度分布

t=101 s和200 s时反应器温度分布云图，如图11所示.在整个时间段内，反应器的最高温度为1 770 K，低于还原阶段的1 880 K，而从101 s到200 s这段时间，反应器内高温区域变大.这是因为通入CO2后，CO2质量流量比N2的质量流量大，且CO2气体的比热容大于N2的比热容，故在一定的热流密度下，氧化阶段的流场温度低于还原阶段时的温度.此外，在整个氧化阶段，反应是放热的，所以随着时间的增加，反应器整体升温.

图11 不同时刻壁面和计算域内温度分布云图

反应器壁面轴向温度在不同时刻的温度分布，如图12所示.从图12可知，随着轴向距离的增加，壁面温度先增高后下降.因为气体被壁面不断加热，随着流动距离温度变高，在临近出口时，由于与入口低温气体(1 625 K)换热，所以温度出现了下降，由于反应放热，所以随着时间增加，壁面温度升高.从图13可知，随着时间的增加，多孔介质体积平均温度不断增加，进一步说明了氧化阶段发生的反应是放热的.

图12 反应器壁面轴向温度分布图13 多孔介质体积平均温度随时间的变化

(2)氧空位浓度δ与CO产率

从图14可看出，随着轴向距离增加，氧空位浓度δ逐渐降低，且随着时间的增加，δ也在减少，特别是x=0 m处，δ下降最快，在t=200 s时多孔介质区域δ的分布云图，如图15所示.可见在x=0 m处，δ约为0.03，体积平均δ约为0.020，与还原阶段的初始值较接近，这是因为氧化阶段，CO2高温下发生热分解，释放出大量O2，使得此时的氧气分压大于平衡压力，CeO2将O2吸收来填补自身的氧空位，所以δ出现下降趋势，且δ越大，氧气平衡压力就越小，CeO2对O2的吸收速率更高，因此在x=0 m处δ下降最快.

图14 轴向氧空位浓度δ分布情况(不同时刻)图15 t=200s时δ分布云图

从图16可看出，随着x的增加，CO质量分数先增加后减少，且随着时间的增加，CO质量分数不断下降，这是因为当通入CO2时，CO2热分解释放的O2被二氧化铈CeO2吸收用以填补氧空位，进而促进CO2的热分解，所以从x=0 m到x=0.15 m处，CO质量分数不断升高；而CO质量分数在x=0.15 m之后出现下降趋势，是因为这段区域的氧空位浓度δ较小，对O2的吸收能力有限，导致CO2热分解产生的O2没有被完全吸收，与CO再结合，导致CO质量分数下降.

CO摩尔流率随时间的变化曲线，如图17所示.随时间增加，CO摩尔流率逐渐下降，这是因为随着反应的进行，氧空位浓度δ逐渐减少，对CO2热分解反应的促进表少，导致CO产率逐渐减少.

图16 轴向CO质量分数分布 图17 出口CO的摩尔流率

3 结 论

本文以Bader设计的反应器为基础，建立了用于二氧化铈小温差下氧化还原循环分解CO2的太阳能高温反应器的三维瞬态模型，考虑了流动过程和非化学计量反应的相互作用及其对温度、压力、组分浓度和燃料产率的影响.研究发现，反应器热流密度集中分布在靠近光圈的一侧，热量损失取决于吸收腔壁面温度和材料保温性能，应尽量减少辐射能量直接作用于吸收腔壁面；反应器温度与反应阶段密切相关，还原阶段温度上升，氧化阶段温度下降；还原阶段，二氧化铈释放氧气，氧空位浓度δ受氧气的抑制作用，在气体入口处最大，出口处最小；氧化阶段，二氧化碳在高温下热分解，二氧化铈吸收分解后的氧原子，促进了CO2的热分解作用；氧气释放量和二氧化碳热分解作用与氧空位浓度、气体流率、反应时间和温度密切相关，在太阳能反应器设计阶段需综合考虑上述因素.