110 kN 核熱發(fā)動機(jī)推力室非平衡流動傳熱數(shù)值模擬研究

李子亮，徐 凱

(北京航天動力研究所，北京100076)

1 引言

隨著深空探測技術(shù)的不斷發(fā)展，核熱推進(jìn)的優(yōu)勢逐漸顯現(xiàn)，與化學(xué)推進(jìn)相比，核熱推進(jìn)具有比沖更高、工作時間更長且不依賴太陽能的優(yōu)點[1-2]。 核熱推進(jìn)的核心是核熱火箭發(fā)動機(jī)，早在上世紀(jì)五十年代中期，國外就已開始核熱火箭發(fā)動機(jī)的研制工作，具有代表性的是美國930 kN 推力的NERVA 發(fā)動機(jī)和蘇聯(lián)35 kN 推力的RD—0410 發(fā)動機(jī)[3-4]。 我國核熱推進(jìn)技術(shù)研究起步晚，基礎(chǔ)較為薄弱，但反應(yīng)堆工程技術(shù)基礎(chǔ)比較雄厚[3]。 核熱發(fā)動機(jī)發(fā)展的主要趨勢是模塊化和小型化，目前美俄重點發(fā)展30 ～100 kN、比沖900 s以上的核熱火箭發(fā)動機(jī)，將多個核熱火箭發(fā)動機(jī)捆綁使用，以滿足不同任務(wù)需求[5-7]。

核熱火箭發(fā)動機(jī)中的氫通過反應(yīng)堆芯被加熱至高溫并經(jīng)噴管膨脹做功產(chǎn)生推力，國外核熱火箭發(fā)動機(jī)研制中普遍將核反應(yīng)堆的設(shè)計作為重點，而對堆芯下游高溫氫在推力室的流動傳熱過程關(guān)注較少[8]。 對于以氫為推進(jìn)劑的核熱火箭發(fā)動機(jī)，氫在流經(jīng)推力室收縮擴(kuò)張噴管時，氣體的壓力和溫度不斷變化，并伴隨有限速率的離解反應(yīng)，因此核熱火箭發(fā)動機(jī)推力室流動傳熱的計算需要考慮氫氣非平衡離解的影響。 非平衡流動廣泛用于化學(xué)推進(jìn)火箭發(fā)動機(jī)推力室流場的計算，目前核熱發(fā)動機(jī)內(nèi)氫的流動研究主要集中在低壓高溫非平衡流動二維計算和高溫高壓平衡流動的一維計算[9-10]，而對推力室內(nèi)高溫高壓氫非平衡流動對發(fā)動機(jī)性能參數(shù)的影響研究不充分。 鑒于此，本文以110 kN 核熱火箭發(fā)動機(jī)推力室為研究對象，通過數(shù)值模擬方法對推力室的流動傳熱進(jìn)行建模仿真，研究堆芯下游推力室內(nèi)高溫高壓氫的非平衡離解作用對發(fā)動機(jī)性能和再生冷卻傳熱的影響。

2 模型與計算方法

2.1 幾何建模與網(wǎng)格劃分

根據(jù)110 kN 核熱火箭發(fā)動機(jī)設(shè)計參數(shù)建立推力室身部和噴管延長段三維幾何模型，其中堆芯出口直徑為416 mm，噴管面積比300，身部采用銑槽式再生冷卻通道且數(shù)目可變，具體參數(shù)如表1 所示。

表1 推力室設(shè)計參數(shù)Table 1 Design parameters of thrust chamber

為節(jié)省計算成本，根據(jù)模型的周向?qū)ΨQ性，選擇夾角為1.2°的單個通道的推力室模型進(jìn)行計算。 采用ICEM 軟件對模型進(jìn)行網(wǎng)格劃分，計算域網(wǎng)格總數(shù)為654.974 萬，其中以六面體網(wǎng)格為主，對與流體接觸的壁面附近網(wǎng)格進(jìn)行加密，保證壁面y＋≈1，如圖1 所示。

圖1 推力室?guī)缀文Ｐ团c網(wǎng)格劃分Fig.1 Geometric model and meshing of thrust chamber

2.2 控制方程

氣體在推力室的流動過程滿足粘性N-S 方程，遵循質(zhì)量守恒、動量守恒和能量守恒，氣體控制方程如式(1)～(3)所示。

式中，μ 為分子粘度，qi為熱流量，σij為粘性張量，cp為比熱容，ρ 為密度，T 為溫度，p 為壓力。

固體傳熱控制方程如式(4)所示。

式中，Q 為熱源，k 為導(dǎo)熱系數(shù)。

2.3 湍流模型

SST k-ω 湍流模型是由k-ε 模型和k-ω 模型組成的復(fù)合模型，并充分集成了二者的優(yōu)點。SST k-ω 湍流模型用湍流計算時，近壁面低雷諾數(shù)區(qū)采用k-ω 模型，遠(yuǎn)壁面完全湍流區(qū)采用k-ε模型，在計算受限空間流場時，SST k-ω 模型較其他雙方程模型具有顯著的優(yōu)越性，因此，采用SST k-ω 模型對推力室流場進(jìn)行計算。

SST k-ω 湍流模型的輸運(yùn)方程如式(5)～(9)所示。

SST k-ω 模型湍流粘度表達(dá)式如式(10) ～(12)所示。

式中，τij為雷諾應(yīng)力張量，Ω 為渦度張量，μ為層流粘度，μt為湍流粘度，y 為到壁面的距離，模型常數(shù)a1＝0.31， β*＝0.09， σk1＝0.85， σω1＝0.5， σk2＝1.0， σω2＝0.856。

2.4 化學(xué)反應(yīng)模型

采用有限速率反應(yīng)機(jī)理計算氫在推力室噴管內(nèi)的非平衡流動。 氫離解的可逆反應(yīng)式為式(13)。

式中，M 代表H 或H2，ΔE 代表反應(yīng)吸收或釋放的能量。 氫的離解反應(yīng)采用有限速率/渦耗散模型計算。

2.5 介質(zhì)物性

經(jīng)堆芯換熱后的高溫氫氣滿足理想氣體狀態(tài)方程，比熱為溫度的多項式函數(shù)，粘度滿足蘇薩蘭公式，導(dǎo)熱率滿足動力學(xué)理論。 冷卻通道內(nèi)超臨界氫的密度、比熱、導(dǎo)熱率和粘度隨壓力和溫度變化，擬合為壓力和溫度的多項式，再生冷卻傳熱內(nèi)壁面與肋的材質(zhì)為鋯銅，外壁材質(zhì)為電鑄鎳，將其比熱、導(dǎo)熱率擬合為溫度的函數(shù)。

2.6 邊界條件

將計算域邊界設(shè)定為以下4 類邊界條件:

1)入口邊界條件:高溫氫氣的入口為壓力入口邊界，根據(jù)已知堆芯出口參數(shù)給定氫的入口壓力、溫度和組分含量；超臨界氫的入口為質(zhì)量流量入口，給定入口質(zhì)量流量、壓力和溫度，并給定相應(yīng)湍流的湍動強(qiáng)度和水力直徑。

2)出口邊界條件:氫氣和超臨界氫的出口均為壓力出口，給定出口壓力、溫度和相應(yīng)的湍動強(qiáng)度與水力直徑。

3)壁面邊界條件:內(nèi)壁面為流固耦合壁面，給定粗糙度為4 μm，外壁面為非耦合壁面，定義其為光滑絕熱壁面。

4)對稱邊界條件:計算域中過軸線的燃?xì)狻⑹冶诤屠鋮s通道兩側(cè)面為對稱面。 無離解態(tài)是指氫分子不發(fā)生離解的狀態(tài)；平衡離解凍結(jié)態(tài)是指在推力室入口氫氣離解反應(yīng)已達(dá)到平衡，噴管內(nèi)的流動過程中氣體組分不再發(fā)生變化；無熱沉積非平衡離解態(tài)是指在不考慮熱沉積的條件下氫氣在噴管內(nèi)非平衡流動的狀態(tài)；熱沉積下非平衡離解態(tài)是指考慮熱沉積的條件下氫氣在噴管內(nèi)非平衡流動的狀態(tài)。 在計算中通過定義入口初始組分和有限化學(xué)反應(yīng)機(jī)理進(jìn)行計算。 在計算熱沉積條件下的非平衡流動時，將5 MPa 和2750 K 下離解平衡態(tài)的組分作為推力室入口的初始組分，假設(shè)反應(yīng)堆中子在推力室的熱沉積能量均勻分布，在計算中作為體熱源。 推力室計算域主要參數(shù)如表2 所示。

表2 計算域主要參數(shù)Table 2 Parameters of thrust chamber

2.7 計算求解

將流固界面作為耦合內(nèi)壁面，對燃?xì)夂屠鋮s劑的流動與壁面的傳熱進(jìn)行同步耦合求解，采用DO 模型計算輻射傳熱過程；流體項采用二階迎風(fēng)格式進(jìn)行離散，采用SIMPLE 算法處理壓力和速度的關(guān)系，采用24 個CUP 進(jìn)行并行計算求解，能量方程解的收斂精度為10-6，連續(xù)性和動量方程收斂精度為10-3。

3 結(jié)果與討論

高溫氫在推力室噴管發(fā)生非平衡離解反應(yīng)時計算域主要參數(shù)分布如圖2 所示。 可以看出推力室冷卻通道內(nèi)流體域的壓力最大且溫度最低；高溫高壓氫經(jīng)噴管膨脹做功，軸向壓力和溫度不斷降低，速度不斷升高，并在出口達(dá)到最大值。

圖2 非平衡離解態(tài)下推力室計算域內(nèi)主要參數(shù)分布Fig.2 Main parameters distribution in the calculation domain of thrust chamber under non-equilibrium dissociation

圖3 和圖4 分別為非平衡離解態(tài)的氫分子和氫原子在推力室噴管內(nèi)的質(zhì)量分?jǐn)?shù)分布。 高溫高壓氫在收縮擴(kuò)張過程中，氫原子不斷重組為氫分子，氫分子的質(zhì)量分?jǐn)?shù)逐漸增加，其中在噴管收斂段氣體組分質(zhì)量分?jǐn)?shù)變化最大。 在相同的軸向位置，由于推力室身部再生冷卻的作用，身部壁面附近氣體溫度較低，有利于氫原子的重組，從而使壁面附近的氫分子質(zhì)量分?jǐn)?shù)高于壁面遠(yuǎn)端流場中的質(zhì)量分?jǐn)?shù)，如圖3 所示。

表3 為氫在無離解態(tài)、平衡離解凍結(jié)態(tài)、無熱沉積非平衡離解態(tài)和熱沉積下非平衡離解態(tài)的推力室再生冷卻壓降、溫升和比沖計算結(jié)果，從表中可以看出，在4 種狀態(tài)下，再生冷卻氫的壓降均為0.47 MPa，而溫升和比沖存在差異。 非平衡離解態(tài)下的溫升較平衡離解凍結(jié)態(tài)和無離解態(tài)稍高，這主要是因為噴管內(nèi)非平衡態(tài)的氫隨著軸向溫度和壓力的降低，壁面附近已離解的氫原子重新組成氫分子，并放出能量，對再生冷卻傳熱起到促進(jìn)作用。 無離解態(tài)和平衡離解凍結(jié)態(tài)下氫的溫升接近，平衡離解凍結(jié)態(tài)的比沖要高于無離解態(tài)，這是因為平衡離解凍結(jié)態(tài)的氫中含有一定量的氫原子，氣體密度更低，相同條件下其比沖更高；非平衡離解態(tài)的溫升和比沖要高于平衡離解凍結(jié)態(tài)和無離解態(tài)，這是因為非平衡離解態(tài)的氫原子在噴管收縮擴(kuò)張段重組為氫分子(圖3 和圖4)，該過程放出的熱量一部分使壁面?zhèn)鳠嵩鰪?qiáng)，溫升變大，另一部分熱量轉(zhuǎn)化為動能使噴管出口比沖增加。由于推力室上的中子熱沉積能量很小，使得無熱沉積非平衡離解態(tài)和熱沉積下非平衡離解態(tài)的溫升與比沖十分接近，因此，熱沉積對噴管比沖和再生冷卻傳熱的影響可以忽略。

圖3 非平衡離解態(tài)下的氫氣質(zhì)量分?jǐn)?shù)分布Fig.3 Mass fraction distribution of hydrogen under non-equilibrium dissociation

圖4 非平衡離解態(tài)下的氫原子質(zhì)量分?jǐn)?shù)分布Fig.4 Mass fraction distribution of hydrogen atom under non-equilibrium dissociation

表3 不同狀態(tài)下推力室主要技術(shù)指標(biāo)數(shù)值計算結(jié)果Table 3 Calculation results of main technical indicators of thrust chamber under different conditions

圖5 和圖6 為不同離解態(tài)下推力室身部燃?xì)獗诿鏌崃髅芏群蜏囟鹊妮S向分布曲線。 從圖中可以看出，熱流密度與壁面溫度的峰值(57 MW/m2與522 K)出現(xiàn)在喉部(x ＝0)上游附近，且燃?xì)獗诿鏈囟染陀诓牧系臏囟壬舷蕖錇榉瞧胶怆x解態(tài)時，喉部壁面熱流與溫度均高于無離解態(tài)與平衡離解凍結(jié)態(tài)，此外，當(dāng)氫為無離解態(tài)和平衡離解凍結(jié)態(tài)時，壁面熱流密度與溫度十分接近。

圖5 不同離解態(tài)下燃?xì)獗诿鏌崃髅芏容S向分布曲線Fig.5 Axial distribution curve of heat flow density on thrust gas wall under different dissociation state

圖6 不同離解態(tài)下推力室燃?xì)獗诿鏈囟容S向分布曲線Fig.6 Axial distribution curve of temperature on thrust gas wall under different dissociation state

4 結(jié)論

1)110 kN 核熱發(fā)動機(jī)推力室流動傳熱仿真中應(yīng)考慮非平衡有限速率離解反應(yīng)對再生冷卻傳熱和發(fā)動機(jī)比沖的影響，該計算方法更加合理。

2)高溫離解的氫原子在噴管內(nèi)會發(fā)生重組并放出熱量，非平衡離解態(tài)下計算得到的溫升、比沖、喉部壁面熱流密度與溫度均高于無離解態(tài)和平衡離解凍結(jié)態(tài)的相應(yīng)值。

3)反應(yīng)堆中子在推力室熱沉積對再生冷卻傳熱影響很小，在計算中可忽略不計。