劉明皓，張玉龍，王成龍，張大林，秋穗正，李 毅，尹莎莎，劉 航，楊紅發(fā)

(1.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室，四川 成都 610041； 2.西安交通大學(xué) 能源與動力工程學(xué)院，陜西 西安 710049)

熔鹽堆作為第4代核反應(yīng)堆中唯一一種液態(tài)燃料堆型，已得到全世界范圍內(nèi)的廣泛關(guān)注[1-3]。福島核事故后，人們對核反應(yīng)堆系統(tǒng)提出了更高的非能動安全要求。為提升熔鹽堆非能動安全性，同時使系統(tǒng)滿足小型化、模塊化的設(shè)計需求，研究者結(jié)合高溫?zé)峁芗夹g(shù)提出了熔鹽堆熱管式非能動余熱排出系統(tǒng)(HP-PRHRS)概念設(shè)計，相關(guān)理論和實驗研究已先行展開，初步驗證了高溫?zé)峁苡糜谌埯}堆系統(tǒng)的可行性[4-8]。HP-PRHRS主要由卸料罐、排熱煙囪和高溫?zé)峁芙M成。當(dāng)反應(yīng)堆發(fā)生一回路大破口、失流等事故時，壓力容器中燃料鹽溫度迅速上升使冷凍閥熔斷，冷凍閥開啟，燃料鹽依靠重力作用快速下泄到卸料罐中，插在卸料罐中的高溫?zé)峁苎杆賳訉?dǎo)出其熱量。HP-PRHRS運(yùn)行由3個自然循環(huán)(對流)耦合完成：1) 卸料罐內(nèi)燃料鹽與高溫?zé)峁荛g的自然對流傳熱；2) 高溫?zé)峁軆?nèi)部工質(zhì)的氣液兩相自然循環(huán)；3) 排熱煙囪內(nèi)空氣的自然循環(huán)。HP-PRHRS的運(yùn)行不需額外的外界驅(qū)動力，僅依靠系統(tǒng)自身的自然循環(huán)(對流)即可實現(xiàn)燃料鹽有效冷卻，使系統(tǒng)具備良好的非能動安全特性。同時系統(tǒng)結(jié)構(gòu)得到大幅簡化，有效避免了中間環(huán)節(jié)故障，對未來熔鹽堆模塊化、小型化設(shè)計非常有利。

本文基于HP-PRHRS結(jié)構(gòu)和運(yùn)行特點，建立一套較為完整的數(shù)學(xué)物理模型，耦合熔鹽堆堆芯物理熱工模型、高溫?zé)峁苣Ｐ秃头悄軇佑酂崤懦鱿到y(tǒng)模型等，開發(fā)熔鹽堆HP-PRHRS分析程序PRAC，并采用MSRE基準(zhǔn)題和瞬態(tài)實驗數(shù)據(jù)對程序進(jìn)行驗證。

1 主要數(shù)學(xué)物理模型

1.1 堆芯物理熱工模型

堆芯物理熱工模型采用點堆模型與燃料鹽流動換熱模型耦合，將堆芯功率分解為時間和空間的函數(shù)，采用反應(yīng)堆中子動力學(xué)模型計算反應(yīng)堆瞬態(tài)功率，考慮了燃料多普勒效應(yīng)、熔鹽密度等反應(yīng)性反饋和燃料鹽流動帶來的緩發(fā)中子先驅(qū)核遷移機(jī)制影響。方程組如下：

(1)

(2)

(3)

燃料鹽在石墨通道內(nèi)流動，熱量由熔鹽帶出堆芯，同時燃料鹽同石墨間存在傳熱，進(jìn)而得出下列傳熱方程：

(4)

(5)

Pout(t)=w(t)cpf(Tout-Tin)

(6)

其中：Mf、Mg分別為燃料鹽和石墨的質(zhì)量，kg；cpf、cpg分別為燃料鹽和石墨的比定壓熱容，J·K-1·kg-1；Pf、Pg分別為核釋熱分別進(jìn)入燃料鹽和石墨的功率，W；Pout為燃料鹽流動傳出堆芯的功率，W；h為石墨通道與燃料鹽之間的對流換熱系數(shù)，W·K-1·m-2；A為石墨通道與燃料鹽的對流面積，m2；Tg、Tf分別為石墨和燃料鹽的平均溫度，K；Tout、Tin分別為堆芯出口和進(jìn)口溫度，K。

燃料鹽溫度和石墨溫度對反應(yīng)性的溫度反饋可按下式計算：

ρr,T(t)=αf(Tf-Tf0)+αg(Tg-Tg0)

(7)

其中：ρr，T(t)為溫度反應(yīng)性反饋；αf為燃料鹽溫度反饋系數(shù)，K-1；αg為石墨溫度反饋系數(shù)，K-1；Tg0、Tf0分別為初始時刻石墨和燃料鹽的平均溫度，K。

1.2 高溫?zé)峁苣Ｐ?/h3>

圖1 高溫?zé)峁茉韴DFig.1 Schematic of high temperature heat pipe

熱管是一較為復(fù)雜的系統(tǒng)，原理圖如圖1所示。熱管依靠由吸液芯毛細(xì)力驅(qū)動的氣液兩相自然循環(huán)完成熱量傳遞，其工作過程包含了多種不同的傳熱形式和連續(xù)的相變過程，因此國際上尚未有準(zhǔn)確獲得熱管內(nèi)部瞬態(tài)特性分析解的理論模型。而實際應(yīng)用關(guān)注的是高溫?zé)峁艿恼w宏觀性能，主要是熱管傳熱過程中的功率和溫度等關(guān)鍵參數(shù)?；谏鲜隼砟?，本文采用熱管熱阻網(wǎng)絡(luò)模型，進(jìn)行高溫?zé)峁艿那蠼夥治觯蠓嵘龜?shù)值穩(wěn)定性和計算效率。基于熱管傳熱流程，可將熱管的傳熱過程以熱阻網(wǎng)絡(luò)的形式表示。經(jīng)合理簡化后，最終的高溫?zé)峁軣嶙杈W(wǎng)絡(luò)模型如圖2所示。圖2中：熱阻1和2分別為蒸發(fā)段管壁徑向熱阻和吸液芯徑向熱阻；熱阻3和4分別為蒸發(fā)段和冷凝段的液環(huán)熱阻；熱阻5和6分別為冷凝段管壁徑向熱阻和吸液芯徑向熱阻；熱阻7和8分別為絕熱段管壁軸向熱阻和吸液芯軸向熱阻；Q為輸入功率，W；T∞為環(huán)境溫度，K。

圖2 熱管熱阻網(wǎng)絡(luò)模型Fig.2 Network model of high temperature heat pipe

對于網(wǎng)絡(luò)中的單個熱阻控制體，有如下控制方程：

(8)

其中：δi為熱阻導(dǎo)熱厚度，m；αi為熱擴(kuò)散系數(shù)，m2·s-1；Ti，1和Ti，2為熱阻兩側(cè)邊界溫度，K；Ti為熱阻中心溫度，K。

進(jìn)而聯(lián)立各熱阻單元，可獲得如下熱管熱阻網(wǎng)絡(luò)模型控制方程，通過求解可獲得熱管各區(qū)溫度，進(jìn)而計算出熱流密度、溫度變化率、熱管效率等參數(shù)。

(9)

εn(n+1)-2)Tn+ε(n+1)nTn+1)

n=2,3,4,5

(10)

(11)

(12)

(13)

(14)

(15)

(16)

其中：ki為熱導(dǎo)率，W·m-1·K-1；h∞，c為冷凝段與環(huán)境換熱系數(shù)，W·K-1·m-2；Sc為冷卻換熱面積，m2；T∞，c為環(huán)境冷卻溫度，K。

熱管的正常工作依賴其內(nèi)部工質(zhì)的兩相循環(huán)，因此需補(bǔ)充工質(zhì)循環(huán)判別模型[9]，Ψ>1表明熱管建立內(nèi)部工質(zhì)循環(huán)。

(17)

其中：Ψ為循環(huán)判別無量綱數(shù)；μv為蒸汽動力黏度，Pa·s；hfg為工質(zhì)汽化潛熱，J·kg-1；ρv為蒸汽工質(zhì)密度，kg·m-3；Rv為蒸汽腔當(dāng)量半徑，m；Le、La、Lc分別為熱管蒸發(fā)段、絕熱段、冷凝段長度，m；δ為吸液芯厚度，m；D為熱管直徑，m；THP為熱管運(yùn)行溫度，K；keff為吸液芯有效熱導(dǎo)率，W·m-1·K，考慮吸液芯材料熱導(dǎo)率kw和吸液芯內(nèi)工質(zhì)熱導(dǎo)率kl，keff按照下式計算：

(18)

其中，φ為吸液芯孔隙率。

然而在實際運(yùn)行中，熱管的傳熱能力受其傳熱極限的限制。當(dāng)對熱管的輸入功率超過熱管傳熱極限時，熱管內(nèi)部的循環(huán)便會遭到破壞，導(dǎo)致熱管失效。傳熱極限與熱管尺寸、工質(zhì)、吸液芯結(jié)構(gòu)、工作溫度等因素有關(guān)[10-11]。因此在熱管的計算分析中，需結(jié)合熱管的傳熱極限(主要是高溫?zé)峁?，尤其是堿金屬熱管，一般為毛細(xì)極限、聲速極限和攜帶極限等)，對熱管傳熱功率進(jìn)行判別校對，以確定熱管處于正常工作狀態(tài)。

1.3 非能動余熱排出系統(tǒng)模型

非能動余熱排出系統(tǒng)模型分為卸料罐內(nèi)燃料鹽流動換熱和卸料罐外排熱煙囪內(nèi)空氣流動換熱兩部分，二者通過高溫?zé)峁苣Ｐ蛯崿F(xiàn)參數(shù)交互。

1) 卸料罐內(nèi)燃料鹽流動換熱

卸料罐內(nèi)燃料鹽控制方程為：

(19)

其中：m為燃料鹽質(zhì)量，kg；Tsalt為燃料鹽溫度，K；τ為系統(tǒng)運(yùn)行時間，s；cp為燃料鹽比定壓熱容，J·K-1·kg-1；Qd為衰變熱功率，W；Qc為系統(tǒng)冷卻功率，W。

對于燃料鹽的衰變熱功率，若其衰變曲線未知，可采用Shure修正預(yù)測模型進(jìn)行燃料鹽衰變熱功率計算[12]。

(20)

其中：Q0為停堆前反應(yīng)堆的運(yùn)行功率，W；Qd為燃料鹽排出時間τ后的衰變熱功率，W；τ0為停堆前反應(yīng)堆已連續(xù)運(yùn)行時間，s；a和b為時間相關(guān)系數(shù)。

卸料罐內(nèi)燃料鹽自然循環(huán)。罐內(nèi)熱管蒸發(fā)段與燃料鹽間的自然對流傳熱可采用下列模型計算。

對于單根獨(dú)立的熱管[6]：

(21)

其中：Hhp為熱管高度，m；Hmax為最大液位高度，m。

對于熱管豎列管束，熱管間距P與管徑D比值為4.33≤P/D≤8.67時[7]：

Nu=3.710 9Ra0.116 06

(22)

為使高溫?zé)峁芴幱诹己玫墓ぷ鳡顟B(tài)，熱管傾斜插入，罐外熱管冷凝段高于罐內(nèi)蒸發(fā)段。對于傾斜管，浮升力(重力)可分解為平行和垂直于表面的兩方向作用力，因此圓管的長度和管徑對于自然對流都有影響[13]。因此對于傾斜圓管的自然對流傳熱，其特征長度需采用傾斜角進(jìn)行修正：

(23)

其中：Lc為傾斜管特征長度，m；L為圓管長度，m；d為圓管直徑，m；θ為圓管傾斜角。當(dāng)傾角θ=0°時即為水平管，此時Lc=d；當(dāng)θ=90°時即為豎直管，此時Lc=L。

2) 卸料罐外排熱煙囪內(nèi)空氣流動換熱

空氣由排熱煙囪底部流入，經(jīng)過卸料罐區(qū)域?qū)峁芾淠芜M(jìn)行冷卻，然后從頂部出口流出。由于空氣的流速遠(yuǎn)低于聲速，且空氣壓縮性影響很小，因此采用一維不可壓縮模型進(jìn)行空氣的流動和換熱計算?？諝饬鲃訐Q熱整體控制方程如下。

連續(xù)性方程：

wi=wair

(24)

動量守恒方程：

(25)

(26)

其中：wair、wi分別為空氣流量和第i個空氣流動換熱控制體內(nèi)流量，kg·s-1；Δpi為排熱煙囪內(nèi)第i個空氣流動換熱控制體內(nèi)的壓降，Pa；Δpstack為排熱煙囪內(nèi)總壓降，Pa；Δp0為空氣從排熱煙囪出口到入口高度所對應(yīng)的環(huán)境壓降，Pa；fi為傳熱區(qū)摩擦阻力系數(shù)；Kj為局部阻力系數(shù)；ρi為第i個控制體內(nèi)空氣密度；Ai為第i個控制體換熱面積；g為重力加速度；Δzi為高度差，m；De為當(dāng)量直徑，m。

能量守恒方程：

(27)

其中：Vi為排熱煙囪內(nèi)第i個空氣流動換熱控制體的容積，m3；Si為控制體i的換熱面積，m2；qi為控制體i換熱面上的熱流密度，W·m-2；hi-1、hi分別為控制體i的進(jìn)、出口焓，J·kg-1。

空氣與熱管冷凝段間通過對流傳遞熱量。當(dāng)空氣流動較弱、流速較低時，可用Corcione模型進(jìn)行空氣外掠豎列管束的自然對流換熱計算[14]。

Nu=Ra0.235(0.292ln((P/d)0.4·

N-0.2)+0.447)

(28)

其中：5×102≤Ra≤5×105，P/d≤10-lgRa；N為豎列管束中的管數(shù)量。

此時，空氣流動可視為外掠光滑圓管，阻力計算可采用如下關(guān)系式：

(29)

其中：Δp為管外空氣流動的壓力損失，Pa；Gmax為最大質(zhì)量流速，kg·m-2·s-1；ξ為考慮管束排列方式的修正系數(shù)；gc為重力加速度，m·s-2。

2 PRAC程序開發(fā)

2.1 程序結(jié)構(gòu)

HP-PRHRS分析程序PRAC采用標(biāo)準(zhǔn)FORTRAN 90程序設(shè)計語言編寫?；谀K化設(shè)計思路，根據(jù)功能進(jìn)行模塊劃分，使每個模塊間具有較高的獨(dú)立性。PRAC程序的模塊結(jié)構(gòu)如圖3所示。其中堆芯物理熱工模塊由中子動力學(xué)模塊和流體動力學(xué)模塊耦合；余熱排出系統(tǒng)模塊由高溫?zé)峁苣K、卸料罐模塊和排熱煙囪模塊3部分耦合；輔助模塊、物性模塊和數(shù)值計算模塊為公用模塊，可根據(jù)需要受到其他模塊調(diào)用。

2.2 程序流程

PRAC程序流程圖如圖4所示。程序邏輯與HP-PRHRS實際運(yùn)行流程一致，主要包括堆芯計算和非能動余熱排出系統(tǒng)計算兩大部分。首先依據(jù)輸入?yún)?shù)開展熔鹽堆堆芯的熱工物理耦合計算，獲得溫度、功率、流量和反應(yīng)性反饋等堆芯參數(shù)。在計算過程中程序進(jìn)行判斷，如果冷凍閥開啟，燃料鹽進(jìn)入余熱排出系統(tǒng)，余熱排出系統(tǒng)隨即開始運(yùn)作，然后開始非能動余熱排出系統(tǒng)計算，獲得相關(guān)系統(tǒng)參數(shù)。相對地，如果無需冷凍閥開啟，在堆芯模塊計算至設(shè)定時間時程序便會終止，余熱排出系統(tǒng)計算模塊不會啟動。

2.3 數(shù)值算法

PRAC程序在運(yùn)行求解中會同時涉及點堆中子動力學(xué)模型和系統(tǒng)熱工水力模型，結(jié)合其余輔助模型，構(gòu)成一套封閉的方程組。由于不同模型間時間量級差距較大，構(gòu)成的方程組具有很大的剛性。針對這一問題，程序采用吉爾(Gear)算法進(jìn)行方程組求解，確保數(shù)值計算的精確穩(wěn)定。對于熱管熱阻網(wǎng)絡(luò)模型建立的一階常微分方程組，則采用高精度單步算法龍格庫塔法(R-K法)，選取區(qū)間上若干點的斜率進(jìn)行加權(quán)平均，通過基于泰勒級數(shù)展開的待定系數(shù)法實現(xiàn)微分方程組的步進(jìn)求解。

圖3 PRAC程序模塊結(jié)構(gòu)Fig.3 Module structure of PRAC code

3 程序及模型驗證

3.1 堆芯物理熱工模塊驗證

本文采用MSRE的啟泵基準(zhǔn)題和停泵基準(zhǔn)題數(shù)據(jù)來驗證其堆芯物理熱工模塊的可靠性和準(zhǔn)確性[15]。MSRE堆芯主要參數(shù)列于表1。MSRE在啟泵、停泵變流量的過程中，通過調(diào)節(jié)控制棒插入的深度來保持堆芯臨界，并記錄控制棒的位置，通過計算控制棒移動所引入的反應(yīng)性就可得出熔鹽流量變化所導(dǎo)致的堆芯反應(yīng)性的變化。

MSRE啟泵和停泵基準(zhǔn)題的計算結(jié)果如圖5所示。由圖5a可看出：啟泵時反應(yīng)性降低，通過控制棒提升來補(bǔ)償損失的反應(yīng)性，10 s后流量達(dá)到穩(wěn)定，控制棒停止運(yùn)動，控制棒補(bǔ)償?shù)姆磻?yīng)性與熔鹽流動造成的反應(yīng)性損失平衡。由圖5b可看出：熔鹽流動造成的反應(yīng)性損失與實驗結(jié)果吻合較好；停泵時熔鹽流量減少，導(dǎo)致堆芯反應(yīng)性增加，控制棒插入堆芯來減少多余反應(yīng)性，20 s后回路中的流量降為0，且最終控制棒減少的反應(yīng)性與流量減少所增加的反應(yīng)性相平衡。兩個基準(zhǔn)題的反應(yīng)性計算值均與實驗數(shù)據(jù)符合較好，證明了程序與模型的準(zhǔn)確性和可靠性。

3.2 余熱排出系統(tǒng)模塊驗證

余熱排出系統(tǒng)模塊采用前期研究獲得的瞬態(tài)實驗數(shù)據(jù)進(jìn)行對比驗證[8]。瞬態(tài)實驗在熔鹽堆HP-PRHRS實驗系統(tǒng)上進(jìn)行，如圖6所示。瞬態(tài)實驗系統(tǒng)和高溫?zé)峁苤饕獏?shù)列于表2。卸料罐總高為0.92 m，內(nèi)徑為0.6 m，側(cè)面安裝6根高溫鉀熱管。卸料罐放置于排熱煙囪內(nèi)，煙囪高度在0～4 m內(nèi)可調(diào)。

圖4 PRAC程序流程圖Fig.4 Flow chart of PRAC code

表1 MSRE基準(zhǔn)題堆芯參數(shù)Table 1 Core parameter in MSRE benchmark

圖5 MSRE啟泵(a)和停泵(b)基準(zhǔn)題計算結(jié)果Fig.5 Benchmark calculation result of MSRE pump strat-up (a) and pump stop (b)

本文針對系統(tǒng)瞬態(tài)實驗中的不同工況進(jìn)行了計算，將程序計算結(jié)果與實驗數(shù)據(jù)進(jìn)行對比，如圖7所示。對比實驗值與PRAC計算值，在不同實驗工況下PRAC計算值與實驗值均吻合良好。PRAC燃料鹽溫度計算值與實驗值相對偏差小于1.2%，熱管壁溫計算相對偏差小于2.8%。PRAC計算結(jié)果能準(zhǔn)確反映出系統(tǒng)冷卻過程中氟鹽溫度的變化情況和熱管壁溫的變化情況，證明了程序與模型的準(zhǔn)確性和可靠性。

圖6 熔鹽堆HP-PRHRS實驗系統(tǒng)Fig.6 HP-PRHRS experimental system for MSR

表2 瞬態(tài)實驗系統(tǒng)和高溫?zé)峁艿闹饕獏?shù)Table 2 Main parameter of transient experimental system and high temperature heat pipe

圖7 瞬態(tài)實驗燃料鹽溫度和熱管壁溫對比Fig.7 Comparison salt temperature and heat pipe temperature for transient experiment

4 結(jié)論

本文根據(jù)熔鹽堆堆芯和HP-PRHRS的結(jié)構(gòu)及運(yùn)行特點建立了一套合理完整的系統(tǒng)模型，主要由堆芯物理熱工模型、高溫?zé)峁苣Ｐ秃头悄軇佑酂崤懦鱿到y(tǒng)模型3部分耦合。采用模塊化編程設(shè)計，開發(fā)了熔鹽堆HP-PRHRS分析程序PRAC。采用MSRE基準(zhǔn)題和瞬態(tài)實驗數(shù)據(jù)進(jìn)行對比驗證，PRAC計算結(jié)果與基準(zhǔn)值吻合良好，且與燃料鹽溫度實驗值的相對偏差小于1.2%，與熱管壁溫實驗值的相對偏差小于2.8%，驗證了模型與程序的合理性與準(zhǔn)確性。本文模型和程序能為后續(xù)開展熔鹽堆HP-PRHRS的深入設(shè)計提供模型和軟件基礎(chǔ)。