离子能量分析器（IEA）也称为阻滞势分析器，作为一种就位探测等离子体能量的重要工具，被广泛应用于电离层探测卫星中^[1-5]。由于电离层离子能量较低，在eV量级，目前在地面通过实验还无法定量地对仪器的测量特性进行研究。仿真软件因为能够仿真出稳定的低能量离子源，而且可以观察到离子在仪器内部的运动情况，进而对一些实验现象进行分析，成为了一种研究离子能量分析器的有力工具。

离子能量分析器不仅能用于测量离子速度，也可测量离子密度和温度，此前的研究多针对离子温度和离子密度的测量，对离子速度测量的研究比较少。Goldan^[6]，Troy^[7]，Knudsen^[8]在2000年之前先后研究了不均匀的电场分布、有限的栅网透过率和收集极面积、有限的栅网半径和厚度引起的离子密度和离子温度的测量误差，不涉及离子速度；Chao于2000年研究了离子在离子能量分析器中的运动^[9]，于2003年利用仿真程序研究了栅网对离子温度和速度测量产生的影响^[10]。囿于当时的计算资源，他只建立了一个缩小版的模型；Klenzing于2009年利用ANSYS分析了不同栅网的对齐方式和栅网层数对离子温度和速度测量的影响^[11]。但在他的模型中将阻滞栅网以外的其他栅网理想化，且第二层阻滞栅网的网孔密度比第一层减小了一半；冯宇波等人于2010年利用SIMION研究了等离子体鞘层和电场畸变引起的测量误差^[12]。他在建模时将栅网设为无限大，因此忽略了撞击到器壁上而损失掉的离子的影响；Fisher^[13]于2016年使用支持向量回归的方法分析了离子能量分析器的测量数据，但是只考虑了等离子体鞘层引起的测量误差。

上述的研究工作还存在以下问题：建立的模型比较理想化；对各种误差因素考虑得不够全面；只考虑了离子均值速度这一点的误差，没有对整个能量分布的误差进行分析。基于这些问题，本文利用仿真软件COMSOL对离子能量分析器的测量特性开展了进一步的研究：仿真模型尽可能还原了离子能量分析器的真实结构；通过离子透过率仿真分析确定了一种较优的仪器设计方案；仿真分析了电场畸变、等离子鞘层、栅网对齐方式和离子温度等多种误差因素对离子能量分布测量结果的影响。这些工作对离子能量分析器的结构优化，实验测试和误差分析与校正具有很好的指导意义。

1.   离子能量分析器工作原理

离子能量分析器的典型结构如图1所示^[14]，图中G1～G4为四层栅网，C为收集极。外部离子通过栅网G1进入仪器，栅网G1接仪器地，以屏蔽内部所加电压对外部离子的影响。离子进入仪器后，遇到两层阻滞栅网G2和G3，阻滞栅网加正的扫描电压，通过栅网之间形成的电场对不同能量的离子进行筛选，只有高于电场能量的离子能够通过阻滞栅网到达G4。栅网G4相对于仪器地加负偏压，用来阻止前面进入的电子到达收集极，也可抑制收集极受高能离子激发的二次电子的逸出。穿过G4的离子最终到达收集极形成电流。

图  1  离子能量分析器典型结构

Figure  1.  Typical structure of ion energy analyzer (IEA)

下载: 全尺寸图片 幻灯片

如果离子能量分析器的轴向平行于z轴放置，则收集极电流I与阻滞栅网所加扫描电压$U_{\rm{RV}}$满足以下关系

式中：$q$为离子所带电荷量；$A$为收集极面积；$N$为离子密度；${v_{x}},{v_{y}},{v_{z}}$为离子速度沿三个坐标轴的分量；$\chi ({{v}},{U_{{\rm{RV}}}})$表示速度为${{v}}$的离子在阻滞电压为${U_{{\rm{RV}}}}$时的透过率；${\cal D}({{v}})$为离子速度分布函数。注意对${v_x}$和${v_y}$的积分是在整个实数域内，而对${v_{\textit{z}}}$的积分只包含了${v_{\textit{z}}} \geqslant {\rm{0}}$的部分，排除了沿−z轴运动即远离仪器的离子。

理想情况下，空间中的离子速度满足麦克斯韦分布，${v_x}$和${v_y}$是在整个实数域积分，满足归一化条件，三重积分化简为一重积分。离子沿z轴的速度分布是相对于航天器运动的一维麦克斯韦分布函数

式中：${v_{{\rm{th}}}} = \sqrt {2kT/m}$称为最概然速率，与离子质量m、玻耳兹曼常数k和离子温度T有关；v_z0为离子群相对于探测卫星的速度，以卫星的速度占主导。

假设离子能量分析器对离子能量具有理想的选择性，即大于电场能量的离子都能达到收集极，而低于电场能量的离子均被阻止，则$\chi ({{v}},{U_{{\rm{RV}}}})$是一个理想的阶梯函数，可写为以下形式

式中：$\chi$为栅网的透过率。对于单层栅网，透过率等于栅网的光学透过率，即网孔面积占栅网总面积的百分比；对于多层栅网，透过率取决于栅网的对齐方式，最大透过率等于单层栅网的光学透过率，最小透过率等于单层栅网光学透过率的乘积。

将式（2）和式（3）代入式（1），同时假设进入离子能量分析器的离子只有一种，可得

式中：$\kappa = ({v_{{\textit{z}}0}} - \sqrt {2q{U_{{\rm{RV}}}}/m} )/{v_{{\rm{th}}}}$；${\rm{erf}}( \cdot )$为误差函数，是以e为底的指数函数积分得到的。

对式（4）进行微分，即可得到离子能量分布函数

实际测量时，通过记录每个扫描电压点对应的收集极电流，得到一组离散的I-V数据。按照式（5）的形式对数据进行拟合，得到I-V曲线，对曲线进行微分，即为待测的离子能量分布。图2是理想情况下离子能量分析器的测量曲线。

图  2  理想情况下离子能量分析器测量曲线

Figure  2.  Ideal IEA measurement curve

下载: 全尺寸图片 幻灯片

2.   仿真工具及计算模型

COMSOL全称COMSOL Multiphysics，是一款多物理场仿真软件，能够对多物理场耦合问题进行建模求解。对于此次需要研究的问题来说，进入仪器的带电离子的运动会改变电场，电场的改变又反过来影响到离子的运动。COMSOL专门提供了粒子-场相互作用接口来对此类问题进行求解，这保证了仿真结果的准确性。

仿真中选择COMSOL版本为5.4，因为计算过程中需要大量的计算资源，使用COMSOL的客户端-服务器模式进行操作，把繁重的计算任务交给服务器完成，客户端只用来下达命令。使用的服务器有128 G内存和两颗英特尔至强E5-2660系列处理器，操作系统为Ubuntu 18.04。

模型的网格剖分影响着求解结果的精度。仿真中采用了多种网格剖分方式，对模型分块进行了剖分。栅网作为模型中的关键部件，在它的上下两侧加密了网格。最终生成的网格质量为0.83，比较高的网格质量，也能保证求解过程尽快地收敛。仿真同时对模型的默认求解器进行了修改，重新设置了求解的时间步长，增加了停止条件，使所有的离子不再运动时计算过程停止。

为了进一步加快求解速度和节省计算资源，根据模型的旋转对称性，在不产生新的误差的前提下，使用COMSOL提供的周期性边界条件^[15]，对模型进行了一步简化。对于网孔为正方形的栅网只建模其中的四分之一，对于网孔为正六边形的栅网只建模其中的六分之一。模型简化过程如图3所示。

图  3  模型简化过程

Figure  3.  Model simplification processes

下载: 全尺寸图片 幻灯片

通过对国际参考电离层IRI-2016^[16-17]导出的电离层探测卫星轨道区域内的数据进行计算，确定了模型的一些尺寸。栅网网孔孔径应小于两倍的德拜半径^[18]，取为1 mm。考虑到栅网的透过率要尽可能高和实际的工程条件，栅网丝径定为0.1 mm，厚度同样为0.1 mm。进入离子能量分析器的离子沿x轴和y轴有一定的横向速度，为了尽量避免离子在到达收集极之前打到内壁上，通过计算确定仪器的内壁为40 mm^[19]。电离层探测卫星轨道区域内，原子氧离子O^＋占绝大多数^[20]，仿真中只考虑O^＋。

3.   离子能量分析器结构仿真分析

式（3）假设透过率曲线是一个理想的阶梯函数，即离子能量分析器对离子能量具有理想的选择性。事实上，选择的非理想性是造成测量误差的主要因素。为了能找到一种透过率曲线最接近阶梯函数的设计方案，以减小测量误差，我们提出了三种不同的仪器设计方案，并对这三种方案进行了离子透过率仿真分析。

方案A：模型结构如图1所示。栅网网孔为正方形，孔径1 mm，丝径0.1 mm，栅网距离4 mm，仪器内径40 mm。

方案B：将方案A的栅网网孔改为正六边形，其他保持不变。

方案C：将方案A的栅网距离改为7 mm，其他保持不变。

仿真中离子能量分析器外壳、栅网G1和收集极接地，栅网G4相对于外壳加−12 V偏压，保持阻滞栅网的电压为5 V不变，离子沿x轴和y轴的入射能量为0，沿z轴的入射能量从4 eV增加到7 eV。仿真开始时在模型开口处入射10 000个离子，入射位置在开口处随机分布。考虑到不同的入射位置可能会对结果产生影响，为每次仿真设置了相同的随机数种子，保证每次仿真随机的位置相同。得到的三种方案的离子透过率仿真曲线如图4所示。图中的阶梯函数为理想情况下的透过率曲线。理论最大透过率为单层栅网的透过率，即${\left( {1/1.1} \right)^2} = 82.64{\text{%}}$。理论最小透过率为四层栅网透过率的乘积，即${\left( {82.64{\text{%}} } \right)^4}$。

图  4  三种设计方案的离子透过率曲线

Figure  4.  Ion transmission curve of three schemes

下载: 全尺寸图片 幻灯片

从透过率大小来看，方案A、方案B和方案C在低能量段的离子透过率曲线相似。在高能量段，方案C的透过率要比方案A以及方案B大。这可以通过栅网产生的畸变电场对离子的聚焦作用解释。增大栅网距离后的电场畸变减小，小的电场畸变足以影响能量较低的离子的运动，方案C表现出和方案A以及方案B相似的透过率。随着离子能量增大，小的电场畸变逐渐不能够影响到离子的运动，方案C的透过率要比方案A以及方案B增大得快。从透过率曲线的变化来看，方案A和方案B的波动较小，更接近一个阶梯函数。为了能够区分方案A和方案B的优劣，我们做了进一步的仿真。

在进一步的仿真中，不仅赋予离子沿z轴4～7 eV变化的能量，也赋予它们沿x轴和y轴T=1000 K时的随机速度，速度分布满足麦克斯韦分布，得到的新的离子透过率曲线如图5所示。可以看到方案B的离子透过率曲线更接近于理想的阶梯函数。方案A在离子以一定的夹角入射时，透过率的波动太大。因此，我们最终选择了方案B。从这两次的结果还可以看出，离子能量在接近电场能量时会出现一个极大值，在高于电场能量10%左右会出现一个凹陷，这是由电场畸变改变了离子的运动轨迹引起的。

图  5  方案A和方案B的离子透过率曲线

Figure  5.  Ion transmission curve for scheme A and scheme B

下载: 全尺寸图片 幻灯片

4.   离子能量分析器误差仿真分析

电场畸变、等离子体鞘层、栅网对齐方式和离子温度是造成测量误差的主要因素。我们通过COMSOL计算离子在仪器中的运动轨迹，分析了离子受电场畸变的影响；通过改变阻滞栅网的对齐方式，分析了栅网对齐的影响；设置仪器为负电位，分析了等离子体鞘层对结果的影响；最后，通过测量不同离子温度下的I-V曲线，将测量得到的离子能量分布和输入的能量分布进行比较，综合分析了所有误差因素对测量结果的影响。

4.1   电场畸变的影响

图6是离子能量分析器内部的电场图。理想情况下的内部电场是沿z轴的均匀电场，而实际上在栅网网孔处存在严重的电场畸变。这种畸变会使低于电场能量的离子通过栅网，也会对靠近栅网的离子施加一个横向力，横向力会改变离子的运动轨迹，对离子产生聚焦的作用。

图  6  离子能量分析器内部电场图

Figure  6.  Internal electric field of IEA

下载: 全尺寸图片 幻灯片

图7是不同入射能量下的离子在仪器内部的运动轨迹。阻滞栅网加5 V固定电压，当入射能量为4.6 eV时，所有的离子都不能通过栅网。当入射能量增加到4.8 eV时，已经有一部分离子能够通过栅网。这说明电场畸变造成低于电场能量的离子也有可能通过栅网，解释了为什么图5中入射能量低于5 eV时透过率已经不为0。

图  7  不同入射能量下的离子运动轨迹

Figure  7.  Ion trajectories at different initial energies

下载: 全尺寸图片 幻灯片

图7（c）和图7（d）可以解释图4中入射能量为电场能量的10%左右时，离子透过率曲线出现凹陷的原因。栅网网孔处的电场畸变对离子施加的横向力使离子在运动过程中形成一个焦点，这个焦点随着离子入射能量的增加向后移动。离子在通过焦点后发散，当入射能量为电场能量的10%左右时，发散的离子有较大的概率撞到第三层栅网。因此，离子的透过率曲线出现凹陷。之后，随着离子能量的增加，离子的发散程度减小，有更多的离子到达收集极，透过率曲线开始上升。

4.2   栅网对齐的影响

由于栅网的网孔只有1 mm，而半径有40 mm，网孔比较稠密，在安装中很难保证每个网孔严格对齐。图8中展示了栅网的两种对齐方式：图8（a）的前一层栅网和后一层栅网网孔保持严格对齐，图8（b）的后一层栅网网丝的交界处在前一层栅网网孔的中心。可以想象，如果网孔严格对齐，离子在通过前一层栅网后有较大的概率能通过后一层栅网。而如果栅网未对齐，通过前一层栅网的离子将有极大的概率撞到后一层栅网。我们按照这两种对齐情况，改变两层阻滞栅网的对齐方式，研究了栅网对齐对测量结果的影响。

图  8  栅网对齐的两种情况

Figure  8.  Two cases of grid alignment

下载: 全尺寸图片 幻灯片

图9是两种对齐方式下的离子透过率曲线。可以看出，未对齐栅网的透过率不仅偏低，而且曲线已经严重偏离了阶梯函数形式，这会严重影响测量结果。这指导我们在安装离子能量分析器时，要尽可能地保证栅网对齐。

图  9  两种对齐方式下的离子透过率曲线

Figure  9.  Ion transmission curve in two alignments

下载: 全尺寸图片 幻灯片

4.3   等离子体鞘层的影响

等离子体在与离子能量分析器壁接触时，会形成非中性的电荷过渡区，称为等离子体鞘层^[21]。由于电子运动速度一般大于离子运动速度，电子会率先到达器壁而使器壁相对于等离子体呈现负电位，负电位的大小在电子动能的量级。电离层离子动能约为5 eV，因此，可以通过设置仪器电位为−2 V来研究等离子体鞘层对测量结果的影响。

图10为仪器电位为−2 V和0 V两种情况下，T=1 000 K时的离子透过率曲线。可以看出，离子能量在低于电场能量时，鞘层使离子透过率升高。而在略高于电场能量时，鞘层使透过率降低，透过率曲线凹陷处对应的入射能量变低。当入射能量增加到6.6 eV之后，鞘层几乎不再对透过率产生影响。相比于仪器电位为0 V的情况，仪器电位为−2 V时相当于离子在进入离子能量分析器前抬高了一定的能量，因此在入射能量低于电场能量时，鞘层使离子透过率升高。抬高的这部分能量使离子撞击到第三层栅网时的入射能量减小，因此透过率曲线凹陷处对应的入射能量变低。

图  10  不同仪器电位下的离子透过率曲线

Figure  10.  Ion transmission curve at different instrument potentials

下载: 全尺寸图片 幻灯片

图11为仪器电位为−2 V和0 V两种情况下，仿真得到的I-V曲线。电离层探测卫星的运行速度为7600 m/s左右，对应氧离子能量均值为4.79 eV。由于等离子体鞘层使能量略高于电场能量的离子透过率降低，能量低于电场能量的离子透过率升高。如果阻滞栅网加0～10 V的扫描电压，则在扫描电压低端，鞘层的存在会使收集极电流降低。在扫描电压的中间区域，鞘层会使收集极电流升高。相比于仪器电位为0 V的情况，仪器电位为−2 V时的I-V曲线向右下方偏移，使离子均值速度的拟合结果增大。图11中仪器电位为0 V时，离子均值速度的拟合结果为7626 m/s。仪器电位为−2 V时的拟合结果为7699 m/s。

图  11  不同仪器电位下的I-V曲线

Figure  11.  I-V curve at different instrument potentials

下载: 全尺寸图片 幻灯片

4.4   测量误差的综合分析

最后通过改变离子温度，对电场畸变和等离子体鞘层的影响做一个综合分析。由麦克斯韦速度分布函数式（2）可知，离子温度会影响离子的速度分布，进而影响离子进入离子能量分析器的角度。当离子温度很低时，离子沿x轴和y轴的释放速度很小，使离子的发散角也很小，离子基本沿z轴进入仪器。随着温度升高，离子的发散角将增大。由IRI-2016^[17]导出的数据可知，电离层探测卫星轨道区域内的离子温度主要在600～1800 K之间变化。我们取离子温度为500～2000 K范围内，以500 K为步进长度的四种不同离子温度，计算了不同温度下测量结果和理论值的误差。

图12是500 K和2000 K两种离子温度下，仿真得到的I-V曲线和微分得到的离子能量分布图。当离子温度比较低时，离子基本沿z轴进入仪器。沿z轴入射的离子透过率介于最大透过率和最小透过率之间，使收集极的电流比理论值要高，如图12（a）所示。随着温度升高，离子束的发散角增大使离子透过率低于最小透过率（这与图9中未对齐的栅网离子透过率偏低的原因相同），收集极收集到了比理论值低的电流，如图12（b）所示。

图  12  不同温度下的I-V曲线和离子能量分布

Figure  12.  I-V curve and ion energy distribution at different temperatures

下载: 全尺寸图片 幻灯片

Chao^[11]，Klenzing^[12]和冯宇波^[13]的论文中只对均值速度的误差进行了计算，这里我们使用均方根误差（RMSE）来计算整个离子能量分布曲线的误差，使用的计算公式如下

式中：$v_z^{\left( {10} \right)}$表示离子能量为10 eV时对应的离子速度；${\cal D}{({v_z})_{{\rm{fit}}}}$为拟合得到的离子能量分布；${\cal D}{({v_z})_{{\rm{input}}}}$为输入的离子能量分布；$\left| \cdot \right|$为绝对值符号。

表1中列出了不同温度下计算得到的均方根误差和拟合的离子均值速度。从表中可以看出，离子均值速度最大相对误差为2.77%，离子能量分布最大均方根误差为0.73 eV。这说明，虽然电场畸变、等离子体鞘层和离子温度会对测量结果产生影响，但这些影响并不大，尤其是在离子温度比较低时。通过仿真结果来看，我们的设计方案可以较为准确地测量出离子能量分布。

表  1  不同温度下离子速度拟合结果

Table  1.  Ion velocity fitting results at different temperatures

temperature/K	average velocity/（m·s−1）	relative error/%	${\varepsilon _{\rm{RMSE}}}$/eV
500	7593	−0.09	0.05
1000	7626	0.34	0.11
1500	7657	0.75	0.07
2000	7811	2.77	0.73

下载: 导出CSV 

| 显示表格

5.   结　论

由于电离层离子能量较低，地面产生这个能量段的稳定离子源比较困难，目前通过实验对离子能量分析器的测量特性进行研究存在许多问题。本文在前人研究的基础上，从离子能量分析器测量离子能谱的角度出发，使用多物理场仿真软件COMSOL对仪器的测量特性进行了更进一步的研究。

首先介绍了离子能量分析器的工作原理，对收集极电流公式进行了详细的推导，这有助于分析不同情况下I-V曲线的变化。为了得到最优的仪器结构设计方案，提出了三种不同的设计方案。之后，通过仿真分析，选择出了一种透过率曲线最接近理想阶梯函数的方案作为最终的设计参考。理论上，这种方案测量结果的误差最小。最后，分析了电场畸变、等离子体鞘层、栅网对齐方式和离子温度对测量结果的影响。通过比较在多种误差因素影响下的测量结果和理论值的差距，得出我们的方案在电离层探测卫星轨道区域的离子温度变化范围内，可以较为准确地测量离子能量分布。

后续将利用建立好的地面模拟原子氧离子源，对真实的离子能量分析器开展各种实验测试，以观测仪器在实际环境中的性能，并与仿真中的结果做对比分析。