林驍雄,溫 正,陶家生

（中國空間技術(shù)研究院 通信衛(wèi)星事業(yè)部，北京 100094）

基于非結(jié)構(gòu)網(wǎng)格的離子推力器羽流數(shù)值模擬方法研究

林驍雄,溫 正,陶家生

（中國空間技術(shù)研究院 通信衛(wèi)星事業(yè)部，北京 100094）

為適應(yīng)復(fù)雜的航天器外形，文章采用基于非結(jié)構(gòu)網(wǎng)格的數(shù)值模擬方法研究離子推力器真空羽流的基本特征。首先，介紹了非結(jié)構(gòu)網(wǎng)格劃分的方法，并分析了建模區(qū)域生成各種拓?fù)潢P(guān)系的特點；然后，提出了粒子穿越網(wǎng)格的具體算法，包括有利于減少計算量的結(jié)合矩陣重標(biāo)號的帶寬壓縮算法；最后，開展了模擬結(jié)果和試驗結(jié)果的對比研究，驗證了以上數(shù)值模擬方法的可行性。

離子電推力器；羽流；數(shù)值模擬；非結(jié)構(gòu)網(wǎng)格

0 引言

采用高比沖的電推進(jìn)系統(tǒng)（如霍爾推力器、離子推力器）來實現(xiàn)星箭分離后的衛(wèi)星變軌、入軌后的位置保持等任務(wù)，具有良好應(yīng)用前景[1-2]。但是，電推進(jìn)系統(tǒng)的應(yīng)用中存在羽流污染問題。羽流中由于離子和中性原子碰撞所引起的低能量電荷交換（CEX），使離子速度降低，較易受到羽流自洽電場的影響，進(jìn)而改變運(yùn)動方向，撞擊到太陽電池陣以及航天器表面等，產(chǎn)生污染。國內(nèi)外學(xué)者均對如何減小羽流對航天器的污染影響進(jìn)行過研究[3-8]。目前，數(shù)值模擬研究多采用二維軸對稱網(wǎng)格或三維結(jié)構(gòu)網(wǎng)格，較少采用非結(jié)構(gòu)網(wǎng)格[9]。

結(jié)構(gòu)網(wǎng)格中每一內(nèi)部節(jié)點周圍的網(wǎng)格單元數(shù)是完全相同的，并能通過任一節(jié)點立即找到與其鄰近的節(jié)點。因此，結(jié)構(gòu)網(wǎng)格用在推力器局部羽流特征分析時有較大優(yōu)勢。然而，多個電推力器的應(yīng)用，加上復(fù)雜的航天器外形，使得航天器周圍環(huán)境越來越復(fù)雜。對于多推力器陣列的羽流分析[10]、推力器羽流形成的航天器周圍等離子體環(huán)境分析[11]等就不太適宜用結(jié)構(gòu)網(wǎng)格的模擬。而非結(jié)構(gòu)網(wǎng)格中內(nèi)部節(jié)點周圍的單元數(shù)是不固定的，更容易離散外形復(fù)雜的區(qū)域，對任意復(fù)雜外形的高度適體性是其他計算網(wǎng)格所無法比擬的。因此，本文采用基于非結(jié)構(gòu)網(wǎng)格的方法模擬離子推力器羽流，并與試驗結(jié)果對比，驗證這一方法的可行性。

1 物理及數(shù)值模型

1.1 羽流物理模型

羽流等離子體中存在著大量的帶電粒子，每個粒子同時受到其他粒子的作用，如洛倫茲力和碰撞等。本文主要對推力器噴射的羽流粒子在電場、磁場和碰撞作用下的運(yùn)動狀況進(jìn)行數(shù)值模擬。對于電場的作用，需要計算等離子體運(yùn)動產(chǎn)生的電勢，并與外加電場產(chǎn)生的靜態(tài)電勢進(jìn)行疊加以得到電勢分布。

電場作用方程為

式中：ε0為真空介電常數(shù)；ρ為空間電荷密度；E為電場強(qiáng)度。

引入電位Φ，由于E=-?Φ，得到泊松（Poisson）方程為

通過求解此方程可以得到空間電勢，進(jìn)而得到電場強(qiáng)度。

電磁力作用方程為

式中：F為電磁力；q為粒子電荷；v為速度；B為磁感應(yīng)強(qiáng)度。

牛頓運(yùn)動方程為

式中：P為動量，有P=mv，m為粒子質(zhì)量。

離子電推力器產(chǎn)生的磁場在其外部很微弱，因此可以忽略磁場的影響[14]。當(dāng)外部磁場較弱時，帶電粒子的運(yùn)動方程簡化為

本文沒有對電推力器羽流的碰撞模型、入口分布等問題進(jìn)行詳細(xì)討論，具體分析可見文獻(xiàn)[3-7]。

1.2 計算參數(shù)

本文以某衛(wèi)星的LIPS-200氙離子推力器為研究對象。該推力器為雙柵極系統(tǒng)，柵極表面向外凸出。推力器的主要工作參數(shù)如表1所示。模擬區(qū)域的邊界采用固定電勢邊界。

表1 離子推力器工作參數(shù)Table 1 Ion thruster parameters

表1（續(xù)）

2 模擬計算

2.1 計算區(qū)域和計算網(wǎng)格

圖1為衛(wèi)星上離子推力器的安裝位置示意圖，有2個推力器固定在矢量推力機(jī)構(gòu)上。

圖1 衛(wèi)星上推力器安裝位置Fig.1 Installation location of thrusters on satellite

計算區(qū)域為包含推力器出口和中和器出口的一個立方體區(qū)域，其網(wǎng)格采用Gmsh軟件通過腳本語言生成，如圖2所示。推力器加速柵極曲面的弧度為15°。建立網(wǎng)格的目的：一方面是便于粒子碰撞時的采樣，另一方面便于應(yīng)用 PIC（Particle in Cell）方法求解空間電勢。結(jié)構(gòu)網(wǎng)格可以方便建立泊松方程的差分格式，但是對于復(fù)雜邊界，比如極小的中和器出口和較細(xì)的天線結(jié)構(gòu)，會遇到網(wǎng)格生成程序無法收斂的情況；另外在結(jié)構(gòu)網(wǎng)格中確定粒子所在網(wǎng)格單元需要花費(fèi)較大計算量。非結(jié)構(gòu)網(wǎng)格能自適應(yīng)生成符合要求的貼體網(wǎng)格，而且根據(jù)網(wǎng)格面片和體積單元的鏈接關(guān)系可以方便確定粒子所在的網(wǎng)格單元。因此，本文采用了非結(jié)構(gòu)網(wǎng)格，并在電勢計算中將體積網(wǎng)格統(tǒng)一選作四面體網(wǎng)格。選取四面體網(wǎng)格的原因是一方面它可以模擬各種幾何形狀，另一方面在四面體網(wǎng)格中可以非常方便地實施粒子穿越網(wǎng)格算法。這個算法是保證在非結(jié)構(gòu)網(wǎng)格中粒子運(yùn)動定位的關(guān)鍵步驟，如果選用其他類型的非結(jié)構(gòu)網(wǎng)格可能導(dǎo)致這個關(guān)鍵步驟無法實施，具體分析見2.3節(jié)。

圖2 模擬區(qū)域和推力器及中和器局部網(wǎng)格Fig.2 Simulation domain and local grids around thruster and neutralizer

2.2 網(wǎng)格預(yù)處理

網(wǎng)格預(yù)處理主要用來生成非結(jié)構(gòu)網(wǎng)格間各種層次關(guān)系，而層次關(guān)系是指求解區(qū)域（Group）、四面體單元（Cell）、三角形表面（Face）、棱邊（Edge）、節(jié)點（Node）等拓?fù)湓氐逆溄雨P(guān)系。在羽流的數(shù)值模擬中主要包含以下幾種：

1）在縮減矩陣帶寬中需要用到計算區(qū)域內(nèi)邊到節(jié)點的映射Edge2Node。

2）在求解電勢的有限元計算中需要用到四面體單元到節(jié)點的映射Cell2Node、推力器表面三角形單元到節(jié)點的映射 ThrFace2Node、航天器表面三角形單元到節(jié)點的映射ScFace2Node，以及邊界表面三角形單元到節(jié)點的映射 BdFace2Node等。其中，ThrFace指推力器表面三角形單元、ScFace指航天器表面三角形單元、BdFace指模擬區(qū)域邊界表面三角形單元。這些映射已經(jīng)由Gmsh生成網(wǎng)格時提供。

3）在粒子穿越網(wǎng)格算法中需要用到四面體單元到三角形單元的映射Cell2Face、三角形單元到四面體單元的映射Face2Cell、表面三角形單元局部標(biāo)號到全局三角形單元標(biāo)號的映射 ThrFace2Face、ScFace2Face、BdFace2Face。這些拓?fù)潢P(guān)系數(shù)組在結(jié)果后處理過程中也會用到，提前生成這些拓?fù)潢P(guān)系數(shù)組就可以在粒子運(yùn)動過程中直接查表，極大提高了計算速度。拓?fù)潢P(guān)系如圖3所示。

圖3 網(wǎng)格間的拓?fù)潢P(guān)系Fig.3 Topological relationship among grids

2.3 粒子穿越網(wǎng)格算法

粒子穿越網(wǎng)格是非結(jié)構(gòu)網(wǎng)格中模擬粒子運(yùn)動的關(guān)鍵步驟之一。在結(jié)構(gòu)網(wǎng)格中，可以借助網(wǎng)格下標(biāo)直接索引粒子將要穿越的相鄰網(wǎng)格，而在非結(jié)構(gòu)網(wǎng)格中只能憑借網(wǎng)格層級關(guān)系進(jìn)行搜索。

粒子穿越網(wǎng)格表面的過程和表面的特征關(guān)系密切。如果粒子穿越的是計算區(qū)域內(nèi)用于計算的輔助網(wǎng)格表面（如圖4中的粒子3），將不對粒子的運(yùn)動產(chǎn)生任何影響；如果粒子穿越的是刪除表面（如圖4中的粒子4），則粒子運(yùn)動到此表面后，在維護(hù)粒子的列表中將被刪除；如果粒子穿越的是反彈表面（如圖4中的粒子5），則粒子運(yùn)動到此表面之后，將和壁面發(fā)生碰撞，反彈之后仍位于計算區(qū)域內(nèi)。本文還提出了一種粒子更新算法，如圖5所示。該更新算法主要根據(jù)粒子當(dāng)前位置、粒子在模擬步長內(nèi)移動的距離和當(dāng)前單元邊界的關(guān)系進(jìn)行判斷，粒子在模擬步長內(nèi)移動的距離通過公式(5)得到的速度對時間積分而得到。

圖4 粒子穿越網(wǎng)格模型Fig.4 Model of particle traversing through the meshs

圖5 粒子運(yùn)動更新算法Fig.5 Update algorithm of particle motion

2.4 Poisson方程的求解

采用有限元方法求解 Poisson方程（式(2)），形成的代數(shù)方程為Ax=b，其中A為系數(shù)矩陣，x為待求解各節(jié)點的電勢，b為結(jié)合邊界條件形成的矢量。最終形成的系數(shù)矩陣A為對稱和稀疏的矩陣，其中稀疏矩陣是指矩陣中絕大部分元素皆為0。本文采用帶狀壓縮存儲的方法存儲矩陣A中的上三角部分。為了使矩陣重排序后的帶寬及輪廓較小，借助于RCM算法[12]尋找這樣的一個排序數(shù)組。RCM示例網(wǎng)格如圖6中左圖，圖中4、5為內(nèi)部節(jié)點，1、2、3、6、7、8為邊界節(jié)點，將各個節(jié)點和邊形成的無向圖定義為圖G，算法過程簡單概括為：

1）找出圖G的邊界頂點r，并賦予x1（在此示例中r取為節(jié)點1）；

2）以頂點r為根，生成圖G的層次結(jié)構(gòu){x1,x2, …,xn}（如圖中虛線區(qū)分的3個層，括號中的數(shù)值為節(jié)點在該層中的度）；

3）對i=1, 2, …,n，找出圖G的層次結(jié)構(gòu){x1,x2, …,xn}中頂點xi的所有未編號相鄰頂點，并按其度的增加順序編號；

4）將上述頂點編號按逆序排列，得到RCM排序數(shù)組（如圖6中右圖所示）。

圖6 RCM算法示例網(wǎng)格Fig.6 Grids of RCM algorithm

使用RCM前、后的非零矩陣元素及半帶寬如圖7所示。

圖7 使用RCM算法前、后的非零矩陣元素及半帶寬Fig.7 Non zero elements and semi-bandwidth before and after using RCM algorithm

可以看到，原來矩陣的半帶寬為 5，變換后的半帶寬為 4。由于本計算為示例，半帶寬并沒有明顯減少。但是在實際的應(yīng)用中，因為內(nèi)部節(jié)點數(shù)遠(yuǎn)遠(yuǎn)大于邊界節(jié)點數(shù)，半帶寬可以減少到原來的 1/5甚至更多。采取這種方法形成矩陣可以極大減少方程的運(yùn)算次數(shù)，而且在用稀疏帶狀矩陣進(jìn)行Cholesky分解時引入的非零元素都在帶寬及輪廓范圍內(nèi)，不會過多占用額外內(nèi)存，可以提高運(yùn)算效率。

3 計算結(jié)果及對比分析

圖8給出了計算所得的束流離子數(shù)密度分布，基本分布在軸線±15°附近，而邊緣附近的鋸齒狀圖案是由于粒子數(shù)偏少導(dǎo)致的局部統(tǒng)計漲落。

圖8 計算所得離子數(shù)密度分布Fig.8 Simulation results of ion density distribution

圖9給出的是本文模擬計算結(jié)果與文獻(xiàn)[13]試驗測量結(jié)果的比較?？梢婋x子數(shù)密度在不同徑向距離處其值都隨軸向距離的增加而減小；在靠近推力器中心線位置處，數(shù)密度值的變化幅度較大，在遠(yuǎn)離推力器中心線位置處，變化幅度較小。在徑向距離140 mm的曲線和徑向距離180 mm的曲線上均出現(xiàn)了明顯偏離主束流特征的數(shù)值，說明這些位置點已經(jīng)不在主束流區(qū)。模擬結(jié)果和試驗結(jié)果的趨勢與數(shù)值范圍基本一致。

圖9 離子數(shù)密度的軸向分布的模擬結(jié)果和試驗結(jié)果對比Fig.9 Comparison of axial profiles of ion density between simulation and experimental results

圖10給出的是在整星狀態(tài)下計算得到的羽流電勢分布，左側(cè)為衛(wèi)星，右側(cè)為太陽電池陣。對于這種復(fù)雜結(jié)構(gòu)組合形成的計算區(qū)域中，非結(jié)構(gòu)網(wǎng)格方法更顯示其優(yōu)勢。

圖10 羽流電勢Fig.10 Plume potential

4 結(jié)束語

本文介紹了采用三維非結(jié)構(gòu)網(wǎng)格對離子推力器羽流進(jìn)行數(shù)值模擬的方法。模擬結(jié)果與試驗測量數(shù)據(jù)進(jìn)行了比較，數(shù)據(jù)基本一致，證明了基于非結(jié)構(gòu)網(wǎng)格方法進(jìn)行羽流模擬的可行性。本方法在求解邊界情況更加復(fù)雜的航天器構(gòu)型的周圍等離子體環(huán)境時會更有優(yōu)勢。

（References）

[1]周志成, 高軍.全電推進(jìn) GEO衛(wèi)星平臺發(fā)展研究[J].航天器工程, 2015, 24(2): 1-6 ZHOU Z C, GAO J.Development approach to all-electric propulsion GEO satellite platform[J].Spacecraft Engineering, 2015, 24(2): 1-6

[2]胡照, 王敏, 袁俊剛.國外全電推進(jìn)衛(wèi)星平臺的發(fā)展及啟示[J].航天器環(huán)境工程, 2015, 32(5): 566-570 HU Z, WANG M, YUAN J G.A review of the development of all-electric propulsion platform in the world[J].Spacecraft Environment Engineering, 2015, 32(5): 566-570

[3]孫安邦, 毛根旺, 陳茂林, 等.離子推力器羽流特性的粒子模擬[J].強(qiáng)激光與粒子束, 2010, 22(2): 401-405 SUN A B, MAO G W, CHEN M L, et al.Particle simulation of ion thruster’s plume characteristics[J].High Power Laser and Particle Beams, 2010, 22(2): 401-405

[4]李娟, 楚豫川, 曹勇.離子推力器羽流場模擬以及Mo~＋CEX 沉積分析[J].推進(jìn)技術(shù), 2012, 33(1): 131-137 LI J, CHU Y C, CAO Y.Numerical simulations of ion thruster plume contamination interactions on spacecraft[J].Journal of Propulsion Technology, 2012, 33(1): 131-137

[5]任軍學(xué), 顧左, 郭寧, 等.離子發(fā)動機(jī)羽流特性的數(shù)值模擬[J].航空動力學(xué)報, 2013, 28(6): 1372-1379 REN J X, GU Z, GUO N, et al.Numerical simulation of characteristics of ion thruster plume[J].Journal of Aerospace Power, 2013, 28(6): 1372-1379

[6]任軍學(xué), 王艷, 仇釬, 等.離子發(fā)動機(jī)羽流二維軸對稱數(shù)值模型與驗證[J].北京航空航天大學(xué)學(xué)報, 2011, 37(12): 1498-1503 REN J X, WANG Y, QIU Q, et al.Validation of two-dimensional axi-symmetric plume model for ion thruster[J].Journal of Beijing University of Aeronautics and Astronautics, 2011, 37(12): 1498-1503

[7]HU Y, WANG J J.Electron kinetic characteristics in plasma plumes: fully kinetic simulations[C]//53rdAIAA Aerospace Sciences Meeting.Kissimmee, Florida, 2015: 1-10

[8]王虹玥, 賀碧蛟, 蔡國飆.穩(wěn)態(tài)等離子體推力器羽流仿真及其對微波的衰減和相移作用[J].推進(jìn)技術(shù), 2011, 32(6): 839-844 WANG H Y, HE B J, CAI G B.Numerical simulation of the stationary plasma thruster plume and its effects on the microwave attenuation and phase shifts[J].Journal of Propulsion Technology, 2011, 32(6): 839-844

[9]王學(xué)德, 伍貽兆, 夏健.三維非結(jié)構(gòu)網(wǎng)格 DSMC方法的實現(xiàn)及其應(yīng)用[J].計算力學(xué)學(xué)報, 2007, 24(2): 231-235 WANG X D, WU Y Z, XIA J.Implementation of 3D unstructured DSMC method and its application[J].Chinese Journal of Computational Mechanics, 2007, 24(2): 231-235

[10]FOSTER J E, PATTERSON M, PENCIL E, et al.Neutralizer characterization of a next multi-thruster array with electrostatic probes[C]//42ndAIAA/ASME/ SAE/ASEE Joint Propulsion Conference.Sacramento, California, 2006: 1-21

[11]楊集, 陳賢祥, 周杰, 等.伸桿對星載電場探測儀的影響研究[J].電子與信息學(xué)報, 2009, 31(4): 1010-1012 YANG J, CHEN X X, ZHOU J, et al.Study of the influence of booms on spaceborne electric field sensors[J].Journal of Electronics and Information Technology, 2009, 31(4): 1010-1012

[12]于二青, 王春江, 趙金城.矩陣重排序算法在結(jié)構(gòu)分析快速求解中的應(yīng)用[J].空間結(jié)構(gòu), 2010, 16(1): 45-50 YU E Q, WANG C J, ZHAO J C.Application of improved RCM algorithm in fast solution of structural analysis[J].Spatial Structures, 2010, 16(1): 45-50

[13]張尊, 湯海濱.氙氣離子推力器束流等離子體特征參數(shù)的 Langmuir單探針診斷[J].高電壓技術(shù), 2013, 39(7): 1602-1608 ZHANG Z, TANG H B.Single Langmuir probe diagnostics for characteristic parameters of the beam plasma in Xe ion thruster[J].High Voltage Engineering, 2013, 39(7): 1602-1608

[14]任軍學(xué), 李娟, 仇釬, 等.離子發(fā)動機(jī)交換電荷離子返流的粒子模擬[J].強(qiáng)激光與粒子束, 2011, 23(7): 1929-1934 REN J X, LI J, QIU Q, et al.Particle simulation of charge-exchange ions back flow in ion thruster plume[J].High Power Laser and Particle Beams, 2011, 23(7): 1929-1934

（編輯：肖福根）

A method for numerical simulation of ion thruster's plume based on unstructured mesh

LIN Xiaoxiong, WEN Zheng, TAO Jiasheng
(Institute of Telecommunication Satellite, China Academy of Space Technology, Beijing 100094, China)

A method of particle simulation based on unstructured mesh is proposed in this paper to reveal the basic characteristics of ion thruster’s plume with complex spacecraft configurations.The method of grid discretization and the topological relationship among the generation of grids are discussed.An algorithm of particle motion is developed to push the particles in grids, then the method of matrix permutation based on the RCM algorithm is used to reduce the bandwidth of the matrix.Finally, a simulation example is compared with experiment results to validate this method.

ion electric thruster; plume; numerical simulation; unstructured mesh

V439.2

：A

：1673-1379(2016)05-0484-06

10.3969/j.issn.1673-1379.2016.05.005

林驍雄（1990—），男，碩士研究生，研究方向為電推進(jìn)系統(tǒng)總體設(shè)計。E-mail: linxiaoxiong_1990@163.com。

2016-03-16；

：2016-09-03