熔鹽堆熱管式非能動余熱排出系統(tǒng)建模及程序開發(fā)

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

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

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

本文基于HP-PRHRS結構和運行特點，建立一套較為完整的數(shù)學物理模型，耦合熔鹽堆堆芯物理熱工模型、高溫熱管模型和非能動余熱排出系統(tǒng)模型等，開發(fā)熔鹽堆HP-PRHRS分析程序PRAC，并采用MSRE基準題和瞬態(tài)實驗數(shù)據(jù)對程序進行驗證。

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

1.1 堆芯物理熱工模型

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

(1)

(2)

(3)

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

(4)

(5)

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

(6)

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

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

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

(7)

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

1.2 高溫熱管模型

圖1 高溫熱管原理圖Fig.1 Schematic of high temperature heat pipe

熱管是一較為復雜的系統(tǒng)，原理圖如圖1所示。熱管依靠由吸液芯毛細力驅動的氣液兩相自然循環(huán)完成熱量傳遞，其工作過程包含了多種不同的傳熱形式和連續(xù)的相變過程，因此國際上尚未有準確獲得熱管內部瞬態(tài)特性分析解的理論模型。而實際應用關注的是高溫熱管的整體宏觀性能，主要是熱管傳熱過程中的功率和溫度等關鍵參數(shù)。基于上述理念，本文采用熱管熱阻網絡模型，進行高溫熱管的求解分析，大幅提升數(shù)值穩(wěn)定性和計算效率。基于熱管傳熱流程，可將熱管的傳熱過程以熱阻網絡的形式表示。經合理簡化后，最終的高溫熱管熱阻網絡模型如圖2所示。圖2中：熱阻1和2分別為蒸發(fā)段管壁徑向熱阻和吸液芯徑向熱阻；熱阻3和4分別為蒸發(fā)段和冷凝段的液環(huán)熱阻；熱阻5和6分別為冷凝段管壁徑向熱阻和吸液芯徑向熱阻；熱阻7和8分別為絕熱段管壁軸向熱阻和吸液芯軸向熱阻；Q為輸入功率，W；T∞為環(huán)境溫度，K。

圖2 熱管熱阻網絡模型Fig.2 Network model of high temperature heat pipe

對于網絡中的單個熱阻控制體，有如下控制方程：

(8)

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

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

(9)

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

n=2,3,4,5

(10)

(11)

(12)

(13)

(14)

(15)

(16)

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

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

(17)

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

(18)

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

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

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

非能動余熱排出系統(tǒng)模型分為卸料罐內燃料鹽流動換熱和卸料罐外排熱煙囪內空氣流動換熱兩部分，二者通過高溫熱管模型實現(xiàn)參數(shù)交互。

1) 卸料罐內燃料鹽流動換熱

卸料罐內燃料鹽控制方程為：

(19)

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

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

(20)

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

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

對于單根獨立的熱管[6]：

(21)

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

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

Nu=3.710 9Ra0.116 06

(22)

為使高溫熱管處于良好的工作狀態(tài)，熱管傾斜插入，罐外熱管冷凝段高于罐內蒸發(fā)段。對于傾斜管，浮升力(重力)可分解為平行和垂直于表面的兩方向作用力，因此圓管的長度和管徑對于自然對流都有影響[13]。因此對于傾斜圓管的自然對流傳熱，其特征長度需采用傾斜角進行修正：

(23)

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

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

空氣由排熱煙囪底部流入，經過卸料罐區(qū)域對熱管冷凝段進行冷卻，然后從頂部出口流出。由于空氣的流速遠低于聲速，且空氣壓縮性影響很小，因此采用一維不可壓縮模型進行空氣的流動和換熱計算。空氣流動換熱整體控制方程如下。

連續(xù)性方程：

wi=wair

(24)

動量守恒方程：

(25)

(26)

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

能量守恒方程：

(27)

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

空氣與熱管冷凝段間通過對流傳遞熱量。當空氣流動較弱、流速較低時，可用Corcione模型進行空氣外掠豎列管束的自然對流換熱計算[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ù)量。

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

(29)

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

2 PRAC程序開發(fā)

2.1 程序結構

HP-PRHRS分析程序PRAC采用標準FORTRAN 90程序設計語言編寫。基于模塊化設計思路，根據(jù)功能進行模塊劃分，使每個模塊間具有較高的獨立性。PRAC程序的模塊結構如圖3所示。其中堆芯物理熱工模塊由中子動力學模塊和流體動力學模塊耦合；余熱排出系統(tǒng)模塊由高溫熱管模塊、卸料罐模塊和排熱煙囪模塊3部分耦合；輔助模塊、物性模塊和數(shù)值計算模塊為公用模塊，可根據(jù)需要受到其他模塊調用。

2.2 程序流程

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

2.3 數(shù)值算法

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

圖3 PRAC程序模塊結構Fig.3 Module structure of PRAC code

3 程序及模型驗證

3.1 堆芯物理熱工模塊驗證

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

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

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

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

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

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

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

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

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

表2 瞬態(tài)實驗系統(tǒng)和高溫熱管的主要參數(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 結論

本文根據(jù)熔鹽堆堆芯和HP-PRHRS的結構及運行特點建立了一套合理完整的系統(tǒng)模型，主要由堆芯物理熱工模型、高溫熱管模型和非能動余熱排出系統(tǒng)模型3部分耦合。采用模塊化編程設計，開發(fā)了熔鹽堆HP-PRHRS分析程序PRAC。采用MSRE基準題和瞬態(tài)實驗數(shù)據(jù)進行對比驗證，PRAC計算結果與基準值吻合良好，且與燃料鹽溫度實驗值的相對偏差小于1.2%，與熱管壁溫實驗值的相對偏差小于2.8%，驗證了模型與程序的合理性與準確性。本文模型和程序能為后續(xù)開展熔鹽堆HP-PRHRS的深入設計提供模型和軟件基礎。