基于PIC-DSMC 算法的嵌套霍爾推進器仿真性能研究

譚 睿，于 博，王平陽

（1.上海交通大學 機械與動力工程學院，上海 200240；2.上海空間推進研究所 電推進事業部，上海 201112）

0 引言

霍爾推進器是一種技術成熟的電推進器，適用于各種各樣的太空任務［1-2］。近年來，隨著航天技術的進一步發展，對大推力、高比沖的霍爾推力器的需求日益提升［3-4］。然而傳統單通道霍爾推進器受到自身幾何條件的局限，其尺寸和運行功率受到限制［5］。將推進器進行并聯集成，可以利用現有型號的推進器，但是存在尺寸過大且推力密度低的問題［6］。采用嵌套形式的霍爾推進器是合理的解決方案。美國國家航天局（NASA）和格林研究中心（Green Research Center）設計研制的第一個10 kW雙流道嵌套霍爾推進器X2 取得了點火成功和相關數據測量［7］。基于嵌套推進器X2 和多個傳統霍爾推進器型號的實驗數據，研制的三通道嵌套霍爾推進器X3 目前正處于地面實驗階段［8-9］。

嵌套霍爾推進器以傳統單通道霍爾推進器為原型，將多個通道嵌套在同一推進器上，嵌套霍爾推進器能夠實現大功率運行，以及任何通道之間的組合［10-11］。與此同時，由于地面試驗受到成本、設備等因素的制約，最優運行工況往往難以找到，而探索不同的配置和參數，找到最佳運行工況是優化推進器結構的重要手段［12］。數值模擬在很大程度上解決了這個問題，同時也彌補了地面試驗成本高、有背壓等諸多不利因素，成為推進器可行性驗證和設計優化的重要方法［13］。由于目前國內嵌套霍爾推進器還處于初步發展階段［14］，數值模擬結果對嵌套推進器型號設計和實驗平臺搭建起到重要的作用。

PIC（Particle-in-Cell）和流體方法是確定所有帶電物種空間分布的2 種主要數值方法［15-16］。PIC方法相對于流體方法具有較高的計算精度，因為該模型采用基本方程沒有過度簡化，跟蹤真實的物理過程，包括粒子傳輸、粒子碰撞和電場計算［17］。靜電場和磁場通過求解泊松方程和麥克斯韋方程來計算［18］。本文所采用的幾何和磁場結構參考了NASA-457M，NASA-400M 以 及X2 等多種型號［19-20］。仿真羽流場分布結果對于磁場設計改進，幾何結構進一步優化和實驗設備搭建具有重要的參考意義。

1 PIC-DSMC 算法概述

1.1 PIC 算法

PIC 方法是等離子體物理中粒子模擬的最重要方法之一。國內外也有很多學者采用了有限粒子法對等離子體流動問題進行模擬［21，22］。在PIC 模擬中，通常求解在網格點上進行，網格內的仿真粒子基于特定原則被分配到節點上。一般情況下對平面區域，主要應用面積權重法，而三維時則選擇體積權重法。PIC 方法中帶電粒子的電荷分配和仿真流程如圖1 所示。

圖1 PIC 方法中帶電粒子的電荷分配和PIC 算法仿真流程Fig.1 Charge assignment of charged particles to nodes and simulation program of a single time step in the PIC method

電荷分配給網格節點的電量為

式中：qi為第i個節點分配到的電量；q為電荷電量。

通過PIC 算法進行模擬時的流程如圖1（b）所示：1）設有大量帶電粒子，對其位置和速度參數進行初始化，基于統計平均的方法計算求解域中的電磁場；2）獲得每個粒子在該時刻電磁場作用下所受到的電場力和洛倫茲力；3）接著計算出其加速度和速度；4）循環迭代，達到穩定狀態后對粒子統計平均輸出結果。

1.2 DSMC 算法

DSMC（Direct-Simulation-of-Monte-Carlo-Collision）是為了求解等離子體粒子時間碰撞而引入的算法［23-24］。DSMC 方法是基于如下幾方面假設建立的［25］。1）假設在所有計算網格內均發生二元碰撞；2）分子平均間距遠大于分子維度，僅在碰撞瞬間考慮分子間作用力；3）碰撞后分子的運動狀態只可基于隨機抽樣方法計算出。

1.3 PIC-DSMC 算法

鑒于電推力器羽流屬于超音速稀薄等離子體［26］，其實際流動過程極其復雜。單一的PIC（不考慮粒子間的動量和電荷交換碰撞）或DSMC 法（不考慮電磁場作用）都不能有效模擬其流動問題［27］。20 世紀90 年代，OH 等［28-29］結合上述2 種方法，提出了PIC-DSMC 混合方法，并運用該法對SPT-100、NASA GRC 實驗室SPT 等離子體羽流穩態和非穩態進行流動仿真，并取得令人滿意的效果。仿真中DSMC 模擬與PIC 模擬相對獨立，前者主要用于計算粒子的運動及碰撞，而后者主要用于模擬等離子體自洽電場及其對帶電粒子的加速作用，兩者涉及的粒子加速、排序、運動等能在一個時間步內同時完成，使得PIC-DSMC 混合法編程過程易于實現。此外，DSMC 仿真在單元進行，而PIC 仿真則在網格點上進行［29］。

PIC-DSMC 算法流程如圖2 所示。

2 幾何條件和邊界條件

嵌套霍爾推進器有內外2 個通道，共有3 種工作模式，推進工質為氙氣。與傳統單通道霍爾推進器相比，嵌套推進器工作范圍更廣，推進力更高，能適應未來航天發展多任務化的需求。仿真共采用20 組算例，研究在內外流道采用不同質量流量，以及磁場強度條件下，對推進器羽流流場分布乃至推進器整體推進性能的影響。推進器幾何結構如圖3 所示。

為了探究在不同工況下嵌套霍爾推進器的推進性能，總共采用20 組工況分別研究內外通道在一系列質量流量與磁場強度下的羽流場特征。質量流量和磁場強度的選擇參考了NASA-457M，NASA-400M 以及X2 等多種型號，并通過初步仿真篩選出潛在的最佳工況組合。20 組工況見表1。

表1 自變量為質量流量和磁場強度的最佳工況陣列Tab.1 Optimal working condition array with independent variables of the magnetic field strength and mass flow rate

為提高程序計算效率，如圖4 所示，對求解域進行劃分而得到7 個子區域，其中各子區域的網格長度為等比級數，各網格含4 個子網格，進行劃分后總共獲得336 00 個矩形計算網格。計算網格的最小尺度為0.1 mm×0.1 mm，而嵌套霍爾推進器羽流出口處的碰撞的平均自由程約為1 m，德拜長度約為5×10-5m［12］，因此模型基本滿足等離子體二元碰撞假設。粒子運動時間步長取10-8s，滿足等離子體頻率限制和碰撞過程仿真要求。在流場達到穩態后，取100 次抽樣均值當做為一組結果，共輸出50次，最后再取均值。仿真中總共使用了240 000 個仿真粒子（不包含背壓粒子）。

圖4 計算網格與邊界條件Fig.4 Computational grid and boundary conditions

本次仿真過程中，內通道內外半徑分別為51.5、84.5 mm；外通道內外半徑分別為151.1、198.5 mm。模擬區域采用二維旋轉對稱模型，旋轉對稱軸虛線表示如圖4（b）所示，對稱軸處采用反射邊界條件，徑向電場強度等于0；固體壁面均采用反射邊界條件，電勢取0；計算區域長度取1 000 mm，計算區域半徑取600 mm，邊界面均認為是與真空交界面，采用出流邊界條件，且出垂直于邊界面的電場強度取0。回流區域長度取300 mm。

3 仿真結果及分析

3.1 羽流場結果云圖

本次模擬結果所展示云圖均為Case1 工況條件下結果云圖，如圖5 所示，圖5 中，橫軸表示沿著推力器軸線方向的距離，m；縱軸表示距離推力器中心的距離，m。

圖5 Case1 工況條件下結果云圖（長度單位：m）Fig.5 Result contours under Case1 working condition（m）

3.2 羽流發散角

嵌套霍爾推進器羽流發散半角為25°～33°。從發散角總體變化趨勢來看。隨著推進器質量流量減小，推進器羽流發散角會減小。在保持質量流量恒定時，內通道磁場強度適當增大有利于羽流發散角減小，外流道磁場強度增大則會導致羽流發散角增大。因此，一定范圍內提高內通道的質量流量和磁場強度有利于約束羽流發散角，提高推進效率，從而改善推進器整體性能。羽流發散角隨質量流量和最大磁場強度的變化如圖6 所示。

圖6 不同工況下的羽流發散角Fig.6 Plume divergence angles under different working conditions

3.3 推力和比沖

推進力和比沖受質量流量、磁場強度的影響見表2。

表2 推進力和比沖受質量流量、磁場強度的影響Tab.2 Effects of the mass flow and magnetic field strength on the propulsion and specific impulse

在相同磁場強度下，推進力與質量流率成正相關。在相同質量流量下，推進力與磁場強度成正相關。推進力大小為1.5～2.5 N，比沖為2 300～3 200 s。

4 結束語

推進力、比沖和羽流發散角是衡量推進器推進性能的重要物理參數。推進力可以衡量推進器在短時間內的最大工作負荷。比沖和羽流發散角可以衡量推進工質利用效率。本次數值模擬針對嵌套霍爾推進器羽流進行數值仿真，得到雙通道質量流率和最大磁場強度對推進性能的影響。

本次模擬推進力大小為1.5～2.5 N，與磁場強度和質量流率成正相關。推進器比沖為2 300～3 200 s，與質量流量和磁場強度成正相關。質量流率越大，經過電離而產生的正價粒子也越多，推進力越大。同時，在特定條件下找到合適的最大磁場強度能夠實現推進器的高效率運行。

本次模擬半羽流發散半角為25°～33°。嵌套霍爾推進器羽流發散角高于一般單通道霍爾推進器，這是由嵌套霍爾推進器特殊的幾何結構所決定的。隨著質量流率的減小，羽流發散角也會隨之減小。在保持質量流量恒定時，內通道磁場適當增大有利于羽流發散角減小，外通道磁場強度增大則會導致羽流發散角增大。因此在設計嵌套霍爾推進器時，一定范圍內提高內通道的質量流量和磁場強度有利于提升推進效率，改善推進性能。

基于此次仿真結果，后續將進一步增大內通道質量流量和磁場強度進行羽流仿真，從而達到進一步約束羽流發散角的目的。同時，在該型號嵌套霍爾推進器完成點火實驗并進行性能評估后，會進一步與本次仿真結果和國外X2 等現有型號進行性能對比，從而對推力器結構進行進一步的優化。