微型电子回旋共振离子推力器(ECRIT)具有无阴极烧蚀、结构简单、易于启动等特点, 可用于小型航天器的飞行控制^[1]或材料处理^[2]等方面. 目前日本已经研究出直径2 cm的μ1型ECRIT, 并于2014年分别用在50 kg HODOYOSHI-4和70 kg PROCYON的微小卫星上^[3,4]. 国内西北工业大学已经研究出直径10和2 cm的ECRIT样机, 对其离子源进行了束流引出实验、朗缪尔探针诊断实验和数值模拟等研究工作^[5-9].

ECR中和器是ECRIT的关键部件, 其作用是中和离子束, 避免推进器和航天器表面电荷堆积, 是限制ECR推力器使用寿命的关键部件^[10], 需要从数值模拟和实验两方面开展细致的研究. 但是实验研究^[7,11,12]难以观察到ECR中和器内部的等离子体瞬时变化过程, 这就需要借助数值模拟方法, 揭示中和器内等离子体参数、内部的碰撞过程以及壁面电流损失随磁路结构的变化规律. 文献[13]利用三维PIC/MCC模拟研究了微型ECR中和器的电子引出机理. 其结果表明E × B是电子引出的重要参数, 认为方位角电势和径向磁场对电子的引出至关重要. 文献[14]对微型ECR中和器进行了三维PIC/MCC模拟. 模拟结果得到的孔板上电流密度分布与实验测量结果比较吻合. 此外, 数值结果表明等离子体静电场对电子的引出起主导作用, 虽然静电场阻碍电子引出, 但是引出电子的轨迹是由鞘层和静磁场共同作用下形成的. 文献[15]通过三维PIC/MCC模拟研究分析了孔板形状和磁场结构对微型ECR中和器的电子的引出效率的影响. 模拟结果表明, 孔板形状对放电特性和电子引出效率没有显著影响. 但在一种新型磁场结构中, 电子引出效率达到了1.5倍, 这一改善是由于引出孔板下游表面上电子回流和放电室内壁上电子损失的减少.

ECR中和器内部磁路和ECR区控制着等离子体分布特征, 腔体长度影响引出孔附近磁力线特征, 是控制电子引出效率和ECR等离子体性能的重要因素. 然而国内还少见这方面的研究, 因此开展不同磁路结构和腔体长度对等离子体参数的影响研究对研究分析中和器的性能十分重要. 与上述文献[13-15]对比, 这些国外学者主要研究中和器的电子EXB引出机制和外加线圈磁场对电子引出效率的影响, 计算中未考虑收集板电压和引出的电子电流. 然而本文考虑了电子引出以及收集板电压的影响. 模拟结果揭示了微型ECR中和器磁路和腔体长度对等离子体分布和引出电子电流的影响, 以及带电粒子损失和输运的变化规律. 另外, 本文研究不同磁路结构下ECR区位置变化和几何结构对中和器引出的影响. 模拟结果表明, ECR区位于天线上游时, 引出孔附近的电势阱会更大, 从而阻碍了电子引出. 除了ECR区与天线位置的的影响, 引出孔前的磁场强度和磁力线分布对电子引出同样有重大影响.

为此本文采取PIC/MCC方法对微型ECRIT中和器等离子体分布及电子引出进行模拟研究. 鉴于在真实的介电常数和离子质量条件下开展三维数值模拟, 计算量大、耗时长并且对计算机性能要求苛刻. 例如文献[16]可知, 一个三维PIC/MCC计算离子源的算例在4核3.0 GHz处理器和12 G内存工作站上需要一周时间才能得到结果, 而本文的二维轴对称模型算例在普通台式机上仅需要40 h. 另外, 国外文献[13-15]采用三维PIC/MCC方法是因为他们主要研究电子E × B引出机制, 这需要三维计算才能准确体现. 然而, 本文主要研究不同磁路结构和腔体长度对电子引出的影响, 采用二维轴对称PIC/MCC模型能展现其影响机制, 并节省计算时间. 本文第2节描述数值计算模型, 第3节给出不同磁路和腔体的模拟结果, 第4节讨论分析磁场结构和腔体长度对电子引出影响的机制, 最后是结论. 本文研究将为揭示并理解微型ECR中和器内的物理过程和性能优化提供指导.

微型ECR中和器主要组件为腔体、环形天线、磁轭、内外环形永磁体和引出孔板等, 如如图1所示. 频率为4.2 GHz的微波通过环形天线馈入中和器, 氙气通过底部磁轭上的进气孔被注入放电室, 电子在ECR区被微波振荡持续加热成高能电子并与原子碰撞, 从而生成等离子体. 中和器单独实验时, 通过对收集板施加正电压将电子引出; 而当中和器与离子源协同工作时, 凭借离子羽流与中和器之间的电势差引出电子.

如图1所示, 定义H₁和H₂分别为外磁环和内磁环高度, W₁和W₂分别为外磁环和内磁环宽度. L₁和L₂定义为表示天线-内磁环和天线-屏的间距^[7]. 本文中不同磁路模拟时保持L₁和L₂值相同.

本文的中和器模型能简化为二维轴对称模型, 模拟区域以图1中该平面对称边界上的点划线为对称轴, 轴向和径向的正方由图中z, r矢量确定. 简化为二维模型的理由有如下几点. 一是由于中和器和离子源内部结构相似, 根据文献[8, 17]可认为中和器内的等离子体分布同样近似轴对称分布, 静磁场和高频电场近似为轴对称. 二是本文中和器模拟的电子密度分布结果与文献[16]三维离子源模拟的电子密度分布结果对比, ECR区都在天线下游时, 等离子体密度分布大致一样, 这也说明了本文二维轴对称计算模拟结果的合理性. 三是本文在柱坐标系下, 二维位置空间2D (r, z)和三维速度空间3D (V_r, V_θ, V_z)中进行粒子模拟. 保留V_θ方向速度, 也考虑了径向力V_θ × B的存在以及它对总能量有贡献, 这使得模拟结果更接近三维计算.

由文献[13]可知, 三维计算结果同样需要在二维轴对称截面表示电子引出, 故采用二维模拟能加快计算速度, 节省时间和成本, 并且计算结果处理方便. 并且二维的计算网格尺寸可以更小(本文采用的网格尺寸为0.025 mm), 而国外三维模拟采用0.2 mm的网格尺寸, 这限制等离子体放电的定性特性. 综上所述, 三维模拟结果能更准确、真实反映中和器引出机制, 但其由于三维计算量大、耗时长, 对计算机性能要求高, 并且其网格尺寸限定了功率吸收和实验无法匹配. 如果想进一步提升三维计算精度和速度, 需采用更细的网格、更多的宏粒子和并行计算方法, 目前实现这个目标有一定困难. 而本文研究目标主要是磁场构型和腔体长度对中和器束流引出大小的影响, 简化后的二维轴对称模型能满足要求, 并且大大提高计算速度, 节省计算时间.

将微型ECR中和器内部的电磁场主要由以下三个模块组成: 静电场、静磁场和微波高频电场^[9]. 静电场通过有限体积法离散泊松方程$\nabla \cdot \nabla \varphi = $$ - {\rho }/{{\varepsilon }_{0}}$, 生成矩阵方程, 然后选用SuperLU库^[18]直接求解矩阵得到电势. 电场为电势场的负梯度, 离散条件下采用梯度构造方法^[19]计算电场. 对称边界和羽流区上边界为第二类条件为${\partial \varphi }/{\partial r}=0$, 其余边界为第一类条件. 除了收集板电势为不同电压值, 其余内部边界电势均为0, 电势边界条件如图2所示. 静磁场采用COMSOL计算得出, 作为背景场; 微波高频电场采用COMSOL计算, 预先得到高频场的幅角和相位, 计算时可忽略等离子体电流, 因为在较低的微波功率下等离子体电流对微波电场影响可忽略不计. PIC程序计算时, 令每个时间步上高频场的瞬时值与静电场叠加, 用于电子运动推进; 离子只靠静电场推动, 是由于离子质量较大而忽略高频电场的作用.

中和器以氙气为工质, 把氙原子当作背景粒子, 均匀分布在计算域中. 考虑的粒子间碰撞类型如下所示:

a) (弹性碰撞)$ \mathrm{e}+\mathrm{X}\mathrm{e}\to \mathrm{e}+\mathrm{X}\mathrm{e} $,

b) (激发碰撞)$ \mathrm{e}+\mathrm{X}\mathrm{e}\to \mathrm{e}+{\mathrm{X}\mathrm{e}}^{*} $,

c) (电离碰撞)$ \mathrm{e}+\mathrm{X}\mathrm{e}\to \mathrm{e}+{\mathrm{X}\mathrm{e}}^{+}+\mathrm{e} $,

d) (电荷交换碰撞)$ {\mathrm{X}\mathrm{e}}^{+}+\mathrm{X}\mathrm{e}\to \mathrm{X}\mathrm{e}+{\mathrm{X}\mathrm{e}}^{+} $,

e) (库伦碰撞)$ {\mathrm{X}\mathrm{e}}^{+}+{\mathrm{X}\mathrm{e}}^{+}\to {\mathrm{X}\mathrm{e}}^{+}+{\mathrm{X}\mathrm{e}}^{+} $, $\mathrm{e}+ $$ {\mathrm{X}\mathrm{e}}^{+}\to \mathrm{e}+{\mathrm{X}\mathrm{e}}^{+}$, $ \mathrm{e}+\mathrm{e}\to \mathrm{e}+\mathrm{e} $.

碰撞截面采用Lxcat^[20]的截面数据, 上述碰撞处理和粒子运动方法详见文献[9, 21, 22]. 图3为计算流程图. 图3中E_EM表面高频电场瞬时幅值, 随着计算时间t发生变化, 即$ {E}_{\mathrm{E}\mathrm{M}}={E}_{0} {cos}\left(2{\text{π}}ft+\phi \right) $, 其中E₀表示COMSOL计算得到的高频电场幅值, f表示微波频率, φ表示高频电场的相位角; B_ST表示COMSOL计算得到的静磁场强度.

在表1三种磁路下进行微型ECR中和器电子引出模拟. 利用COMSOL软件仿真得到三种磁路下中和器内部与羽流区磁场分布, 且ECR区与天线位置关系如图4所示. 随着H₁和W₁增加, ECR区逐渐远离天线并靠近屏栅, 同时引出孔板附近的磁通密度也逐渐增强.

中和器的输入微波频功率和氙气流量分别为1 W, 0.3 sccm. 根据公式${p}_{0}=\dfrac{4{k}_{\mathrm{B}}{T}_{\mathrm{g}}{Q}_{0}}{{v}_{\mathrm{g}}{A}_{\mathrm{g}}}$和$n=\dfrac{{p}_{0}}{{k}_{\mathrm{B}}{T}_{\mathrm{g}}}$估算得到的中和器内部的中性原子压强和密度分别为0.537 Pa和$ 1.3\times {10}^{20}\;{\mathrm{m}}^{-3} $. 羽流区的中性原子密度设置为内部原子密度的0.1%. 其中$ {k}_{\mathrm{B}} $为玻尔兹曼常数, $ {T}_{\mathrm{g}} $为气体温度(等于室温298 K), $ {Q}_{0}=4.48\times {10}^{17}Q $(原子/s, Q单位为sccm), $ {v}_{\mathrm{g}} $为氙原子的热速度(对应298 K时约为218 m/s), $ {A}_{\mathrm{g}} $为开孔面积(约为18.85 mm²).

在真实的介电常数和离子质量条件下, 模拟时间步长Δt和网格步长Δx分别取为5 × 10^–12 s和2.5 × 10^–5 m^[9]. 若以离子在推力器中迁移速度表示扰动在中和器内的传播速度, 则以扰动在离子源内的传播时间为模拟总时长. 2 cm ECRIT中和器内离子迁移速度为10⁴—10⁵ m/s, 则离子源内等离子体演化特征时间为10^–6—10^–7 s. 当计算时间步长在250—300万步以内的粒子数相对误差变化均低于0.1%时, 计算达到稳态, 模拟总时间达到1.5 × 10^–5 s. 设定初始电子和离子温度分别为2 eV和0.5 eV, 初始等离子体密度为$ 1\times {10}^{16}\;{\mathrm{m}}^{-3} $.

对不同磁路下中和器模拟时, 保持其他几何参数一致, 即L₁ = 3.9 mm, L₂ = 0.7 mm, 其模拟结果与实验结果^[7]如表2所列. 从表2可以看出, 不同磁路下中和器在引出相同大小1 mA的电子电流所需收集板电压不同. 模拟中采用与实验相同的收集板电压, 模拟结果中的收集电流与实验结果电流的相对误差小于15%.

不同磁路的电子密度分布的计算结果如图5所示. 从图5中可以看出, 结构1, 2和3的等离子体集中分布在ECR区附近^[9,14], 并且电子密度峰值分别是$ 1.8\times {10}^{17}\;{\mathrm{m}}^{-3} $, $ 1.2\times {10}^{17}\;{\mathrm{m}}^{-3} $, $ 1.3\times {10}^{17}\;{\mathrm{m}}^{-3} $. 三种结构的等离子体密度峰值区面积大小为: 结构1 > 结构3 > 结构2. 结构2的电子密度峰值区域最小, 是因为ECR区与天线相交, 从而影响在磁镜场运动的电子从ECR区中获能, 所以峰值密度区域最小. 结构1的等离子体密度峰值位于天线上游, 因为ECR区位于较小的磁镜区, 电子被ECR区加热后能有效地在磁镜内发生电离碰撞, 所以内外磁环间的等离子体密度较高, 而迁移扩散到引出孔的过程中由于天线阻碍或损失在其他壁面, 导致引出孔附近的等离子体密度较低. 结构3的ECR区面积大于结构2, 并且ECR区位于天线下游磁镜区, 所以结构3的峰值大于结构2. 但部分未电离的高能电子容易扩散并被收集板电势引出, 所以等离子体密度低于结构1.

不同磁路下孔板内侧附近的电子密度范围分别是$ 8\times {10}^{14}—3\times {10}^{15}\;{\mathrm{m}}^{-3} $, $1.8\times {10}^{15}—4\times {10}^{15}\;{\mathrm{m}}^{-3}$和$ 2\times {10}^{15}—4.2\times {10}^{15}\;{\mathrm{m}}^{-3} $. 孔板外侧电子羽流的密度范围是$ 2\times {10}^{14}—6\times {10}^{14}\;{\mathrm{m}}^{-3} $. 图5中羽流区黄色虚线区域为引出电子, 但由于模拟的引出电流较小, 没有文献[23]中羽流区的电子密度明显.

如图6所示, 不同磁路结构的电势分布和等离子体密度分布基本一致, 结构1, 2和3内部峰值电势出在ECR区附近, 分别为15.8 V, 10.0 V和9.5 V. 而且从图6中可以明显看到引出孔附近存在电势阱.

图7为结构2(L₁ = 3.9 mm)引出不同电流的模拟结果与实验结果曲线. 模拟得到的引出电子电流的I-V曲线与实验结果不同. 这是因为当ECR中和器的引出电流较低时, 电离速率较低, 故实际上参与电离的中性粒子数目较少. 然而, 当ECR中和器的引出电流较高时, 被引出的电子越多, 通过离子轰击中和器壁上的复合形成的中性粒子越多, 使得中性粒子的密度增加. 因此, 本文根据实验电引出电子电流的大小, 调整不同收集板电压下背景中性粒子的密度, 使模拟得到的I-V曲线接近实验曲线. 在未来的PIC/MCC求解中, 将考虑中性粒子的动态调整. 相同电压下, 模拟得到的电子电流与实验结果误差均小于15%, 这意味着通过调节中性原子的密度, 可以使模拟电流可以与实验结果相吻合.

根据文献[7], 选用性能优异的结构2进行不同腔体长度模拟研究, 即, L₂尺寸与上文一致, L₁增加0.4 mm. 图8为相同计算条件下, 结构2的腔体加长0.4 mm后的模拟结果. 其孔板内侧的电子密度和离子密度范围分别是$2\times {10}^{15}—4.5\times $$ {10}^{15}\;{\mathrm{m}}^{-3}$和$ 2\times {10}^{15}—5\times {10}^{15}\;{\mathrm{m}}^{-3} $. 孔板外侧的电子羽流的密度和离子密度范围分别是 $2\times {10}^{14}— $$ 6\times {10}^{14}\;{\mathrm{m}}^{-3}$和$ 2\times {10}^{14}—8\times {10}^{14}\;{\mathrm{m}}^{-3} $. 图8(a)中羽流区黄色虚线区域为引出电子, 但由于模拟的引出电流较小, 所以羽流区的电子密度在$ 1\times {10}^{14}\;{\mathrm{m}}^{-3} $. 与原腔体尺寸相比, 同样大小的收集板电压15 V能引出1.98 mA电子电流. 实验中功率流量条件相同下, 引出2 mA电流需要收集板电压为13.8 V, 证实了模拟结果可靠.

首先分析不同磁路下电子引出性能差异的原因. 中和器的孔板附近的鞘层阻碍电子引出, 即电子需要一定能量穿越孔板上游的电势阱才能被引出. 图9表示了孔中心轴线上放电室到羽流区电势变化, 即图6中引出孔中心黄色直线上的电势变化. 不同磁路下到达引出孔板下游的电子需要不同能量才能越过电势阱, 如图9所示其电子能量大小为: 结构1 > 结构2和3. 这解释了为什么结构1的电子引出性能最差. ECR区位于天线上游时会影响电势阱大小, 从而影响电子引出. 然而, 结构2和3的电势阱大小相似, 这表明当ECR区位于下游或与天线切割时, 可忽略电势阱对电子引出性能的影响. 所以需要进一步从磁场构型上分析结构2和3电子引出性能差异的原因.

根据磁场分布结果(图4), 中和器内电子引出的轨迹可分为两类: 第一类是电子沿磁力线运动, 第二类是电子横越磁力线运动. 比于第二类运动, 电子更易沿磁力线运动到收集板^{[24, 25]}. 由图4可见, 结构1收集到的电子主要是靠横越磁力线, 结构2和3收集到的电子主要是沿磁力线运动和少量电子横越磁力线运动. 这也解释了结构2和3电子引出性能优于结构1. 在磁化等离子体中, 当电子回旋频率$ {\omega }_{\mathrm{c}} $远大于碰撞频率ν时, 垂直于磁场的经典扩散系数为$ {D}_{\perp } $, 可以写成${D}_{\perp }=\dfrac{{k}_{\mathrm{B}}T}{m\nu }\dfrac{1}{{\omega }_{\mathrm{c}}^{2}{\nu }^{-2}}\approx $$ \dfrac{{mk}_{\mathrm{B}}T\nu }{{q}^{2}}\dfrac{1}{{B}^{2}}$. 由扩散公式可知$ {D}_{\perp } $与B平方成反比, 说明引出孔附近磁场越强, 电子横越磁场到羽流区的可能性越低. 由图4可知结构3的引出孔板附近磁场强度大于结构2, 所以结构2内部电子横越磁力线的扩散速度大于结构3. 这解释了结构2电子引出性能优于结构3.

为了证实增加沿磁力线运动的电子数目, 对结构2进行了不同腔体长度的影响研究. 图10为相同磁路结构下不同腔体长度时引出孔板上各表面损失的电子电流密度. 结合图4和图10, 发现增加腔体长度后, 孔通道表面上电子损失减小, 这表明减少了第二类磁力线越过引出孔. 孔板内侧表面上电子损失减小, 这表明更多电子被约束在磁镜中, 从而减少了孔板上的电子损失, 并且. 模拟结果与实验结果都表明相同磁路下, 增加腔体长度能增加引出的电子电流. 根据上述分析有两种原因: 一是引出孔板表面的电子损失减小; 二是增加腔体长度减小了引出孔前的磁场强度从而增大了电子横越磁力线的扩散速度.

本文采用二维轴对称PIC/MCC模型, 计算了不同磁路和腔体下微型ECRIT中和器的等离子体分布和电子引出过程, 得到结论如下:

1) ECR区与天线位置的不同, 会影响等离子体分布以及电势分布. 模拟结果发现ECR区位于天线上游时, 引出孔附近的电势阱会更大, 从而阻碍了电子引出, 使得收集板需要更高电压才能引出相同大小的电子电流.

2) 除了ECR区与天线位置的影响, 引出孔前的磁场强度和磁力线分布对电子引出同样有重大影响. 从模拟与实验结果得知, 应尽量保证引出孔附近磁力线平行于轴线, 并且适当的增加腔体避免2电子损失在引出孔板, 有助于提高电子引出.

本文PIC/MCC模型存在一定局限性和误差, 未考虑中性原子的流动. 因为真实情况下中和器内部的原子密度是不均匀的, 并且腔体长度增加时, 原子密度会降低. 未来, 首先需要思考如何进行实验测量中性原子分布, 其次模型考虑中性原子流动后, 相互对比本文模拟结果才能进一步定量评估误差大小, 以及三维模拟下能更真实地反映中性原子从引出孔泄露到真空环境中, 从而得到更准确的中性原子密度. 而且二维轴对称模型下将引出孔等效会带来一定误差. 由于是二维模拟, 所以引出孔未能未能表示出六个孔数. 本文是将引出孔等效为圆环, 模拟得到的引出束流需进行换算, 即乘以六个引出孔面积与圆环面积之比. 虽然本文中也提及了引出电子束流的模拟结果与实验结果相差15%以内, 但如果需要进一步提高计算结果准确性, 需要更高的计算条件以及考虑中性原子流动和进行三维模拟.