王 波， 张 根， 张文杰， 祁 超， 严 辉

(1.北京工业大学材料与制造学部新型功能材料教育部重点实验室， 北京 100124;2.北京工业大学材料与制造学部固体微结构与性能北京市重点实验室， 北京 100124)

受控热核聚变一直是社会各界持续关注的能源危机解决办法. 通过受控核聚变可以获得无穷尽的清洁能源来推动社会可持续发展一直是各国科学家追求的最终目标[1]. 但是，要将热核聚变反应堆实现工业应用或商业应用，仍旧存在很多有待攻克的难题. 其中关键问题之一就是选择一种合适的面向等离子体材料(plasma-facing material，PFM)[2-3]. 为实现对PFM的性能检测，实验室模拟聚变环境起着非常重要的作用. 在各种实验模拟中，直线等离子体装置(linear plasma device，LPD)能够产生受约束的强等离子体流，具有用途广泛、可靠性强、操作灵活和成本较低的优势，因此各种LPD已被用来研究等离子体与材料的相互作用(plasma-material interaction，PMI)相关过程. 现有较成熟的LPD设备有STEP[4]、PSI-2[5]、NAGDIS-Ⅱ[6]、PISCES-A/B[7]和试验- PSI/MAGNUM-PSI[8]，其参数如表1所示.

表1 各类辐照装置的等离子体参数比较[4-8]

ECR直线等离子体装置，是以微波电子回旋共振(microwave electron cyclotron resonance,MWECR)技术为核心，能够在低气压(0.1～1.0 Pa)下获得较高的等离子体密度[9]. 由北京工业大学薄膜材料实验室搭建的ECR离子源装置配合有永磁体单元和线圈，设计简单、运行可靠并且可以调控共振面位于放电腔室中的位置[10]. ECR装置的基本原理是当电子回旋频率与微波输入频率同为2.45 GHz时就会在磁感应强度为875 Gs的位置点产生电子回旋共振作用，电子获得能量，从而使腔体内的工作气体电离率提高. 但是在具体的实验过程中，共振点如果接近微波窗口，产生的热量会导致微波窗口处的石英玻璃破裂，影响装置的运行稳定性，因此需要考虑共振点的位置. 同时磁场具有调控粒子轨迹的作用，离子通量和等离子体密度等参数强烈地依赖于磁场环境，所以对磁场的优化对于实现高通量辐照实验的稳定运行具有非常重要的意义.

鉴于此，本文采用有限元方法进行仿真计算并通过模拟计算的结果与实际磁场环境的测量结果相对比，得到共振点位置随电磁场电流不同的分布情况，以及对共振点处粒子出射到达样品台的束流密度进行了分析.

1 仿真建模

1.1 仿真建模与参数设定

本模型采用磁场分布与带电粒子追踪，求解腔体内磁场空间分布以及带电粒子在瞬态下的运动轨迹. 图1(a)为本实验室ECR装置腔体示意图. 腔体由微波窗口、励磁线圈、样品台以及永磁体组成. 在讨论ECR装置腔体内磁场分布时，首先采用二维轴对称模型等效替代三维模型计算腔体内磁场分布的情况，其理由有：1) 本实验室中ECR装置结构呈轴对称，根据文献[11-12]对静磁场和线圈磁场模拟结果都为近似轴对称分布； 2) 由文献[13]可知，三维束流分布结果同样需要在二维轴对称的截面表示粒子的出射，采用二维轴对称模拟磁场分布可以加快计算速度，并且可以使网格尺寸划分得更小，使结果更加准确，对于计算的结果后处理更方便. 如图1(b)所示，构建了上半部分为7.0 cm×24.0 cm与下部分为12.5 cm×30.0 cm的矩形区域作为圆柱求解域的半个切面. 矩形的求解域设定为真空. 永磁体材料设定为钕铁硼，同时励磁线圈材料设定为金属铜. 考虑到为了模拟与实验室ECR装置一致的条件，设定永磁体为高度2 cm、半径3 cm的圆柱形，剩余磁化强度设置为1.3 T，磁场方向沿z轴向上. 同理，设定励磁线圈截面积为14 cm×12 cm，线圈金属铜管的圆截面直径为8 mm，匝数为220. 实验装置中使用的励磁电源为MP15200D，最大稳定输出电流为200 A，所以在仿真中给定线圈电流为100～200 A. 仿真区域内的总磁场B相当于永磁铁的磁场Bmagnet和电磁铁磁场Bcoil的矢量和[11，14]，对于线圈磁场Bcoil，安培定律决定磁矢势A的角向分量

(1)

式中：μr为相对磁导率；μ0为真空磁导率；J为电流密度.

电流密度Jφ有角分量

Jφ=NI/S

(2)

式中：N为螺线管的匝数；I为螺线管电流；S为螺线管的横截面积.

磁场Bcoil可以表示为

(3)

对于静磁场Bmagnet，有

(4)

式中:μr为相对磁导率；μ0为真空磁导率；Vmagnet为永磁铁的磁势.

图1 ECR装置示意图与仿真模型对照Fig.1 Comparison of ECR device schematic diagram and simulation model

1.2 网格划分

在绘制完仿真模型以及设定好参数后，需要进行网格划分. 在有限元分析中，网格划分的质量与数目直接关系到计算规模的大小和预计算结果的精确度. 考虑到二维模型计算规模较小，计算速度快，为了提高计算精度，将图2空气域网格部分设定为较细化，磁场线圈和永磁体网格部分设定为极细化.

图2 网格划分Fig.2 Diagram of computational mesh

2 仿真结果与分析

2.1 腔体磁场分布模拟与实验对比

在后处理中添加一维绘图组，手动选择以样品台表面为0点的中轴线，调节励磁线圈电流为100～200 A，得到腔室轴线处的磁感应强度，如图3所示. 可以看出在励磁线圈与永磁体构成的叠加磁场中，距离样品台越近，磁感应强度越大，在距离样品台表面14 cm处磁感应强度衰减到最小点为0.008 T，在距离样品台表面高度为42 cm上下两侧磁感应强度呈对称分布.

图3 不同电流下纵向磁感应强度Fig.3 Longitudinal magnetic induction intensity at different currents

实验过程中，使用高斯计沿腔轴线测量出不同高度的磁感应强度. 选取励磁线圈电流120、140、160 A实验数据与模拟数据对比. 图4中模拟与实验结果增减趋势大致相同，模拟数据整体稍大于实验数据，主要原因如下[15].

1) 在模拟中，简化了模型，且只考虑了电磁场的情况，而没考虑其他的加载. 例如：实际情况为了冷却，在铜管中有循环水冷. 会对电磁场产生一定影响.

2) 在实际测量中，高斯计测量时数据跳动，测量存在误差.

3) 励磁线圈铜管采用线包的方式，缠绕精度不够.

图4 不同线圈电流为时磁场强度模拟与实验对比Fig.4 Simulation of magnetic field strength for different coil currents compared with experiments

2.2 共振点位置

在实验过程中，测得共振点位置距离微波窗口非常近，所以在放电过程中微波窗口会不断受到等离子体的轰击[16]. 等离子体的轰击对微波窗口的影响主要是大功率时将使微波窗口的温度急剧升高，会在一瞬间导致微波窗口受热不均而发生爆裂的现象. 需要确保只有一个共振点在腔室内部的基础上尽可能地使共振点远离微波窗口，因此，调节共振点的位置对于保护微波窗口和微波系统至关重要.

如图5所示，在结果的后处理中定义775、875、975 Gs等值线，绿色线条位置即为共振层的位置. 以微波窗口为零点，通过测量可以得到中轴线上下两共振点到样品台的距离，分别用La、Lb表示，微波窗口到样品台的距离为L0，则上共振点到微波窗口的距离a=La-L0，下共振点到微波窗口的距离b=L0-Lb. 表2列出了励磁线圈电流为100～200 A时上下共振点距微波窗口的距离. 在100 A时上下2个共振点都在微波窗口之下，120～200 A时随着电流增大，上下共振点的位置也随之增大. 对调节励磁线圈组的位置以及电流的大小提供了一个可靠的依据.

图5 共振点位置示意Fig.5 Esonance point positions

ECR装置中的励磁线圈磁场与永磁体磁场可以构成磁镜效应，当带电粒子向磁场变强的方向运动时，带电粒子由于磁矩的不变性会使垂直速度变大，同时其总动能由于磁场不对带电粒子做功而维持不变，因此平行速度会相应地减小. 当平行速度为零，会受磁场的反射向相反的方向运动，当运动至与之前磁感应强度大小值相同的地方时，带电粒子又会被此处的磁场反射，周而复始，实现了对带电粒子在磁场中的约束，同理，若平行速度较大的带电粒子运动到磁感应强度最大处时，其平行速度不能减小到零，那么这些带电粒子就能够穿越磁场，成为逃逸粒子[17-18]. 如图6所示，在励磁线圈磁场和永磁体单元磁场分别存在最大磁感应强度位置点Bmax1、Bmax2，并且在两者之间存在最小磁感应强度位置点Bmin，当磁镜比Rc=Bmax/Bmin越大时，磁镜对粒子的约束性能越强. 因此通过模拟研究了本实验室ECR直线等离子体装置磁场对束流密度的影响.

表2 不同线圈磁场电流的共振点在轴线处位置

图6 ECR装置磁镜示意Fig.6 Schematic diagram of ECR device magnetic mirror

2.3 磁场电流对粒子束流密度的影响

对二维模型进行轴对称旋转以分析带电粒子在三维磁场中的分布情况，选择共振层为发射面，为了减小计算量，本文采用了一定的近似，因对束流的贡献主要来自速度垂直向下的离子，因此计算中只重点考虑了垂直向下的离子. 粒子在发射面呈高斯分布，选择每次释放的粒子数N为100个，粒子的初始动能为70 eV. 由图7可以看出，线圈电流为160 A时，通过微波窗口观察，实验中等离子体聚集情况基本与模拟结果一致，均呈梨形分布. 为进一步了解束流密度在三维磁场中的分布情况，选择每次释放的粒子数N为1×106个，其他初始条件不变，在距样品台表面处0、1、2、3、4、5、10、15、20 mm处添加截面，在截面中心径向位置划分51个2 mm×2 mm的正方形截面单元，截面单元编号及位置具体分布如图8(a)所示，观察通过各个截面的束流密度径向分布情况，其束流密度分布如图8(b)所示，在各截面处的束流密度均呈对称分布，且束流密度最大值均在中心位置处，连接各截面中心位置处，发现在距样品台表面20～5 mm处，束流密度增长缓慢，在5 mm至样品台表面处增长迅速，其样品台表面中心位置处束流密度可达2.2×108粒子数/m2，边缘位置的粒子受磁场约束能力较小，成为逃逸粒子.

图8 截面位置及各位置在160 A电流下的束流密度Fig.8 Section position and beam density at current 160 A at each position

在结果后处理中，样品台表面所在空间截面添加庞加莱截面，即对粒子连续运动的轨迹将其横截，得到轨迹穿过截面的情况以及粒子的数量N，其磁场电流为160 A时样品台表面粒子分布情况如图9所示，红色点为沉积到样品台表面处的粒子，聚集形状均呈圆形，外侧旋状红色点则为漂移到样品台外的粒子，通过测量得出通过样品台表面处束斑的面积S，则有样品台截面处单次束流密度为N/S. 由上文得励磁线圈在100 A时，上下共振点都在微波窗口之下，所以只取线圈电流由120 A增大到200 A时对束流密度的影响，如表3所示，束斑面积随励磁线圈电流的增大变化不大，其粒子数随励磁线圈电流的增大而增大，导致其束流密度逐渐增大. 为进一步了解励磁线圈电流对束流密度的影响，取样品台处如图10(a)所示的11个2 mm×2 mm的正方形截面单元，且6号位置为中心位置. 计算出通过这11个位置的粒子数得到瞬时束流密度，其具体分布情况如图10(b)所示，在样品台表面处，束流密度呈正态分布，相同位置处束流密度随着电流增大而增大，中心位置处束流密度随电流变化较为明显. 其原因是粒子在非均匀磁场的背景下受洛伦兹力的影响沿着磁感线方向运动[19]. 随着线圈电流的增大，共振层的位置远离微波窗口，即粒子的发射层靠近样品台，粒子行程变短，束流密度变大. 在永磁体单元与线圈磁场中间的磁场强度最小值大约为114 Gs，因此在单线圈磁场中心处的磁镜比大约为13.1，而在永磁体单元中心处的磁镜比大约为25.7，符合轴向磁镜比趋向样品台方向越来越大的趋势，显著地提升了磁镜对粒子的约束能力[20].

图9 样品台表面所在空间截面粒子分布情况Fig.9 Distribution of particles in the space section where the surface of the sample table

表3 样品台表面处粒子参数

图10 截面单元位置及各位置在不同电流下的束流密度Fig.10 Position of the section element and the beam density at different currents at each position

3 结论

ECR离子源工作稳定，结构相对简单，工作寿命长，并且在低气压下产生高密度、高电离度和体积均匀的等离子体. 虽然由于微波窗口靠近共振点在大功率时容易产生过大热应力从而破碎，且装置对粒子的约束能力较弱，但是可以通过改变其磁场结构来改善该缺陷.

1) 通过有限元计算模拟，腔体内磁场分布结果和实验数据比较符合，可以调整共振点的位置，在确保只有一个共振点在腔室内部的基础上尽可能地使共振点远离微波窗口.

2) 通过模拟带电粒子在三维磁场中的分布，解释了放电区形状的构成原理，实现励磁线圈与永磁体叠加磁场对粒子约束能力大小的可控，并给出了励磁电流对束流密度的影响关系.

3) 带电粒子轨迹模拟结果与实验中近似梨形等离子体束的分布一致，到达样品台的离子数随着励磁线圈电流的增大而增大，而束斑的面积的变化不大，因而其束流密度随励磁电流的增大而增大.