衛星尾跡帶電數值模擬研究

秦曉剛，李得天，湯道坦，柳 青，陳益峰

(蘭州空間技術物理研究所，真空低溫技術與物理重點實驗室，甘肅蘭州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.