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

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

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

反應(yīng)堆熱工系統(tǒng)分析程序是開展反應(yīng)堆設(shè)計(jì)與安全評(píng)價(jià)的重要工具。針對(duì)不同類型的堆型，國內(nèi)外研究機(jī)構(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)堆的熱工計(jì)算，應(yīng)用范圍有限。針對(duì)氦氙氣冷反應(yīng)堆，李楊柳等[10]開發(fā)了用于環(huán)形與圓形流道的反應(yīng)堆單通道分析程序，但無法滿足系統(tǒng)層面的熱工分析。美國愛達(dá)荷國家實(shí)驗(yàn)室開發(fā)了RELAP5-3D/ATHENA程序[11]，可拓展用于超臨界CO2、He-Xe等工質(zhì)的高溫氣冷堆熱工系統(tǒng)分析，但該程序?qū)鴥?nèi)不開放。為更好地對(duì)氦氙氣冷空間核反應(yīng)堆系統(tǒng)開展熱工分析與安全評(píng)價(jià)，開發(fā)相應(yīng)熱工系統(tǒng)分析程序十分必要。

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

1 RELAP5/MOD4.0概述

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

圖1 RELAP5/MOD4.0程序框架

2 程序功能拓展

2.1 He-Xe物性計(jì)算模塊

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

(1)

λ=λrTb

(2)

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

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

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

(3)

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

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

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

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

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

圖3 DB公式計(jì)算值與實(shí)驗(yàn)值對(duì)比

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

圖4 實(shí)驗(yàn)裝置示意圖

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

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

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

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

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

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

圖6 RELAP5程序Tw計(jì)算值與實(shí)驗(yàn)值的對(duì)比

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

圖7 RELAP程序Tb計(jì)算值與實(shí)驗(yàn)值的對(duì)比

表5 拓展后的RELAP5程序計(jì)算He-Xe流動(dòng)壓降的誤差

4 結(jié)論

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