脈沖等離子體推力器特氟隆燒蝕過程數值模擬①

謝澤華，周 進

(國防科技大學高超聲速沖壓發(fā)動機技術重點實驗室，長沙 410073)

0 引言

脈沖等離子體推力器(Pulsed Plasma Thruster，PPT)具有結構簡單、易于集成、控制方便、響應快速、能產生小而精確的推力脈沖、所需星上電源功率小等優(yōu)點，可廣泛用于微小衛(wèi)星的姿態(tài)控制和軌道保持等任務，是微小衛(wèi)星推進技術研究的熱點和重要方向之一［1］。PPT利用脈沖放電燒蝕或加熱推進劑，并使其電離形成等離子體，等離子體在電磁力和氣動壓力作用下加速噴出推力器，產生脈沖推力。根據所使用的推進劑，PPT有多種不同類型，其中最常見的是使用固體聚四氟乙烯(俗稱特氟隆)的燒蝕型脈沖等離子體推力器。在數十年的研究與應用過程中，PPT一直面臨著效率較低的問題［2］。作為推力器工作的重要階段，推進劑燒蝕對推力效率具有直接重大的影響。為了掌握PPT的工作機理，提升推力效率，有必要對推進劑燒蝕過程開展細致深入的研究。目前，已經發(fā)展了一些模型用于PPT中特氟隆推進劑的燒蝕計算。Brito等［3］提出了一個零維模型，將燒蝕質量表達為放電參數和推力器幾何參數的函數，適用于推力器的性能預估。Turchi等［4］在一維準穩(wěn)態(tài)高磁雷諾數高磁壓等離子體流動假設下，得到了用放電電流表示的燒蝕質量流率的表達式，可用于PPT工作過程的一維簡單分析。Thomas［5］和 Henrikson［6］在對 PPT 的研究過程中，采用過相似的簡化燒蝕模型，將燒蝕質量流率表示為熱流密度與特氟隆的蒸發(fā)潛熱之商，這一模型并不能得到推進劑溫度分布的信息，適合于只關心推力器總體性能和等離子體流動過程的研究場合。Mikellides［7］將特氟隆燒蝕視為一個物質升華的過程，假設在特氟隆表面存在一層飽和蒸氣，表面壓力由平衡蒸氣壓曲線計算，表面溫度根據熱傳導方程計算，燒蝕質量流率根據等離子體流動的動量方程由流動入口處的壓力梯度確定，Lin［8］和 Antonsen［9］則在各自研究中，利用平衡蒸氣壓曲線直接將燒蝕質量流率表示為特氟隆表面溫度的函數。同樣利用平衡蒸氣壓曲線，Keidar和Boyd［10］建立了一個動理學模型，假設推進劑表面存在一個動理學非平衡層和一個流體力學非平衡層，推進劑燒蝕質量流率根據推進劑表面溫度和等離子體密度、溫度計算。這些模型被應用到PPT的數值研究中，并取得了大量的研究成果。然而，由于它們并沒有反映出特氟隆燒蝕的具體過程，當前研究者們對PPT中推進劑燒蝕特性的認識還較有限。在將特氟隆用作熱防護材料的研究過程中，Clark［11］提出了第一個較為全面考慮真實物理化學過程的一維兩相燒蝕模型，將特氟隆劃分為晶體和無定形體兩種狀態(tài)。Holzknecht［12］建立了一個相似模型，考慮了熱膨脹效應和高分子化合物的生成，模型計算復雜，但與Clark的實驗結果吻合很好。Arai［13］簡化了燒蝕過程，在 Clark模型的基礎上，考慮了特氟隆對輻射的吸收和透射的影響，證明了兩相假設對精確預測特氟隆溫度分布和熱傳遞屬性有著關鍵作用。Clark、Holzknecht和Arai對特氟隆燒蝕防護做了大量具有重要意義的研究工作，他們建立的模型代表著目前對特氟隆燒蝕過程最詳細的分析。PPT的工作過程呈現出明顯的多維特性，而且推進劑燒蝕過程與等離子體流動傳熱過程聯系緊密相互影響。

為深入探究PPT的工作特點和性能規(guī)律，本文以Clark等人的工作為基礎，建立了特氟隆的二維兩相燒蝕模型，在與等離子體流動傳熱過程耦合計算的情況下，考察了PPT中特氟隆的燒蝕過程及對推力效率的影響。

1 計算模型

1．1 燒蝕傳熱方程

在低于600 K的溫度下，特氟隆是白色晶狀長鏈聚合物，溫度達到600 K時發(fā)生相變，成為具有極高粘性的無定形體，同時高聚物分子開始解聚、分解，生成單體小分子產物［13］。特氟隆燒蝕過程如圖1所示，在晶體和無定形體兩個區(qū)域分別考慮傳熱過程，給出導熱微分方程如下:

其中，ρ、T、k分別表示特氟隆的密度、溫度和導熱系數;下標c表示晶體;下標a表示無定形體;Qp表示單位體積的特氟隆分解產生的熱量，用Arrhenius反應定律表示為

式中 Ap為頻率因子;Hp為熱解熱;E為活化能。

圖1 燒蝕過程示意圖Fig．1 diagram of the ablation process

1．2 燒蝕傳熱邊界條件

在燒蝕表面上，凈傳入的熱流密度由外部傳入熱流密度、表面輻射損失及燒蝕質量帶走的能量確定，邊界條件可表示為

其中，下標s表示燒蝕表面;δ為當地表面切線方向與y軸的夾角;hs為燒蝕產物的比焓。

推進劑燒蝕發(fā)生在燒蝕表面附近非常薄的微米量級厚度的區(qū)域內，燒蝕表面的后退速度和燒蝕質量通量可根據特氟隆分解速率由下式計算:

式中 ρ0為特氟隆的參考密度;s和θ分別為燒蝕表面的后退距離和無定形體區(qū)域的長度。

在相界面上，推進劑的溫度等于相變溫度Tm。相界面的運動速度可根據熱流平衡條件求出，對其積分可計算出相界面的位置。

式中 ρm為相界面處的平均密度;hm為相變潛熱。

其余表面采用絕熱邊界條件。

1．3 等離子體流動與傳熱

大量的計算結果表明［7］，熱傳導是PPT中等離子體向推進劑傳遞熱量的主要方式。因此，忽略對流和輻射傳熱，假設傳給推進劑表面的熱流密度由表面附近的重粒子和電子的導熱系數及等離子體的溫度梯度確定。

等離子體的流動與傳熱屬性采用磁流體力學方法進行描述，控制方程組為

式中 ρ為密度;U為速度;B為磁場強度;p為壓強;et為總比能;νe為磁擴散率;τ=為粘性應力張量;q為熱流密度矢量;μ0為真空磁導率。

2 數值方法

2．1 坐標變換

推進劑燒蝕過程中，燒蝕表面和相界面的位置是隨時間變化的。為計算方便，作如下坐標變換:

經過變換后，新的坐標χ在無定形體區(qū)域的取值范圍和ξ在晶體區(qū)域的取值范圍均為［1，0］。在新的坐標系(χ，y)中，燒蝕傳熱方程(2)變換后的形式為

其中

燒蝕傳熱方程(1)在坐標系(ξ，y)中具有與式(10)相似的形式。此外，原坐標系中的邊界條件也需進行坐標變換，這里對變換結果予以省略。

2．2 離散求解

本文采用有限差分方法，對燒蝕傳熱方程進行求解。在特氟隆燒蝕過程中，無定形體區(qū)域的厚度非常小，為避免顯示格式計算時因穩(wěn)定性限制帶來的計算時間難以承受的問題，這里采用隱式格式進行計算。

在均勻網格上對時間導數項采用一階前向差分，空間導數項采用二階中心差分，對內熱源項進行泰勒級數展開，變換后的控制方程可離散成如下九點形式:

式中 a1，a2，…，a10為系數。

直接求解式(12)，需對一個大型稀疏矩陣作矩陣運算，難度很大。本文采用加權-對稱超松弛迭代算法(WSSOR)進行求解。對磁流體力學方程組引入雙曲型散度清除方法，以避免偽磁場散度對數值求解的影響，對無粘通量采用M-AUSMPW+格式進行差分，對粘性項采用中心差分，整個控制方程組采用基于隱式LU-SGS格式的雙時間步方法進行求解。

3 驗證算例

3．1 導熱過程

考慮單相物體內的一維非穩(wěn)態(tài)導熱問題。假設物性參數恒定，初始溫度為T0，當熱流密度q保持恒定時物體的溫度分布存在解析解［14］:

式中 α為熱擴散率。

采用二維程序進行一維計算，計算結果見圖2。圖2中，3個時刻的數值解與解析解的溫度分布曲線均高度重合，表明程序能對非穩(wěn)態(tài)熱傳導過程進行準確模擬。

圖2 導熱過程中的溫度分布Fig．2 Temperature distribution in the heat conduction process

3．2 相變過程

考慮相界面運動和相變潛熱對溫度分布的影響。假設特氟隆在初始時刻處于晶狀固相，初始溫度為相變溫度，t＞0時在其前端表面施加一個熱流密度保持為q的外部熱源，不考慮特氟隆的分解吸熱，假設相變后的區(qū)域內溫度分布是坐標的二次方函數，則在物性參數保持恒定的條件下，可得到溫度分布近似解［14］:

其中，系數μ和相界面的位置X由式(15)計算:

圖3給出了考慮相變過程的溫度分布曲線，圖3中的數值解和近似解較接近，說明程序能較好地對特氟隆的相變過程進行模擬。

圖3 相變過程中的溫度分布Fig．3 Temperature distribution in the phase transition process

3．3 燒蝕過程

假設燒蝕表面附近的溫度呈線性分布，特氟隆密度保持恒定，Kemp［15］給出了穩(wěn)態(tài)燒蝕情況下，由燒蝕表面溫度Ts確定燒蝕質量通量和燒蝕表面后退速度的近似計算公式:

式中 Bp為活化溫度;hs－h(huán)－∞表示推進劑燒蝕前后的焓差。

利用編寫的程序計算特氟隆在恒定熱流密度作用下燒蝕過程達到穩(wěn)定狀態(tài)時的燒蝕特征參數，將計算結果與據式(16)、式(17)得到的近似值進行比較，結果見圖4。圖4中的數值解與近似解符合較好，表明程序能較準確地計算出特氟隆的燒蝕特性。

圖4 穩(wěn)態(tài)燒蝕質量通量和燒蝕表面后退速度Fig．4 Steady-state ablated mass flux and surface recession velocity

4 計算結果及分析

4．1 計算條件

本文對常見的平行板電極PPT的典型工況進行模擬，PPT的結構參數及電參數見表1。

表1 PPT結構參數及電參數Table 1 Configuration and electric parameters of PPT

數值模擬的計算區(qū)域分為2部分，即特氟隆燒蝕傳熱區(qū)A和由放電通道B1及推力器出口下游延伸區(qū)B2組成的等離子體流動傳熱區(qū)，如圖5所示。A區(qū)特氟隆傳熱距離L0取為20 μm，B1區(qū)的等離子體流動入口參數由燒蝕過程和放電電流確定，電極表面處理為等溫滑移理想導體邊界，延伸區(qū)B2是自由流動空間。

圖5 計算區(qū)域Fig．5 Computational region

燒蝕過程數值模擬中，特氟隆的密度、比熱容和導熱系數都取為溫度的函數，初始溫度取為300 K，計算所需的物性參數取自Clark的文獻［11］。

4．2 結果分析

PPT的放電過程可等效為一個RLC電路放電過程，在放電產生的強脈沖電流作用下，推進劑被加熱分解，產生燒蝕。推進劑表面溫度關于推進劑中心位置近似呈對稱分布，圖6為推力器的放電電流波形和y方向上距離電極0．05、0．10、0．20、0．50 倍電極間距處的表面溫度變化曲線。放電開始時，放電電流急速增大，推進劑表面被迅速加熱到分解溫度以上。在放電電流達到脈沖峰值時，推進劑表面溫度在高熱流密度作用下繼續(xù)升高，在放電電流達到峰值后約0．5 μs時，推進劑表面溫度達到最大值。之后，隨著電流振蕩衰減而緩慢震蕩回落。PPT工作時電極附近的電流密度最大，放電對等離子體的歐姆加熱作用最強，傳導到推進劑表面的熱流密度也就最大，這使得靠近電極位置的推進劑表面溫度明顯高于其余位置，這也意味著這個位置的推進劑燒蝕速率更大，與實驗結果相吻合。此外，由于電極附近等離子體中的電流密度相對更接近于放電回路中的電流密度，靠近電極位置的推進劑表面溫度變化具有更明顯的放電電流波動跟隨性。

圖6 放電電流和特氟隆表面溫度Fig．6 Current waveform and surface temperature

圖7給出了特氟隆燒蝕速率和燒蝕質量隨時間的變化關系，對比圖6可看出，在PPT放電的前1/4個振蕩周期內，放電電流從0增大到最大值，而燒蝕速率和燒蝕質量卻不大。在PPT放電的第二個1/4振蕩周期內，推進劑燒蝕速率在放電電流達到峰值后，約0．5 μs時達到最大值，與此時的推進劑表面溫度達到最高相對應。在這個1/4振蕩周期內燒蝕的推進劑質量快速增加，占總燒蝕質量的比重超過95%。之后，隨著推進劑表面溫度下降，燒蝕速率迅速減小，燒蝕質量增加相當緩慢。

作為一種電磁推力器，PPT產生推力主要依靠等離子體受到的洛倫茲力J×B的加速作用，大的放電電流意味著高電流密度和強感應磁場，也就對應著大的加速力。PPT中推進劑燒蝕速率峰值時刻滯后于放電電流峰值時刻，大部分的燒蝕質量產生于放電電流達到最大值后不斷減小的時間段內，這既不利于推進劑燒蝕產物的電離，也不利于等離子體的加速。這種推進劑燒蝕過程和推力器放電過程之間的匹配失衡，會極大地降低推進劑的使用效率和推力器的能量利用效率，導致推力效率低下。

圖8給出了電極附近的推進劑在多個時刻點時的內部溫度分布，圖9為推進劑燒蝕過程中無定形體區(qū)域特氟隆的體積隨時間的變化關系。

圖7 特氟隆燒蝕速率和燒蝕質量Fig．7 Ablation rate and mass of Teflon

圖8 y=0．05 h處的Teflon內部溫度分布Fig．8 Temperature distribution in Teflon at y=0．05 h

圖9 無定形體特氟隆體積Fig．9 Volume of Teflon at amorphous phase

在放電開始的最初幾微秒內，推進劑表面溫度迅速上升，隨著大量推進劑相變之后開始分解，燒蝕質量迅速增加。在PPT放電半個振蕩周期之后，盡管有熱量持續(xù)不斷地傳遞給推進劑，使得發(fā)生相變的推進劑不斷增多，推進劑遠端溫度不斷升高，但由于熱流密度減小，推進劑表面溫度逐漸降低，推進劑燒蝕速率迅速下降到極低的水平，使得這一部分傳入能量并沒有被有效利用。此外，放電電流反向以后的電磁加速作用有所減弱，甚至計算發(fā)現部分區(qū)域的洛倫茲力會使等離子體運動減速，導致前期燒蝕產物長時間滯留于放電通道內?？梢?，PPT放電后期的能量利用效率較低，減小電流振蕩，有助于提高推力效率。

Spanjers等［16］在實驗中觀察到在PPT放電結束后有大量特征尺度在1～100 μm范圍的特氟隆微粒以低速噴出，這種粒子發(fā)射現象損耗了大量推進劑，由此產生的推力卻可忽略。粒子發(fā)射現象的產生原因可能是由于推進劑表面附近具有很大的溫度梯度及超過特氟隆相變溫度的高溫，使得處于表面附近的推進劑在因相變分解導致結構強度降低的情況下受熱應力和表面之后的推進劑受熱分解產生氣體的擠壓作用而被撕裂和沖刷帶走。Spanjers等人是在放電開始100 μs后拍攝到粒子發(fā)射圖像的，而在放電的前10 μs并沒有找到這一現象存在的證據，在圖9中成為無定形體的特氟隆體積在放電過程中不斷增大，這意味著隨后可能會有不少推進劑因粒子發(fā)射而被耗費掉。發(fā)生相變的推進劑體積由熱流密度決定，在放電后期迅速關斷放電電流和對推進劑進行冷卻，可減緩無定形體特氟隆的體積增長，甚至使部分無定形體特氟隆再結晶成為晶狀相。因此，這有助于減小因粒子發(fā)射帶來的推進劑損失，從而提高推力效率。

5 結論

(1)靠近電極位置的推進劑表面溫度更高，燒蝕速率更大。

(2)推進劑燒蝕速率峰值時刻滯后于放電電流峰值時刻，大部分的燒蝕質量產生于PPT放電的第2個1/4振蕩周期內，燒蝕過程和放電過程之間的匹配失衡降低了推力效率。

(3)PPT放電后期的能量利用效率較低，減小電流振蕩，有助于提高推力效率。

(4)在放電后期，迅速關斷放電電流和對推進劑進行冷卻，有助于減小因粒子發(fā)射帶來的推進劑損失，從而提高推力效率。

［1］ Burton R L．Pulsed plasma thruster［J］．Journal of Propulsion and Power，1998，14(5)．

［2］ Spanjers G G，Mcfall K A，Gulczinski III F S，et al．Investigation of propellant inefficiencies in a pulsed plasma thruster［R］．AIAA 96-2723．

［3］ Brito C M，Elaskar S A，Brito H H，et al．Zero-dimensional model for preliminary design of ablative pulsed plasma Teflon thrusters［J］．Journal of Propulsion and Power，2004，20(6)．

［4］ Turchi P J，Mikellides I G，Mikellides P G，et al．Optimization of pulsed plasma thrusters for microsatellite propulsion［R］．AIAA 99-2301．

［5］ Thomas H D．Numerical simulation of pulsed plasma thrusters［D］．The University of Tennessee，2000．

［6］ Henrikson E M．An experimental and theoretical study towards performance improvements of the ablation fed pulsed plasma thruster［D］．Arizona State University，2010．

［7］ Mikellides Y G．Theoretical modeling and optimization of ablation-fed pulsed plasma thrusters［D］．Ohio State University，1999．

［8］ Lin G，Karniadakis G E．High-order modeling of micro-pulsed plasma thrusters［R］．AIAA 2002-2872．

［9］ Antonsen E L，Burton R L，Reed G A，et al．Effects of postpulse surface temperature on micro-pulsed plasma thruster operation［R］．AIAA 2004-3462．

［10］ Keidar M，Fan J，Boyd I D．Vaporization of heated materials into discharge plasmas［J］．Journal of Applied Physics，2001，89(6)．

［11］ Clark B L．An experimental and analytical investigation of Teflon abltion heat transfer parameters by the method of non-linear estimation［D］．Cornell University，1971．

［12］ Holzknecht B．An analytical model of the transient ablation of polytetrafluoroethylene layers［J］．International Journal of Heat and Mass Transfer，1977，20．

［13］ Arai N．Transient ablation of Teflon in intense radiative and convective environments［J］．AIAA Journal，1979，17(6)．

［14］ Gatsonis N A，Juric D，Stechmann D P，et al．Numerical analysis of Teflon ablation in pulsed plasma thrusters［R］．AIAA 2007-5227．

［15］ Kemp N H．Surface recession rate of an ablating polymer．［J］．AIAA Journal，1968，6(9)．

［16］ Spanjers G G，Lotspeich J S，McFall K A，et al．Propellant losses because of particulate emission in a pulsed plasma phruster［J］．Journal of Propulsion and Power，1998，14(4)．