基于非結構網格的離子推力器羽流數值模擬方法研究

林驍雄,溫 正,陶家生

（中國空間技術研究院 通信衛星事業部，北京 100094）

基于非結構網格的離子推力器羽流數值模擬方法研究

林驍雄,溫 正,陶家生

（中國空間技術研究院 通信衛星事業部，北京 100094）

為適應復雜的航天器外形，文章采用基于非結構網格的數值模擬方法研究離子推力器真空羽流的基本特征。首先，介紹了非結構網格劃分的方法，并分析了建模區域生成各種拓撲關系的特點；然后，提出了粒子穿越網格的具體算法，包括有利于減少計算量的結合矩陣重標號的帶寬壓縮算法；最后，開展了模擬結果和試驗結果的對比研究，驗證了以上數值模擬方法的可行性。

離子電推力器；羽流；數值模擬；非結構網格

0 引言

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

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

1 物理及數值模型

1.1 羽流物理模型

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

電場作用方程為

式中：ε0為真空介電常數；ρ為空間電荷密度；E為電場強度。

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

通過求解此方程可以得到空間電勢，進而得到電場強度。

電磁力作用方程為

式中：F為電磁力；q為粒子電荷；v為速度；B為磁感應強度。

牛頓運動方程為

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

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

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

1.2 計算參數

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

表1 離子推力器工作參數Table 1 Ion thruster parameters

表1（續）

2 模擬計算

2.1 計算區域和計算網格

圖1為衛星上離子推力器的安裝位置示意圖，有2個推力器固定在矢量推力機構上。

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

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

圖2 模擬區域和推力器及中和器局部網格Fig.2 Simulation domain and local grids around thruster and neutralizer

2.2 網格預處理

網格預處理主要用來生成非結構網格間各種層次關系，而層次關系是指求解區域（Group）、四面體單元（Cell）、三角形表面（Face）、棱邊（Edge）、節點（Node）等拓撲元素的鏈接關系。在羽流的數值模擬中主要包含以下幾種：

1）在縮減矩陣帶寬中需要用到計算區域內邊到節點的映射Edge2Node。

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

3）在粒子穿越網格算法中需要用到四面體單元到三角形單元的映射Cell2Face、三角形單元到四面體單元的映射Face2Cell、表面三角形單元局部標號到全局三角形單元標號的映射 ThrFace2Face、ScFace2Face、BdFace2Face。這些拓撲關系數組在結果后處理過程中也會用到，提前生成這些拓撲關系數組就可以在粒子運動過程中直接查表，極大提高了計算速度。拓撲關系如圖3所示。

圖3 網格間的拓撲關系Fig.3 Topological relationship among grids

2.3 粒子穿越網格算法

粒子穿越網格是非結構網格中模擬粒子運動的關鍵步驟之一。在結構網格中，可以借助網格下標直接索引粒子將要穿越的相鄰網格，而在非結構網格中只能憑借網格層級關系進行搜索。

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

圖4 粒子穿越網格模型Fig.4 Model of particle traversing through the meshs

圖5 粒子運動更新算法Fig.5 Update algorithm of particle motion

2.4 Poisson方程的求解

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

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

2）以頂點r為根，生成圖G的層次結構{x1,x2, …,xn}（如圖中虛線區分的3個層，括號中的數值為節點在該層中的度）；

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

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

圖6 RCM算法示例網格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。由于本計算為示例，半帶寬并沒有明顯減少。但是在實際的應用中，因為內部節點數遠遠大于邊界節點數，半帶寬可以減少到原來的 1/5甚至更多。采取這種方法形成矩陣可以極大減少方程的運算次數，而且在用稀疏帶狀矩陣進行Cholesky分解時引入的非零元素都在帶寬及輪廓范圍內，不會過多占用額外內存，可以提高運算效率。

3 計算結果及對比分析

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

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

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

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

圖10給出的是在整星狀態下計算得到的羽流電勢分布，左側為衛星，右側為太陽電池陣。對于這種復雜結構組合形成的計算區域中，非結構網格方法更顯示其優勢。

圖10 羽流電勢Fig.10 Plume potential

4 結束語

本文介紹了采用三維非結構網格對離子推力器羽流進行數值模擬的方法。模擬結果與試驗測量數據進行了比較，數據基本一致，證明了基于非結構網格方法進行羽流模擬的可行性。本方法在求解邊界情況更加復雜的航天器構型的周圍等離子體環境時會更有優勢。

（References）

[1]周志成, 高軍.全電推進 GEO衛星平臺發展研究[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].航天器環境工程, 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].強激光與粒子束, 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].推進技術, 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]任軍學, 顧左, 郭寧, 等.離子發動機羽流特性的數值模擬[J].航空動力學報, 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]任軍學, 王艷, 仇釬, 等.離子發動機羽流二維軸對稱數值模型與驗證[J].北京航空航天大學學報, 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]王虹玥, 賀碧蛟, 蔡國飆.穩態等離子體推力器羽流仿真及其對微波的衰減和相移作用[J].推進技術, 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]王學德, 伍貽兆, 夏健.三維非結構網格 DSMC方法的實現及其應用[J].計算力學學報, 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].電子與信息學報, 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]于二青, 王春江, 趙金城.矩陣重排序算法在結構分析快速求解中的應用[J].空間結構, 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]張尊, 湯海濱.氙氣離子推力器束流等離子體特征參數的 Langmuir單探針診斷[J].高電壓技術, 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]任軍學, 李娟, 仇釬, 等.離子發動機交換電荷離子返流的粒子模擬[J].強激光與粒子束, 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—），男，碩士研究生，研究方向為電推進系統總體設計。E-mail: linxiaoxiong_1990@163.com。

2016-03-16；

：2016-09-03