范益朋，夏廣慶*，鹿 暢，孫 斌，丁 亮

（1. 大連理工大學 工業(yè)裝備結構分析國家重點實驗室； 2. 大連理工大學 遼寧省空天飛行器前沿技術重點實驗室：大連 116024； 3. 北京衛(wèi)星環(huán)境工程研究所，北京 100094）

0 引言

離子液體推力器（ionic liquid thruster, ILT）是在微小衛(wèi)星的需求下迅速發(fā)展的一種新型微牛級微電推力器，其利用靜電噴射原理，具有結構精簡、體積小、重量輕、功率低、推力精度高等特點[1]，可用作微納衛(wèi)星的動力系統(tǒng)，或作為分布式推力器應用于大型航天器[2]。作為推進劑的離子液體在室溫下通常呈液態(tài)，且表面張力較小，電導率較高，故發(fā)射時無需為推進劑離解提供輔助部件；同時，離子液體黏度較高，能提供較大的流阻，適用于毛細作用被動式供給，無須壓力供給系統(tǒng)；此外，離子液體可交替發(fā)射陰、陽離子，實現(xiàn)電荷自中和，無須配置中和器[3-4]。因此，離子液體推力器的應用前景極為廣泛。

目前關于離子液體推力器的仿真研究主要是針對從離子液體中發(fā)射離子的過程，即關于泰勒錐形成過程與泰勒錐破碎發(fā)射離子的過程。例如：Borner 等利用分子動力學（molecular dynamics,MD）仿真模擬了不同電場強度和推進劑質量流量下發(fā)射電流的大小以及各種離子所占比例[5-6]；程世豪通過COMSOL 仿真軟件研究了發(fā)射極陣列密度、電極間距、發(fā)射極內徑對離子發(fā)射所需閾值電壓的影響[7]；Kim 將有限元法（FEM）與有限差分法（FDM）結合到MD 仿真中，分析了引出極孔徑、發(fā)射電壓、發(fā)射極溫度和潤濕性等因素對發(fā)射電流大小、電流密度和離子速度的影響[8]。但截至目前關于離子液體推力器的研究并未涉及推力器結構對整體性能影響的多工況粒子網(wǎng)格（particle in cell, PIC）仿真分析。

本文著重研究離子液體推力器發(fā)射陣列幾何結構參數(shù)對推力器性能的影響：首先基于PIC 法建立單個發(fā)射極的二維仿真模型；而后分析電極間距與引出極孔徑改變引起的推力器推力和角向效率變化規(guī)律，以期為離子液體推力器的優(yōu)化設計提供參考。

1 仿真模型

1.1 ILT 工作原理

離子液體推力器主要由發(fā)射極、引出極和推進劑儲箱等組成，單個發(fā)射極結構如圖1 所示。其基于靜電噴射原理，通過在發(fā)射極和引出極間施加發(fā)射電壓，使發(fā)射極尖端產生強電場，發(fā)射極尖端的離子液體受靜電場力和表面張力的共同作用變形成泰勒錐[9]（Taylor cone），隨后泰勒錐尖端破碎噴射出帶電液滴以及溶劑化和未溶劑化的離子，并在電場中加速噴出。部分推力器在引出極后增添了加速極，實現(xiàn)對引出離子的進一步加速，以提高推力器推力和比沖[10]。

由于單個發(fā)射極產生的推力非常低，通常在十幾nN 量級，難以滿足實際應用需求，因此需要將發(fā)射極陣列化，以獲得較大的發(fā)射電流和推力[11]。

1.2 仿真模型建立

陣列化后的發(fā)射極有多種類型：按陣列形式可分為一維陣列式和二維陣列式；按推進劑供給形式可分為毛細管型、外部浸潤型和多孔基質型。其中，毛細管型發(fā)射極泰勒錐形成位置和離子發(fā)射位置相較于其他兩種類型發(fā)射極更容易確定，且較為穩(wěn)定，適合通過仿真探索一些規(guī)律性問題。離子液體推力器常見的工作模式是錐-噴流（cone-jet）模式，在此模式下，泰勒錐的母線可能是直線也可能是曲線[12]，并且隨著使用液體的電導率的減小，噴流的區(qū)域逐漸向圓錐底部延伸[7]。因本文主要針對推力器宏觀性能進行研究，故此處考慮母線為直線形式的標準泰勒錐模型，對單個毛細管型發(fā)射極建立了二維軸對稱模型，如圖2 所示，具體結構參數(shù)及工況如表1 所示。

圖2 單個毛細管發(fā)射極的仿真模型Fig. 2 Simulation model of a single capillary emitter

考慮工作在陽極工況下的發(fā)射極模型，其發(fā)射電壓設置為3 kV，采用1-乙基-3-甲基咪唑四氟硼酸鹽（1-ethyl-3-methylimidazolium tetrafluoroborate,EMI-BF4）作為推進工質。發(fā)射極工作在純離子模式（purely ionic regime, PIR）下，即發(fā)射物中只包含單體（EMI+）、二聚體（[EMI-BF4]EMI+）、三聚體（[EMI-BF4]2EMI+）三種離子。離子總發(fā)射電流320 nA，各種離子所占比例取于Romero-Sanz 的實驗結果[13]，如表2 所示。

表2 各種離子所占電流組分Table 2 Proportions of various ions in the current

上述模型基于PIC 方法建立。PIC 方法作為一種計算機數(shù)值模擬方法，通過追蹤單個粒子的運動，對單個微觀粒子的運動情況進行數(shù)值描述，再對組成物體的大量微觀粒子進行統(tǒng)計平均，由此可得到宏觀物體的物質特性和運動規(guī)律，可用于推力器性能研究。

如圖3 所示，PIC 方法需要通過網(wǎng)格劃分的方式將力場（一般為電場、磁場）離散化，同時將時間變量離散為一個個時間步長，分析粒子在當前時刻所受離散化力場的作用力，便可求得下一個時間步粒子位置與速度；進而再與外力場耦合，解算出新的力場；以此循環(huán)，直至系統(tǒng)達到穩(wěn)態(tài)。

圖3 PIC 方法流程圖Fig. 3 Flow chart of PIC method

其耦合過程的核心思想是將網(wǎng)格內粒子的電量、粒子運動產生的電流等參數(shù)按一定的權重（如面積權重法、體積權重法等）分配至網(wǎng)格節(jié)點上，再計算網(wǎng)格節(jié)點上的力場參數(shù)[14]。

由于離子液體推力器為靜電式電推力器，發(fā)射場中只包含靜電場，故在PIC 算法中求解電場時不必求解完整的Maxwell 方程組，只需求解Poisson方程[15]

式中E為電場矢量。

仿真模型中，除發(fā)射極與引出極外，其余邊界條件均設置為第二類邊界條件（諾依曼邊界條件），即

對于帶電粒子在靜電場中的運動，一般可描述為

式中：m為帶電粒子質量；v表示粒子速度矢量；q為粒子所帶電量；dt為時間步長。

最后，將上述各式離散化[16]，以獲得各網(wǎng)格節(jié)點上的數(shù)值，并基于網(wǎng)格節(jié)點進行PIC 仿真計算。

由于離子發(fā)射所需電場強度的閾值約在109V/m量級[9]，電勢在發(fā)射極與引出極之間的變化梯度極大，因此為保證仿真結果的準確性，兼顧到仿真時長，將該模型網(wǎng)格大小控制在1×10-6m。

1.3 推力器性能計算方法

本文主要研究內容為電極間距與引出極孔徑改變引起的推力器性能變化，即以初始模型為基礎，分別改變模型電極間距和引出極孔徑，通過仿真，獲得二者改變引起的推力器推力及角向效率變化曲線。

推力器推力T為[17]

式中：η為推力器總效率；Iemit為發(fā)射電流；Uemit為發(fā)射電壓；q/m為羽流中包含粒子的平均荷質比。

而仿真時，離子液體推力器推力能夠通過對所有粒子的速度和質量進行統(tǒng)計獲得，即

式中：n為通過推力統(tǒng)計截面的離子數(shù)量；mi為第i個離子的質量；viz表示第i個離子沿軸向的速度。

推力器角向效率ηangle為羽流軸向動能Ekz與總動能Ek之比，即

其中，Ekz和Ek可分別表示為：

此外，還可通過編寫仿真程序獲得推進劑質量流量，進而計算出推力器比沖。推進劑質量為

再根據(jù)式(7)和式(11)即可計算出推力器比沖

式中g為重力加速度，取9.8 kg/(m·s-2)。

收集電流Ico計算方法與質量流量計算方法類似，可得

式中qi為第i個離子所帶電量。

2 結果與討論

2.1 初始仿真結果

在PIC 仿真達到穩(wěn)態(tài)后對結果進行輸出，可獲得離子數(shù)密度分布、離子速度分布圖像，以及推力、比沖曲線等。

圖4 為穩(wěn)態(tài)后的羽流離子數(shù)密度分布圖像。圖中羽流整體形貌呈圓錐狀發(fā)散，具有一定的束流散射半角。在發(fā)射極尖端附近離子數(shù)密度最大，達1019/m3量級；經電場加速后離子擴散到引出極右側的較大區(qū)域，此時軸線附近離子數(shù)密度降至1015/m3量級，且數(shù)密度從軸線沿徑向向外呈現(xiàn)先減后增再減的趨勢，羽流邊緣處離子數(shù)密度僅達1014/m3量級。圖5 為組成羽流的三種離子的數(shù)密度分布圖。可見：雖然由于各種離子的質量及發(fā)射數(shù)量不同，三種離子數(shù)密度分布存在一定差異，但離子整體分布規(guī)律不變，即：離子數(shù)密度在發(fā)射極尖端位置最高，其次在靠近軸線附近相對較高，整體分布呈沿徑向從軸線向外先減后增再減的規(guī)律。

圖4 離子數(shù)密度分布Fig. 4 Distributions of ion number density

圖5 三種離子數(shù)密度分布Fig. 5 Distributions of ion number density of EMI+, [EMI-BF4]EMI+ and [EMI-BF4]2EMI+

圖6 為電場強度分布圖。從圖中可以看出：引出極左側的電場強度整體上比引出極右側高很多，在發(fā)射極尖端及引出極左側尖端，電場強度更是比其他區(qū)域大了幾個數(shù)量級。這表明離子能在通過引出極前即能夠將勢能轉化為動能，實現(xiàn)較為完全的加速。如圖7 所示，離子在通過引出極前整體平均速度已加速至約40 000 m/s；而電場矢量主要從發(fā)射極指向引出極，離子在電場中加速后，除了獲得軸向的速度外，還將獲得徑向的速度。

圖6 電場強度分布Fig. 6 Electric field intensity distributions

圖7 離子速度分布Fig. 7 Ion velocity distributions

整體而言，離子在電場中加速后獲得的速度以軸向速度為主，如圖8 所示，其平均軸向速度可達約45 000 m/s；相比于軸向速度，徑向速度最高僅約22 000 m/s。從圖9 的電場強度分布圖可以看出，軸向電場強度整體大于徑向電場強度，因此離子在電場中加速后能獲得高于徑向速度的軸向速度。而羽流邊緣附近的離子受電場加速后，其軸向速度雖低于羽流內部區(qū)域的離子速度，但徑向速度比羽流內部區(qū)域離子的要高，因此羽流邊緣離子的總速度相較于軸線附近的離子速度更高。

圖8 離子軸向和徑向速度分布Fig. 8 Ion velocity distributions in axial and radial directions

圖9 電場強度的軸向和徑向分布Fig. 9 Electric field intensity distributions in axial and radial directions

另外由于三種離子比荷不同，其在電場中的加速效應也不一樣（如圖10 所示），到達推力統(tǒng)計截面的時間會有差異。由圖可見：單體離子比荷最高，在電場中加速后速度也最大，可達約60 000 m/s，能夠最早到達推力統(tǒng)計截面；其次是二聚體，電場能將其加速至35 000 m/s 以上；最后是三聚體，受電場加速后速度約在30 000 m/s，最晚到達推力統(tǒng)計截面。故推力曲線呈現(xiàn)圖11 所示的明顯的階梯狀上升現(xiàn)象。

圖10 三種離子速度分布Fig. 10 Ion velocity distributions of EMI+, [EMI-BF4]EMI+and [EMI-BF4]2EMI+

圖11 推力和比沖曲線Fig. 11 Curve of thrust and specific impulse

而離子液體推力器被證明在發(fā)射極陣列化后，其單個發(fā)射極的工作狀況幾乎不受影響，因此可近乎成比例地提升推力[18]。如圖11 所示，單個發(fā)射極的推力約為34 nN，考慮在1 cm2內集成480 個發(fā)射極，則其推力可達約16 μN/cm2，比沖保持在約4100 s。

2.2 電極間距影響

在初始模型基礎上，保持其他參數(shù)不變，將電極間距由80 μm 逐漸增加至300 μm，獲得平均推力及角向效率隨電極間距變化曲線，如圖12 所示?？梢姡涸陔姌O間距增加初期，推力增幅相對較大；隨著電極間距的增大，推力增幅減?。辉陔姌O間距約280 μm 時推力達到最大值，約34.7 nN，隨后出現(xiàn)減小的趨勢；而角向效率隨電極間距單調增加。

圖12 平均推力及角向效率隨電極間距變化曲線Fig. 12 Variations of average thrust and angular efficiency against the electrode spacing

在確保未出現(xiàn)羽流被引出柵極攔截的情況下，對該現(xiàn)象產生的原因進行了分析，如圖13 所示。由于電場矢量主要由發(fā)射極指向引出極，在增加電極間距后，整體上電場矢量與軸線的夾角減小，即離子通過電場加速后軸向速度占總速度的比例增加，羽流的準直性有所提升，因此角向效率隨電極間距的增加而不斷增大；但電極間距的增加同樣會導致電場強度減小，這將導致離子加速后的速度不斷減小，因此隨著電極間距的不斷增加，離子經電場加速后獲得沿軸向的動能所占比例不斷增加，但離子總動能逐漸減小，在這兩種因素的耦合作用下，離子軸向速度先增加后減小，在280 μm 時獲得最大平均軸向速度，反映在推力上即為推力出現(xiàn)先增加后減小的趨勢，在280 μm 時獲得最大推力。

圖13 不同電極間距的離子軸向速度Fig. 13 Axial velocity of ions with different electrode spacing

2.3 引出極孔徑影響

通過改變初始模型引出極孔徑，獲得平均推力及角向效率隨引出極孔徑的變化曲線，如圖14所示。

圖14 平均推力及角向效率隨引出極孔徑變化曲線Fig. 14 Variations of average thrust and angular efficiency against the extractor aperture

隨著引出極孔徑由初始模型的300 μm 逐漸減小至80 μm，推力器推力逐漸增加，而角向效率卻不斷下降；當引出極孔徑＜120 μm 時，推力曲線與角向效率曲線出現(xiàn)驟變。這是由于引出極孔徑減小至120 μm 以下時出現(xiàn)了大量離子被引出柵極攔截的情況，造成明顯的推力損失與柵極侵蝕，嚴重影響推力器壽命，故應盡量避免出現(xiàn)該情況。

研究離子未被引出柵極攔截時引出極孔徑變化對平均推力和角向效率的影響，結果如圖15 和圖16 所示。在此情況下，引出極孔徑的減小使得電場強度顯著增加，從而使離子在電場中加速后獲得的軸向動能與徑向動能均有提升，而其徑向速度增幅相較于軸向速度更為顯著，即離子軸向動能占總動能的比例隨引出極孔徑的減小而逐漸減小。因此，在出現(xiàn)羽流被柵極攔截的情況之前，其角向效率隨引出極孔徑的減小而逐漸減?。欢S向動能占比雖然隨引出極孔徑的減小不斷減小，但總動能不斷增加，且增幅顯著，故推力整體上呈現(xiàn)隨引出極孔徑的減小而不斷增加的趨勢。

圖15 不同引出極孔徑的離子軸向速度Fig. 15 Axial velocity of ions with different extractor apertures

圖16 不同引出極孔徑的離子徑向速度Fig. 16 Radial velocity of ions with different extractorapertures

3 結論

本文基于PIC 方法建立了離子液體推力器單個毛細管型發(fā)射極的仿真模型，獲得了單個發(fā)射極的束流分布特性和性能參數(shù)曲線，隨后著重研究了電極間距和引出極孔徑對推力器推力和角向效率的影響，獲得了發(fā)射極平均推力和角向效率隨電極間距及引出極孔徑的變化規(guī)律。仿真結果表明：

1）推力器角向效率與電極間距成正比，推力則受電場強度和電場方向的共同影響，隨電極間距的增加先增大后減?。?/p>

2）在出現(xiàn)離子被引出柵極攔截的情況前，推力器角向效率與引出極孔徑成正比，推力與引出極孔徑成反比，但過小的引出極孔徑會造成柵極侵蝕，嚴重影響推力器性能和壽命；

3）在上百個發(fā)射極實現(xiàn)陣列化之后，發(fā)射場幾何結構的改變引起的推力變化預計達μN 量級。

由于本研究未考慮電場變化引起發(fā)射離子數(shù)量的變化和不同種類離子所占比例的變化，因此推力曲線與實驗存在一定差異，但推力器角向效率變化規(guī)律不受發(fā)射離子數(shù)量和比例影響。在后續(xù)研究中將考慮電場變化引起發(fā)射離子數(shù)量和比例的變化，耦合多種影響推力器性能的因素，使仿真結果更貼合實際。