针-板空气间隙流注放电起始过程的三维PIC/MCC仿真研究∗

李晗蔚 孙安邦 张幸 姚聪伟 常正实 张冠军

(西安交通大学,电力设备电气绝缘国家重点实验室,西安 710049)

1 引 言

流注放电(streamer discharge)多呈现出明亮、丝状通道结构,通常在施加电场大于介质本身击穿电场的区域形成,且在其头部存在强烈的电场增强效应,具有向低于击穿电场的未电离区域传播的特性[1−3].作为高电压领域气体放电的一种常见初始形式,流注放电在自然界中大气放电、特高压输电线路长空间间隙放电以及工业应用中都扮演着非常重要的角色.

流注放电理论于20世纪40年代初由Raether,Loeb和Meek提出[4,5],长期以来,为了解流注放电的物理机理,国内外学者进行了大量的实验与数值仿真研究.在实验研究方面,最具有代表性的为荷兰埃因霍温理工大学应用物理系,采用棱镜法等对流注光电特性进行了大量实验研究,包括其三维结构、传播速度、流注尺寸等[6,7].然而,由于流注放电的产生与发展过程非常快,在大气压下一般为纳秒量级,且实验主要集中在毫米至厘米尺度放电间隙内进行,因此实验手段的精确测量具有较大难度,当前流注物理特性的研究仍主要以数值仿真为主.

流注放电的数值仿真模型主要分为三类[8]:流体模型、粒子模型和流体粒子混合模型.流体模型从宏观角度描述等离子体的放电过程,计算效率较高.如,Kulikovsky[9]首次将光电离加入到流注放电的二维轴对称漂移-扩散近似流体模型中,研究了两平板电极空气间隙的流注放电特性,得到流注传播行为主要受头部电荷层宽度影响的结论;Luque和Ebert[10]引入自适应可变网格技术,采用三维流体模型研究了空气中多个流注在压强、氮氧混合比等因素影响下的运动行为,发现相邻流注互相作用的机制主要是电荷的排斥和非局域光电离造成的相互吸引.但是,流体模型的使用范围具有很大局限性,也很难描述放电过程的随机行为(比如分叉等).基于第一性原理建立的粒子模型从微观角度直接跟踪放电过程中带电粒子的运动过程,进而统计得到其宏观参数的变化.如,西安交通大学李永东教授课题组采用粒子模型研究了大气中光电离对流注形成和传播的影响,获得了逃逸电子产生过程的物理图像以及流注的分叉结构[11];我们与Teunissen博士和Ebert教授合作,开发了一套流注放电三维粒子仿真模型,研究了均匀电场下自然界中光电离、背景电离(自由电子在强电场作用下从O−2解吸附)对流注形成的影响[12−14];粒子模型的主要缺点在于计算量相对较大.结合粒子和流体模型各自的优点,流体-粒子混合模型在电场梯度较小区域应用流体模型,在流注外围及其他电场梯度大的区域应用粒子模型[15],然而,由于存在粒子流体的过渡区域,常规的并行计算以及可变网格等新技术较难实现,总体计算量依然比较大.

纵观空气中流注放电的数值研究现状,宏观流体模型仍是流注数值仿真的常用手段,且大多局限于二维仿真,仅适用于对单个流注的传播特性进行研究.本文侧重于对针-板电极结构中流注放电开展三维粒子模拟研究,研究价值在于:首先,针-板电极结构产生不均匀电场,电场在间隙中呈不均匀发散状,影响流注演化特性与分叉结构;第二,针型电极附近产生高强度电场区域,具有较高的电子数密度,流注放电易于从针电极处产生;第三,实际放电大多是在不均匀电场结构中产生,而现有数值研究大多局限于二维仿真或均匀电场的仿真,亟需对以针-板电极结构类的非均匀电场下流注放电过程开展三维流注放电研究,反映真实放电形貌和结构特征.

2 仿真物理模型介绍

2.1 三维PIC/MCC数值模型

本文采用我们与Ebert教授课题组合作开发的三维PIC/MCC模型对针-板电极放电起始过程开展数值研究.该模型前期已用于研究均匀电场空气中背景电离对流注形成的影响[12−14].PIC/MCC数值模拟的流程如图1所示.

为了减少计算时间,本模型中采用的模拟环境和基本假设如下:

1)模拟气体为干燥空气,主要为氮氧混合气体组成;

2)由于离子的运动速度远小于电子,在所关注的较短的放电时间内 (大气压下ns时间尺度),离子假设为静止;

3)整个模拟过程处于弱电离状态,只考虑电子-中性粒子之间的碰撞,忽略电子-电子碰撞以及电子-离子之间的库仑碰撞,中性气体假设为背景,参与电子之间的碰撞行为,而本身的密度不做改变;电子与中性气体的碰撞界面由Lxcat数据库获得[16];

4)忽略电荷运动产生的自洽磁场,通过求解Poisson方程更新模拟区域的电场变化.

图1 PIC/MCC模拟流程图Fig.1.Flow chart of PIC/MCC simulation.

2.2 电子运动与碰撞

2.2.1 电子运动

电子的运动遵循牛顿-洛伦兹定律,不考虑磁场作用,运动方程为:

式中,me为电子质量,q为电子电荷,v和x分别为电子速度及其位置矢量,E为电场.

采用Verlet算法描述电子运动,具体计算方法为:

式中,vk,xk为电子运动k个时间步长后的速度和位置;ai为k个时间步长后的电子加速度.

2.2.2 粒子碰撞

模型中,考虑粒子之间的弹性碰撞和非弹性碰撞(例如电子激发、电离、吸附等),主要反应方程如下:

弹性碰撞

电离碰撞

激发碰撞

电子吸附

电子在不同碰撞截面积σi与背景气体分子发生碰撞时,其碰撞频率可以表示为

式中:N为背景气体分子密度;v为电子相对速度.

在Δt时间内,电子与中性气体碰撞概率可表示为

2.3 光电离(photoionization)

流注从正极性针型电极起始和发展过程中,其传播方向和外加电场方向相同,与电子漂移方向相反,这意味着外部电子需要从其流注前方反向运动到流注通道内部.在空气中,光电离是产生外部自由电子的主要机制:空气中的N2受激发之后,处于激发态的N2分子退激发为基态N2分子时释放出的光子与O2分子发生碰撞电离,从而释放出自由电子[8].虽然光电离机制被广泛采纳并且应用已久,然而其内部具体的物理过程还不甚清楚.本模型采用Zheleznyak等提出的光电离模型,计算流注放电过程中的光电离过程:

定义湮灭因子ηq:

式中:p为压力;pq为湮灭压力,pq=30 mbar.

当p=pq时,一半激发态的分子会由于与其他分子的碰撞而退激发,并不会对外释放出光电子.在整个模型中,平均每单位光电离会产生η个光电子,定义η=ηE·ηq,ηE是一个与电场强度有关的因子[18].

2.4 粒子数值模拟加速技术

三维粒子模型的计算量极大,因此采取并行计算的方法,加速电势及电场的计算速度.同时,也为了在计算过程中保持计算准确度,使用了自适应可变网格(adaptive mesh re fi nement,AMR)、粒子可变权重和并行计算技术.

2.4.1 自适应可变网格

根据电子的密度梯度和电场大小对计算区域的网格大小进行调整.当区域内电子数密度超过设定阈值时,区域内网格将自动加密处理,新网格大小为原网格的一半;电子数密度低于阈值时,网格将自动粗化处理.图2是多重自适应网格示意图.在电场强度高、粒子数密度梯度较大的区域内,网格划分的较细;在电场较低、粒子数密度梯度较低的区域内,网格划分的相对较粗.

图2 多重自适应网格Fig.2.Multi-scope adaptive grid.

不同尺寸网格数量、最大和最小网格尺寸决定着模拟的计算精度和计算时间.本文计算中,允许网格细化次数为5,在经过5次细化后,达到最小网格,不再继续进行细化.网格进行细化时的最小电子数密度为1×1011m3;每隔2.0×10−11s判断网格是否需要细化.

2.4.2 粒子可变权重

放电过程中,计算域内粒子数目多,随放电的发展,粒子数目也急剧增长.为了控制模拟速度并减小计算量,对计算粒子权重进行可变调节,如图3所示,每个计算粒子可以分解为两个计算粒子,或与其他计算粒子合并成为一个计算粒子.

图3 计算粒子的分裂与合并Fig.3.The splitting and merging of calculated particles.

详细的权重调整方法以及权重变化后粒子的速度和位置分布处理等,可参见文献[19,20].

3 结果与讨论

本文主要探讨电压幅值、气体组分和电极尖端形状尺寸三个参数对流注放电的影响,着重关注放电过程中形成的初始电离云、流注分叉结构和电场强度差异等.

3.1 仿真初始条件

3.1.1 计算域和电极

本模型中,计算域为图4(a)所示的(5×5×5)mm3立方体区域,电极尖端正对平面为地电极.

图4 计算域示意图及电极尖端放大图Fig.4.The schematic of the computational domain and the enlarged electrode tip.

计算域的选择主要出于以下考虑:

1)在立方体计算域中,电极尖端处于立方体正中心处,与各平面距离相同,流注向空间各方向发展概率均等;若采用长方体区域,电极尖端与各平面距离不等,导致空间中电场强度分布存在较大差别,可能也会影响流注的发展趋势;

2)在一个相对较小的区域内,较低电压即可实现流注放电.

针型电极由球体、圆柱与圆锥模拟而成,如图4(b)所示,以在电极表面放置点电荷的形式描述外施加电压.

3.1.2 初始参数设置

模型中,模拟区域内背景电离密度为1012m−3,等离子体在整个计算域空间内随机分布.为加快流注放电仿真速度,在电极尖端正上方0.05 mm处人为放置一定数量的电子和正离子(宏观为中性),带电粒子数密度呈高斯分布,半径为50µm,包含3×106个电子并呈现为高斯分布.

模拟中,气体压强为1 bar,温度293 K,使用N2和O2的混合气体模拟空气,比例为N2:O2=80%:20%.

3.2 电压幅值对流注放电的影响

选择电压参数分别为:3.0,3.5,4.0和4.5 kV,分别在此四个电压幅值下进行模拟.

3.2.1 初始电离云

流注放电开始阶段,在针尖附近观察到形成球状类似于电晕的初始电离云(inception cloud),这与之前实验中观测到的现象是一致的[6].如图5所示,电离云大小随施加电压的增大而增大.

图5 不同电压幅值下初始电离云图像 (a)3.0 kV,1.0 ns;(b)3.5 kV,0.8 ns;(c)4.0 kV,0.6 ns;(d)4.5 kV,0.4 nsFig.5.Photos of inception clouds under different voltage amplitudes:(a)3.0 kV,1.0 ns;(b)3.5 kV,0.8 ns;(c)4.0 kV,0.6 ns;(d)4.5 kV,0.4 ns.

假定初始电离云是一个半径为R的球体,其表面近似等电位U0,电极间隙的击穿场强Ec,则其最大半径可近似估算为R0=V0/Ec.因此,电压越高,初始电离云半径越大,当电离云发展到最大半径R0位置时,单个流注通道开始出现并向前发展,在此类低场强区域,流注的传播主要依靠头部空间电荷产生的较强自洽电场(可达到3Ec—5Ec)的作用.

图6给出了不同电压下产生初始电离云时的电场强度分布.所有图像均是x-z轴平面,在y轴方向取计算区域的中心位置截断取得的,即y=2.5 mm处.后文中所有电场强度图像均是按照此方法得到.

可以发现:各电压条件下,针型电极表面存在一个较薄的高场强区域,即为上述所描述的电场增强效应;在人为放置种子电子的影响下,电压较低(3.0 kV及3.5 kV)时,种子电子处整体电场强度较高,而当电压升高时(4.0 kV及4.5 kV),仅其头部成为高电场强度区域.因此,高斯种子电子的放置会改变电极尖端附近的电场强度分布.

图6 初始电离云电场强度截面图 (a)3.0 kV,1.0 ns;(b)3.5 kV,0.8 ns;(c)4.0 kV,0.6 ns;(d)4.5 kV,0.4 nsFig.6.Cross Section view of inception clouds’electrical field strength:(a)3.0 kV,1.0 ns;(b)3.5 kV,0.8 ns;(c)4.0 kV,0.6 ns;(d)4.5 kV,0.4 ns.

3.2.2 流注形态结构

图7 不同电压幅值下流注放电电子数密度演化过程 (a)3.0 kV;(b)3.5 kV;(c)4.0 kV;(d)4.5 kVFig.7.Electron density evolution of streamer discharge under different voltage amplitudes:(a)3.0 kV;(b)3.5 kV;(c)4.0 kV;(d)4.5 kV.

图7是流注放电中电子数密度随时间演化过程的模拟结果,为了便于比较,所有流注图像的观测角度均相同.因不同电压下的流注放电持续时间不同,图中进行比较的时间选取略有不同.

常用者为头孢哌酮、头孢哌酮舒巴坦、头孢克肟，菌素对肠杆菌科细菌等革兰阴性杆菌具有强大抗菌作用，头孢他啶和头孢哌酮除肠杆菌科细菌外对铜绿假单胞菌亦具高度抗菌活性。

可以看出,在流注的发展过程中,易产生分叉形成类似树枝状结构.不同电压幅值作用下,流注发展形态有很大差别.电压越高,流注分叉时间越早,其分支结构越多.

电压为3.0 kV时,流注形态与其他电压下流注形态有较大差别.第一,如3.2.1节所述,几乎未形成初始电离云;第二,其他电压工况下流注分叉,形成多流注通道;而此工况下形成单流注通道,直到5.4 ns时流注开始分叉、产生多流注通道的迹象.此外,在此电压工况下,流注通道内电子密度比其他更高电压下略高,可能是由于流注发展过程中通道相对集中所导致的.随放电发展,现有流注通道将会形成更多的流注通道.

3.3 气体组分对流注放电的影响

根据文献,N2-O2气体混合物中,N2与O2的比例直接影响流注放电的结构特征.本节在外加4.0 kV电压下,确定另外两种工况下N2和O2的比例为:0:100和99.9:0.1.前一种情况等效为纯净O2,后一种情况等效为N2中掺杂少量O2.

Nijdam等[21]在实验中观察了不同O2浓度下流注放电的发展过程,发现在100%O2中几乎看不到初始电离云的存在;N2体积分数越大,流注分叉越多,且分支越细.

3.3.1 初始电离云

在本文的仿真条件下,流注放电在不同气体组分下形成的初始电离云形貌如图8所示:初始电离云都具有较高的电子数密度,随着O2浓度提高,电离云半径增大,且单独流注通道形成时间提早.根据流注理论,光电离是空气中产生二次电子崩的主要原因,而二次电子崩是流注向前发展的重要因素;在本模型中同时考虑了电子的吸附与解吸附过程.因此,O2所占体积越高,光电离和电子解吸附过程越容易提供自由电子,利于正流注的向前发展.

图8 不同气体组分影响下的初始电离云图像 (a)0.1%O2,0.8 ns;(b)20%O2,0.6 ns;(c)100%O2,0.4 nsFig.8.Comparison of inception clouds with different gas components:(a)0.1%O2,0.8 ns;(b)20%O2,0.6 ns;(c)100%O2,0.4 ns.

3.3.2 流注形态结构

三个工况下流注的发展形态如图9所示.

图9 不同气体组分流注发展图像 (a)0.1%O2;(b)20%O2;(c)100%O2Fig.9.Photos of the development of streamer discharge with different gas components:(a)0.1%O2;(b)20%O2;(c)100%O2.

图10 不同气体组分流注放电电场强度截面图 (a)0.1%O2;(b)20%O2;(c)100%O2Fig.10.Sectional views of electrical field strength of streamer discharge with different gas components:(a)0.1%O2;(b)20%O2;(c)100%O2.

可以看出,在0.1%O2浓度工况中,流注的电子数密度明显大于高O2浓度工况下流注的电子数密度;在3.2 ns内,形成的分支较少,前期近似为单流注向板电极发展,在2.2 ns时可见明显流注分支结构;流注表面整体较粗糙,类似于存在许多细小凸起.随O2浓度的增高,20%O2工况下,流注分叉数目,尤其是向两侧分叉数目逐渐变多;100%O2工况下,除形成向上的主流注通道外,还向四周形成了细小的流注放电通道.

图10给出了三种工况下流注发展时的场强变化图.流注放电初期,尖端电极表面存在高场强区域薄层;随着流注发展,薄层消失;此后,流注头部电场强度较高.初期流注可能向空间各个方向发展,后期由于流注通道之间的相互抑制,部分流注由于电场强度过低或少量自由电子注入而熄灭,最终产生随机形态的分叉结构.

在0.1%O2工况时,流注头部始终存在大片高场强区域,如图10(a)所示,由于光电离项的影响,该工况下流注发展较慢,流注通道没有完全形成;随着流注继续向前发展,会在通道内部产生大量电荷,电场强度相应降低.

3.4 电极形状尺寸对流注放电的影响

上文模拟设置中,电极尖端球体半径为50µm,本节考虑了另外两种工况为尖端球体半径5µm和250µm,如图11所示.

图11 不同尺寸电极尖端示意图 (a)5µm;(b)250µmFig.11. Photos of different electrode tip sizes:(a)5µm;(b)250µm.

3.4.1 初始电离云

三种工况下,初始电离云图像及其电场强度截面图分别如图12和图13所示.

图12 不同电极尖端尺寸工况下初始电离云图像 (a)5µm,0.4 ns;(b)50µm,0.6 ns;(c)250µm,1.2 nsFig.12.Photos of inception clouds with different electrode tip sizes:(a)5µm,0.4 ns;(b)50µm,0.6 ns;(c)250µm,1.2 ns.

图13 不同电极尖端尺寸形成的初始电离云电场强度截面图 (a)5µm,0.4 ns;(b)50µm,0.6 ns;(c)250µm,1.2 nsFig.13.Sectional views of electrical field strength of inception clouds with different electrode tip sizes:(a)5 µm,0.4 ns;(b)50µm,0.6 ns;(c)250µm,1.2 ns.

可见,尖端半径为250µm的工况下,几乎没有形成初始电离云,且电极尖端处电场强度较低.在5µm和50µm电极尺寸尖端均有密度较高的电场强度区域.在尖端半径较大时,低电场强度会导致碰撞电离减少,放电不足以形成高电子数密度聚集的初始电离云.

3.4.2 流注形态结构

尖端电极尺寸影响下流注的发展形貌如图14所示.

图14 不同电极尖端流注放电发展图 (a)5µm;(b)50µm;(c)250µmFig.14.Photos of streamer discharge with different electrode tip sizes:(a)5 µm;(b)50 µm;(c)250 µm.

对于电极尖端半径为5µm的工况,与原先50µm时流注形态对比,流注后期发展无明显区别.这是由于流注分叉且独立通道形成后,空间电荷产生的自洽电场成为影响流注发展的主要因素,因此后期流注发展与电极尖端形状尺寸关系不大.尖端半径为250µm的工况下,流注放电形态与前两种工况出现了明显区别,未形成明显的初始电离云,且随放电时间发展,基本呈现单流注形式,并在5 ns左右产生分叉.可以看出,针型电极尺寸可显著影响流注发展过程分叉的时间和结构.图15给出了5 ns后,250µm半径流注发展的图像.

综合图15和图14(c)可知,由于电极尺寸较大,电极上方始终存在一片高电子数密度区域;然而随时间发展,流注后期仍会产生分叉结构,与其他工况下的流注发展无异.这同样是由于后期流注的发展主要取决于空间电荷的自洽电场,与电极附近的电场强度高低关系不大.

图15 尖端球体半径为250µm的电极5—8 ns流注图像Fig.15.Photos of streamer discharge from 5 ns to 8 ns with 250µm radius of electrode tip.

4 结 论

本文针对针-板电极结构空气中流注放电现象,采用三维粒子仿真模型,确定计算域和相关计算参数,调整电压幅值、气体组分和电极尖端形状尺寸参数并进行相关仿真,分析各不同工况下的初始电离云及后续流注发展形态结构,得出以下结论.

1)放电初始阶段,电极尖端附近形成近似为球体的初始电离云,具有很高的电子数密度.不同工况下的初始电离云半径不同,其发展为多个单独流注通道的时间也不同.后期流注的形态发展方向具有很大的随机性.

2)随电压幅值增大,电离云的半径越大,电压工况为4.0 kV和4.5 kV时流注半径无明显差异;流注放电所形成的树枝状分叉结构也随着电压的增高变多.

3)N2-O2气体混合物的比例显著影响流注放电形貌.随N2-O2混合物中O2浓度的增大,初始电离云的半径逐步增大,分叉结构和数目明显变化.

4)改变电极尖端形状尺寸对放电初期影响较大,而在放电后期,电极尖端的改变对流注发展无明显影响,这主要是因为流注后期发展过程中自洽电场起主导作用.

[1]Nijdam S 2011Experimental Investigations on the Physics of Streamers(Eindhoven:Technische Universiteit Eindhoven)pp2–4

[2]Nijdam S,van Veldhuizen E,Peter B 2012Plasma Chemistry&Catalysis in Gases&Liquids(Germany:WILEY-ICH)pp1–44

[3]Zhang C,Tarasenko V F,Shao T,Beloplotov D V,Lomaev M I,Wang R,Sorokin D A,Yan P 2015Phys.Plasmas22 033511

[4]Raether H 1939Z.Phys.112 464

[5]Loeb L B,Meek J M 1940J.Appl.Phys.11 438

[6]Briels T M P,van Veldhuizen E M,Ebert U 2008J.Phys.D:Appl.Phys.41 234008

[7]Nijdam S,Moerman J S,Briels T M P,van Veldhuizen E M,Ebert U 2008Appl.Phys.Lett.92 101502

[8]PengQ J2012Ph.D.Dissertation(Chongqing:Chongqing University)(in Chinese)[彭庆军 2012博士学位论文(重庆:重庆大学)]

[9]Kulikovsky A A 1997J.Phys.D:Appl.Phys.30 441

[10]Luque A,Ebert U 2014New J.Phys.16 013039

[11]Li Y D,Wang R P,Zhang Q G,Zhou Y,Wang H G,Liu C L 2011IEEE Trans.Plasma Sci.39 2226

[12]Sun A B,Teunissen J,Ebert U 2013Geophys.Res.Lett.40 2417

[13]Sun A B,Teunissen J,Ebert U 2014J.Phys.D:Appl.Phys.47 445205

[14]Teunissen J,Sun A B,Ebert U 2014J.Phys.D:Appl.Phys.47 365203

[15]Li C,Ebert U,Hundsdorfer W 2010J.Comput.Phys.229 200

[16]Lxcat,http://lxcat net laplace univ-tlse fr[2017-3-22]

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

[18]Teunissen J,Ebert U 2016Plasma Sources Sci.Technol.25 044005

[19]Sun A B,Becker M M,Loffhagen D 2016Comput.Phys.Commun.206 35

[20]Teunissen J,Ebert U 2014J.Comput.Phys.259 318

[21]Nijdam S,van de Wetering F M J H,Blanc R,van Veldhuizen E M,Ebert U 2010J.Phys.D:Appl.Phys.43 145204