推力器真空羽流熱效應(yīng)計(jì)算模型修正及 誤差分析

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

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

0 引言

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

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

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

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

1 推力器羽流流場(chǎng)計(jì)算方法

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

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

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

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

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

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

圖1 點(diǎn)源法計(jì)算流程Fig.1 The flow chart of source flow method

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

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

2.1 Knudsen 模型分析

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

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

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

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

表1 漫反射系數(shù)α 變化對(duì)羽流效應(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 時(shí)，α隨表面溫度T的變化如圖3所示。由圖可見，α的變化范圍是比較寬的，對(duì)于分子量比較小的氣體（H2），α的數(shù)值分布總體上＜0.4，并且隨溫度的升高而減??；對(duì)于分子量比較大的氣體（H2O、N2、CO2、CO），α的數(shù)值分布總體上＞0.7。而一般情況下，羽流與航天器表面發(fā)生撞擊時(shí)，存在能量交換的漫反射的發(fā)生占比最大，少量分子會(huì)發(fā)生鏡面反射和吸附在航天器表面。因此，在以往航天器羽流熱效應(yīng)分析中，常選用完全漫反射模型進(jìn)行羽流熱效應(yīng)計(jì)算。

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

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

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

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

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

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

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

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

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

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

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

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

3.2 模型修正

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

（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]張建華, 蔡國(guó)飆.用Simons 法計(jì)算真空羽流[J].推進(jìn)技術(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]沈青.稀薄氣體動(dòng)力學(xué)[M].北京∶國(guó)防工業(yè)出版社, 2003∶228-230

[7]李志輝, 李中華, 楊東升, 等.衛(wèi)星姿控發(fā)動(dòng)機(jī)混合物羽流場(chǎng)分區(qū)耦合計(jì)算研究[J].空氣動(dòng)力學(xué)學(xué)報(bào), 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é)報(bào), 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ì)算[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 推力器羽流試驗(yàn)數(shù)據(jù)建立羽流場(chǎng)數(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