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

譚 睿，于 博，王平陽

（1.上海交通大學(xué) 機械與動力工程學(xué)院，上海 200240；2.上?？臻g推進研究所 電推進事業(yè)部，上海 201112）

0 引言

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

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

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

1 PIC-DSMC 算法概述

1.1 PIC 算法

PIC 方法是等離子體物理中粒子模擬的最重要方法之一。國內(nèi)外也有很多學(xué)者采用了有限粒子法對等離子體流動問題進行模擬［21，22］。在PIC 模擬中，通常求解在網(wǎng)格點上進行，網(wǎng)格內(nèi)的仿真粒子基于特定原則被分配到節(jié)點上。一般情況下對平面區(qū)域，主要應(yīng)用面積權(quán)重法，而三維時則選擇體積權(quán)重法。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

電荷分配給網(wǎng)格節(jié)點的電量為

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

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

1.2 DSMC 算法

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

1.3 PIC-DSMC 算法

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

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

2 幾何條件和邊界條件

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

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

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

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

圖4 計算網(wǎng)格與邊界條件Fig.4 Computational grid and boundary conditions

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

3 仿真結(jié)果及分析

3.1 羽流場結(jié)果云圖

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

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

3.2 羽流發(fā)散角

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

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

3.3 推力和比沖

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

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

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

4 結(jié)束語

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

本次模擬推進力大小為1.5～2.5 N，與磁場強度和質(zhì)量流率成正相關(guān)。推進器比沖為2 300～3 200 s，與質(zhì)量流量和磁場強度成正相關(guān)。質(zhì)量流率越大，經(jīng)過電離而產(chǎn)生的正價粒子也越多，推進力越大。同時，在特定條件下找到合適的最大磁場強度能夠?qū)崿F(xiàn)推進器的高效率運行。

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

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