房玉良，秦 浩，王成龍，蘇光輝，田文喜，秋穗正

(西安交通大學(xué) 動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

核熱推進(jìn)技術(shù)是利用反應(yīng)堆產(chǎn)生熱能將推進(jìn)劑加熱到高溫狀態(tài)，高溫工質(zhì)從噴管膨脹產(chǎn)生高速氣流形成推力的新型推進(jìn)技術(shù)，具有比沖高、推力大等優(yōu)點(diǎn)[1-2]。核熱推進(jìn)系統(tǒng)采用氫氣作為推進(jìn)劑并將其加熱到2 500 K以上，產(chǎn)生與化學(xué)推進(jìn)相媲美的推力，比沖可達(dá)到化學(xué)推進(jìn)的兩倍。核熱推進(jìn)技術(shù)研究起源于20世紀(jì)50年代，經(jīng)過(guò)美、俄(蘇)等60多年的研究，固態(tài)堆芯核熱推進(jìn)技術(shù)逐漸成熟。美國(guó)宇航局在載人火星探測(cè)報(bào)告DRA 5.0中明確提出將核熱推進(jìn)火箭作為太空擺渡車的首選方案[3]。固態(tài)堆芯核熱推進(jìn)反應(yīng)堆冷卻通道可分為軸向型、徑向型，其中軸向型通道具結(jié)構(gòu)簡(jiǎn)單、易于加工制造等優(yōu)點(diǎn)使其應(yīng)用更為廣泛，如NERVA型堆芯、CERMET型堆芯等[4]。但由于軸向型通道具有很大的長(zhǎng)徑比，加上堆芯的熱流密度、溫度很高，因此通道內(nèi)的物性變化很大，換熱特性隨之變化。為避免換熱特性變差威脅堆芯熱工安全，有必要對(duì)大長(zhǎng)徑比通道內(nèi)冷卻劑流動(dòng)換熱特性進(jìn)行研究。

基于美俄開(kāi)展的核熱推進(jìn)計(jì)劃，研究人員開(kāi)展了大量的氫氣流動(dòng)換熱特性理論和實(shí)驗(yàn)研究，獲取了豐富的研究資料。McCarthy等[5]、Taylor[6-7]、Petukhov等[8]、Perkins[9]、McEligot等[10]開(kāi)展了涵蓋反應(yīng)堆熱工參數(shù)的氫氣流動(dòng)換熱實(shí)驗(yàn)研究，獲取了相關(guān)的實(shí)驗(yàn)關(guān)聯(lián)式。Hess與Kunz[11]對(duì)超臨界氫氣流動(dòng)換熱特性開(kāi)展理論分析并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比，絕大多數(shù)數(shù)據(jù)誤差在±20%之內(nèi)。Sajan[12]針對(duì)典型核熱推進(jìn)反應(yīng)堆堆芯單通道進(jìn)行了參數(shù)敏感性分析。受限于當(dāng)時(shí)實(shí)驗(yàn)研究條件的落后和核熱推進(jìn)計(jì)劃的終止，氫氣在反應(yīng)堆通道內(nèi)熱工水力行為的研究仍存在著一定的缺陷，此外，高溫、高流速、高熱流密度等特殊環(huán)境下的氫氣流動(dòng)換熱特性的研究資料較少，相關(guān)數(shù)值計(jì)算模型還不完善。本研究采用數(shù)值計(jì)算方法開(kāi)展高溫、高流速氫氣在圓管內(nèi)流動(dòng)換熱特性研究，并驗(yàn)證模型及算法計(jì)算的正確可行性。

1 數(shù)學(xué)物理模型與計(jì)算方法

1.1 控制方程

本研究基于圓管流道采用二維軸對(duì)稱模型進(jìn)行數(shù)值模擬，模型假設(shè)如下：穩(wěn)態(tài)無(wú)內(nèi)熱源，忽略輻射換熱、化學(xué)反應(yīng)等。質(zhì)量、動(dòng)量、能量守恒方程如下：

(1)

(2)

(3)

div(λ·grad(T))+Φdissipation

(4)

式中：ρ為密度，kg·m-3；u、v分別為軸向速度、徑向速度，m/s；x、r分別為軸向位置、徑向位置，m；p為壓力，Pa；T為溫度，K；μ為動(dòng)力黏度，Pa·s；cp為比定壓熱容，J/(kg·K)；λ為導(dǎo)熱系數(shù)，W/(m·K)；Φdissipation為黏性耗散項(xiàng)。

1.2 數(shù)值計(jì)算方法

本研究采用商用CFD軟件ANSYS FLUENT，求解器為基于壓力的耦合方法，同時(shí)求解質(zhì)量、動(dòng)量方程后再進(jìn)行能量、湍流方程等標(biāo)量方程的求解，收斂速度更快，更適合用于控制方程存在較強(qiáng)依賴關(guān)系的高速領(lǐng)域?？臻g離散采用二階迎風(fēng)格式。湍流模型采用SSTk-ω模型，并考慮壓縮效應(yīng)和黏性加熱效應(yīng)。

1.3 邊界條件

本研究計(jì)算模型幾何結(jié)構(gòu)如圖1所示，其中圓管內(nèi)徑2.95 mm，粗糙度0.44 μm，入口段長(zhǎng)度41.25 mm，加熱段長(zhǎng)度228.60 mm。不計(jì)算壁面固體域。計(jì)算域邊界條件如下。1) 入口：質(zhì)量流量、溫度；2) 出口：壓力；3) 入口段壁面：絕熱、無(wú)滑移壁面、給定壁面粗糙度；4) 加熱段壁面：第2類邊界條件、無(wú)滑移壁面、給定壁面粗糙度。

圖1 計(jì)算模型幾何結(jié)構(gòu)Fig.1 Structure of calculation model

1.4 網(wǎng)格無(wú)關(guān)性驗(yàn)證

以質(zhì)量流量m=0.870 6 g/s、入口流體溫度Tb,in=315.56 K、出口壓力pout,abs=147.33 kPa、壁面熱流密度q=1 849 864 W/m2為基礎(chǔ)工況進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，本研究共給出了5套網(wǎng)格，網(wǎng)格信息列于表1。

表1 網(wǎng)格無(wú)關(guān)性驗(yàn)證Table 1 Grid independency verification

以出口主流總溫Tb,total,out、進(jìn)出口壓差Δp作為參考量，網(wǎng)格無(wú)關(guān)性驗(yàn)證結(jié)果表明，在節(jié)省計(jì)算資源和保證計(jì)算精度的條件下，選用網(wǎng)格36×1 800即可滿足需要。

2 氫氣物性模型

氫氣屬于小分子量(摩爾質(zhì)量Mmole=2.015 9 g/mol)輕氣體，存在兩種自旋異構(gòu)體形態(tài)：正氫(Ortho-hydrogen)和仲氫(Para-hydrogen)，兩者物理性質(zhì)有所差異，化學(xué)性質(zhì)相同并可相互轉(zhuǎn)化。在室溫環(huán)境中，處于平衡狀態(tài)的氫氣稱為標(biāo)準(zhǔn)氫(normal hydrogen)，其中正氫與仲氫體積之比為3∶1。標(biāo)準(zhǔn)氫的臨界壓力為1.296 4 MPa，臨界溫度為33.145 K，臨界密度為31.262 kg/m3。本研究采用標(biāo)準(zhǔn)氫作為冷卻工質(zhì)進(jìn)行數(shù)值計(jì)算，在計(jì)算工況范圍內(nèi)，氫氣物性變化很大，物性變化對(duì)流動(dòng)換熱影響很大，因此本研究采用變物性模型。

2.1 狀態(tài)方程

計(jì)算工況范圍溫度遠(yuǎn)超過(guò)氫氣的臨界溫度、壓力低于臨界壓力，氫氣作為實(shí)際氣體考慮，狀態(tài)方程采用Aungier-Redlich-Kwong(ARK)模型：

(5)

2.2 熱物理與輸運(yùn)性質(zhì)

氣體分子隨溫度的升高裂解成更簡(jiǎn)單的分子或原子，這種現(xiàn)象稱為熱解或熱離(thermal dissociation or atomization)。在足夠高的溫度下，氣體分子可全部熱解成原子。在更高的溫度下，氣體還會(huì)發(fā)生電離過(guò)程，電子從原子中分離逃逸形成等離子體。

常壓狀態(tài)下，氫氣H2的熱解溫度約為1 500 K[13-14]，在溫度達(dá)到5 000 K以上時(shí)，H2分子幾乎全部熱解變成原子氫氣體。在核熱推進(jìn)反應(yīng)堆中，堆芯溫度很高且伴有各類放射性射線和高能帶電粒子，這導(dǎo)致了氫氣分子的電離和熱解[15]。

研究結(jié)果表明，溫度、壓力以及堆芯中能量沉積均對(duì)氫氣熱解有重要影響，圖2為不同溫度、壓力下的氫氣熱解份額[13-14]，處于熱解平衡狀態(tài)下的分子氫與原子氫熱物性與輸運(yùn)性質(zhì)發(fā)生了很大的變化[16]。物性變化使冷卻劑流動(dòng)換熱特性發(fā)生變化，欲真實(shí)模擬氫氣在管內(nèi)熱工水力行為，須考慮高溫氫氣熱解對(duì)物性的影響。

圖2 氫氣熱解份額Fig.2 Degree of thermal dissociation of hydrogen

流體的熱物理與輸運(yùn)特性取決于溫度和密度(或壓力)，為簡(jiǎn)化計(jì)算過(guò)程、節(jié)省計(jì)算資源，通常擬合為溫度和壓力的多項(xiàng)式形式。由于工況計(jì)算范圍內(nèi)物性受溫度變化影響程度遠(yuǎn)大于壓力，本研究忽略壓力變化對(duì)物性的影響，將物性簡(jiǎn)化為溫度的單值函數(shù)。NIST(National Institute of Standards and Technology)給出了氫氣物性計(jì)算數(shù)據(jù)庫(kù)，但溫度適用范圍上限僅為1 000 K且未考慮高溫?zé)峤鈱?duì)物性的影響，無(wú)法滿足計(jì)算需要。在相應(yīng)高溫條件下，氫工質(zhì)為熱解平衡狀態(tài)下的氫原子和氫分子的混合物。本文根據(jù)NASA(National Aeronautics and Space Administration)和NBS(National Bureau of Standards)推薦的熱平衡條件下氫氣物性數(shù)據(jù)[17]擬合了多項(xiàng)式函數(shù)，更適合本文研究。

(6)

式中：Y代表cp、λ、μ等變量；N為系數(shù)。

熱物理與輸運(yùn)特性列于表2。

表2 氫氣低壓條件下的熱物理與輸運(yùn)特性Table 2 Thermo-physical and transport property of low pressure hydrogen

3 模型驗(yàn)證與換熱特性分析

3.1 驗(yàn)證實(shí)驗(yàn)

Taylor[6]研究了氫氣、氦氣在鎢金屬圓管內(nèi)流動(dòng)換熱特性，其中熱流密度由實(shí)驗(yàn)段電阻加熱提供，實(shí)驗(yàn)段結(jié)構(gòu)與尺寸參數(shù)如1.3節(jié)所述。實(shí)驗(yàn)范圍雷諾數(shù)Re=7 600～39 500，實(shí)驗(yàn)段入口壓力0.27～0.70 MPa，Tw/Tb最大值為5.6。

Taylor實(shí)驗(yàn)采用電阻式加熱方式，即在實(shí)驗(yàn)段兩側(cè)直接加載電流，由于鎢金屬圓管材料的電阻率隨溫度的升高而增大，實(shí)驗(yàn)中軸向上的熱流密度并不均勻。由式(7)可知，加熱段軸向熱流密度分布實(shí)際為壁面溫度的函數(shù)。

(7)

本研究選取驗(yàn)證工況列于表3[6]。圖3為驗(yàn)證工況熱流密度qe和壁溫Tw分布[6]。

表3 驗(yàn)證工況Table 3 Verification case

由圖3可知，壁溫分布與熱流密度分布基本一致。在較低熱流密度時(shí)，如Run1、3、7，實(shí)驗(yàn)段前半部分熱流密度、壁溫分布近似線性，在較高熱流密度時(shí)迅速上升達(dá)到峰值。如Run19、20、21，質(zhì)量流量基本相同，而隨功率輸入的增加，兩個(gè)因素導(dǎo)致加熱段入口部分壁面溫升較大。首先，Tw/Tb增加伴隨著傳熱系數(shù)降低，進(jìn)一步提高了壁面溫度。其次，在較高溫度下鎢電阻率增加放大了提高Tw/Tb的效果。連接部件傳導(dǎo)損耗使加熱段入口處和出口處的軸向溫度梯度較大。

3.2 結(jié)果對(duì)比

本研究基于Taylor實(shí)驗(yàn)進(jìn)行驗(yàn)證分析，研究湍流模型對(duì)壁溫分布的影響。Taylor實(shí)驗(yàn)結(jié)果處理中考慮了輻射散熱、軸向?qū)岬?利用熱平衡計(jì)算得到了等效熱流密度分布，本研究將得到的熱流密度分布作為邊界條件輸入，對(duì)比了Standard(標(biāo)準(zhǔn))k-ε、Realizable(可實(shí)現(xiàn))k-ε以及SSTk-ω等3種湍流模型的適用性，計(jì)算值與實(shí)驗(yàn)工況值對(duì)比結(jié)果如圖4所示，其中實(shí)驗(yàn)工況值的誤差帶為±10%。結(jié)果表明，SSTk-ω模型的適用性較其他兩種模型要好得多。這主要是在研究范圍內(nèi)近壁面處溫度、速度、物性等變化非常劇烈，湍流動(dòng)能、比耗散率等很大，而SSTk-ω模型直接穿過(guò)黏性底層應(yīng)用于壁面[18]，在近壁面處的處理要優(yōu)于k-ε模型，因此計(jì)算值更接近于實(shí)驗(yàn)值。在實(shí)驗(yàn)加熱段的進(jìn)出口處，計(jì)算值與實(shí)驗(yàn)值相差較大，主要原因是計(jì)算工況未考慮軸向?qū)帷?shí)驗(yàn)加熱段有連接部件傳導(dǎo)損耗。此外，實(shí)驗(yàn)中還存在測(cè)量位置誤差、測(cè)量?jī)x器誤差等。本研究認(rèn)為，采用SSTk-ω湍流模型計(jì)算高溫、高流速氫氣圓管內(nèi)流動(dòng)換熱特性是合理可靠的。

3.3 換熱特性分析

本文以質(zhì)量流量m=0.870 6 g/s、入口流體溫度Tb,in=315.56 K、出口壓力pout,abs=147.33 kPa、壁面熱流密度q=1 849 864 W/m2為基礎(chǔ)工況，進(jìn)行熱流密度加熱條件下變物性可壓縮高溫氫氣流動(dòng)換熱特性分析。通過(guò)改變邊界條件，研究質(zhì)量流量和壁面熱流密度等熱工參數(shù)對(duì)換熱特性的敏感性影響。

1) 流場(chǎng)與溫度場(chǎng)分布

本研究對(duì)基礎(chǔ)工況進(jìn)行流道內(nèi)的流場(chǎng)與溫度場(chǎng)分布情況分析。由于流道直徑較小，氫氣密度較低，因此即使入口質(zhì)量流量很低，入口的速度量級(jí)也相當(dāng)可觀。氫氣流場(chǎng)與溫度場(chǎng)分布示于圖5。

圖3 驗(yàn)證工況熱流密度和壁溫分布Fig.3 Heat flux and wall temperature distributions of verification cases

圖4 不同湍流模型計(jì)算值與實(shí)驗(yàn)值對(duì)比結(jié)果Fig.4 Comparison of calculation values and experiment data under different turbulent models

圖5 氫氣流場(chǎng)與溫度場(chǎng)分布Fig.5 Distributions of flow and temperature fields of hydrogen

由圖5可知，隨軸向壓力的降低、溫度的升高，流道內(nèi)的氫氣膨脹密度降低、速度升高，并在靠近出口位置處馬赫數(shù)Ma>1，實(shí)現(xiàn)了跨聲速流動(dòng)。由于入口段效應(yīng)，壁溫與主流溫度之比也出現(xiàn)先增大后減小的趨勢(shì)。此外，黏度、比熱、導(dǎo)熱系數(shù)等物性隨軸向主流溫度的升高而升高，軸向雷諾數(shù)Re、努塞爾數(shù)Nu均降低，即管內(nèi)氫氣流動(dòng)換熱變差。近壁面的氫氣密度較小、黏度很大形成黏滯力，阻礙了徑向傳熱并加劇了壁溫升高，進(jìn)一步造成壁面與主流溫差增大形成正反饋，因此，導(dǎo)致氫氣換熱特性變差。軸向不同位置處的徑向速度與徑向溫度分布如圖6所示，出口位置處的近壁面與中心線溫差可達(dá)500 K以上。

圖6 軸向不同位置處速度與溫度徑向分布Fig.6 Radial distributions of velocity and temperature at different axial positions

圖7 熱工參數(shù)對(duì)流動(dòng)換熱特性的敏感性分析Fig.7 Sensitive analysis of flow and heat transfer characteristics under different thermal parameters

2) 參數(shù)敏感性分析

為探究熱工參數(shù)對(duì)氫氣在圓管內(nèi)的流動(dòng)換熱影響情況，將基礎(chǔ)工況的入口流量、壁面熱流密度分別變化±10%、±20%，進(jìn)行了參數(shù)敏感性分析。其中摩擦因子計(jì)算如下。

Δp=pin-pout=Δpa+Δpf

(8)

(9)

(10)

式中：Δp壓降，Pa；G為質(zhì)量流速，kg/(m2·s)；f為摩擦因子；下標(biāo)in、out、a、f、av分別為入口、出口、加速值、摩擦值、平均值。

圖7為熱工參數(shù)對(duì)流動(dòng)換熱特性的敏感性分析。結(jié)果表明，隨入口流量的增大，Nu增大、f減小，即增大入口流量起到了強(qiáng)化換熱的效果。入口流量每增大10%，Nu增大10%～11%、f減小約4%。隨熱流密度的增加，Nu與f均呈下降趨勢(shì)，即增大熱流密度對(duì)氫氣在圓管內(nèi)換熱不利。熱流密度每增大10%，Nu減小3%～4%、f減小不足1%。

當(dāng)傳熱壁面與流動(dòng)氣體之間存在較大溫差時(shí)，流動(dòng)截面與流動(dòng)方向上物性變化很大，管內(nèi)氣體傳熱不再符合一般的定常物性傳熱關(guān)系式。近壁面處的氣體溫度近似等于壁溫，徑向上的物性隨氣體溫度變化而變化，而影響氣體密度的壓力在垂直于流動(dòng)的橫截面上通常是恒定的。因此，考慮變物性、適合于定壁溫或定熱流條件的氣體流動(dòng)換熱關(guān)系式采用溫度、熱入口段修正，即：

(11)

文獻(xiàn)中常用的適用于氫氣流動(dòng)換熱的實(shí)驗(yàn)關(guān)聯(lián)式列于表4[19]。

為研究不同實(shí)驗(yàn)關(guān)聯(lián)式的適用性，本研究將計(jì)算工況結(jié)果與表4的實(shí)驗(yàn)關(guān)聯(lián)式展開(kāi)對(duì)比，計(jì)算工況Nu[20-21]為：

(12)

表4 氫氣流動(dòng)換熱實(shí)驗(yàn)關(guān)聯(lián)式Table 4 Heat transfer correlation of hydrogen

式中：Pr為普朗特?cái)?shù)；h為對(duì)流換熱系數(shù)，W/(m2·K)；下標(biāo)w、b表示壁面、主流。

圖8為計(jì)算值與換熱模型對(duì)比。由圖8可知，計(jì)算值與實(shí)驗(yàn)關(guān)聯(lián)式相對(duì)誤差均在±6%以內(nèi)，本研究范圍內(nèi)計(jì)算工況符合實(shí)驗(yàn)關(guān)聯(lián)式模型的。此外，考慮了溫度修正的Miller-Taylor實(shí)驗(yàn)關(guān)聯(lián)式與計(jì)算值偏差最小，因此本研究推薦該關(guān)聯(lián)式評(píng)價(jià)圓管內(nèi)變物性、高溫、高流速氣體流動(dòng)換熱特性。

圖8 計(jì)算值與換熱模型對(duì)比Fig.8 Comparison of calculation value and result of heat transfer model

4 小結(jié)

核熱推進(jìn)技術(shù)是未來(lái)深空探測(cè)、載人航天最有競(jìng)爭(zhēng)力的技術(shù)之一，反應(yīng)堆熱工設(shè)計(jì)又是最關(guān)鍵問(wèn)題。本文采用數(shù)值計(jì)算方法開(kāi)展高溫、高流速氫氣圓管內(nèi)流動(dòng)換熱特性研究，可得出以下結(jié)論。

1) 采用壓力基耦合算法、SSTk-ω湍流模型以及物性模型研究高溫、高流速氫氣圓管內(nèi)流動(dòng)換熱特性合理可行，壁溫計(jì)算結(jié)果與實(shí)驗(yàn)值符合較好，驗(yàn)證了數(shù)值計(jì)算模型的選擇是正確的。

2) 由于軸向壁溫升高，近壁面氫氣黏度變大、密度變小，徑向傳熱較差使換熱變差，進(jìn)一步提升了壁溫，形成正反饋，因此熱工設(shè)計(jì)時(shí)要著重考慮變物性對(duì)換熱特性的影響。

3) 通過(guò)進(jìn)行入口流量、壁面熱流密度等熱工參數(shù)敏感性分析發(fā)現(xiàn)，隨入口流量的增加，換熱效果增強(qiáng)，隨熱流密度的增加換熱效果變差。研究范圍內(nèi)，本文推薦修正Miller-Taylor實(shí)驗(yàn)關(guān)聯(lián)式評(píng)價(jià)圓管內(nèi)變物性、高溫、高流速氣體流動(dòng)換熱特性。