大气压甲烷针-板放电等离子体中粒子密度和反应路径的数值模拟∗

赵曰峰1) 王超1)2) 王伟宗3) 李莉4) 孙昊2) 邵涛2) 潘杰1)

1)(山东师范大学物理与电子科学学院,济南 250014)

2)(中国科学院电工研究所,北京 100190)

3)(比利时安特卫普大学化学学院,安特卫普 2610,比利时)

4)(国家电网济南供电公司,济南 250012)

(2017年10月10日收到;2018年2月11日收到修改稿)

1 引 言

甲烷作为天然气(气田气、油田气与页岩气等)、沼气和矿井气的主要成分,不仅是制造氢气、炭黑、甲醛、乙炔和一氧化碳等物质的基本原料,而且是民用与工业加热的重要燃料,在供热取暖、交通运输、燃气发电、合成气制取和农作物施肥等领域得到了广泛的关注与应用[1,2].随着我国原油资源的日益紧缺,采用重油加氢技术处理高硫、高残炭、高金属的劣质渣油,将原油中占比较高的重油转化为轻质油,从而优化能源结构并降低污染物排放,已经成为国民经济和社会发展中迫切需要解决的关键科技问题[3].传统的重油加氢技术由于在高温高压条件下使用临氢环境和催化剂,因而存在能耗高、轻油收率低和催化剂易结焦失活等缺点.千伏级高压电源驱动的针-板放电能够在宽广的气压范围内产生高化学活性与非平衡性的甲烷放电等离子体,适宜开展放电等离子体甲烷转化的科学研究和工业应用[4−12].由甲烷针-板放电与重油加氢反应耦合形成的甲烷转化重油加氢技术,可以实现重油高效加氢并增产高附加值低碳烯烃,具有广阔的应用前景和重要的研究意义.

比利时安特卫普大学Bogaerts课题组是较早对放电等离子体甲烷转化进行系统研究的科研团队.2013年,该课题组整合实验诊断数据与数值模拟结果,分析甲烷氮气混合气体介质阻挡放电(dielectric barrier discharge,DBD)等离子体中氮气浓度对甲烷转化制氢的产额、机理和动力学过程的影响[13].2015年,该课题组与荷兰埃因霍温大学van Dijk等[3]开展合作,探讨甲烷-氧气和甲烷-二氧化碳混合气体DBD等离子体中甲烷转化产生合成气以及产生高分子碳氢化合物和含氧化合物的主导路径.2017年,美国德雷克塞尔大学Fridman等[1]使用纳秒脉冲DBD和大气压辉光放电(atmospheric pressure glow discharge,APGD)等离子体活化甲烷,将气态甲烷分子直接并入液态碳氢化合物燃料,从而实现了非平衡放电等离子体作用下的甲烷液化.上述研究表明,非平衡放电等离子体技术为甲烷转化提供了新的发展方向,但在已经公开发表的放电等离子体甲烷转化工作中大多采用DBD和APGD的放电形式,鲜有针对大气压甲烷针-板放电等离子体的实验诊断和数值模拟研究报道.

电子温度和粒子密度的时空演化特性与反应路径贡献是非平衡等离子体及其应用研究方向的长期学术热点[14−21].对大气压甲烷针-板放电等离子体而言,由于相应的原位实时诊断方法尚不完善并且所适用的粒子种类有限,数值模拟成为探索上述问题的必要研究方法.天津大学Wang等[22]使用零维全局模型,研究等离子体能量和放电气隙间距对大气压DBD微等离子体反应器中甲烷转化产生高分子碳氢化合物的影响.美国德克萨斯大学奥斯汀分校Raja和Levco[23]结合全局模型与流体模型,研究在正极性和负极性触发电压作用下,甲烷流注放电中击穿机理、等离子体组分和活性粒子生成的共性与区别.然而,现有的数值模拟通常采用零维全局模型或一维流体模型[24,25],缺乏与实验诊断结果相符合的二维模拟以及对大气压甲烷针-板放电等离子体中活性粒子的产生、输运、时空分布和反应路径等问题的全面分析与合理阐释.

本文对负直流高压电源激励下的大气压甲烷针-板放电等离子体建立二维流体模型并进行实验验证,开展数值模拟研究,呈现电场强度、电子温度和粒子密度的空间与轴向分布,计算产生带电粒子和中性粒子的反应贡献,探讨生成各种粒子的关键路径,为调控大气压甲烷针-板放电等离子体的产物选择性、提高甲烷转化重油加氢技术的能量利用效率,以及优化放电等离子体甲烷转化与重油加氢反应的耦合时间和作用方式等问题提供依据.

2 模型描述

大气压甲烷针-板放电等离子体二维流体模型的几何结构如图1所示,其中直径为2 mm的圆柱金属针(材料为钨)电极位于平行金属板(材料为铜)电极上方,针电极的长度为3 mm,针尖的曲率半径为1 mm,针电极尖端距离板电极表面的距离为2 mm,电压振幅为5 kV的负直流高压电源经过阻值为10 kΩ的限流电阻接针电极,板电极接地.求解区域大小为10 mm×5 mm.

图1 大气压甲烷针-板放电等离子体二维流体模型的几何结构Fig.1.Geometry structure of two-dimensional fluid model of methane needle-plane discharge plasma at atmospheric pressure.

工作气体甲烷在求解区域内设为均匀分布,气体压强设为1 atm(1 atm=101325 Pa),气体温度设为300 K.除基态甲烷分子CH4以外,数值模型中包含的21种带电和中性粒子有e,CH,CH2,CH3,C2H2,C2H4,C2H5,C2H6,C3H8,H和H2,粒子种类的选择参照文献[3,13].甲烷放电等离子体中的等离子体物理化学反应及其反应速率系数如表1所列,其中电子碰撞反应R1—R16的反应速率系数由碰撞截面计算,R1,R2,R3,R8和R13的碰撞截面取自Bolsig+软件,R4—R7,R9—R12和R14—R16的碰撞截面取自LXcat数据库,其他反应的反应速率系数取自参考文献[26—29].

表1 甲烷放电等离子体中的等离子体物理化学反应及其反应速率系数Table 1.Plasma physical and chemical reactions and the corresponding reaction rate coefficients in methane discharge plasma.

流体模型中的基本控制方程主要包括连续性方程、泊松方程和电子能量平衡方程[14,23].甲烷放电等离子体中各种粒子的动态行为满足如下连续性方程:

式中nk为第k种粒子的粒子数密度;Γ为粒子通量;S为粒子源项.

在连续性方程中,带电粒子通量为Γk=µknkE−Dk∇nk;中性粒子通量仅考虑扩散项,表示为Γk=−Dk∇nk.其中E为电场强度,µ和D分别为粒子的漂移率和扩散系数.

电场强度E由泊松方程描述如下:

式中φ为电位;ε0为真空介电常数;εr为气体的相对介电常数;e为电子电荷;qk为第k种粒子的电量.

电子温度Te由电子能量平衡方程描述如下:

式中KB为玻尔兹曼常数;下标e表示电子;ψe为电子能量通量;me为电子质量;ne为电子数密度;M为气体的平均分子质量;vel为弹性碰撞频率;Tg为气体温度;ri和∆εi分别为第i次反应的速率和由非弹性碰撞引起的能量损失.

在针电极和板电极处,电子通量Γe和电子能量通量ψe满足如下边界条件:

式中n为向外单位向量;re为反射系数;νe,th为电子热速度;γk为二次电子发射系数;Γk为第k种粒子的粒子通量;Γt为热发射通量;ϵk为第k种粒子的热离子能量;ϵ为平均热离子能量.

使用Comsol 5.0软件进行网格剖分和数值运算,在求解区域内使用轴对称结构和不均匀三角形网格,剖分结果如图2所示,其中生成的网格数量约为2.3×105,总自由度个数约为3.7×106.此外,初始电子密度设为1.0×1016m−3,初始正离子密度总和设为1.0×1016m−3,中性粒子密度的初始值设为1.0×1011m−3,电子温度的初始值设为3 eV.需要指出的是,对于大气压针-板放电等离子体,上述带电粒子密度和电子温度的初始值设置很高.

图2 求解区域内使用轴对称结构和不均匀三角形网格的剖分结果Fig.2. Subdivision results in calculation domain using axial symmetry structure and asymmetry triangle mesh.

为了验证建立的二维流体模型的正确性,使用Stark展宽法估算大气压甲烷针-板放电等离子体中的电子密度[30].需要指出的是,对于大气压针-板放电,Stark展宽法估算等离子体中的电子密度存在一定误差.图3为在二维流体模型中使用的放电条件下,开展实验测量得到的针电极下方区域的Hα谱线,测量时光纤垂直于放电通道放置,与放电等离子体反应器外边沿的距离为2 cm,将图中半高全宽区域的Stark展宽代入可得相应位置的电子密度数量级为1021m−3.

图3 Hα谱线Fig.3.Hαline.

图4 给出了外施电源的作用时间为6 ns时刻电子密度的空间分布的数值模拟结果,表明大气压甲烷针-板放电等离子体具有接近“等离子体子弹”的空间分布特性.从图中可以看出,在针电极尖端的下方区域内,电子密度最大值超过1.0×1020m−3,表明数值模拟结果与实验诊断数据基本符合,因此本文呈现的大气压甲烷针-板放电等离子体中电场强度、电子温度和粒子密度的空间与轴向分布以及反应路径贡献均为6 ns时刻的数值模拟结果.

图4 电子密度的空间分布Fig.4.Spatial distribution of electron densities.

3 数值模拟结果与讨论

图5 (a)电场强度和(b)电子温度的空间分布Fig.5.Spatial distributions of(a)electric intensity and(b)electron temperature.

图5 为大气压甲烷针-板放电等离子体中的电场强度和电子温度的空间分布,图5(a)中的箭头表示电场强度矢量的方向.如图所示,在半球形针电极头部的外侧边缘(与针电极外侧的距离小于0.1 mm)区域,电场强度和电子温度均达到最大值,分别高于5×106V/m和4.5 eV.已有研究指出,大气压甲烷DBD等离子体中电子密度和电子温度最大值分别约为5.5×1019m−3和3 eV[13],与之相比,大气压甲烷针-板放电等离子体有更高的电子密度和电子温度最大值.此外,在由针电极尖端指向板电极中心位置的轴向方向,电场强度和电子温度的变化趋势均为在与针电极尖端的距离约为0.1—0.6 mm的区域内先减小后增大、在与针电极尖端的距离大于0.6 mm的区域内再逐渐减小.上述电场强度和电子温度的数值模拟结果表明,在负直流高压电源的激励和空间电荷的作用下,大气压甲烷针-板放电的等离子体鞘层不明显,在靠近阴极的区域内能够形成阴极暗区和负辉区,在远离阴极的区域内则主要由正柱区构成[31].

图6 带电粒子密度的空间分布Fig.6.Spatial distributions of charged particle densities,including(a)

图7 带电粒子密度的轴向分布Fig.7.Axial distributions of charged particle densities,including

图8和图9分别为数值模拟结果中密度较高的6种中性粒子(CH,CH2,CH3,C2H4,C2H5和H)密度的空间分布和轴向分布.从图中可以看出,甲基CH3和氢原子H是大气压甲烷针-板放电等离子体中粒子密度最高的两种中性粒子,其粒子密度最大值达到1020m−3数量级.并且,CH3和H的粒子密度空间分布和轴向分布几乎相同.与甲烷-二氧化碳混合气体DBD等离子体的数值模拟结果相比,甲烷针-板放电等离子体中的CH3密度更高,与H密度较为接近[3].亚甲基CH2、乙烯C2H4和乙基C2H5是另外3种粒子密度较高的中性粒子,其粒子密度分布在靠近阴极的区域内明显不同,最大值分别处于1019,1018和1017m−3数量级,而在正柱区内则呈现均匀且相像的分布规律,粒子密度分别约为5.6×1015,7.2×1015和3.7×1015m−3.在与针电极尖端的距离约为0.1 mm处,次甲基CH具有与C2H5接近的粒子密度最大值(约为2.1×1017m−3),但在正柱区内CH的粒子密度降低至1012m−3数量级.此外,CH粒子的轴向分布与其他中性粒子有显著区别,而与带电粒子更为相似,是本文模拟结果中仅有的在与针电极尖端距离约为0.1—0.6 mm的区域内粒子密度先减小后增大的中性粒子.

图8 中性粒子(a)CH,(b)CH2,(c)CH3,(d)C2H4,(e)C2H5和(f)H密度的空间分布Fig.8.Spatial distributions of neutral particle densities,including(a)CH,(b)CH2,(c)CH3,(d)C2H4,(e)C2H5and(f)H.

图9 中性粒子CH,CH2,CH3,C2H4,C2H5和H密度的轴向分布Fig.9.Axial distributions of neutral particle densities,including CH,CH2,CH3,C2H4,C2H5and H.

图10 为大气压甲烷针-板放电等离子体中带电粒子的反应路径,其中实线表示生成箭头所指向的带电粒子的惟一反应或主要反应,虚线表示生成箭头所指向的带电粒子的次要反应,红色和黑色的线段分别表示是否有工作气体甲烷参与的反应,线段起点和线段中注明的粒子表示参与反应的粒子,线段终点的粒子表示该反应生成的粒子,线段中注明的百分比表示该反应产额占生成该粒子的所有反应总产额的百分比,括号内的序号为该反应在表1中的反应序号.如图10所示,电子与基态甲烷分子CH4间的电子碰撞电离反应R5和R4是生成的主要路径.分别与CH4发生分子碰撞解离反应R18和R17生成大部分的电子与氢气分子H2间的电子碰撞电离反应R10生成发生分子碰撞解离反应R20生成分别与CH4和C2H4分子发生反应R21和R23生成少量的与乙烷分子C2H6间的反应R19也生成少量的由于电子碰撞电离反应R5和R4是生成的主导反应,而电场强度和电子温度为在与针电极尖端距离约为0.1—0.6 mm的区域内先减小后增大,因此和的粒子密度在该区域内也呈现先减小后增大的变化规律.

图10 带电粒子的反应路径Fig.10.Reaction pathway of charged particles.

图11 碳氢化合物分子的反应路径Fig.11.Reaction pathway of hydrocarbon molecules.

生成碳氢化合物分子的反应路径如图11所示,其中蓝色和黑色的线段分别表示是否有电子参与的反应,其他的设置与图10相同.图11表明,电子与基态CH4分子间的电子碰撞分解反应R6,R7和R8分别是生成甲基分子CH3、亚甲基分子CH2和次甲基分子CH的主要路径.CH2与CH4分子间的反应R30也生成CH3,该反应产额占CH3分子总产额的12.68%.CH与CH4分子发生反应R29生成乙基分子C2H5.CH3分别与CH3和C2H5分子发生反应R24和R27生成乙烷分子C2H6和丙烷分子C3H8.尽管有多条反应路径可以产生乙烯分子C2H4,但CH2与CH4分子间的反应R31是生成C2H4的主导反应.电子与C2H4分子间发生的电子碰撞分解反应R14是生成乙炔分子C2H2的关键路径.

产生氢原子H的反应路径有R5,R6,R8,R9,R12,R20,R32和R33,其中电子与CH4间的电子碰撞分解反应R6的反应产额比其他反应的总产额高105数量级,是生成H的主导反应.反应R7,R11,R14,R18,R19,R21,R22,R23,R28,R31,R34和R35产生氢分子H2,其中电子与CH4间的电子碰撞分解反应R7和CH2与CH4发生的反应R31的产额各占H2总产额的52.15%和47.85%,是生成H2的主要路径.

4 结 论

对由负直流高压电源驱动的大气压甲烷针-板放电等离子体建立二维流体模型并进行实验验证,数值模拟研究其中的粒子密度和反应路径.模拟结果表明,大气压甲烷针-板放电能够形成阴极暗区、负辉区和正柱区.粒子密度的最大值都超过1018m−3数量级.和密度在沿轴向达到最大值后逐渐减小.和的粒子密度轴向演化与电场强度和电子温度的轴向演化接近.的粒子密度最大值出现在针电极尖端下方约0.1 mm处,在其他区域内的粒子密度分布较为均匀.CH3和H的粒子密度最大值达到1020m−3数量级,且空间分布和轴向分布几乎相同.CH2,C2H4和C2H5的粒子密度分布在靠近阴极的区域内明显不同而在正柱区内较为相像.CH的粒子密度最大值与C2H5接近,但在正柱区内粒子密度降低至1012m−3数量级.电子与CH4间的电子碰撞电离是生成和的主要路径.分别与CH4发生分子碰撞解离生成电子与H2间的电子碰撞电离生成发生分子碰撞解离生成电子与CH4间的电子碰撞分解是生成CH3,CH2,CH和H的主导反应.CH与CH4反应生成C2H5.CH3与CH3和C2H5反应分别生成C2H6和C3H8.CH2与CH4和电子与C2H4发生的反应分别是生成乙烯C2H4和乙炔C2H2的关键路径.电子与CH4间的电子碰撞分解反应和CH2与CH4发生的反应的产额各占H2总产额的52.15%和47.85%,是生成H2的主要路径.

[1]Liu C,Chernets I,Ji H,Smith J,Rabinovich A,Dobrynin D,Fridman A 2017IEEE Trans.Plasma Sci.45 683

[2]Kang H,Lee D,Kim K,Jo S,Pyun S,Song Y,Yu S 2016Fuel Process.Technol.148 209

[3]Bie C D,van Dijk J,Bogaerts A 2015J.Phys.Chem.C119 22331

[4]Xu Y,Zhang X,Yang C,Zhang Y,Yin Y 2016Plasma Sci.Technol.18 1012

[5]Wang C,Zhang Z,Cui H,Xia W,Xia W 2017Chin.Phys.B26 085207

[6]Liu J L,Park H W,Chung W J,Park D W 2016Plasma Chem.Plasma Proc.36 437

[7]Zhang Z B,Wu Y,Jia M,Song H M,Sun Z Z,Li Y H 2017Chin.Phys.B26 065204

[8]Wang W,Snoeckx R,Zhang X,Cha M S,Bogaerts A 2018J.Phys.Chem.CDOI:10.1021/acs.jpcc.7b10619

[9]Bi Z H,Hong Y,Lei G J,Wang S,Wang Y N,Liu D P 2017Chin.Phys.B26 075203

[10]Zhang D Z,Wang Y H,Wang D Z 2017Chin.Phys.B26 065206

[11]Shao T,Wang R X,Zhang C,Yan P 2018High Voltage3 14

[12]Gao Y,Zhang S,Liu F,Wang R X,Wang T L,Shao T 2017Trans.China Electrotechnical Soc.32 61(in Chinese)[高远,张帅,刘峰,王瑞雪,汪铁林,邵涛 2017电工技术学报32 61]

[13]Snoeckx R,Setareh M,Aerts R,Simon P,Maghari A,Bogaerts A 2013Int.J.Hydrogen Energy38 16098

[14]Pan J,Li L 2015J.Phys.D:Appl.Phys.48 055204

[15]Sun A B,Li H W,Xu P,Zhang G J 2017Acta Phys.Sin.66 195101(in Chinese)[孙安邦,李晗蔚,许鹏,张冠军2017物理学报66 195101]

[16]Pan J,Li L,Wang Y,Xiu X,Wang C,Song Y 2016Plasma Sci.Technol.18 1081

[17]Niu Z T,Zhang C,Ma Y F,Wang R X,Chen G Y,Yan P,Shao T 2015Acta Phys.Sin.64 195204(in Chinese)[牛宗涛,章程,马云飞,王瑞雪,陈根永,严萍,邵涛2015物理学报64 195204]

[18]Pan J,Li L,Chen B,Song Y,Zhao Y,Xiu X 2016Eur.Phys.J.D70 136

[19]Babaeva N Y,Zhang C,Qiu J,Hou X,Tarasenko V F,Shao T 2017Plasma Sources Sci.Technol.26 085008

[20]Yang D P,Li S Y,Jiang Y F,Chen A M,Jin M X 2017Acta Phys.Sin.66 115201(in Chinese)[杨大鹏,李苏宇,姜远飞,陈安民,金明星2017物理学报66 115201]

[21]Yang W B,Zhou J N,Li B C,Xing T W 2017Acta Phys.Sin.66 095201(in Chinese)[杨文斌,周江宁,李斌成,邢廷文2017物理学报66 095201]

[22]Wang B,Yan W,Ge W,Duan X 2013Chem.Eng.J.234 354

[23]Levko D,Raja L L 2017Plasma Sources Sci.Technol.26 035003

[24]Yin Z Q,Wang Y,Zhang P P,Zhang Q,Li X C 2016Chin.Phys.B25 125203

[25]Wang Q,Yu X L,Wang D Z 2017Chin.Phys.B26 035201

[26]Herrebout D,Bogaerts A,Yan M,Gijbels R,Goedheer W,Dekempeneer E 2001J.Appl.Phys.90 570

[27]Lefkowitz J K,Guo P,Rousso A,Ju Y 2015Phil.Trans.R.Soc.373 20140333

[28]Adamovich I V,Li T,Lempert W R 2015Phil.Trans.R.Soc.373 20140336

[29]Takana H,Nishiyama H 2014Plasma Sources Sci.Technol.23 034001

[30]Nikiforov A Y,Leys C,Gonzalez M A,Walsh J L 2015Plasma Sources Sci.Technol.24 034001

[31]Yao C W,Ma H C,Chang Z S,Li P,Mu H B,Zhang G J 2017Acta Phys.Sin.66 025203(in Chinese)[姚聪伟,马恒驰,常正实,李平,穆海宝,张冠军 2017物理学报 66 025203]