圓管內氦氙混合氣體與超臨界二氧化碳換熱特性對比分析

寧可為, 劉凱, 孫汝雷, 趙富龍, 游爾勝, 余霖, 譚思超

(1.哈爾濱工程大學 核科學與技術學院, 黑龍江 哈爾濱 150001; 2.哈爾濱工程大學 黑龍江省核動力裝置性能與設備重點實驗室, 黑龍江 哈爾濱 150001; 3.中國核動力研究設計院, 四川 成都 610041)

陸上長時供能是小型核反應堆未來的應用方向之一。出于任務場景對體積、重量的限制,傳統壓水堆已無法滿足要求,新型冷卻劑工質得到了廣泛的關注,以氦氙混合氣體為代表的氣體冷卻劑與超臨界二氧化碳為代表的超臨界工質具有壓縮性能好、熱物性穩定、傳熱能力優秀等特點,適合與能量轉換系統結合,構成運行安全、輕小靈活的核動力系統。

冷卻劑換熱特性對反應堆運行性能具有重要影響。近年來,諸多學者以空間應用為背景開展了氦氙混合氣體物性及流動換熱的實驗及數值模擬研究,Tournier等[1]基于查普曼-恩斯科格(Chapman-Enskog)輸運理論及赫希菲爾德(Hirschfelder)氣體混合方法提出了用于預測二元惰性混合氣體物性的半經驗公式;Vitovsky等[2]針對空間反應堆緊密排布的棒束堆芯結構開展了圓形及準三角通道流動換熱實驗,研究得到了氦氙混合氣體的流動結構,并證明適合采用質量流量平均溫度作為定性溫度;Dragunov等[3]對不同通道形式下的層流、湍流及變物性條件下的流動換熱關系式進行了對比,其推薦采用別祖霍夫(Petukhov)關系式用于氦氙混合氣體流動換熱計算,該結論與秦浩模擬計算所得相同,同時,Qin等[4]進一步驗證發現,未經修正的原始Petukhov關系式計算結果預測精度較高。

超臨界二氧化碳作為性能優異的新型冷卻劑工質,其跨臨界及超臨界態下的流動換熱特性得到國內外學者的廣泛研究。Liu等[5]以TID(tube in duct)型燃料組件為基礎,完成了300 MWt超臨界二氧化碳反應堆概念設計;劉新新[6]采取實驗測量和數值模擬方法研究了超臨界二氧化碳在加熱直管和螺旋管的換熱特性和換熱惡化機理;Song等[7]實驗對比了超臨界二氧化碳在不同管徑豎直管內的換熱特性,研究發現超臨界二氧化碳具有優秀的熱力學性質,但在不同的流動條件下,其換熱規律差異較大。

綜上所述,目前的研究對象基本局限于單一工質,對于不同冷卻劑工質換熱特性的對比分析不夠充分。本文將采用數值計算方法對堆芯圓形冷卻劑通道的流動換熱特性進行分析,研究不同參數下氦氙混合氣體及超臨界二氧化碳對堆芯流動換熱特性的影響機理,討論二者作為工質時的反應堆安全特性。

1 幾何物理模型及網格劃分

1.1 幾何模型

陸上多用途的小型核反應堆要求體積高度緊湊,因此,提高燃料比例、組件緊密排布的堆芯設計方案應運而生,六棱柱燃料元件即解決思路之一。中空六棱柱燃料元件以燃料作為基體,圓形冷卻劑孔道在基體內部呈三角排布,該設計能夠增加燃料體積分數,提供更高的熱慣性,有助于抵御異常運行工況下的功率激增。用空氣直接冷卻的Tory-ⅡC型沖壓發動機已對該燃料型式進行了成功實驗;Pope[8]及Liu 等[5]采用該型燃料元件進行了超臨界二氧化碳反應堆概念設計,結果表明堆芯物理及流動換熱特性均滿足反應堆運行要求。本文計算采用7冷卻劑孔道中空六棱柱燃料元件,截面六邊形對邊距14.2 mm,冷卻劑孔道直徑4 mm,燃料元件長0.69 m,具體型式如圖1所示。

圖1 計算模型及簡化Fig.1 Calculation model and simplification

對于該燃料組件型式,由于關注重點在于流體域,因此對模型進行適當簡化。Lu等[9]通過計算發現,中空六棱柱型的燃料元件與圓環型的燃料元件在流體域部分的流動參數特性誤差最大在1.45%左右,因此對流體域單獨建模,得到直徑4 mm的圓形冷卻劑通道。簡化后的模型可以增大網格密度,增加網格的精細度,以提高計算的精確度。以該模型為基礎,開展氦氙混合氣體與超臨界二氧化碳流動換熱對比分析。

1.2 物理模型

本文基于計算流體力學(computational fluid dynamics, CFD)軟件進行堆芯流動換熱數值模擬計算。其中,分離求解器是基于壓力的求解器,按順序逐一的求解每個方程,主要用于不可壓流動和低速微可壓流動,對于冷卻劑最高速度超過0.3 Ma的情況,應當視為不可壓流動,采用耦合流模型;而當冷卻劑最高速度不超過0.3 Ma時,可視為低速微可壓流動,采用分離流模型。

對于氦氙混合氣體,以美“普羅米修斯”反應堆為參考,其冷卻劑最高速度約40 m/s,因此計算采用分離流模型。Meng[10]以實驗值為參考,對多種氦氙混合氣體流動的計算模型進行了對比,結果發現相比于SST(shear stress transport) 模型和V2F模型,RKE模型符合程度較好,計算精度較高。k-ε模型通常適用于熱傳遞的工業型應用,并可有效平衡穩定性、計算量和精度。氦氙混合氣體在堆芯中的流動同時存在徑向與縱向熱交換,因此選擇k-ε模型。其中標準k-ε模型和可實現k-ε模型均為高雷諾數模型,可實現的k-ε兩方程模型可以為全y+壁面處理帶來更多靈活性,因此本文計算基于可實現k-ε兩方程模型開展。

對于超臨界二氧化碳,大量學者已開展直管、螺旋管內的傳熱特性實驗與數值模擬研究,He[11]以實驗值和直接數值模擬(DNS)計算結果為參考,提出k-ε或V2F模型具有較高的預測精確性。由于標準的k-ε方程只適用于離開壁面一定距離的高Re數湍流區域,壁面附近低Re數區域須由壁面函數法來修正相關的參數[12]。本文研究中,工質流動的Re基本均大于104,對比其結論得出的工況,選擇k-ε模型、全y+壁面處理完成計算。

1.3 網格劃分及無關性驗證

為確保計算結果精度不受網格影響,且節省計算時間,在計算前首先完成兩模型的網格無關性驗證。為了將近壁面流體計算更精確,將流體域進行了近壁面網格加密處理,采用了棱柱層網格在流體域近壁面處繪制邊界層網格。網格無關性驗證結果見圖2。由圖2可知,對于氦氙混合氣體模型,不同網格尺寸對于溫度參數的影響微乎其微,但對于對流換熱系數的影響較大,低網格數量下,入口段效應體現并不明顯,這將對工質換熱對比產生不利影響。網格尺寸在2 mm處與0.3 mm處相對誤差為1.12%,在可接受范圍內。因此,為提高計算效率,網格尺寸確定為2 mm,y+=0.96,網格數量為1 218 434。

表1 網格無關性驗證參數及結果Table 1 Grid independence verification parameters and results

圖2 計算模型網格無關性驗證Fig.2 Verification of grid independence of computational model

2 計算結果及對比分析

2.1 工質物性對比

作為二元混合氣體,氦氙混合氣體物性與混合比例直接相關。混合氣體的熱物性采用基于維里系數的稀有氣體修正公式計算得出[1]。不同混合比例對堆內通道流動換熱具有顯著影響,余霖[15]通過研究發現,單冷卻劑通道活性區域內,30～40 g/mol為氣體流動換熱最佳的摩爾質量比例區間,在2 MPa下,31.5 g/mol和40 g/mol氦氙混合氣體物性隨溫度變化如圖3所示。可以看出,不同混合比例并不影響物性的變化趨勢,其中,較小的混合比例具有較小的密度、較大的定壓比熱與導熱能力。對于動力粘度而言,不同混合比例的影響微乎其微。

圖3 氦氙混合氣體及超臨界二氧化碳物性隨溫度變化Fig.3 Variation of physical properties of He-Xe mixture and S-CO2 over temperature

作為一元工質,超臨界二氧化碳物性僅受溫度、壓力的影響,其超臨界狀態使得物性變化趨勢與常規流體存在較大差異,該特點主要表現為物性參數在擬臨界區附近存在劇烈變化。以10、17 MPa為例,定壓比熱在臨界點附近存在一局部極大值,而在亞臨界向超臨界跨越的過程中,密度、動力粘度出現陡降,在擬臨界區之外,超臨界二氧化碳物性幾乎不存在顯著變化。壓力對物性的影響同樣出現在臨界點附近,跨過擬臨界區后,壓力提升對物性參數的影響較小。

通過圖3中2種冷卻劑工質物性對比可以看出,在溫度大于400 K后,超臨界二氧化碳的密度高于氦氙混合氣體,且其具有2倍于氦氙混合氣體的定壓比熱,這使得超臨界二氧化碳具有較強的載熱能力;氦氙混合氣體導熱系數高于超臨界二氧化碳,但就數值而言,只有約0.15～0.25 W/(m·K),二者導熱能力均較差。

2.2 流動換熱特性對比

為對比二者流動換熱特性,明確不同冷卻劑的熱工安全性能,針對飛行器閉式反應堆堆芯開展不同運行工況下的反應堆熱工安全特性對比。以局部對流換熱系數與平均對流換熱系數作為主要對比參數,從熱流密度、溫度、速度等方面對氦氙混合氣體與超臨界二氧化碳進行對比計算。其中,為了獲得相對較優的換熱能力,氦氙混合氣體相對分子質量選擇31.5 g/mol,該值與美“普羅米修斯”方案保持一致;超臨界二氧化碳壓力的選擇需考慮能量轉換系統特性[16-17],以17 MPa作為總壓進行計算。

2.2.1 加熱功率對換熱系數的影響

為了研究不同加熱功率對換熱系數的影響,在入口溫度、流速相同的條件下,改變加熱功率,得到2種冷卻劑管內流動局部對流換熱系數。入口溫度及流速分別為1 200 K、10 m/s,對加熱功率為5、20、80、140 kW/m2下的工況展開計算。其中,氦氙混合氣體出口壓力2 MPa,超臨界二氧化碳取17 MPa,加熱方式均取余弦加熱。

取點時,沿軸向等距取10個點,導出每點處對應截面的流體平均溫度與壁面上的流體溫度,二者之差即為該位置處的溫度差。計算得到熱流密度數據,根據牛頓冷卻公式得到傳熱系數。圖4(a)、(b)為不同加熱功率下,2種冷卻劑工質各自的對流換熱系數沿軸向的分布。

圖4 不同加熱功率下的氦氙混合氣體及超臨界二氧化碳對流換熱系數沿軸向分布Fig.4 Axial distribution of convective heat transfer coefficients of He-Xe mixture and S-CO2 under different heating power

由圖4可以看出,對于氦氙混合氣體和超臨界二氧化碳,局部對流換熱系數均呈現入口段先下降后上升、湍流充分發展段基本不變、出口位置下降的趨勢,這主要受入口段效應和出口冷卻劑物性變化的影響。入口段熱邊界層存在從零開始生長直至充分發展的過程,后由于邊界層中出現湍流,其擾動與混合作用使對流換熱系數出現波動,并最終趨于穩定;出口段區域,冷卻劑溫度升高、壓力降低,兩者共同作用導致密度加劇降低,冷卻劑載熱能力相對降低,此時,變物性效應導致冷卻劑與固體壁面溫差相對增加,對流換熱系數降低。

整體而言,加熱功率增大會導致對流傳熱系數先增大后維持,低功率下加熱量增加對于對流換熱系數的影響較大;當處于較高功率水平后,功率再次提升所帶來的影響主要反映在入口段。這主要受加熱壁面與流體溫差的影響。低功率下的冷卻劑攜帶熱量的能力充足,因此,功率增加時,固體仍能夠有效地加熱流體;功率到達一定程度時,冷卻劑無法充分帶走壁面熱量,導致氣體和固體間的溫差增大,此時對流傳熱系數基本維持不變,甚至出現下降趨勢。當功率增大到一定值時,由于冷卻劑攜帶熱量有限,入口段效應帶來的氣體溫度上升幅度減小,固體溫度與流體溫差減小,入口段帶來的影響效果減弱。

由圖4可以看出,對于2種工質,在中間段湍流充分發展區域,對流換熱系數均保持基本穩定,不出現較大波動。因此,在該區段取點并作平均,將其作為整段管道的平均對流換熱系數。本文取0.2、0.28、0.35、0.4、0.5 m處的局部對流換熱系數,取5個點的平均值得到整管平均對流換熱系數。

通過將沿軸向位置處的傳熱系數平均,得到不同加熱功率下的平均對流換熱系數如圖5所示。

圖5 氦氙混合氣體與超臨界二氧化碳對流換熱系數隨加熱功率的變化Fig.5 Variation of convective heat transfer coefficient of He-Xe mixture and S-CO2 with heating power

通過對比可以發現,在5～100 kW/m2的加熱功率區間內,超臨界二氧化碳對流換熱系數始終高于氦氙混合氣體,該現象反映了二者換熱能力的區別。

在該工況區間下,超臨界二氧化碳的換熱能力顯著高于氦氙混合氣體,這主要由流體導熱系數、密度、普朗特數等物性參數所決定。

2種工質的對流換熱系數對加熱功率均不十分敏感,在中低功率水平時,超臨界二氧化碳對流換熱系數變化較低,隨加熱功率增加,對流換熱系數小幅上升;與之對比,氦氙混合氣體的對流換熱系數在小加熱功率下的增加幅度較大,而在大功率下仍維持在較高的水平。該現象主要反映了2種不同工質載熱能力的差別。功率水平較高時,超臨界二氧化碳載熱的能力不足以帶走壁面加熱量,故其換熱溫差相對增加值大于氦氙混合氣體。

2.2.2 入口速度對換熱系數的影響

為了研究不同入口速度對換熱系數的影響,在入口溫度、加熱功率相同的條件下,改變入口流速。控制壁面恒熱流密度為80 kW/m2,定入口溫度為1 200 K,入口流速分別為1、5、10、15、20 m/s。其中,考慮到模型在低雷諾數下存在較大的計算誤差及堆內實際應用條件,舍去氦氙混合氣體1 m/s、超臨界二氧化碳20 m/s的工況。

圖6(a)、(b)為不同入口速度下,2種冷卻劑工質各自的對流換熱系數沿軸向的分布。由圖6可知,當入口速度變化時,局部對流換熱系數分布趨勢基本保持一致,即存在入口段、充分發展段和出口段。入口段對流換熱系數先下降后上升,充分發展段基本保持不變,而出口段對流換熱系數大幅降低。2種冷卻劑局部對流換熱系數沿軸向分布趨勢的差異主要表現在充分發展段,氦氙混合氣體在軸向距離中間至后部的位置對流換熱系數基本不變或稍有提升,而超臨界二氧化碳則隨軸向距離增加而始終保持下降趨勢。

圖6 不同入口速度下的氦氙混合氣體及超臨界二氧化碳對流換熱系數沿軸向分布Fig.6 Axial distribution of convective heat transfer coefficients of He-Xe mixture and S-CO2 under different inlet velocity

不同入口速度下的平均對流換熱系數如圖7所示。由圖7可知,隨著入口速度增加,2種冷卻劑工質對流換熱系數均增大,對流換熱效果增強。而在入口速度大于1 m/s后,超臨界二氧化碳對流換熱系數與其隨流速增長的斜率均顯著高于氦氙混合氣體。從數值上看,當流速從5 m/s增長至20 m/s時,氦氙混合氣體對流換熱系數僅從620 W/(m2·K)增長至1 013 W/(m2·K),且在15 m/s后,對流換熱系數隨入口速度的增幅放緩,通過增大流速以改善換熱的效果有限;而超臨界二氧化碳的對流換熱系數隨入口流速近似線性增加,在15 m/s時達到3 060 W/(m2·K)。因此,在反應堆設計時,氦氙混合氣體應盡量選擇較大流速以實現較優的換熱效果;而超臨界二氧化碳可適當放寬對流速的要求以滿足系統運行的要求。

圖7 氦氙混合氣體與超臨界二氧化碳對流換熱系數隨入口速度的變化Fig.7 Variation of convective heat transfer coefficient of He-Xe mixture and S-CO2 with inlet velocity

2.2.3 入口溫度對換熱系數的影響

為了研究不同入口溫度對換熱系數的影響,在入口流速為10 m/s、加熱功率為80 kW/m2條件下,改變入口溫度,分別為600、900、1 200、1 500 K。

圖8(a)、(b)為不同入口速度下,2種冷卻劑工質各自的對流換熱系數沿軸向的分布。由圖8可知,當入口溫度發生變化時,氦氙混合氣體的局部對流換熱系數變化較為劇烈,而超臨界二氧化碳趨勢與之前相似,仍有較為明顯的入口段、充分發展段和出口段。

圖8 不同入口溫度下的氦氙混合氣體及超臨界二氧化碳對流換熱系數沿軸向分布Fig.8 Axial distribution of convective heat transfer coefficients of He-Xe mixture and S-CO2 under different inlet temperature

不同入口溫度下的平均對流換熱系數如圖9所示。由圖9可知,隨入口溫度的增加,工質對流換熱系數呈不同趨勢。氦氙混合氣體隨入口溫度增加,管內平均換熱系數先增加后下降,在1 000～1 100 K左右有換熱系數峰值,約800 K/(m2·k),而超臨界二氧化碳管內平均換熱系數則隨入口溫度的增加持續下降,對于所計算工況,在600 K時有換熱系數最大值,且在入口溫度600～800 K內,換熱系數下降速度大于800～1 200 K。該現象出現的原因主要與工質物性隨溫度變化有關,對于2種冷卻劑工質,當入口溫度增大,導熱系數增加,工質流動的雷諾數降低,而氦氙混合氣體普朗特數Pr略有增加,超臨界二氧化碳Pr有所減小,綜合效果為局部換熱系數的降低,但Pr變化的差異導致兩者曲線呈現不同形態。

圖9 氦氙混合氣體與超臨界二氧化碳對流換熱系數隨入口溫度的變化Fig.9 Variation of convective heat transfer coefficient of He-Xe mixture and S-CO2 with inlet temperature

兩者對比可知,超臨界二氧化碳在較低入口溫度下的換熱效果良好,而氦氙混合氣體的換熱系數隨溫度先上升后下降,在1 000 K左右達到峰值;而在數值上,超臨界二氧化碳對流換熱系數始終高于氦氙混合氣體,對流換熱能力較強。

通過以上因素的分析可知,對于2種冷卻劑工質而言,為了達到較優的堆內換熱效果,在設計時應選擇高加熱功率、大流速,現有超臨界二氧化碳反應堆概念設計中[5,18],堆內冷卻劑流速取值在5～7 m/s,氦氙混合氣體冷卻劑流速約20 m/s; 入口溫度的選擇,氦氙混合氣體適合選擇在1 000～1 200 K,超臨界二氧化碳則適合選擇較低的入口溫度,現有概念設計在550～650 K左右。

3 對流換熱模型修正

盡管以牛頓冷卻公式為基礎計算出的對流換熱系數能夠對冷卻劑管內流動的對流換熱現象做出較好描述,但該參數無法對流動過程做出準確描述,并且基于溫差的計算過程使其討論過程中存在諸多的不確定性。以無量綱參數為基礎、根據相似理論發展而來的傳熱關系式是研究傳熱問題的重要手段。其典型代表即D-B關系式[19]與Petukhov關系式[20],D-B關系式(dittus-boelter correlation)為:

Nu=0.023ReaPrb

(1)

式中加熱流體時a=0.8,b=0.4。

原始型式的別祖霍夫(Petukhov)關系式如下:

(2)

ξ=(1.82lgRe-1.64)-2

K1(ξ)=1+3.4ξ,K2(Pr)=11.7+1.8Pr-1/3

D-B關系式以大量實驗值為基礎進行擬合,適合于流體普朗特數Pr≥0.6、流體與壁面具有中等溫度的換熱情況;Petukhov關系式基于常物性結合邊界層理論推導而來,能夠適應Pr≥0.5及變物性的情況,為簡化計算,K1、K2也可分別取常數1.07、12.7,此可視為原始關系式的修正式。

以式(1)～(2)為基礎,添加物性修正項的思想能夠使之適合不同的實驗或數值模擬條件。諸多研究表明,基礎物性的差距將導致流動換熱特性的差異,因此,有必要對氦氙混合氣體及超臨界二氧化碳的傳熱過程開展差異化分析。

對于冷卻劑工質,物性變化、流道結構是提出不同型式流動換熱關系式的主要考量因素。本研究立足小尺寸圓管內流動,對不同雷諾數下的氦氙混合氣體及超臨界二氧化碳在定常加熱情況下的流動進行廣泛計算,將計算所得數據點與現有公式進行比較,得到適合單一工質的較優流動換熱關系式。計算過程采用進出口平均溫度作為定性溫度,在Nu的計算中,對流換熱系數計算方法與第2節所述相同。

對于兩型冷卻劑工質,選擇5個具有代表性的經驗關系式,開展Nu的計算誤差比較。由圖10可知,當Re在1×104～5×104的范圍內時,Petukhov關系式的預測結果偏差小于20%,而采用常數修正的關系式預測精度相對更高。D-B關系式的一般型式與計算值偏差較大,但其偏離趨勢近似為線性,因此考慮對其系數進行調整,使之符合計算結果。圖10中,實心點為氦氙混合氣體;空心點為超臨界二氧化碳。

圖10 經驗關系式與模擬結果的Nu數比較Fig.10 Comparison of Nu number between correlations and simulation results

對于氦氙混合氣體,采用D-B關系式、調整系數進行修正的結果如圖11(a)所示,圖中實線表示的擬合曲線為:

圖11 氦氙混合氣體與超臨界二氧化碳D-B型公式擬合Fig.11 D-B type formula fitting of He-Xe mixture and S-CO2

Nu=0.34Re0.5Pr0.4

(3)

式中:Re∈[3×103,105];Pr∈[0.15,0.5]。

當Re>104時,修正后曲線與計算數據的符合程度良好,預測精度能夠達到95%以上。其原因與流動結構有關,在該計算模型尺寸下,Re在7×103～104內,Nu/Pr0.4隨Re變化的斜率出現明顯改變,因此,該段可被認為是氦氙混合氣體管內流動的轉捩區,當Re<7×103時,流動狀態為層流,與之對應,Re>104對應湍流的充分發展段。

超臨界二氧化碳的D-B關系式修正如圖11(b)所示,公式為:

Nu=0.1Re0.65Pr0.4

(4)

式中:Re∈[3×103,1.5×105];Pr>0.5。

需要指出,對于超臨界二氧化碳,研究指出浮升力在對流換熱過程中具有不可忽略的影響,而重要的判斷依據之一即Bo數[21]:

(5)

當Bo滿足式(5)時,浮升力的影響即可忽略不計。

本文的計算中,考慮了重力與流動方向相同、相反、垂直3種情況。計算結果表明,不論對于何種工況,Bo?5.6×10-7。對于系統壓力為17 MPa下的4 mm直徑管內流動,浮升力對換熱過程施加的影響微乎其微。因此,圖11(b)中未對3類計算數據加以區分。

對比2種冷卻劑的計算數據及修正關系式可知,氦氙混合氣體與超臨界二氧化碳均可通過修改經驗關系式中的系數使得流動符合D-B關系式的一般型式;而在湍流充分發展段,無論是原始型式還是修正型式,Petukhov關系式均能夠較好地描述流動換熱過程。需要指出,Petukhov關系式的傳統應用范圍要求Pr>0.5,在計算中可以發現,由于超臨界二氧化碳Pr約為0.74,符合應用條件,因此其計算誤差相對較小;氦氙混合氣體作為低Pr工質,Pr≈0.25,超出了公式設定的應用范圍,因此計算誤差相對較大,若使用Petukhov關系式則需對Pr做出修正。

4 結論

1)在高加熱功率、大入口流速時,氦氙混合氣體及超臨界二氧化碳具有換熱系數極大值,在氦氙混合氣體入口溫度1 000～1 200 K、超臨界二氧化碳600 K時有最大對流換熱系數,該參數適合作為反應堆設計時的參考。

2)以D-B關系式的一般型式修正了氦氙混合氣體與超臨界二氧化碳換熱模型,在高Re區域(Re>104),擬合模型與計算值符合良好。