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

李子亮，徐 凱

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

1 引言

隨著深空探測技術(shù)的不斷發(fā)展，核熱推進(jìn)的優(yōu)勢逐漸顯現(xiàn)，與化學(xué)推進(jìn)相比，核熱推進(jìn)具有比沖更高、工作時間更長且不依賴太陽能的優(yōu)點(diǎn)[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ā)展的主要趨勢是模塊化和小型化，目前美俄重點(diǎn)發(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è)計(jì)作為重點(diǎn)，而對堆芯下游高溫氫在推力室的流動傳熱過程關(guān)注較少[8]。 對于以氫為推進(jìn)劑的核熱火箭發(fā)動機(jī)，氫在流經(jīng)推力室收縮擴(kuò)張噴管時，氣體的壓力和溫度不斷變化，并伴隨有限速率的離解反應(yīng)，因此核熱火箭發(fā)動機(jī)推力室流動傳熱的計(jì)算需要考慮氫氣非平衡離解的影響。 非平衡流動廣泛用于化學(xué)推進(jìn)火箭發(fā)動機(jī)推力室流場的計(jì)算，目前核熱發(fā)動機(jī)內(nèi)氫的流動研究主要集中在低壓高溫非平衡流動二維計(jì)算和高溫高壓平衡流動的一維計(jì)算[9-10]，而對推力室內(nèi)高溫高壓氫非平衡流動對發(fā)動機(jī)性能參數(shù)的影響研究不充分。 鑒于此，本文以110 kN 核熱火箭發(fā)動機(jī)推力室為研究對象，通過數(shù)值模擬方法對推力室的流動傳熱進(jìn)行建模仿真，研究堆芯下游推力室內(nèi)高溫高壓氫的非平衡離解作用對發(fā)動機(jī)性能和再生冷卻傳熱的影響。

2 模型與計(jì)算方法

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

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

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

為節(jié)省計(jì)算成本，根據(jù)模型的周向?qū)ΨQ性，選擇夾角為1.2°的單個通道的推力室模型進(jìn)行計(jì)算。 采用ICEM 軟件對模型進(jìn)行網(wǎng)格劃分，計(jì)算域網(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)點(diǎn)。SST k-ω 湍流模型用湍流計(jì)算時，近壁面低雷諾數(shù)區(qū)采用k-ω 模型，遠(yuǎn)壁面完全湍流區(qū)采用k-ε模型，在計(jì)算受限空間流場時，SST k-ω 模型較其他雙方程模型具有顯著的優(yōu)越性，因此，采用SST k-ω 模型對推力室流場進(jìn)行計(jì)算。

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ī)理計(jì)算氫在推力室噴管內(nèi)的非平衡流動。 氫離解的可逆反應(yīng)式為式(13)。

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

2.5 介質(zhì)物性

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

2.6 邊界條件

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

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

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

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

4)對稱邊界條件:計(jì)算域中過軸線的燃?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)。 在計(jì)算中通過定義入口初始組分和有限化學(xué)反應(yīng)機(jī)理進(jìn)行計(jì)算。 在計(jì)算熱沉積條件下的非平衡流動時，將5 MPa 和2750 K 下離解平衡態(tài)的組分作為推力室入口的初始組分，假設(shè)反應(yīng)堆中子在推力室的熱沉積能量均勻分布，在計(jì)算中作為體熱源。 推力室計(jì)算域主要參數(shù)如表2 所示。

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

2.7 計(jì)算求解

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

3 結(jié)果與討論

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

圖2 非平衡離解態(tài)下推力室計(jì)算域內(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)的推力室再生冷卻壓降、溫升和比沖計(jì)算結(jié)果，從表中可以看出，在4 種狀態(tài)下，再生冷卻氫的壓降均為0.47 MPa，而溫升和比沖存在差異。 非平衡離解態(tài)下的溫升較平衡離解凍結(jié)態(tài)和無離解態(tài)稍高，這主要是因?yàn)閲姽軆?nèi)非平衡態(tài)的氫隨著軸向溫度和壓力的降低，壁面附近已離解的氫原子重新組成氫分子，并放出能量，對再生冷卻傳熱起到促進(jìn)作用。 無離解態(tài)和平衡離解凍結(jié)態(tài)下氫的溫升接近，平衡離解凍結(jié)態(tài)的比沖要高于無離解態(tài)，這是因?yàn)槠胶怆x解凍結(jié)態(tài)的氫中含有一定量的氫原子，氣體密度更低，相同條件下其比沖更高；非平衡離解態(tài)的溫升和比沖要高于平衡離解凍結(jié)態(tài)和無離解態(tài)，這是因?yàn)榉瞧胶怆x解態(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ù)值計(jì)算結(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ī)比沖的影響，該計(jì)算方法更加合理。

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

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