趙芃沛，朱 鐳

(陜西空天動力研究院有限公司，西安 710103)

火箭發(fā)動機(jī)噴管羽流特征對于發(fā)動機(jī)推力大小、紅外輻射特性[1-3]以及噪聲特性都有直接影響[4-6]。地面條件下羽流流動可以用連續(xù)介質(zhì)模型來描述，用速度、密度、壓力、溫度以及位置和時間等宏觀量建立流場信息的聯(lián)系。真空羽流場與地面羽流場差別較大。燃?xì)鈴膰姽車姵龊?，在近似真空的環(huán)境中無干擾地迅速膨脹，呈球狀擴(kuò)散，擴(kuò)散區(qū)域大大增加，密度急劇下降，流動狀態(tài)也經(jīng)歷了從連續(xù)流區(qū)至過渡流區(qū)，最終到達(dá)自由分子流區(qū)的過程。此過程將經(jīng)歷連續(xù)流、過渡流、稀薄流等流動狀態(tài)，并伴隨非平衡化學(xué)反應(yīng)和多種激波產(chǎn)生等，使羽流流動狀態(tài)的描述以及相關(guān)問題的研究變得極為復(fù)雜[7-9]。

目前關(guān)于真空條件下噴管羽流主要采用DSMC[7]方法和DSMC/CFD混合方法，但目前并沒有成熟的商業(yè)軟件采用這兩種方法來計算羽流，考慮到本文最終關(guān)注的是羽流外邊界輻射特性對部件的影響，考慮時間和經(jīng)費(fèi)成本投入，論文采用了CFD方法來進(jìn)行計算，可以快速得到初步特性來進(jìn)行定性分析，在一定范圍內(nèi)可作為設(shè)計的初步參考。

1 計算模型與方法

1.1 計算方法和校驗

本文采用的控制方程為雷諾平均N-S方程(RANS)，該方程的矢量形式為：

流體粘性系數(shù)由Sutherland公式確定：

式中:μ0=1.716×10-5NS/m2，T0=121.56 K,μt為湍流粘性系數(shù)，由湍流計算模型給出。

湍流模型采用Menter的SSTk-ω模型，其流動方程如下：

式中：Gk表示湍流的動能，Gω為ω方程，Γk和Γω分別代表k與ω的有效擴(kuò)散項,Yk和Yω分別代表k與ω的發(fā)散項。Dω代表正交發(fā)散項。Sk和Sω用戶自定義。

采用標(biāo)準(zhǔn)壁面函數(shù)法，標(biāo)準(zhǔn)壁面函數(shù)法假定：壁面附近粘性底層以外的地區(qū)，量綱速度服從對數(shù)率分布。

流場求解利用有限體積法離散，定常雷諾平均N-S方程，各方程聯(lián)立耦合隱式求解，對流項為二階迎風(fēng)格式，粘性項采用二階中心差分格式，其他為一階迎風(fēng)格式。計算結(jié)束條件為所有殘差下降3個數(shù)量級且不再變化、關(guān)鍵截面流量不再變化表示計算結(jié)果收斂。

采用以上模型針對文獻(xiàn)[10]中數(shù)據(jù)進(jìn)行了校驗計算。該文獻(xiàn)采用DSMC方法對真空羽流進(jìn)行了數(shù)值模擬。將該方法給出的數(shù)據(jù)與本文采用CFD方法的結(jié)果進(jìn)行了對比，如圖1。

可以看出，在噴管內(nèi)DSMC方法與CFD方法吻合度較好，因為這一區(qū)域滿足連續(xù)性假設(shè)；而在距離噴管出口較遠(yuǎn)的位置，即過渡流區(qū)以外以及倒流區(qū)的區(qū)域有一定的差別。

1.2 幾何模型和計算區(qū)域

計算區(qū)域和噴管幾何構(gòu)型如圖2所示。左下角為噴管所在位置。在近似真空的環(huán)境條件下，壁面附近的氣流將繞過唇口向管壁上游流動，產(chǎn)生倒流現(xiàn)象，因此將計算域水平方向和豎直方向均設(shè)定為2 000 mm。對壁面進(jìn)行加密，并且對噴管出口處對流向方向進(jìn)行適當(dāng)加密。由于結(jié)構(gòu)網(wǎng)格局部加密會導(dǎo)致整體加密的特點(diǎn)，故對遠(yuǎn)場的節(jié)點(diǎn)分布律進(jìn)行調(diào)整，改為均勻分布，這樣既保證噴管處網(wǎng)格質(zhì)量，又保證遠(yuǎn)場網(wǎng)格均勻分布?？偩W(wǎng)格數(shù)量11萬。

圖1 Ma數(shù)等值線圖

其中發(fā)動機(jī)噴管出口直徑36.3 mm，喉部直徑4.525 mm，出口直徑36.3 mm，幾何尺寸如圖2所示。室壓1.8 MPa，工作在120 km高度。由于120 km高度處已接近真空狀態(tài)，分子平均自由程已經(jīng)接近或者很有可能超過噴管的幾何尺寸，因此要特別注意非連續(xù)效應(yīng)的影響。

圖2 計算區(qū)域和噴管幾何構(gòu)型

由于該噴管為軸對稱噴管，主要研究羽流狀態(tài)，并且計算中并未涉及到其他部件，因此整個計算域可以采用軸對稱構(gòu)型以減小計算量。

其中燃燒室采用質(zhì)量流率入口，給定流率、總溫、總壓、組份；遠(yuǎn)場為壓力出口，給定壓力和組份；壁面均為無滑移壁面。

1.3 計算工況

表1給出了計算工況下推進(jìn)劑的流量，燃料為CH6N2；氧化劑為N2O4。

在室壓為1.8 MPa條件下，根據(jù)室壓及燃料混合比，可以計算燃溫、燃燒產(chǎn)物等。采用平衡流假設(shè)，經(jīng)過熱力學(xué)計算軟件CEA計算得到：燃燒室總溫3 103.7 K。以此作為計算的初始條件之一。

工作高度120 km情況下，經(jīng)過查表得到：大氣壓力：0.002 54 Pa；溫度：360 K。

由于帶化學(xué)反應(yīng)的計算需要耗費(fèi)大量的計算資源，為了簡化計算，本次計算直接給定燃燒產(chǎn)物，而不考慮反應(yīng)過程。

經(jīng)過CEA軟件計算，反應(yīng)主要涉及的產(chǎn)物有11種，分別為CO、CO2、H、HO2、H2、H2O、NO、N2、O、OH、O2。出于簡化計算的考慮，選取質(zhì)量分?jǐn)?shù)大于1%的組份，并對這些組份濃度進(jìn)行歸一化處理，如表2。

表2 燃燒室組份濃度

3 計算結(jié)果

經(jīng)過迭代計算，各項殘差變化趨于平穩(wěn)，可以認(rèn)為計算結(jié)果收斂。若僅從圖3所示的求解結(jié)果看，可以發(fā)現(xiàn)，氣體經(jīng)噴管加速以后流出，在低壓條件下迅速擴(kuò)張，達(dá)到了很高的Ma數(shù)。

圖3 Ma數(shù)等值線圖

由圖4、圖5可以看到，相對于常壓條件，真空條件下尾流擴(kuò)張角度更大，呈現(xiàn)出弓形擴(kuò)張區(qū)域。同時，由于真空條件下極低的壓力，導(dǎo)致了垂直于流向方向上較大的擴(kuò)散?？梢钥吹?，整個計算域的總溫均較高，這與常壓條件下完全不同。

圖4 整體流線圖

圖5 總溫分布等值線圖

圖6、圖7是噴管局部的Ma數(shù)、溫度分布云圖?？梢钥闯鲈趪姽艹隹谛纬奢^弱的斜激波，但波系結(jié)構(gòu)和常壓條件下相比較少。

圖6 噴管附近Ma數(shù)云圖

圖7 噴管附近溫度云圖

由于噴管附近的區(qū)域(不包括水平壁面上方的倒流區(qū))分子數(shù)目較多，符合連續(xù)性假設(shè)，因此可以認(rèn)為對該區(qū)域的求解結(jié)果真實度較高。而過渡流區(qū)以外的區(qū)域需要進(jìn)一步探討。

4 結(jié)論

運(yùn)用CFD方法來進(jìn)行計算真空條件下拉伐爾噴管產(chǎn)生的羽流外邊界輻射特性對部件影響，可以得出真空條件下尾流擴(kuò)張角度更大，呈現(xiàn)出弓形擴(kuò)張區(qū)域。同時，由于真空條件下極低的壓力，導(dǎo)致了垂直于流向方向上較大的擴(kuò)散。整個計算域的總溫均較高，這與常壓條件下完全不同。在噴管出口形成較弱的斜激波，波系結(jié)構(gòu)和常壓條件下相比較少。

由于噴管附近的區(qū)域(不包括水平壁面上方的倒流區(qū))分子數(shù)目較多，符合連續(xù)性假設(shè)，可以認(rèn)為對該區(qū)域的求解結(jié)果真實度較高。