王黎珍，史紀(jì)鑫，鄭世貴

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

2 推力器羽流熱效應(yīng)表面特性模型

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

2.1 Knudsen 模型分析

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

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

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

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

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

2.3 漫反射系數(shù)α 取值范圍分析討論

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

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

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

其中：μ=Mg/Ms，Mg為氣體物質(zhì)的分子量，Ms為固壁材料的分子量；F代表吸收層的覆蓋分?jǐn)?shù)，由固體表面特性和表面溫度Ts決定。

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

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

3 羽流熱效應(yīng)表面特性模型修正

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

3.1 MBB 10 N 推力器試驗條件及結(jié)果

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

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

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

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

3.2 模型修正

在羽流沖擊航天器表面時，吸附現(xiàn)象一般發(fā)生在溫度低于200 K 的敏感器等儀器表面，而對于溫度較高的表面，很少發(fā)生吸附現(xiàn)象。因此，本文考慮的表面溫度較高，計算時設(shè)吸附系數(shù)λ為0。根據(jù)推力器羽流熱試驗工況的設(shè)置情況，可采用Knudsen 完全漫反射模型進行分析，計算結(jié)果見 表4?？梢钥闯?，計算得到的羽流熱流密度值普遍比試驗值大，相對誤差基本在50%以上，表明該模型夸大了羽流沖擊航天器表面產(chǎn)生的熱效應(yīng)，將這些計算結(jié)果應(yīng)用于工程設(shè)計顯然過于保守。

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

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

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

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

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

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

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

4 結(jié)束語

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

（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].推進技術(shù), 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 模型研究衛(wèi)星 羽流污染[J].航天器環(huán)境工程, 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]沈青.稀薄氣體動力學(xué)[M].北京∶國防工業(yè)出版社, 2003∶228-230

[7]李志輝, 李中華, 楊東升, 等.衛(wèi)星姿控發(fā)動機混合物羽流場分區(qū)耦合計算研究[J].空氣動力學(xué)學(xué)報, 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]梁杰, 閻超, 楊彥廣, 等.過渡區(qū)側(cè)向噴流干擾的并行DSMC 數(shù)值模擬研究[J].宇航學(xué)報, 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]周志雄, 魏蔚, 汪榮順.真空下氣-固界面熱適應(yīng)系數(shù)的數(shù)值計算[J].低溫與超導(dǎo), 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 推力器羽流試驗數(shù)據(jù)建立羽流場數(shù)學(xué)模型[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