推力器真空羽流熱效應計算模型修正及 誤差分析

王黎珍，史紀鑫，鄭世貴

（北京空間飛行器總體設計部，北京 100094）

0 引言

在航天器進行變軌及姿態控制調整時，推力器工作可能產生羽流擾動力/力矩、羽流熱和污染效應。隨著航天器安裝的推力器種類和數量越來越多，與羽流相關的各種問題日趨凸顯，已引起了設計部門高度重視，開展了羽流效應尤其是熱效應的大量研究。對羽流熱效應的數值計算至關重要：若分析過于保守將會導致熱防護過設計；若估計不足則可能導致星上設備過熱。因此，深入開展推力器真空羽流熱效應精確評估是合理進行熱控設計的前提和條件。

推力器真空羽流熱效應的計算分兩大步驟：羽流流場的計算和羽流對航天器表面沖擊作用的計算分析。推力器真空羽流流動狀態極為復雜，可分為連續流（羽流核心區）、過渡流和自由分子流（羽流外圍區）等3 種流動狀態[1]。對于連續流區域，通常用CFD 方法數值求解N-S 方程以得到流場[2-3]。對于過渡流和自由分子流區域，采用半經驗的工程模型（點源模型）[4-5]或基于DSMC 算法[6-8]的數值仿真分析得到流場。目前，無論是采用CFD/DSMC的數值求解方法，還是采用半經驗的工程模型求解方法，羽流內外流場的求解精度均在工程可接受的范圍內。

在羽流對航天器表面沖擊作用的計算分析中，羽流沖擊航天器表面的特性模型的選擇是影響羽流熱效應計算結果的關鍵因素。模型給出氣體分子在固壁上的運動方式，反映氣體分子與固壁間相互作用時力和熱的傳遞過程。根據氣體與固壁間相互作用方式的不同，模型分為Knudsen模型、Nocilla 模型和Lord 模型[9]。Knudsen 模型是對經典的Maxwell 模型的改進，認為羽流分子撞擊航天器表面時存在吸附、漫反射和鏡面反射3種運動方式，并用相應的系數表征各運動方式的占比。Nocilla 模型和Lord 模型則以Knudsen 模型為基礎，對羽流分子反射方向和范圍做出一定的限制和假定。工程中一般推薦采用Knudsen 模型計算羽流熱流效應。

本文首先對推力器羽流內外流場進行計算；然后對Knudsen 模型中的關鍵參數進行討論，明確這些參數對羽流熱效應計算的敏感度，并用理論經驗公式對關鍵參數的取值范圍進行初步估算；最后利用MBB 公司10 N 推力器羽流熱流試驗數據，采用遍歷搜索的方法對Knudsen 模型中關鍵參數的取值進行精確評估，使得羽流熱流密度計算值與實驗值間的偏差減小到7%左右。整個修正計算大幅度地提高了羽流熱效應計算的精度和可信度，驗證了所采用的熱流計算模型的正確性和工程可用性。

1 推力器羽流流場計算方法

本文根據燃燒室內的總溫、總壓等燃燒條件，氣體熱物理性質及噴嘴的形狀計算推力器羽流核心區的氣體流場。對于推力器內部的流場，考慮黏性的作用，采用有限體積方法求解N-S 方程得到。對于推力器外部噴嘴附近區域（100 倍喉部半徑內）的流場，可以忽略黏性作用，采用數值求解Euler方程的方法得到。噴嘴邊緣處的流場可以由在噴嘴壁面處的普朗特-邁耶（Prandtl-Mayer）方程展開得到。

本文采用點源模型法對羽流的外流場進行模擬。點源模型法的特點是簡單、易理解、計算量少且又能夠滿足一定的工程精度要求。點源模型認為在噴管外流場中，由于氣流具有高馬赫數與低溫特性，即使氣體分子之間發生碰撞，但它們的相對速度很低，在垂直于羽流軸線方向上由于分子碰撞引起的熱散射很小，可以采用點源模型來描述。具體做法是在流場中定義一個凍結面，并在其上布置一定數目的自由分子點源，而流場中任意一點的流動參數可以視為所有自由分子點源產生的流場在該點的疊加。自由分子點源產生的流場密度ρ遵循輻射衰減律，即

式中：V為羽流擴展極限速度；s*為噴管喉部流率；R為距噴口出口處的距離；θ為偏離羽流軸線角度；θL為P-M 擴展角；f(θ)為偏離羽流軸線角度的函數。

由于流動的基本參量不同，f(θ)在羽流核心區和邊界層膨脹區的表達式不同。f(θ)在羽流核心區的表達式為

其中：γ為氣體的比熱比；θ∞是噴管流動極限偏轉角；θ0是羽流核心區流動極限偏轉角。在邊界層膨脹區，引入系數β，并認為f(θ)在該區是以指數形式衰減，其表達式為

流場中的其他物理量，如壓力P、溫度T和極限速度V可根據一維等熵關系求出。用點源模型計算得到羽流外流場的流程見圖1。

圖1 點源法計算流程Fig.1 The flow chart of source flow method

2 推力器羽流熱效應表面特性模型

羽流流場確定后，羽流沖擊模型的選擇就成為關鍵因素。本文首先對Knudsen 模型的特點及其關鍵參數對羽流熱流密度的敏感度進行分析，然后通過理論經驗公式對Knudsen 模型中關鍵參數的取值范圍進行分析和估算。

2.1 Knudsen 模型分析

圖2 Knudsen 模型中氣體和表面的相互作用Fig.2 The interaction between gas and surface for the Knudsen model

羽流沖擊到星體表面后會發生與固壁的相互作用及能量交換。Knudsen 模型認為羽流分子沖擊到衛星表面會產生被衛星表面吸附、鏡面反射和漫反射3 種運動方式（見圖2），它們的發生占比分別由吸附系數λ、鏡面反射系數τ和漫反射系數α確定。 漫反射系數α又稱壁面熱適應系數，對羽流在星體表面產生的熱流密度影響很大，是真空羽流熱效應數值模擬的關鍵參數。該參數表征了反射分子的溫度在多大程度上“適應”了星體表面的溫度狀況。當羽流分子在固壁上發生Maxwell 反射即完全漫反射時，α=1；完全鏡面反射時，α=0。Knudsen模型假設：若入射的氣體分子質量為m，氣體分子在與壁面作用的過程中，入射流量的λm被表面吸附，τ(1-λ)m被鏡面反射，(1-τ)(1-λ)m被漫反射，即系數α、τ和λ間的關系為α=(1-τ)(1-λ)。

2.2 漫反射系數α 對羽流熱流密度的敏感度分析

為了討論漫反射系數α對羽流熱效應的具體影響，本文進行了以下的數值計算假設：在距離推力器噴口100 cm 處放置接收板，假設氣體與接收板表面只發生吸附和漫反射2 種作用形式。通過計算α從0（氣體被星體表面完全吸附）變到1（氣體被星體表面完全漫反射）時羽流熱流密度的變化，得到α對羽流熱效應的影響關系。計算結果如 表1所示，可以看出，羽流熱流密度隨α的增加而增加。氣體分子被壁面完全吸附（α=0）時熱流密度最小，完全漫反射（α=1）時熱流密度達到最大。隨著α的變化，熱流密度在比較大的范圍內波動，完全漫反射時的羽流熱流密度值是完全吸附狀態的50 倍左右。由此可見，羽流熱流密度對α的變化高度敏感。

表1 漫反射系數α 變化對羽流效應的影響（τ=0）Table 1 The influence of accommodation ratio on the plume effect（τ=0）

2.3 漫反射系數α 取值范圍分析討論

通常情況下，漫反射系數的取值與燃燒產物氣體分子特性、材料表面特性和材料表面溫度相 關[10-11]。經過調研和資料收集，得到了H2、H2O、N2、CO2和CO 等5 種燃燒產物在表面材料為Al和Fe 時，α隨表面溫度T的變化如圖3所示。由圖可見，α的變化范圍是比較寬的，對于分子量比較小的氣體（H2），α的數值分布總體上＜0.4，并且隨溫度的升高而減?。粚τ诜肿恿勘容^大的氣體（H2O、N2、CO2、CO），α的數值分布總體上＞0.7。而一般情況下，羽流與航天器表面發生撞擊時，存在能量交換的漫反射的發生占比最大，少量分子會發生鏡面反射和吸附在航天器表面。因此，在以往航天器羽流熱效應分析中，常選用完全漫反射模型進行羽流熱效應計算。

圖3 常見氣體在常見材料表面的漫反射系數Fig.3 The accommodation ratio of the some gases and surface

Yasar Demirel 等[12]和周志雄等[13]通過對一系列氣體分子在不同溫度（273～1250 K）的金屬表面的漫反射系數進行分析，總結出如下關系：

其中：μ=Mg/Ms，Mg為氣體物質的分子量，Ms為固壁材料的分子量；F代表吸收層的覆蓋分數，由固體表面特性和表面溫度Ts決定。

對于航天器常用的單組元和雙組元推力器，其燃燒產物的摩爾分數和質量分數見表2。無論對于單組元還是雙組元推力器，其燃燒產物中大分子量氣體（NH3、H2O、N2）均約占總質量的99%左右。表2還給出了航天器表面溫度為273 K 時各燃燒產物的漫反射系數，并按照各燃燒產物的質量分數對漫反射系數進行加權平均，計算出單組元推力器和雙組元推力器的綜合漫反射系數分別為0.767 和0.782。

表2 推力器燃燒產物組成及漫反射系數估算Table 2 The evaluation of accommodation ratio for combustion products of thrusters

3 羽流熱效應表面特性模型修正

由前述可知，羽流熱流密度對漫反射系數α的變化非常敏感，而在以往的工程計算中，對α的取值比較保守（即將α取為1），這樣得到的羽流熱流值較大。因此在熱防護設計中會出現過防護，增大質量負荷。為此，本文利用MBB 公司的 10 N推力器羽流試驗數據重點開展對羽流表面特性模型的修正工作，將羽流熱流密度的計算誤差減小到工程允許的范圍內，同時驗證該模型的正確性和可信度。

3.1 MBB 10 N 推力器試驗條件及結果

MBB 公司使用標稱10 N 雙組元推力器[14-15]，在地面高真空試驗艙內進行了羽流熱效應測量試驗。真空試驗艙為直徑1 m、長約2 m 的圓柱容器。推力器安裝在圓柱容器中心線上，且二者的中心線重合。真空試驗艙采用液氦低溫系統和液氮低溫系統進行冷卻，最低溫度可以達到20 K；真空抽氣速度為0.014 m3/s，在推力器工作前艙內真空度可以達到1×10-4Pa。

羽流熱效應試驗重點研究了距噴口不同距離R，不同羽流角θ和不同羽流入射角β情況下的羽流熱流分布，其中θ和β的定義見圖4，試驗結果見表3。由于受真空艙直徑限制，測量時θ的最大有效值僅為40°。

圖4 推力器羽流角θ 和羽流入射角β 定義Fig.4 The definition of plume flow angle θ and the plume incidence angle β

表3 MBB 10 N 推力器不同工況下羽流熱試驗結果Table 3 The experimental plume heating results of MBB 10 N thruster

3.2 模型修正

在羽流沖擊航天器表面時，吸附現象一般發生在溫度低于200 K 的敏感器等儀器表面，而對于溫度較高的表面，很少發生吸附現象。因此，本文考慮的表面溫度較高，計算時設吸附系數λ為0。根據推力器羽流熱試驗工況的設置情況，可采用Knudsen 完全漫反射模型進行分析，計算結果見 表4?？梢钥闯觯嬎愕玫降挠鹆鳠崃髅芏戎灯毡楸仍囼炛荡?，相對誤差基本在50%以上，表明該模型夸大了羽流沖擊航天器表面產生的熱效應，將這些計算結果應用于工程設計顯然過于保守。

為此，我們采用遍歷搜索的算法，通過對Knudsen 模型中的漫反射系數α進行調整，將每次調整后的Knudsen 模型再進行5 個工況的計算分析；并將計算結果與試驗結果的相對誤差進行比對，得到5 組相對誤差的平均值相對于α的變化關系（如圖5所示）。由圖可以看出，隨著α從0 到1 的變化過程，熱流密度的相對誤差先減小、后增大。其中α在0.75～0.90 之間時，相對誤差降到20%以下。這個α的取值范圍與2.3 節利用理論經驗公式討論所得的結果一致。α為0.8 時計算結果的相對誤差最小，約為7.17%，以該α的取值作為模型修正的5 種工況下熱流密度計算結果及相對誤差見表4?？梢钥吹剑斯r2，其他4 個工況的相對誤差均在5%以下。

表4 不同模型計算結果與MBB 試驗結果對比Table 4 The comparision between MBB experimental results and simulation results for different models

圖5 計算結果與試驗結果相對誤差隨漫反射系數 α 的變化Fig.5 The variation of relative error between experimental and simulation results against accomodation ratio α

為了進一步驗證修正后Knudsen 模型的有效性，利用α=0.8 的修正模型分別計算了距10 N 推力器噴口的距離R=55 cm 處、在不同羽流角θ下的羽流熱流密度值。以θ=0°為基準，將計算結果和 試驗結果進行歸一化后再對比，結果見表5和圖6。

表5 數值仿真和MBB 10 N推力器試驗的歸一化結果比較Table 5 The comparision between MBB experimental normalized results and numerical simulation normalized results

可以看出，無論是計算結果還是試驗結果，在R不變的情況下，羽流熱流密度均隨著θ的增大而迅速減小，在趨勢上表現出了非常好的一致性。整體上分析結果與試驗結果吻合得較好，驗證了修正后的Knudsen 模型的有效性。

圖6 推力器不同羽流角時熱流密度計算結果與 試驗結果對比Fig.6 The comparision between numerical simulation results and the experimental results at different plume flow angles

4 結束語

本文分析了Knudsen 模型中關鍵參數對羽流熱效應計算結果的影響；并針對典型單組元和雙組元推力器燃燒產物，討論了模型中漫反射系數的取值范圍；最后利用MBB 10 N推力器地面試驗數據，采用遍歷搜索的算法，對Knudsen 模型進行了修正，使得修正后的羽流熱效應計算結果與試驗結果的相對誤差減小到約7%。相比保守的完全漫反射模型約50%的相對誤差，修正后的模型大幅度地提高了羽流熱效應計算的精度和可信度。

（References）

[1]Bird G A.Molecular gas dynamics and the direct simulation of gas flows[M].Oxford∶Clarendon Press, 1994∶2-4

[2]Gatsonis N A, Nanson R A, Lebeau G L.Simulations of cold-gas nozzle and plume flows and flight data comparison[J].Journal of Spacecraft and Rockets, 2000, 37(1)∶39-48

[3]Ebrahimi H B.Numerical investigation of multi-plume rocket phenomenology, AIAA 1997-2942[R], 1997

[4]張建華, 蔡國飆.用Simons 法計算真空羽流[J].推進技術, 2002, 23(5)∶406-409 Zhang Jianhua, Cai Guobiao.Computation based on the Simons model for vacuum plume[J].Journal of Propulsion Technology, 2002, 23(5)∶406-409

[5]焦子龍, 龐賀偉, 楊東升.利用Simons 模型研究衛星 羽流污染[J].航天器環境工程, 2007, 24(3)∶160-163 Jiao Zilong, Pang Hewei, Yang Dongsheng.Application of Simons model for satellite plume contamination analysis[J].Spacecraft Environment Engineering, 2007, 24(3)∶160-163

[6]沈青.稀薄氣體動力學[M].北京∶國防工業出版社, 2003∶228-230

[7]李志輝, 李中華, 楊東升, 等.衛星姿控發動機混合物羽流場分區耦合計算研究[J].空氣動力學學報, 2012, 30(4)∶483-491 Li Zhihui, Li Zhonghua, Yang Dongsheng, et al.Coupled simulation of mixture plume for attitude-control satellite thruster[J].Acta Aerodynamica Sinica, 2012, 30(4)∶483-491

[8]梁杰, 閻超, 楊彥廣, 等.過渡區側向噴流干擾的并行DSMC 數值模擬研究[J].宇航學報, 2011, 32(5)∶1012-1018 Liang Jie, Yan Chao, Yang Yanguang, et al.Parallel DSMC simulation of lateral jet interaction in rarefied transitional region[J].Journal of Astronautics, 2011, 32(5)∶1012-1018

[9]Andrew D K, Dean C W, Muntz E P.Influence of gas-surface interaction models on predicted performance of a micro-resisitoje, AIAA 2000-2430[R]

[10]Thomas L B, Olmer F G.The accommodation coefficient of He, Ne, A, H2, D2, O2, CO2and Hg on platinum as a function of temperature[J].Journal of the Americal Chemical Society, 1943, 65∶1036-1043

[11]Goodman F O, Wachman H Y.Formula for thermal accommodation coefficient, AD631007[R], 1966

[12]Demirel Y, Saxena S C.Heat transfer in rarefied gas at a gas-solid interface[J].Energy, 1996, 21(2)∶99-103

[13]周志雄, 魏蔚, 汪榮順.真空下氣-固界面熱適應系數的數值計算[J].低溫與超導, 2007, 35(1)∶36-40 Zhou Zhixiong, Wei Wei, Wang Rongshun.The mathematic calculation of thermal accommodation coefficients at gas-solid interface in vacuum[J].Cryogenics and Superconductivity, 2007, 35(1)∶36-40

[14]孫寶祥.利用10N 推力器羽流試驗數據建立羽流場數學模型[J].航天控制, 2005, 23(3)∶26-29 Sun Baoxiang.Modeling plume field using plume test data for 10N thruster[J].Aerospace Control, 2005, 23(3)∶26-29

[15]Trinks H.Evaluation of the exhaust plume of the bipropellant MBB 10N thruster for the INMARSAT2 satellite[C]//AIAA/ASME/SAE/ASEE 25thJoint Propulsion Conference.Monterey, CA, 1989.AIAA 89-2507