2 cm 电子回旋共振离子推力器离子源中磁场对等离子体特性与壁面电流影响的数值模拟*

夏旭 杨涓† 付瑜亮 吴先明 耿海 胡展

1) (西北工业大学 航天学院 西安 710072)

2) (兰州空间技术物理研究所 兰州 730000)

电子回旋共振离子推力器(electron cyclotron resonance ion thruster, ECRIT)离子源内等离子体分布会影响束流引出, 而磁场结构决定的ECR 区与天线的相对位置共同影响了等离子体分布.在鞘层作用下, 等离子体中的离子或电子被加速对壁面产生溅射, 形成壁面离子或电子电流, 造成壁面磨损和等离子体损失, 因此研究壁面电流与等离子体特征十分重要.为此本文建立2 cm ECRIT 的粒子PIC/MCC (particle-in-cell with Monte Carlo collision)仿真模型, 数值模拟研究磁场结构对离子源内等离子体与壁面电流特性的影响.计算表明, 当ECR 区位于天线上游时, 等离子体集中在天线上游和内外磁环间, 栅极前离子密度最低, 故离子源引出束流、磁环端面电流和天线壁面电流较低.ECR 区位于天线下游时, 天线和栅极上游附近的等离子体密度较高, 故离子源引出束流、天线壁面电流和磁环端面电流较高.腔体壁面等离子体分布与电流受磁场影响最小.

1 引 言

ECRIT 具有无阴极烧蚀特点, 除了用于空间飞行器控制[1,2], 还能用于材料处理等方面[3,4].目前国外研究出10 cm 和2 cm ECRIT, 并分别应用于HAYABUSA[2,5]和50 与70 kg 微小卫星[6,7].目前国内外热衷于微小卫星推进装置的研究, 因此开展2 cm ECRIT 研究具有重要意义.离子源是2 cm ECRIT 重要的部件, 其内部物理过程包括电离、碰撞以及鞘层加速带电粒子对壁面轰击, 这些过程决定了等离子体和壁面电流特征.等离子体特征影响推力器引出束流大小, 同时推力器内壁面溅射电流是约束工作寿命和引起等离子体损失的重要因素,而且等离子体和壁面电流存在密切关系并同时受磁场结构和工作参数等的影响, 因此开展不同磁场结构对等离子体和壁面电流特征的影响研究十分重要.由于实验研究[8−10]很难观测到离子源内瞬态物理过程变化, 需要借助数值模拟手段研究离子源内部等离子体过程.虽然PIC 模拟已经揭示了ECRIT 的等离子体特征[11−14], 但2 cm ECRIT 离子源壁面电流数值模拟少有报道.为此本文采用PIC/MCC 方法研究磁场结构对离子源内等离子体与壁面电流特性的影响, 为推力器的性能预估提供参考.

2 物理和计算模型

2.1 物理模型

如图1 所示, 2 cm ECRIT 离子源由圆柱腔体、环形天线、磁轭、内外环形永磁体和双栅极组成, 其中两磁环分布在离子源底部, 其间环向均布8 个进气孔; 4.2 GHz 的微波能量通过环形天线馈入放电室, 在0.15 T 的磁场区域形成ECR 区[15].在ECR 区, 电子绕磁力线的回旋运动和微波电场变化同步产生谐振, 由此微波电场加热电子形成高能电子, 高能电子和原子发生碰撞激发和电离, 从而形成ECR 等离子体, 等离子体中的离子经由双栅极被高速引出.图1 中的H1,W1,H2和W2分别表示外磁环高度、外磁环宽度、内磁环高度和内磁环宽度.L1为天线环形段下表面与屏栅上表面之间的距离,L2为内磁环下表面与天线环形段上表面之间的距离, 本文中不同磁场结构中L1和L2的值相同.定义A, B, C 和D 分别表示天线环形段的内表面, 下表面, 外表面和上表面; E, F 分别表示内磁环端面、外磁环端面; G 表示腔体内表面.

图1 2 cm ECRIT 离子源结构示意图Fig.1.Schematic diagram of the 2 cm ECRIT ion source internal structure.

2.2 计算模型

图2 2 cm ECRIT 离子源的放电形貌Fig.2.Discharge image of 2 cm ECRIT ion source.

图2 为2 cm ECRIT 离子源的气体放电图像.从图2 中可以看出,z轴周围有一个圆形的等离子体区域, 其亮度较均匀, 环形天线的阴影区域是可见的.由图可见, 离子源内的等离子体分布近似轴对称分布.文献[9]也表明2 cm ECRIT 离子源中静磁场轴对称, 高频场接近于轴对称性.基于以上分析, 计算域等离子体参数分布具有轴对称性, 因此将放电室内的三维的环形电离区简化为二维轴对称区域.

2.3 全粒子PIC/MCC 模型

2.3.1 电磁场

本文将ECRIT 离子源内部的电磁场分解为3 个模块: 静电场、静磁场和微波高频电场, 这里忽略高频磁场, 因为其比静磁场小得多.静电场通过离散泊松方程∇·∇φ=−ρ/ε0直接求解得到, 其中φ,ρ和ε0分别为电势, 电荷密度和真空介电常数.除了对称边界条件为其余边界电势均为0 V; 静磁场作为背景场, 由COMSOL 多物理场软件计算, 静磁场的计算参数设置: 永磁体(钐钴合金)的相对磁导率为1.36; 磁轭(软铁)的相对磁导率为4000; 空气的相对磁导率为1.静磁场的求解域为物理模型的所有区域, 边界条件仅需设置外边界为零磁标量势边界.微波高频电场采用COMSOL 计算, 计算时不考虑等离子体电流, 这是由于微波功率较低时等离子体电流对微波电场影响甚微[16].微波电场的计算参数设置: 壁面(不锈钢)的相对介电常数为1, 电导率为 4×106S/m ;天线(铜)的相对介电常数为1, 电导率为 6×107S/m ;空气的相对介电常数为1, 电导率为0 S/m.微波电场的求解域为放电室区域, 边界条件仅需设置微波输入端口边界, 输入端口类型设置为同轴, 输入功率设为2 W.PIC 程序计算时, 每个时间步上高频场的瞬时值与静电场叠加共同推进电子运动; 而离子质量较大, 故高频场与离子作用弱, 认为离子不受高频场的影响, 只有静电场推动离子运动.

2.3.2 粒子碰撞

离子源以氙气为工质, 假设中性原子为背景粒子, 均匀分布在计算域中.考虑的碰撞类型有: 电子与原子的弹性碰撞、激发碰撞和一级电离碰撞,以及带电粒子间的库伦碰撞.电子与中性原子之间的碰撞是通过蒙特卡洛(MCC)方法处理, 碰撞截面采用Szabo[17]的截面数据; 采用直接蒙特卡洛(DSMC)方法来描述带电粒子间的库伦碰撞[18].带电粒子的运动采用Boris“蛙跳”算法求解[17].粒子运动的边界条件: 对称边界时, 粒子在该边界上反射, 其余边界上, 粒子完全吸收.

图3 不同磁场结构下磁场分布: (a) 1 号源; (b) 2 号源; (c) 3 号源Fig.3.Distribution of magnetic flux density inside of the discharge chamber: (a) 1; (b) 2; (c) 3.

2.3.3 磁场结构

如表1 所示, 计算中模拟了3 种不同磁场结构.通过COMSOL 软件计算离子源中磁场强度,其ECR 区如图3 所示.可以看出, 1, 2 和3 号源的ECR 区分别位于天线上游、与天线相交、天线下游.随着外磁环的高度和宽度增大, ECR 区逐渐向屏栅移动并远离天线.

表1 磁场结构参数Table 1.Magnetic circuit structure parameters.

2.4 计算条件

计算条件为氙气流量0.3 sccm (1 sccm =98.3 µg/s), 微波功率2 W, 微波频率4.2 GHz.中性粒子密度当作背景气体代入PIC 程序, 0.3 sccm流量时, 估算得到的中性原子密度范围为1.5×1020m−3[19].

根据文献[12]给出的离子源内等离子体密度为 2×1017m−3, 再假设平均电子温度为10 eV, 则德拜长度约为0.053 mm.保持真实的介电常数下,为避免数值过热, 将网格步长设定为0.025 mm,即使等离子体密度达到5 × 1017m–3也能保证计算精度.

3 结果和讨论

3.1 计算结果的收敛性

计算开始时, 设初始电子和离子温度分别为2,0.5 eV, 初始等离子体密度为 1×1016m−3.图4 记录了不同磁场结构离子源的带电粒子数随计算时间的变化.如图4 所示, 随着计算时间的进行,在50 万步以后, 电子和离子数都逐渐达到了线性稳态.

图4 带电粒子数随计算步的变化Fig.4.Quantities of charged particles versus computation step.

图5 不同磁场结构下等离子体分布: (a)电子密度; (b)离子密度Fig.5.Plasma distribution of ion source with different magnetic circuit: (a) Electron density; (b) ion density.

3.2 等离子体特性分布

不同磁场结构下离子源内部的电子和离子密度分布计算结果如图5 所示.从图中可以看出不同离子源内的电子和离子密度峰值都出现在ECR 区附近, 且电子密度和离子密度峰值分别是3.5×1017m−3/3.9×1017m−3,2.6×1017m−3/3.0×1017m−3, 4.1×1017m−3/4.2×1017m−3.由此可见, 3 种离子源内部的等离子体密度峰值均大于4.2 GHz微波频率的截止密度ncut-off=ε0meω2/e2≈2.1×1017m−3[20,21], 其中ω为微波频率, 说明了ECR等离子体的高密度特征.不同离子源内的等离子体密度峰值区面积大小为: 3 号源 > 1 号源 > 2 号源.2 号源的密度峰值区域最小, 是因为ECR 区被天线切断, 从而影响磁镜中的电子被ECR 区持续加热.1 号源的等离子体密度峰值位于天线上游,并且ECR 区位于较小的磁镜区, 所以天线上游的内外磁环间的等离子体密度较高而靠近屏栅区域的等离子体密度较低.这是因为天线环形段阻碍了等离子体扩散.3 号源的ECR 区面积大于1, 2 号源, 而且ECR 区位于较大的磁镜区, 所以其等离子体密度最高而且靠近栅极前的等离子体高于1,2 号源.栅极前的等离子体密度越高越有利于离子束流的引出, 符合文献[15]的实验结果, 即ECR区位于天线下游的离子源引出效果优于ECR 区位于天线左侧的离子源.

图6 为不同磁场结构下离子源的电子温度分布.1, 2 和3 号源中电子在ECR 区被加热后, 沿磁力线扩散到ECR 区附近的磁镜中, 电子温度峰值分别为13.7 eV, 13.3 eV 和12.7 eV.1 号源中电子温度最高, 因为电子ECR 加热区靠近微波电场幅值高的天线上游[15].与图6 对比, 1 号源的电子温度峰值出现在内外磁环间的ECR 区和天线上表面附近的ECR 区附近, 但离子密度峰值出现在天线上游附近的ECR 区.这是因为天线上表面ECR 区的高能电子未与原子碰撞电离就损失在天线上表面, 而内外磁环间ECR 区的高能电子会与原子发生电离碰撞.2 号源的电子温度峰值靠近在天线下表面ECR 区, 而高密度区偏移了ECR 区.这也是因为未发生电离碰撞的高能电子损失在天线下表面.3 号源中电子被ECR 区加热后, 沿磁力线扩散到ECR 区和栅极上游附近的磁镜中.栅极上游磁镜中的高能电子多于ECR 区附近的高能电子, 这是因为栅极上游的磁镜范围更大, 而且ECR区附近的高能电子可能损失在天线环形段表面或与原子发生电离碰撞.3 号源的高密度区位于ECR区和栅极上游之间的磁镜中.这是ECR 区加热产生的高能电子向栅极扩散中, 部分高能电子与原子发生电离碰撞, 从而ECR 区和栅极上游之间的电子温度较低.部分高能电子和未与原子发生电离碰撞的电子可能在磁镜的约束下又往复穿越ECR区, 被ECR 区加热并约束在栅极上游的磁镜区.2 和3 号源在靠近栅极表面附近有少量高能电子(大于10 eV), 这是因为壁面附近鞘层的欧姆加热效应.

图6 不同磁场结构下离子源的电子温度Fig.6.Electron temperature of ion source with different magnetic circuit.

3.3 壁面电流

图7 不同磁场结构下离子源内各壁面的电流密度统计Fig.7.Current density statistics on each surface of ion source with different magnetic circuit.

4 结 论

本文采用全粒子PIC/MCC 方法模拟了2 cm ECRIT 离子源内部的等离子体与壁面电流.结论如下:

(1)电子被磁镜约束在天线ECR 区附近, 使得等离子体也分布在ECR 区附近.并且ECR 区与天线的相对位置会共同影响离子源的引出束流.当ECR 区位于天线下游时, 更有利于离子引出.从等离子体分布结果中发现栅极前等离子体主要集中在外磁环与轴线之间, 推测当保持孔数不变时, 而减少壁面附近的栅极孔分布范围, 将栅极孔集中分布在轴线中心附近, 能提高引出束流以及推进剂利用效率.

(2) ECR 区位于天线下游时, 等离子体主要分布在天线下游而且栅极前等离子体密度高, 但离子对天线溅射严重.这样的磁场设计虽然能够提高引出的离子束流, 但会使天线工作寿命缩减.不同壁面的电流分布规律为: 天线壁面 > 磁体端面 > 腔体壁面.由此可见离子源天线是壁面电流的主要来源, 有可能是影响工作寿命的主要因素.未来设计推力器时, 应权衡引出性能和寿命来选取合适的磁场结构.