基于三维PIC数值模拟的2 cm ECRIT推力控制计算
胡展1, 杨涓1, 陈茂林1, 于达仁2, 朱悉铭2     
1. 西北工业大学 航天学院, 陕西 西安 710072;
2. 哈尔滨工业大学 能源科学与工程学院, 黑龙江 哈尔滨 150001
摘要: 无拖曳控制是空间引力波探测的关键技术,主要由微型推力器完成。微型电子回旋共振离子推力器(ECRIT)体积小、推力可调,可用于空间引力波的无拖曳控制。基于三维PIC数值模拟方法计算微型2 cm ECRIT的推力控制范围,分析其用于无拖曳控制系统的可行性。首先计算不同栅极孔径下的推力性能和栅极聚焦特性,获得合理栅极结构,再计算栅极电压、栅极前离子密度对推力器性能的影响,获得满足无拖曳控制要求的推力器性能参数范围。结果表明:减小栅极孔径能降低推力,但同时影响栅极聚焦效果;调节栅极前离子密度可大范围调节推力;在给定的栅极结构和栅前离子密度下,存在合适的栅极加速电压区间保证离子的良好聚焦。综合考虑推力性能和栅极聚焦特性,选择屏栅孔径0.6 mm、加速栅孔径0.34 mm的栅极,当栅极前离子密度分别为1×1017,0.7×1017,0.4×1017,0.2×1017 m-3时,通过调节加速电压,可实现5.05~141.44 μN的推力调节。此研究将为分析ECRIT应用于引力波探测的可行性奠定基础。
关键词: 无拖曳控制    电子回旋共振    束流引出    PIC数值模拟    推力计算    

空间引力波探测任务要求采用无拖曳(drag-free)控制技术抵消太空环境对卫星的作用力,从而使卫星平台达到空间实验技术指标要求[1]。日本引力波探测计划DECIGO[2]对相关微推进系统提出要求:推力范围5~100 μN,推力精度小于0.1 μN,推力噪声小于0.1 μN/Hz1/ 2。微型2 cm ECRIT具有推力调节方便、结构简单等优势可应用于引力波探测,受国内外关注。日本Koizumi、Funaki等[3-4]等通过实验手段分析了微型ECRIT的工作参数和栅极孔径变化对推力的影响;Izumi等[5]通过调节微波功率和工质气体流量等参数降低了微型ECRIT的推力范围,同时测量了推力噪声;国内通过调整磁路结构、天线位置和推进剂工质,实现了2 cm ECRIT的推力优化[6-7],同时发现了改变ECRIT推力的有效途径。

本文在相关研究基础上,以2 cm ECRIT的双栅极系统为研究对象,采用三维PIC方法,对栅极的离子引出过程进行三维数值模拟,以获得推力性能参数。计算首先从栅极结构参数入手,根据不同栅极孔径的束流引出计算结果,比较推力性能和栅极聚焦特性,获得合适的栅极孔径。然后固定该孔径不变,计算分析不同工况下的推力调节范围。由于程序无法模拟等离子体的生成过程,而微波输入功率和工质气体流量影响推力的本质是改变栅极引出孔前的离子密度,因此用栅前离子密度参数反映二者的影响。通过调节栅前离子密度和栅极加速电压,获得相应的推力变化范围,从而为分析微型ECRIT用于无拖曳控制系统的可行性奠定基础。

1 ECR离子源结构和栅极计算模型

2 cm ECRIT主要由离子源和中和器两部分组成,其中离子源通过双栅极引出离子束形成推力,中和器引出电子中和离子束保持系统电位平衡。离子源是决定推力器性能的关键部件,其结构如图 1a)所示,主要由内外磁环、磁轭、环形天线、屏栅和加速栅组成,其中栅极结构如图 1b)所示。

图 1 2 cm ECR离子源结构

采用PIC方法进行栅极系统数值模拟时,主要考虑耦合电场分布、电子和离子运动以及电荷交换碰撞过程。其中电场按静电场处理,通过离散泊松方程求解得到电场强度。离子采用粒子方法处理,其运动满足Newton-Lorentz方程,而电子则当做流体处理,密度分布根据玻尔兹曼关系式求解。碰撞方面考虑弹性碰撞和CEX碰撞,采用文献[8]的碰撞截面参数,并用Monte Carlo方法处理。

由于栅极的对称性, 选择如图 2所示的包含2个1/4孔的长方体为计算域。网格步长服从德拜长度λD的限制, 即ΔX < λD, ΔY < λD, ΔZ < λD; 时间步长满足Δt·ωp≤1, 其中ωp是等离子体频率。为了维持PIC粒子模型计算的稳定性及减少粒子运动噪声, 要满足CFL条件[9]

图 2 计算区域

边界条件如图 2所示, ABCD为入口边界面, 等离子体电势为:

(1)

式中:Φsc为屏栅电势;Φw为壁面相对于预鞘的电势降;Φp为等离子体与预鞘的电势差。A′B′C′D′为出口边界面, 当离子离开出口边界时, 认为离子逸出计算域, 删除粒子信息。AA′B′B, AA′D′D, BB′C′C, CC′D′D为对称边界面。当离子运动越过对称面时, 认为离子被镜面反射。对于栅极固壁, 电势为恒定值; 栅极内部电场强度为0 V/m。当离子撞击到栅极表面时, 认为离子被壁面吸收, 删除粒子信息。

在计算中, 记录每个时间步上的离子信息, 统计打到屏栅和加速栅上的离子数目,计算屏栅电流Isc和加速栅电流Iac, 统计轴向电流密度参数计算引出束流Iion, 离子通过率为; 根据离子轴向速度计算推力。获得计算区域的性能参数后, 进而扩展到整个栅极区域得到推力器各性能参数。由于模型中使用2个1/4栅孔计算区域计算整个栅极性能, 而在离子源放电室出口的不同径向位置, 等离子体密度变化较大, 因此在计算中合理的栅前离子密度参数选取很重要。基于离子推力器原理[10], 在引出束流未饱和, 即未达到Child-Langmuir鞘层最大引出束流密度限制前, 引出束流Ib满足

(2)

式中,e, n0, Te, M, As, Ts分别为元电荷、栅前离子密度、电子温度、气体原子质量、栅极面积、栅极离子透过率。因此可根据离子源束流引出实验结果, 依据(2)式反推估算相应的栅前离子密度参数。在微波输入功率0.5~2 W、氙工质流量9.7~29.3 μg/s的工作区间内, 栅前离子密度的变化范围约为0.2×1017~2×1017 m-3

图 3所示, 在离子聚焦加速过程中, 存在不同的聚焦状态。在正常聚焦条件下, 等离子体与屏栅形成的离子鞘面略微弯曲, 发射的离子将全部加速通过栅孔被引出。如果栅极加速电压不足或栅前等离子体密度过高, 引出束流密度偏大, 鞘面变平, 离子聚焦性变差, 为欠聚焦状态。随着引出束流密度的进一步增大, 束流直径大于加速栅孔径, 离子束边缘的离子直接轰击加速栅孔边缘, 导致加速栅截获电流急剧增大, 对应的极限束流称为导流极限。相反的, 如果栅前等离子体密度较低或栅极加速电压过高, 引出束流密度偏低, 鞘面过于弯曲, 会造成离子束的过聚焦。随着引出束流密度的进一步减小, 通过屏栅孔边缘的离子交叉穿过束流中轴线轰击另一侧的加速栅孔边缘, 导致加速栅截获电流急剧增大, 对应的极限束流称为交叉极限[11]。当栅极未能正常聚焦时, 部分高能离子直接撞击到加速栅上, 不仅会使引出的有效束流减小, 还会加快加速栅的腐蚀失效。为避免该现象的发生, 在栅极数值模拟时, 需对栅极的聚焦性能进行分析。定义加速栅截获电流比γ为聚焦性的表征:

(3)

式中:Iac为加速栅截获电流;Ib为引出的离子束电流。加速栅截获电流比参数可直观表征离子轰击带来的栅极腐蚀现象, γ越小表示栅极离子引出聚焦性越好。

图 3 栅极离子引出的不同聚焦状态

为了验证栅极仿真模型的可靠性, 选择现有2 cm ECRIT结构进行模拟, 其屏栅孔直径0.72 mm, 加速栅栅孔直径0.4 mm。在输入微波功率2 W、氙气工质流量29.3 μg/s、栅极加速电压1 850 V工况下, 栅极入口平均等离子体密度取2×1017 m-3, 上游电子温度取10 eV, 下游电子温度取1 eV。在该工况下, 计算引出束流为5.268 mA, 推力为367.7 μN;

实验中测得引出束流为5.3 mA, 采用推力公式(4)根据束流值估算推力器推力为368.6 μN。

(4)

可以看出, 采用该仿真模型计算出的推力性能与实验结果基本一致, 验证了模型和数值方法的准确性。

2 计算结果分析 2.1 不同栅极孔径计算

其他条件不变, 取5种屏栅/加速栅孔径(mm)0.72/0.4, 0.66/0.36, 0.6/0.34, 0.54/0.32, 0.48/0.3进行推力性能计算。参考前期实验参数[6], 取栅前离子密度2×1017m-3为高密度工况, 0.2×1017m-3为低密度工况。加速电压分别取1 850和850 V。

高密度工况下不同栅极孔径离子源计算结果如图 4所示, A, B, C, C, D, E分别表示5种栅极孔径组合(mm)0.48/0.3, 0.54/0.32, 0.6/0.34, 0.66/0.36, 0.72/0.4。可以看出, 在相同的栅极加速电压下, 随着栅极孔径的减小, 栅极透明度降低, 离子通过率不断减小, 引出束流和推力随之降低。不同栅极加速电压下, 加速栅截获电流比γ变化趋势略有不同。在1 850 V加速电压下, 随着栅极孔径的减小, γ变化不大, 栅极离子引出均处于良好聚焦状态; 而在850 V的栅极加速电压下, γ随栅孔减小呈现增长趋势, 随着栅极孔径的减小, 栅极引出达到欠聚焦状态。对于屏栅孔径0.48 mm、加速栅孔径0.3 mm的栅极结构, 在850 V加速电压下, 已接近束流导流极限。这是由于栅前离子密度较高, 随着栅孔减小, 离子有效引出面积减小, 虽然引出束流也有所降低, 但束流引出密度有所增长, 更容易达到欠聚焦状态。

图 4 高密度工况下不同栅极孔径离子源计算结果

这说明过小的栅极孔径不适合高密度工况的束流引出。

低密度工况下不同栅极孔径离子源计算结果如图 5所示, 可以看出, 与高密度工况相似, 随着栅极孔径的减小, 离子通过率和推力也随之降低。而在不同的栅极加速电压下, 栅极的聚焦效果有很大差异。从图 5b)可以看出, 在850 V的栅极加速电压下, 不同栅极孔径的离子引出过程加速栅截获电流比都很低, 均处于正常聚焦状态; 而在1 850 V的栅极加速电压下, 加速栅截获电流比都很高, 这是由于低密度工况下的栅前离子密度很低, 引出束流密度相对较低, 在较高的栅极加速电压下, 栅极离子引出更容易达到束流交叉极限, 加剧加速栅的截获。这也说明, 在较低的离子密度时, 不宜采用过高的栅极加速电压。

图 5 低密度工况下不同栅极孔径离子源计算结果

根据以上计算结果, 可以看出, 虽然减小栅极孔径能有效降低推力范围, 但过小的栅极孔径不适合高密度工况的束流引出。综合考虑孔径对推力性能和聚焦特性的影响, 选择屏栅孔径0.6 mm, 加速栅孔径0.34 mm的栅极结构进行后续仿真计算。考虑到改变离子密度后, 放电室需要一定时间达到电离平衡, 而改变栅极加速电压响应较快, 推力变化及时[12], 因此根据控制变量法, 先固定栅极加速电压不变, 分析栅前离子密度对推力性能的影响, 获得满足5~100 μN推力要求的离子密度范围, 然后从中选择几个密度保持不变, 通过调节栅极电压的方式来实现推力的快速调节。

2.2 不同栅前离子密度计算

取栅极加速电压分别为1 850和850 V, 对屏栅孔径0.6 mm, 加速栅孔径0.34 mm的栅极结构进行不同栅前离子密度条件下的仿真计算, 计算结果如图 6所示。

图 6 不同栅前离子密度下计算性能参数

从图中可以看出, 随着栅前离子密度的不断减小, 推力下降趋势明显, 而离子通过率略有增加。这是由于随着栅前离子密度的减小, 可供引出的离子总数在减小, 束电流强度降低, 鞘层往屏栅孔外发展, 覆盖区域增大, 从而使离子通过率略有增加。从图 6b)可以看出, 在计算的大部分栅前离子密度下, 栅极系统均能保持良好的聚焦效果, 仅在1 850 V栅极加速电压下, 曲线左端加速栅截获电流比出现明显增长。这是由于随着栅前离子密度的减小, 引出束流密度减小, 栅极离子引出更容易达到过聚焦状态。在较低的栅前离子密度下, 过高的栅极加速电压将使离子鞘面更加弯曲, 加剧栅极过聚焦, 离子引出更接近束流交叉极限, 使加速栅截获电流比急剧增长。这也说明随着栅前离子密度的变化, 保持良好聚焦状态的栅极加速电压范围也在改变。

2.3 不同栅极加速电压计算

根据计算结果, 可以看出, 当栅前离子密度为0.2×1017~1×1017 m-3时, 基本可以满足5~100 μN的推力范围, 且当栅前离子密度分别为1×1017, 0.7×1017, 0.4×1017, 0.2×1017 m-3时, 可实现推力的连续调节。因此, 选择上述4种栅前离子密度, 分别进行不同栅极加速电压下仿真计算。在计算中, 栅极加速电压分别取600, 850, 1 200, 1 450和1 850 V。计算结果如图 7所示, 可以看出, 随着栅极加速电压的增加, 推力上升趋势明显, 这是由于增加栅极加速电压不仅能增加离子通过率, 还能提高引出的离子速度, 从而增大推力。比较图 7中加速栅电流比γ计算结果可以看出, 当栅前离子密度较高时, 过低的栅极加速电压会使离子引出达到束流导流极限, 导致γ增加如图 7a)所示; 当栅前离子密度较低时, 过高的栅极加速电压会使离子引出达到束流交叉极限, 导致γ急剧增加如图 7d)所示。在这2种情况下, 栅极聚焦性能均不理想。

图 7 不同栅极加速电压下的推力性能

由上述计算可知, 对于固定的栅极结构和栅前离子密度, 应存在一个合适的栅极加速电压范围以保证离子引出过程的良好聚焦。因此, 根据图 7计算结果, 重新调整栅极加速电压, 获得4个栅前离子密度良好聚焦状态的栅极加速电压范围和对应的推力范围如表 1所示。可以看出, 计算选取的4类栅前离子密度在保证栅极聚焦状态的同时, 基本能够实现推力从5.05~141.44 μN的连续变化。

表 1 不同离子密度下良好聚焦状态的推力范围
栅前离子密度/m-3 栅极加速电压范围/V 推力范围/μN
1×1017 700~3 000 62.34~141.44
0.7×1017 600~2 500 33.92~79.03
0.4×1017 500~1 900 13.61~39.07
0.2×1017 300~1 300 5.05~13.67
3 结论

通过对2 cm ECRIT的栅极数值仿真, 探究了推力器的性能变化规律。获得的主要结论如下:

1) 减小栅极孔径能降低推力, 但同时也会影响栅极系统聚焦效果。过小的栅极孔径不适合高密度工况束流引出。综合考虑, 选取屏栅孔径0.6 mm, 加速栅孔径0.34 mm的栅极结构。

2) 随着栅前离子密度的减小, 尽管栅极离子通过率略有增加, 推力仍会明显降低。

3) 在给定的栅极结构和栅前离子密度下, 存在一个合适的栅极加速电压区间保证引出离子的良好聚焦。在这个区间, 调节栅极加速电压可以实现推力连续调节。

4) 对于屏栅孔径0.6mm、加速栅孔径0.34 mm栅极结构的离子源, 当栅前离子密度分别为1×1017, 0.7×1017, 0.4×1017, 0.2×1017 m-3时, 通过调节栅极加速电压, 可实现5.05~141.44 μN的推力调节。

参考文献
[1] 胡明, 李洪银, 周泽兵. 无拖曳控制技术及其应用[J]. 载人航天, 2013, 19(2): 61-69.
HU M, LI H Y, ZHOU Z B. Drag-Free Control Technology and Its Applications[J]. Manned Spaceflight, 2013, 19(2): 61-69. (in Chinese)
[2] FUNAKI I, NAKAYAMA Y, HORISAWA H, et al. Micro-Thruster Options for the Japanese Space Gravitational Wave Observatory Missions[C]//International Electric Propulsion Conference, 2011
[3] KOIZUMI H, KUNINAKA H. Performance of the Miniature and Low Power Microwave Discharge Ion Engine μ1[C]//46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, 2010
[4] NAKAYAMA Y, FUNAKI I, KUNINAKA H. Sub-Milli-Newton Class Miniature Microwave Ion Thruster[J]. Journal of Propulsion and Power, 2007, 23(2): 495-499. DOI:10.2514/1.21565
[5] IZUMI T, KOIZUMI H, YAMAGIWA Y, et al. Performance of Miniature Microwave Discharge Ion Thruster for Drag-Free Control[C]//48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. Atlanta, Georgia, 2012
[6] 夏旭, 杨涓, 金逸舟, 等. 磁路和天线位置对2 cm电子回旋共振离子推力器性能影响的实验研究[J]. 物理学报, 2019, 68(23): 235202.
XIA Xu, YANG Juan, JIN Yizhou, et al. Experimental Study of Magnetic Circuit and Antenna Position Influence on the Ion Beam Extraction and Coupling Voltage of 2 cm ECRIT[J]. Acta Physics Sinica, 2019, 68(23): 235202. (in Chinese)
[7] 汤明杰, 杨涓, 金逸舟, 等. 微型电子回旋共振离子推力器离子源结构优化实验研究[J]. 物理学报, 2015, 64(21): 215202.
TANG Mingjie, YANG Juan, JIN Yizhou, et al. Experimental Optimization in Ion Source Configuration of a Miniature Electron Cyclotron Resonance Ion Thruster[J]. Acta Physica Sinica, 2015, 64(21): 215202. (in Chinese)
[8] SUN A B, MAO G W, YANG J, et al. Particle Simulation of Three-Grid ECR Ion Thruster Optics and Erosion Prediction[J]. Plasma Science and Technology, 2010, 12(2): 240-247. DOI:10.1088/1009-0630/12/2/21
[9] FARNELL C C, WILLIAMS J D, WILBUR P J. NEXT Ion Optics Simulation via FFX[C]//39th AIAA, ASME, SAE, and ASEE Joint Propulsion Conference and Exhibit, 2003
[10] GOEBEL D M, KATZ I. Fundamentals of Electric Propulsion:Ion and Hall Thrusters[M]. New York: John Wiley & Sons Inc, 2008: 100-141.
[11] 刘金声. 离子束技术及应用[M]. 北京: 国防工业出版社, 1995.
LIU Jinsheng. Ion Beam Technology and Application[M]. Beijing: National Defense Industry Press, 1995. (in Chinese)
[12] 温正, 王敏, 仲小清. 多任务模式电推进技术[J]. 航天器工程, 2014, 23(1): 118-123.
WEN Zheng, WANG Min, ZHONG Xiaoqing. Multitask Mode Electric Propulsion Technologies[J]. Spacecraft Engineering, 2014, 23(1): 118-123. (in Chinese)
Calculation on 2 cm ECRIT Thrust Regulation Based on 3D PIC Numerical Simulation Method
HU Zhan1, YANG Juan1, CHEN Maolin1, YU Daren2, ZHU Ximing2     
1. School of Astronautics, Northwestern Polytechnical University, Xi'an 710072, China;
2. School of Power Science and Engineering, Harbin Institute of Technology, Harbin 150001, China
Abstract: Drag-free control is the key technology for space gravitational wave detection, which will be completed by miniature thruster. 2 cm Electron Cyclotron Resonance Ion Thruster (ECRIT) may be used for the drag-free control. Therefore, the 3D PIC method is used to estimate the thrust regulation of 2 cm ECRIT through simulating the property of ion extracting from the double-grid system. Through analyzing the thrust performance and grid focusing property under different grid apertures, the better grid structure is selected for the subsequent calculation. On this basis, the influence of grid voltage and the ion density on the performance of the thruster is analyzed. The result shows that reducing grid aperture can reduce the thrust, but it also affects the focusing property of the grid. Changing the ion density can adjust the thrust in a wide range. At given grid structure and ion density, there is a suitable grid accelerating voltage range to ensure the best ion focus characteristics. Considering the thrust performance and focusing characteristics, the grid structure with screen and accelerating grid aperture in 0.6 mm and 0.34 mm is selected. When the ion density is 1×1 017, 0.7×1 017, 0.4×1 017 and 0.2×1 017 m-3, the thrust can range from 5.05 to 141.44 μN through adjusting the grid accelerating voltage. This study can provide information for analyzing the possibility of the thruster application.
Keywords: drag-free control    electron cyclotron resonance    beam extraction    PIC numerical simulation    thrust calculation    
西北工业大学主办。
0

文章信息

胡展, 杨涓, 陈茂林, 于达仁, 朱悉铭
HU Zhan, YANG Juan, CHEN Maolin, YU Daren, ZHU Ximing
基于三维PIC数值模拟的2 cm ECRIT推力控制计算
Calculation on 2 cm ECRIT Thrust Regulation Based on 3D PIC Numerical Simulation Method
西北工业大学学报, 2020, 38(4): 733-739.
Journal of Northwestern Polytechnical University, 2020, 38(4): 733-739.

文章历史

收稿日期: 2019-10-16

相关文章

工作空间