周 彪，孫 倩，孫 俊，孫玉良

(清華大學(xué) 核能與新能源技術(shù)研究院，先進反應(yīng)堆工程與安全教育部重點實驗室，北京 100084)

高效可靠的空間電源是空間任務(wù)順利開展的重要保障。隨著未來軍民空間任務(wù)的不斷拓展，以化學(xué)能、太陽能為主的常規(guī)能源將難以滿足空間任務(wù)需求[1]?？臻g核反應(yīng)堆電源具備工作壽命長、功率覆蓋范圍廣、受空間環(huán)境影響小等特點，是未來開展多類別空間任務(wù)的理想電源形式[2]。自20世紀(jì)50年代，美國和俄羅斯就啟動了空間核反應(yīng)堆電源的相關(guān)研究工作，提出了多種技術(shù)路線。當(dāng)系統(tǒng)功率達(dá)到百千瓦乃至兆瓦量級時，采用氣冷堆結(jié)合布雷頓循環(huán)技術(shù)方案的優(yōu)勢更加明顯，主要在于其能量轉(zhuǎn)換效率高、技術(shù)成熟度高、系統(tǒng)比質(zhì)量更小。為滿足空間布雷頓循環(huán)系統(tǒng)緊湊性的設(shè)計要求，降低葉輪機械氣動載荷，同時獲得相對更佳的換熱性能，通常選擇氦氙混合氣體(He-Xe)作為反應(yīng)堆冷卻劑和布雷頓循環(huán)工質(zhì)[3-5]。

反應(yīng)堆熱工系統(tǒng)分析程序是開展反應(yīng)堆設(shè)計與安全評價的重要工具。針對不同類型的堆型，國內(nèi)外研究機構(gòu)開發(fā)了多種熱工系統(tǒng)分析程序，典型輕水堆系統(tǒng)分析程序有RELAP、CATHARE、TRACE等[6]；液態(tài)金屬反應(yīng)堆系統(tǒng)分析程序有SAS4A/SASSY-1、SIMMER、RELAP5-3D/ATHENA等[7]。盡管有THERMIX、TINTE等氣冷堆熱工系統(tǒng)分析程序[8-9]，但這些程序只能用于球床式反應(yīng)堆的熱工計算，應(yīng)用范圍有限。針對氦氙氣冷反應(yīng)堆，李楊柳等[10]開發(fā)了用于環(huán)形與圓形流道的反應(yīng)堆單通道分析程序，但無法滿足系統(tǒng)層面的熱工分析。美國愛達(dá)荷國家實驗室開發(fā)了RELAP5-3D/ATHENA程序[11]，可拓展用于超臨界CO2、He-Xe等工質(zhì)的高溫氣冷堆熱工系統(tǒng)分析，但該程序?qū)鴥?nèi)不開放。為更好地對氦氙氣冷空間核反應(yīng)堆系統(tǒng)開展熱工分析與安全評價，開發(fā)相應(yīng)熱工系統(tǒng)分析程序十分必要。

考慮到RELAP5程序的物理模型和計算方法較為成熟，且廣泛應(yīng)用于核動力系統(tǒng)安全評審，本文選擇對RELAP5/MOD4.0程序進行拓展，基于程序源碼添加He-Xe物性計算模塊與適用的傳熱關(guān)系式，使程序初步具備計算He-Xe流動換熱的功能。

1 RELAP5/MOD4.0概述

RELAP5是由美國愛達(dá)荷國家實驗室與美國核管會協(xié)同開發(fā)的一維瞬態(tài)熱工分析程序，RELAP5/MOD4.0為當(dāng)前最新的可用版本。程序采用自上而下的編程邏輯，如圖1所示，頂層主要由INPUT、TRNCTL、STRIP 3個功能區(qū)組成，各功能區(qū)再分別由相應(yīng)子程序執(zhí)行特定功能。其中INPUT控制輸入卡數(shù)據(jù)讀取與計算初始化；TRNCTL控制水力學(xué)部件與熱構(gòu)件的步進計算與結(jié)果保存；STRIP提供數(shù)據(jù)提取、繪圖等功能。本文基于RELAP5/MOD4.0開發(fā)He-Xe流動換熱計算模塊，編寫He-Xe密度、動力黏度、熱導(dǎo)率、比熱容4個物性計算子程序，在RNEWP與HYDRO路徑文件中添加物性子程序的調(diào)用接口，實現(xiàn)對He-Xe物性的計算。在HTADV路徑下的dittus.for文件中新增適用的傳熱關(guān)系式，以滿足He-Xe對流換熱的計算需求。

圖1 RELAP5/MOD4.0程序框架

2 程序功能拓展

2.1 He-Xe物性計算模塊

RELAP5/MOD4.0程序中不凝結(jié)氣體密度ρ與比定壓熱容cp采用理想氣體狀態(tài)方程進行求解。El-Genk等[5]研究表明，該方法求解稠密He-Xe物性時將引入較大誤差。不凝結(jié)氣體動力黏度μ與熱導(dǎo)率λ分別采用式(1)、(2)進行計算，這些公式均是基于Sutherlands定律演變而來[12]。Sutherlands定律是基于理想氣體分子動力理論與分子間勢能理論提出的，在空氣黏度等物性計算方面具有較好的適用性。

(1)

λ=λrTb

(2)

式中：S為Sutherlands溫度，程序中S=114.0 K；μr、λr分別為參考溫度Tr(Tr=273.13 K)下混合氣體的動力黏度與熱導(dǎo)率；b為與氣體類型相關(guān)的常數(shù)；T為氣體實際溫度。用βr代表混合氣體μr、λr、b，三者的數(shù)值均可由通式(3)計算得到。式(3)中下標(biāo)1表示氙氣，2表示氦氣，x1為氙氣摩爾分?jǐn)?shù)。純氣體相關(guān)常數(shù)在RELAP5用戶手冊[13]中有詳細(xì)說明，其數(shù)值列于表1。

表1 程序中計算物性的常數(shù)

βr=x1βr1+(1-x1)βr2

(3)

表2 混合稀有氣體物性的半經(jīng)驗公式

基于Tournier等的物性公式，本文利用MATLAB開發(fā)了He-Xe物性計算程序PROHX，并在之前的研究[15]中驗證了程序計算的準(zhǔn)確性。為探究式(1)、(2)對He-Xe物性計算的適用性，將其計算值與PROHX程序計算結(jié)果進行對比，如圖2所示，發(fā)現(xiàn)式(1)、(2)計算值與PROHX計算值差異明顯，表明Sutherlands定律不適用于He-Xe物性計算。因此，需在程序中添加新的物性模塊。本研究在源碼中新增了rhohx.for(密度)、viscshx.for(黏度)、thcnhx.for(熱導(dǎo)率)、cpphx.for(比熱容)4個物性子程序文件以及hexeconstant.for庫文件(定義物性計算時涉及的變量與相關(guān)函數(shù))，構(gòu)成完整的氦氙物性計算模塊。上述物性均為溫度、壓力、混合比相關(guān)的函數(shù)。通過輸入卡給定相應(yīng)工況參數(shù)后，在程序初始化和步進計算過程中調(diào)用相應(yīng)物性子程序，可實現(xiàn)對氦氙混合氣體物性的計算。

圖2 氦氙物性計算結(jié)果對比

2.2 He-Xe傳熱關(guān)系式

RELAP5程序中強迫對流換熱計算關(guān)系式采用Dittus-Bolter(DB)公式。為探究DB公式對He-Xe流動換熱計算的適用性，采用DB公式計算了不同Pr(空間堆推薦使用的He-Xe對應(yīng)Pr為0.16

圖3 DB公式計算值與實驗值對比

圖3中實驗數(shù)據(jù)參考Taylor等[4]的圓管加熱實驗，實驗裝置示意圖如圖4所示，其中實驗測試段為直徑D=5.87 mm的圓管，由絕熱段(l1=56D)和均勻熱流加熱段(l2=60D)組成，該實驗探究了包括He-Xe在內(nèi)的多種混合氣體的流動換熱過程。

圖4 實驗裝置示意圖

DB公式為常物性公式，而氣體流經(jīng)加熱流道時，截面徑向位置處的物性差異明顯，此時氣體變物性對流動換熱的影響不可忽略。如表3所列，本研究首先選取4組變物性傳熱公式[16-18]進行適用性分析。式中：Nub為變物性Nu；Nu0為常物性Nu；Tw為壁面溫度；Tb為He-Xe平均溫度；n為變物性修正因子；ζ為中間變量。為覆蓋目前空間堆設(shè)計方案中使用的工質(zhì)范圍，本文選擇4個不同摩爾質(zhì)量混合氣體的實驗組[4]，如表4所列，工況范圍為：0.2130時，Petukhov公式與實驗值符合較好，原因是其未考慮熱入口段的影響(即x/D的影響)。綜合考慮，本研究采用Taylor公式對程序原有傳熱關(guān)系式進行修正。

圖5 現(xiàn)有Nu關(guān)系式計算值與實驗值的對比

表3 現(xiàn)有的Nu關(guān)系式

表4 氦氙混合氣體的實驗參數(shù)

3 程序功能驗證與結(jié)果分析

為驗證拓展后的RELAP5程序?qū)e-Xe流動換熱計算的準(zhǔn)確性，采用拓展后的程序?qū)ι鲜鰧嶒瀃4]的測試段進行建模。為探究程序拓展后的改進效果，引入程序默認(rèn)的DB公式作對比。計算的邊界條件與實驗工況一致，將不同工況下的程序計算值與實驗結(jié)果進行對比，結(jié)果如圖6所示，發(fā)現(xiàn)程序中傳熱關(guān)系式的選擇對壁面溫度Tw的計算值有重要影響。采用DB公式計算得到的加熱段Tw顯著低于實驗結(jié)果；采用Taylor公式得到的計算值與實驗值符合較好；當(dāng)Tw/Tb較大時，在接近熱充分發(fā)展段處的Tw計算值略高于實驗值，此時計算結(jié)果具有一定的保守性。

圖6 RELAP5程序Tw計算值與實驗值的對比

此外，計算了He-Xe平均溫度Tb，如圖7所示。由圖7可看出，采用DB公式與Taylor公式的計算結(jié)果基本重合，且各工況下Tb計算值與實驗值符合較好。這是因為RELAP5求解流體平均溫度是直接通過能量守恒法來實現(xiàn)的，只要邊界條件確定，各控制體的平均溫度則確定，與傳熱關(guān)系式的選擇無關(guān)。在流動阻力方面，計算了兩組工況下第2個壓力測點處的壓力以及從起始加熱點至第2個壓力測點之間的壓降，結(jié)果列于表5。由表5可看出，兩個工況的壓降計算相對誤差均在10%以內(nèi)。通過以上計算，驗證了拓展后的RELAP5程序在計算He-Xe流動換熱功能上的準(zhǔn)確性。

圖7 RELAP程序Tb計算值與實驗值的對比

表5 拓展后的RELAP5程序計算He-Xe流動壓降的誤差

4 結(jié)論

為開發(fā)適用于氦氙氣冷空間堆的熱工系統(tǒng)分析程序，本文對RELAP5/MOD4.0程序開展了適用性分析，并對源碼進行了二次開發(fā)。研究表明：程序默認(rèn)的Sutherlands定律用于He-Xe物性計算時將引入較大誤差；Dittus-Bolter公式對He-Xe對流換熱時的Nu預(yù)測偏高，壁溫計算結(jié)果不保守。在RELAP5源碼實現(xiàn)了He-Xe物性計算模塊和傳熱關(guān)系式的添加，對程序功能進行了拓展。通過與實驗值的對比，發(fā)現(xiàn)拓展后的程序?qū)e-Xe流動換熱模擬結(jié)果與實驗值吻合較好，驗證了程序?qū)e-Xe流動換熱計算的準(zhǔn)確性。在后續(xù)研究中，將開發(fā)布雷頓循環(huán)系統(tǒng)中相關(guān)部件的計算模塊，基于具體空間堆設(shè)計方案，搭建系統(tǒng)計算模型，從而實現(xiàn)對系統(tǒng)層面的瞬態(tài)及穩(wěn)態(tài)計算分析。