卫星尾迹带电数值模拟研究

秦晓刚，李得天，汤道坦，柳 青，陈益峰

(兰州空间技术物理研究所，真空低温技术与物理重点实验室，甘肃兰州730000)

1 引言

在地球外围空间中，存在着大量的等离子体。以环绕地球的等离子体为主的空间区域称为磁层，它是受太阳风和行星际磁场约束、地磁场起控制作用的一个有限空间范围［1］。磁层等离子体的性质一般用其粒子密度和粒子能量来描述，对其影响最显著的因素是空间高度和纬度的变化［2］。其中有些轨道(如低/极轨道)的等离子体环境特点为高密度(1010～1012m－3)和低能量(0.1～0.3 eV)，同时又存在高能电子(1～100 keV)的注入。这些特点给卫星带电带来了新的内容—卫星的尾区带电效应［3］。

当卫星运行在低温度、高密度等离子体环境中时，在其尾迹形成一明显的“航迹”，这是一个不相等的电子和离子耗尽区。由于卫星轨道速度大于离子热速率而小于电子热速率，因此电子可较容易地进入这个区域从而形成一负电位势垒，这就是所谓的“尾迹效应”。它对卫星的明显作用是在尾区介质表面将充电至较高的负电位，此表面电位主要依赖于收集的电子通量与离子通量之比［4］。卫星因“尾迹效应”而形成的表面高压带电是影响中低轨道特别是极轨卫星安全运行的重要原因之一。由于星尾区表面高压带电引起的放电危害小则可影响卫星的正常工作，大则可使整星失效。实验主要通过数值分析的方法，研究卫星尾迹的带电情况。

2 尾迹带电效应的机理

运行在低轨道中的卫星可以看作是一个浸没在低温、高密度等离子体悬浮高压探针。卫星周围的等离子体的平衡状态可以用poisson方程和无碰撞的Vlasov方程表示［5］:

其中 φ是和无穷远处未绕动等离子体比较的电势，fi和fe分别是离子和电子的分布函数。

当卫星的表面电势处于平衡状态时，到达卫星表面的离子电流和电子电流相等，可以表示为:

式中 je和ji分别是环境电子和离子的电流密度。平衡状态时的电子布局可以用下式表示:

其中 k是玻尔兹曼常数，Φw是卫星的电势。通过上面两式并归一化可以得出:

在给定的电子和离子参数情况下，上式表明离子电流决定了卫星的电势。

而在低轨道中，由于Cs＜＜vsc＜＜Ce(其中vsc是卫星的速度，Ce是电子热运动速度)，离子很难到达卫星的尾部，但是电子可以不受限制，即在式(5)中的离子电流将会减少，而电子电流基本没有变化。从而卫星的尾部，电势将会出现一个更负的电势，这时卫星的尾部将会出现所谓的尾区。

当有极光电子注入时，式(5)可以表示为:

从上式可以发现，在高能电子电流注入时，卫星尾部收集的电子电流将会更大，即其尾部将出现比没有高能电子注入时更负的电位。尾迹效应的示意图如下图1所示。

图1 尾迹效应示意图

3 尾迹带电效应的PIC模型

PIC方法对等离子体粒子的模拟就是通过跟踪大量电子和离子在电磁场中的运动来描述等离子体的动力学行为。通常假设大量的带电粒子具有初始位置和速度，再根据边界条件等，对它们统计平均求出等离子体空间的电荷和电流密度分布，主要是通过计算描述电磁场演化的Maxwell方程组［6～8］来实现。

对于空间低能等离子体与卫星相互作用来说，主要是卫星充电产生的电场和空间等离子体及其自洽电场相互作用的过程。针对LEO轨道空间等离子体与航天器相互作用的PIC模拟技术，由于地磁场的影响被忽略不计了。因此，利用PIC方法来描述等离子体的运动和航天器表面充电过程可以采用静电模型。

在静电模型中，等离子体的运动状态变化主要是由于电荷分离产生静电场所引起的。对于静电模型来说，只需要求解泊松方程即可:

从上面的理论分析可以知道，卫星尾区效应主要是由离子速度和卫星运动速度的可比性造成。因此为了分析尾区的带电特性，需要着重考虑航天器相对与等离子体流的运动运动速度。计算模拟时，可以考虑运动的相对性，即离子可以认为是朝着卫星迎风面进入模拟区域，而电子的运动方向是各项同性的，具体的示意图如图2所示。且模拟计算时，一般二维平面的情况就可以反应出尾区带电的物理特性，本章就二维PIC数值模拟程序进行数值方法分析。

图2 尾区效应计算模型示意图

为了提高模拟精度，程序中采用了leap-frog(蛙跳)格式。计算粒子运动时，在整时间点计算粒子位置，而在半时间点计算粒子速度，以该速度代替相邻两整时间点间的平均速度。网格电场的计算也借鉴leapfrog格式的思想，在网格点上计算电势，而在网格中心处计算电场。

利用PIC方法计算等离子体时，对于空间网格的和时间步长的设置要兼顾计算量以及精度的要求。对于模拟空间大小，网格划分，时间步长取值依据模拟空间的等离子体密度，温度等参数确定。

主程序运行都需要终止的平衡条件，程序运行的平衡条件是由卫星表面电位，空间等离子体和电位分布决定的。当程序运行到这些量基本达到一个比较小的波动变化时，即认为系统状态达到了平衡，这时统计所需要的物理量。

4 结果与讨论

计算时，地球轨道等离子体参数选取中采用DEMETER卫星的一组测量数据:ne=ni=109/cm3，电子温度为0.26 eV，离子温度为0.13 eV，卫星速度vs/c=7.5 km/s。卫星表面材料OSR热控膜，其二次电子发射系数为2.5，对应的能量为0.3 keV。同时我们忽略了地磁场和自洽磁场的作用，卫星的运动方向沿－X方向，卫星的尺寸我们选择为1 m，模拟空间为10×10 m，空间网格为0.05 m。

图3所示的为我们采用自编PIC程序计算的卫星本体及周围电势分布。从上图可以看出，卫星尾部电位大约为－1 V左右，而卫星迎风面的电位可以达到0.1 V左右，尾区的扰动区域为约为1～2 m范围。

卫星尾区带电和卫星的尺寸是密切相关的，国外研究结果表明卫星尾区带电和卫星尺寸呈正比关系，即卫星尺寸越大，尾区效应带电越严重。因此，我们也分析不同卫星尺寸的尾区带电情况。图4所示的就是卫星尺寸为2 m时的卫星周围电势分布。从图中可以看出，2 m卫星的尾部表面电位约为－2 V，大约为1 m尺寸卫星尾部电位的2倍左右，结果符合前面的理论分析。

图3 自编PIC程序计算的1 m卫星周围电势分布

图4 自编PIC程序计算的2 m卫星周围电势分布

从上面的分析结果可以看出，当卫星处于低轨时，其尾区表面带电基本都是负几伏的数量级，对卫星基本不能造成危害。但是当卫星运行到极轨时，由于极光电子(能量约为数KeV)的注入，其尾区带电可能就是主要的带电问题。图5所示的为我们在程序中添加高能电子时，卫星的尾区带电效应。由于网格划分的兼容性问题，添加的高能电子能量为0.5 KeV，密度为107m3。从图中可以看出，卫星迎风面的电位和其尾部的电位相差接近200 V。

图5 高能电子注入时卫星尾区带电情况

5 结论

本文在分析卫星尾区带电物理机制的基础上，建立简单的二维计算分析模型，并利用PIC方法编制二维计算程序。运用程序计算不同卫星尺寸的尾区带电情况和高能电子对尾区的带电影响。通过对比计算结果，我们得出卫星尺寸和卫星尾区带电效应是呈正比关系。模拟计算的结论可以为低轨卫星的带电防护提供指导。

［1］WHIPPLE E C.Potentials of Surfaces in Space.Reports on Progress in Physics.Vol.44，1981.1197 ～1250.

［2］徐家鸾，金尚宪，等.离子体物理学.第一版.北京:原子能出版社，1981.

［3］J.Forest，L.Eliasson，A.Hilgers.a new spacecraft plasma simulation software，PicUp3D/SPIS.7th spacecraft charging technology conference，The Netherlands，2001.

［4］沈超，刘振兴.外辐射带动态演化的粒子模拟研究，空间科学学报，2001，21(3):230～237.

［5］黄文耿，古士芬.大功率高频无线电波对赤道E区双流不稳定性的影响－PIC静电粒子模拟研究，空间科学学报，2003，23(2):102～109.

［6］黄文耿，乐贵明，古士芬，等.电离层中SEE现象的静电粒子模拟，空间科学学报，2002，22(2):113～118.

［7］郭俊，陆全明，王水，等.磁场重联的二维粒子模拟研究，空间科学学报，2003，23(4):248～255.

［8］陆全明，郭俊，窦贤康，等.二维粒子模拟方法及在空间物理中的应用，计算物理，2004，21(2):137～142.

［9］曹晋滨、汪学毅、周国成，等.低轨道磁化等离子体中运动航天器等离子体鞘层特性，地球物理学报，2000，43(4):459～463.

［10］Yu.N.Grigoryev，V.A.Vshivkov，M.P.Fedoruk.Numerical"particle－in－cell"methods:theory and applications，The Netherlands:VSP BV.2002.

［11］R.W.Hockney，J.W.Eastwood.Computer Simulation Using Particles.France:Taylor，1988.

［12］D.P.Grote，A.Friedman，S.M.Lund，I.Haber.3D Particle Simulation of Space－Charge－Dominated Beams in HIF Accelerator Experiments.Particle accelerator conference.Canada，1997.

［13］Chung Chen，Jim Bronwning，Steve Meassick.Current collection in a spacecraft wake;laboratory and computer simulations，5th spacecraft charging technology conference.1989.

［14］A.Soubeyran，L.Levy.Numerical simulation of the wake of non equipotential spacecraft in the ionospere.5th spacecraft charging technology conference.1989.

［15］Eric Engwall，Anders Eriksson，Arne Pedersen，Julien Forest，Goetz Paschmann，et al.Wake effects on positively charged spacecraft in flowing tenuous plasma:cluster observations and modeling.8th spacecraft charging technology conference.2004.

［16］I.Katz，M.J.Mandell，G.A.Jongeward，et al.Three－Dimensional computer models of current collected by active spacecraft in low earth orbit.5th spacecraft charging technology conference.1989.

［17］David L.Cooke，M.S.Gussenhoven，David A.Hardy.Polar code simulation of DMSP satellite auroral charging.5th spacecraft charging technology conference.1989.

［18］G.Vannaroni，M.Dobrowolny，F.De Venuto，L.Iess.Current collection by rapidly moving charged bodies in the ionosphere.6th spacecraft charging technology conference.2000.

［19］L.W.Parker.Contribution to Satellite Sheath and Wake Modeling.In:Proc.17th Symp.On Spacecraft/Plasma Interaction ESASP－198，Noorwijk，Netherlands，1983.

［20］R.Biasca，J.Wang.Transient charging of objects in spacecraft wakes:Numerical simulations.AIAA 94－0331，1994.