滑动弧裂解CO2机理

杨 锐, 游滨川, 刘 潇, 杨家龙, 郑洪涛, 马 彪,2

(1.哈尔滨工程大学动力与能源工程学院, 黑龙江哈尔滨 150001; 2.哈尔滨广瀚燃气轮机有限公司, 黑龙江哈尔滨 150001)

引 言

近年来, 随着国家节能减排的开展和能源资源的紧缺, 低热值燃料燃烧技术在工业界与学术界显得越来越重要，其热值低于6.28 MJ/(N·m3).比较常见的低热值燃料有高炉煤气、炼油伴生气和瓦斯气等[1].高炉煤气在低热值燃料中产量非常大, 我国每年有20%的高炉煤气被浪费掉[2], 总量相当于9×106吨标煤[3].但是由于高炉煤气的热值较低和含H量较少, 会造成燃烧不稳定和易熄火问题.为了有效解决上述问题, 有必要引入等离子体, 通过其产生的高能电子与燃料发生碰撞, 触发一系列碰撞激发、解离等反应, 产生小分子燃料, 自由基和离子等活性粒子, 从而加速燃烧的链式反应[4].高炉煤气中含有16%～19%的不可燃组分CO2, 而CO2通过等离子体裂解可以产生活性粒子和可燃物质, 从而改善燃烧性能.

常压下等离子体中最有效和最有应用前景之一的是滑动弧等离子体[5-8], 简称滑动弧.在常压下滑动弧非稳态放电并产生非平衡等离子体.由于它能产生高密度的高活性组分, 因此在裂解CO2应用中引起了研究者们的注意.

为了改善滑动弧的裂解性能, 各国学者通过使用电学测量[9]、光谱测量[10-11]和高速摄影[12]等实验方法研究了滑动弧的物理和化学特性.除了滑动弧特性实验研究之外, CO2裂解的深层机理研究同样至关重要, 各国学者也尝试多种方式对滑动弧裂解进行深入研究.

Kustova等[13-14]对CO2的振动动力学进行了研究, 但是没有考虑CO2裂解过程中等离子化学反应的影响.Rusanov等[15]结合模拟计算和实验, 对等离子体裂解CO2进行了研究, 得到的结果与实验吻合很好.但是并没有考虑电子的演化.Hokazono等[16]从微观的角度对等离子体裂解CO2进行了研究, 建立CO2裂解理论模型, 在密闭容器内, 在不同功率条件下对中性粒子和振动激发态粒子的数密度进行了探究, 但是忽略了对流对CO2裂解的影响.比利时的Aerts等[17]和大连理工的吴帆[18]开发了零维CO2裂解化学反应动力学模型, 但是同样忽略了对流的影响.Pietanza等[19]计算了裂解CO2过程中的电子能量分布函数, 但是对其中涉及的化学反应并没有研究, 因此并不是很全面.

综上所述, 常压下滑动弧裂解CO2应用前景广泛.由于实验研究滑动弧裂解CO2对实验设备精度的要求较为严苛.因此本文通过仿真软件建立滑动弧裂解CO2的反应机理模型，对滑动弧裂解过程中的电子数密度、电子温度、CO2裂解机制进行分析.等离子化学的复杂性导致计算量巨大.首先建立一维准滑动弧模型, 不考虑流场, 通过引入特征冷却频率[20], 以此来计算由横向气流对流引起的等离子活性粒子的损失.此外还考虑详细的CO2等离子体化学机理模型.本文对电子数密度和电子温度进行了计算, 并与实验进行了对比.同时分析了滑动弧裂解CO2的机制, 并对比不同特征冷却频率下CO2转化率的大小.

1 一维滑动弧模型

1.1 一维滑动弧几何模型

如图1所示, 左侧为截取的一段准圆柱滑动电弧,Vg为气体速度,Va为电弧速度,r为准圆柱滑动弧的径向；右侧为准圆柱电弧沿轴向z的截面α.只考虑准圆柱电弧放电区域的径向, 然后将滑动弧简单地描述为轴对称圆柱形几何体中的放电通道.Wang等[20]使用类似的一维滑动弧模型, 得到的结果与实验拟合较好, 证明了这种简化方式的准确性.本文使用的一维几何模型的长度rmax为3 mm.除此以外, 由于滑动弧的特性, 需要考虑对流引起的等离子组分损失, 下面引入特征冷却频率描述这种损失.

图1 一维模拟区域

1.2 对流冷却处理方法

Fridman等[6]研究表明, 相对于气体流动, 滑动弧的滑动存在明显滞后现象.气体速度与电弧速度之差随着气体速度增大而变大, 随着电弧长度变大而减小[21].很明显, 气体与电弧的相对速度Vc=Vg-Va对滑动弧放电行为有很大的影响.相对速度会造成滑动弧放电通道内的等离子组分损失.本文使用Wang等[20]的对流冷却处理方式, 即引入特征冷却频率描述对流引起的损失.描述方程如下

dφ=φFdt

式中, dφ为变量φ因对流损失的值; dt表示时间;F为特征冷却频率, 由以下方程计算得到

F=Vc/Lc

式中,Vc为气体对电弧的相对速度;Lc为特征长度, 由以下方程计算得到

Lc=πrmax/2

由1.1节知rmax=3 mm, 则特征长度Lc=1.5π(mm).相对速度由理论计算和高速相机实验[21-22]获得, 分别为3.9, 7.9, 15.7 m·s-1, 则对应的特征冷却频率为0.8, 1.8, 3.3 kHz.

1.3 数学模型

1.3.1 粒子平衡方程

电子通量如下

中性粒子通量仅由扩散决定

式中,qs为粒子的电荷量,μs为s粒子的迁移率,Ds为s粒子的扩散系数,E为电场.其中电子的迁移率和扩散系数由BOLSIG+计算得到.CO2+的迁移率由文献[23]获得, 扩散系数由Einstein关系计算得到.中性粒子的扩散系数由COMSOL软件自带的数据库得到.

1.3.2 平均电子能量平衡方程

平均电子能量由如下方程计算

式中,με,e为电子能量迁移率,Dε,e为电子能量扩散系数, 由以下方程计算得到

με,e=(5/3)μe

1.3.3 Poisson方程

放电的电场由Poisson方程计算得到

Δ(ε0εrE)=ρq

式中，ε0为真空介电常数;εr为相对介电常数;ρq为电荷密度.

1.4 化学机理

滑动弧裂解CO2的基本过程见图2.电子碰撞CO2分子发生一系列碰撞、解离等反应, 产生O等活性粒子, 并促进链式反应的发生.

图2 滑动弧裂解CO2示意图

根据文献[18]的计算结果, CO2的振动激发态多数退激发为基态CO2, 对CO2生成的贡献率并不大, 本文为节省计算成本, 忽略CO2的振动激发态.本文创建的化学机理模型共包括7种粒子.

表1 模型中包含的粒子

本文模型包括的反应如表2, 3所示.主要包括电子碰撞反应、中性粒子间的反应, 共包括19个反应.

表2 模型中电子碰撞反应

电子碰撞反应的速率常数通过求解电子能量分布函数(electron energy distribution function, EEDF)得到, 相应的方程为

式中,e为电子所带电荷量,m为电子质量,ε为电子能量,σi为i种粒子的碰撞截面,f(ε)为EEDF.

表3中的M表示反应体系中的中性粒子,二体反应速率常数单位是cm-3/s, 三体反应的速率常数单位是cm-6/s.

表3 模型中中性粒子反应

2 结果与分析

本文进行了准稳态分析和瞬态分析, 准稳态分析使用直流25 mA的电流源, 瞬态分析使用交流I=35sin(2·π·50·t)mA的电流源.准稳态分析主要探究特征冷却频率对滑动弧放电特性的影响以及主要粒子和反应随电弧位置的演化.瞬态分析主要探究特征冷却频率对滑动弧放电特性的影响以及上述变量随时间的演化, 并对不同特征冷却频率下的CO2转化率进行研究; 对于瞬态分析, 将在电弧中心对电流半个周期进行研究, 即时间范围为0～10 ms.下文中, 如未做特殊说明, 特征冷却频率默认为0.8 kHz.

2.1 准稳态分析

2.1.1 滑动弧放电特性

在等离子体放电特性中, 电子密度和电子温度很大程度上决定电子碰撞反应速率的快慢.因此本文滑动弧放电特性将从电子数密度ne和电子温度Te两个方面进行分析.

图3(a)为不同特征冷却频率下电子数密度随电弧位置的演化, 图3(b)为不同特征冷却频率下电子温度随电弧位置的演化.电子数密度和电子温度都是在r=0处达到最大, 其中当特征冷却频率为0.8 kHz时, 电子数密度最大值为3.77×1018m-3, 电子温度最大为1.46 eV, 随后沿r逐渐减小, 证明放电是在电弧中心处发生的,即r=0处.一些实验研究表明滑动弧放电的电子数密度在1017～1020m-3之间, 电子温度在1～3 eV 之间[10,31-33], 证明了本文准稳态模拟的准确性.随着对流的增强, 放电产生的电子被向后吹散, 造成了图3(a)呈现的曲线规律, 即在电弧中心处, 随着特征冷却频率的增加，电子数密度减小, 而沿着r增大, 特征冷却频率和电子数密度成正比关系.除此以外, 对流也会造成能量的损失, 使电子从电场中获得的能量减小, 因而随着特征冷却频率的增加, 电子能量减小, 如图3(b)所示.

(a)Electron number density

2.1.2 主要反应随电弧位置的演化

本小节主要分析CO2消耗反应、CO2生成反应、O2消耗反应及O2生成反应的反应速率随电弧位置的变化.

如图4所示, 由于在电弧中心处发生放电, 电子从电场中获得足够高的能量, 电子温度最高, 此时电子碰撞解离反应R7:e+CO2→e+CO+O的反应速率最快, 对CO2裂解的贡献最大.但是随着r增大, 电子温度降低, R7的反应速率迅速降低.活性粒子引发的链式反应中, C+CO2→CO+CO的反应速率最快, 是对CO2裂解起主要贡献的中性粒子反应.不同反应对CO2裂解的贡献通过对反应速率沿r进行积分计算得到, 不同反应的贡献率如表4所示.

(a)Reaction rates of R7 and R11

由表4可知e+CO2→e+CO+O对CO2裂解的贡献率最大, 为90.63%, C+CO2→CO+CO的贡献率为9.26%, 其余的基元反应的贡献率均小于1%.由此可见, 提高CO2转化率的一个方法是提高电场能量, 使用更高的放电电流, 以此增大反应R7的反应速率.

表4 各反应对CO2裂解贡献率

如图5所示, 反应CO+O+M→CO2+M对CO2生成的贡献最大, CO+O2→CO2+O的反应速率较低.由此可见, 提高CO2转化率的一个方法是消耗CO和O, 以此抑制反应CO+O+M→CO2+M.

图5 CO2生成反应速率随位置的演化

如图6所示, R19的反应速率远大于R16的反应速率, 即O2的生成反应在整个计算域占据主导地位, 从而促进活性粒子O2的产生.

图6 O2生成和消耗反应速率随位置的演化

2.1.3 主要粒子随电弧位置的演化

本节将探究CO2、活性助燃粒子随电弧位置的演化.

如图7所示, CO2数密度沿r逐渐增大, 原因是放电发生在r=0处, 此时是反应R7起主要作用, CO2被裂解, 而随着r增大, CO2生成反应R13占据主导地位, 导致CO2数密度增大.

图7 CO2数密度随位置的演化

如图8所示, CO与O的数密度在电弧中心处最大, 随着r增大而减小.如图4, 5所示, 随着r增大, CO与O重新结合成CO2的反应R13占据主导地位, CO和O逐渐消耗, 从而使CO与O数密度随着r增大而减小以及图4中CO2数密度的升高.而如图6所示, O2的生成速率远大于O2的消耗速率, 因而O2数密度增大.CO数密度最大达到1.5×1024m-3, O最大达到1.4×1024m-3, O2最大达到1.46×1021m-3, 由此可知滑动弧裂解CO2会产生大量助燃的活性粒子, 裂解生成可燃的CO，从而提高低热值燃料中的可燃成分浓度, 进而改善低热值气的燃烧性能.

图8 活性粒子数密度随位置的演化

2.2 瞬态分析

2.2.1 滑动弧放电特性

瞬态分析的滑动弧放电特性同样从电子数密度和电子温度两个方面进行分析, 图9为电子数密度和电子温度随时间的演化.

如图9所示, 电子数密度和电子温度在10-4s时开始增大, 证明放电发生在10-4s, 电子开始从电场中获得能量.电子数密度最大达到4.316×1017m-3, 电子温度最大达到1.43 eV, 据2.1.1节所述, 瞬态模拟的电子数密度和电子温度符合实验得到的滑动弧特性, 证明了瞬态模拟的准确性.

图9 电子数密度和电子温度随时间的演化

如图10所示, 对流会将电子和能量吹离电弧中心, 对流越强, 吹离效果越明显, 电子从电场中获得的能量减少, 从而使电子数密度和电子温度随着对流的增加而减小.

图10 不同特征冷却频率下的最大电子数密度和电子温度

2.2.2 主要反应随时间的演化

本小节主要分析CO2消耗反应、CO2生成反应、O2消耗反应及O2生成反应的反应速率随电弧位置的变化.

如图11所示, 10-4s时开始放电, 电子从电场中获得了足够高的能量, 电子碰撞解离反应R7:e+CO2→e+CO+O的反应速率逐渐增加, 与其他基元反应相比反应速率最快, 对CO2裂解的贡献最大.活性粒子引发的链式反应中, R15:C+CO2→CO+CO的反应速率最快, 是对CO2裂解起主要贡献的中性粒子反应.不同反应对CO2裂解的贡献通过对反应速率沿时间进行积分计算得到[18], 不同反应的贡献率如表5所示.

(a)Reaction rates of R7 and R11

由表5可知e+CO2→e+CO+O对CO2裂解的贡献率最大, 为98.43%, C+CO2→CO+CO的贡献率为1.57%, 其余的基元反应的贡献率均小于0.01%.

表5 各反应对CO2裂解贡献率

如图12所示, 反应CO+O+M→CO2+M对CO2生成的贡献最大, CO+O2→CO2+O的反应速率较低.与准稳态模拟得到的规律相同.

图12 CO2生成反应速率随时间的演化

如图13所示, R19的反应速率在0.001 s后快速上升, 远大于R16的反应速率, 从而促进活性粒子O2的产生.

图13 O2生成和消耗反应速率随位置的演化

2.2.3 主要粒子随时间的演化

如图14所示, 在10-4s附近开始放电, CO2发生裂解, CO2数密度下降.由2.2.1节可知电子数密度和电子温度随着特征冷却频率的增加而减小, 因此电子碰撞解离反应R7的反应速率也将随着特征冷却频率的增加而减小, 而由2.2.2节可知, 电子碰撞解离反应R7在整个CO2裂解反应体系中贡献率最大, 因此CO2最终的裂解量随着特征冷却频率的增加而减少.为了更好地衡量CO2转化的程度, 采用变量CO2的转化率Φ, 计算方程为

图14 不同特征冷却频率下CO2数密度随时间的演化

式中,n1为CO2初始数密度，n2为反应体系稳定后CO2数密度.不同特征冷却频率下CO2转化率如图15所示.

图15 不同特征冷却频率下CO2转化率

CO2转化率随着特征冷却频率的增加而减小, 0.8 kHz下的转化率为11.8%, 在实际应用中应选择合适的气体流量, 以获得较高的CO2转化率.

如图16所示, CO和O的数密度曲线几乎重合, 这也说明了在整个CO2裂解反应体系中, 电子碰撞解离反应R7占据主导地位.由图13可知O2生成反应速率在0.001 s后开始上升, 并超过O2消耗反应, O2数密度增大.由图16可知, 滑动弧裂解CO2会产生大量活性粒子和可燃的CO, 会加速燃烧过程中的链式反应, 并从一定程度上增大燃料浓度, 进而改善低热值气的燃烧性能.

图16 活性粒子数密度随时间的演化

3 结论

本文建立了一维滑动弧模型和CO2裂解化学反应动力学模型, 进行准稳态和瞬态CO2裂解反应过程研究, 并分析主要反应途径和各反应的贡献率, 得到以下结论:

(1)滑动弧裂解CO2会产生活性粒子O和O2以及可燃的CO气体.滑动弧裂解CO2会提高低热值燃料可燃成分浓度, 进而改善低热值气的燃烧性能.

(2)CO2裂解反应体系中, 电子碰撞解离反应e+CO2→e+CO+O的贡献率最大,准稳态中贡献率为90.63%, 瞬态中贡献率为98.43%.而CO2生成反应体系中, 反应CO+O+M→CO2+M占据主导地位.因此为提高CO2转化率, 可使用较高的放电电流, 增强电场, 使电子获得更高的能量, 增大e+CO2→e+CO+O的反应速率, 与此同时, 及时转走或消耗CO和O, 抑制CO+O+M→CO2+M的进行.

(3)对流会造成等离子体组分和能量的损失, 将其吹离放电通道, 从而使电子数密度和电子温度降低, 使裂解CO2的主要反应电子碰撞解离反应e+CO2→e+CO+O反应速率下降, 引起CO2转化率下降.因此在实际应用中, 要选择合适的气体流量, 防止气体流量过大引起对流速度增大而CO2转化率减小.

由于计算资源有限, 本文进行了一维模拟, 做出了很多简化, 例如默认计算域的速度是恒定值, 但是在实际应用中反应器内的速度是随着位置和时间变化的, 因此未来的研究可以尝试建立三维模型, 以获得更准确的结果.