液態金屬冷卻快堆子通道分析軟件SACOS-LMR研發與工程應用

王金順,陳榮華,朱昕陽,田家豪,田文喜,秋穗正,蘇光輝

(西安交通大學 核科學與技術學院,陜西 西安 710049)

我國已確立了壓水堆-快堆-聚變堆三步走的核能發展戰略,快堆在這一戰略路線中發揮著承上啟下的關鍵作用,以鈉冷、鉛鉍或鉛冷為代表的液態金屬冷卻快堆則是其中最主流的堆型。在核反應堆設計中,堆芯的熱工水力特性分析是反應堆安全的前提,為提高堆芯的安全性能,要求盡可能精確地計算出堆芯的熱工參數。子通道分析方法目前是反應堆堆芯詳細熱工水力分析的一種相對精確的計算方法。考慮到液態金屬快堆堆芯的結構特點,有必要在壓水堆子通道分析程序的基礎上進行改進。韓國原子能研究所基于COBRA-Ⅳ-Ⅰ和MATRA程序,開發了適用于鈉冷快堆帶繞絲燃料組件的穩態及瞬態分析子通道程序MATRA-LMR,并利用所開發程序對HYPER堆芯燃料棒束進行了熱工水力分析[1-2]。上海交通大學[3]在COBRA的基礎上開發了COBRA-LM,可用于液態金屬反應堆堆芯子通道計算分析,程序采用ORNL19棒束實驗進行了驗證。中國科學技術大學[4]自主開發了鉛基堆子通道分析程序KMC-Sub,程序通過ORNL19棒束實驗和KYLIN臺架的61棒束實驗完成了驗證。博洛尼亞大學Lodi等[5]開發了適用于液態金屬快堆熱工安全分析的子通道程序ANTEO+,可用于堆芯組件的穩態計算分析。

西安交通大學熱工水力研究室在液態金屬快堆熱工水力特性分析方面開展了大量研究,積累了大量的研究經驗和軟件開發基礎[6-12],本文基于自主開發的子通道程序SACOS[13],通過添加液態金屬快堆特有的模型,如繞絲模型、盒間流模型、液態金屬對流換熱模型等,采用SIMPLE算法和交錯網格技術完成子通道計算模型的求解,形成適用于液態金屬快堆的子通道熱工安全分析程序SACOS-LMR;為了驗證程序的功能和計算精度,對德國卡爾斯魯厄核研究中心開展的37棒束鈉冷瞬態實驗[14]的穩態工況和沸騰前瞬態工況進行對比計算;基于驗證后的SACOS-LMR子通道程序對歐洲鉛冷快堆(ALFRED)[15]開展穩態工況和瞬態事故工況下的熱工安全特性分析,以證明SACOS-LMR子通道程序在液態金屬快堆的堆芯設計和熱工水力計算方面的能力。

1 數值方法

1.1 控制方程

在液態金屬冷卻快堆堆芯內不大可能發生沸騰現象,因此本研究的模型中沒有考慮鉛鉍的沸騰導致的兩相流動和換熱。對于單相液態金屬,其可壓縮性很小,可看成是不可壓縮流體,忽略其中的重力、壓力以及動能變化所做的功,并忽略軸向導熱,假設流體無內熱源,并且假設相鄰子通道之間的湍流交混不產生質量交換,則可得到以下熱工水力分析模型,依次為質量守恒方程、能量守恒方程、軸向動量方程和橫向動量方程。

(1)

(2)

(3)

(4)

式中:A為流通面積,m2;ρ為冷卻劑密度,kg/m3;t為時間,s;m為軸向質量流量,kg/s;z為軸向高度,m;wij為橫向質量流速,kg/(m·s);h為焓值,J/kg;w′ij為湍流質量流速,kg/(m·s);h*為i與j通道間隙處焓值,J/kg;GD為熱導因子;λ為熱導率,W/(m·K);s為間隙寬度,m;l為湍流長度,m;T為通道流體溫度,K;Qrad_con為相鄰通道流體徑向導熱,W;u為軸向速度,m/s;p為壓力,Pa;g為重力加速度;α為流動方向與豎直方向夾角,rad;f為摩擦阻力系數;De為子通道i的水力當量直徑,m;Ks為局部形阻系數;Kg為橫向流動阻力系數;下標i和j為通道編號,N為通道i相鄰的子通道總數。圖1為子通道劃分示意圖。

圖1 子通道劃分示意圖Fig.1 Schematic diagram of sub-channel division

1.2 燃料棒導熱模型

SACOS-LMR程序可對燃料棒或加熱棒進行建模,并可考慮棒在徑向、周向的導熱:

(5)

式中:cp為比定壓熱容;λ為熱導率;r為半徑;φ為周向方向上的角度,考慮到周向控制體數量只能和相鄰子通道個數保持一致,因此無法按照周向角度精細劃分,燃料棒內部導熱計算節點劃分如圖2所示;qV為單位體積熱源。燃料棒包殼表面溫度T采用第三類邊界條件計算:

圖2 燃料棒導熱控制體劃分Fig.2 Heat conduction node in fuel rod

(6)

其中:ht為燃料棒壁面換熱系數,W/(m2·K);下標w和f分別為冷卻劑和壁面。

1.3 輔助模型

1) 繞絲模型

為了封閉上述橫向動量方程(式(3)),需要對沿程摩擦阻力系數fi和局部形阻系數Ksi進行計算。液態金屬快堆組件常采用繞絲定位棒束,如圖3所示。繞絲的存在將大大影響冷卻劑的流量和溫度分布,一方面,棒束總壓降增大,另一方面,繞絲改變了流體在包殼表面的流動方向,增大了通道間的橫流交混。因此,繞絲模型在子通道分析中十分重要。Novendstern[16]結合以前所有的實驗方法和實驗數據,提出了適用于六邊形繞絲棒束的摩擦阻力系數計算模型。本研究采用Novendstern模型或UCTD模型[17]來計算fi,而由于局部形阻系數Ksi過于依賴經驗,常通過用戶根據局部形阻特征給定,在此不予考慮。在Novendstern模型中,通過有效摩擦系數feffect來表征繞絲帶來的影響,在光管阻力系數fsmooth上乘以倍增系數M得到帶繞絲燃料棒的有效阻力系數,摩擦阻力系數fi由下式確定:

圖3 繞絲定位的六邊形組件Fig.3 Hexagonal assembly with helical wire-wrap spacer

(7)

(8)

式中:L為燃料元件的長度,m;Vi為通道i的冷卻劑流速,m/s;M為倍增系數,主要由節徑比(P/D)、柵徑比(H/D)以及雷諾數確定。

(9)

Novendstern模型適用范圍為:棒數19～217,P/D范圍1.06～1.42,H/D范圍8.0～96.0,且適用于過渡區和紊流區,雷諾數范圍2 600～105。

UCTD模型是Chen等[17]在CTD模型的基礎上改進得到的,優化了層流雷諾數的預測精度并修正了組件壓降隨棒數變化趨勢異常的問題。相比于Novendstern模型,UCTD模型在層流湍流及過渡流做了更細致的區分和計算,相對更加精確。

(10)

(11)

式中,CfiL和CfiT為湍流常數,與燃料棒和繞絲的幾何參數有關。

2) 換熱系數模型

國內外關于液態金屬冷卻的棒束通道內的流動換熱系數模型研究得較多,其中 Mikityuk模型[18]的適用范圍更廣,平均偏差低于5%,且開發年代較近,模型如下:

Nu=0.047(1-e-3.8(P/D-1))(Pe0.77+250)

(12)

式中,Pe為佩克萊數,表示對流與擴散之比。Mikityuk模型適用范圍為1.1

此外,程序中還內置了Kazami和Carelli開發的Westinghouse換熱關系式[19],即:

(13)

3) 湍流交混模型

相鄰子通道間冷卻劑在流動過程中存在著橫向質量、動量和能量的交換,這種橫向的交換統稱為交混。由于交混的存在,各通道內的冷卻劑流量沿軸向將不斷發生變化,熱通道內冷卻劑的溫度和焓會有所降低。因此,準確模擬冷卻劑交混是子通道分析的關鍵之一。交混主要分為湍流交混、橫向流動、流動散射和流動掃掠4種形式。橫向流動是由于子通道間的橫向壓差導致的定向流動,通過式(4)右邊第一項橫向壓差作用力項考慮。流動散射和流動掃掠分別是由定位件引起的非定向交混和定向橫流,現有研究還不能夠精細描述機理規律,因此在SACOS-LMR軟件中將定位件的影響疊加在湍流交混效應中考慮。湍流交混本質上是子通道間流體脈動時自然渦團擴散引起的非定向交混。方程(2)～(4)中的湍流交混項w′ij為:

(14)

β=0.007 479Re-0.1

(15)

4) 液態金屬對流換熱模型

本文開發的液態金屬子通道分析程序適用于液態鈉、鉛鉍和鉛作為冷卻劑的快堆。世界經濟合作組織核能署(OECD/NEA)在2015年出版的鉛鉍手冊中對前人的一些實驗結果進行了整理和分析,總結出一套液態鉛鉍合金的物性模型[21],如表1所列。Fink等[22]在1995年發布的鈉物性報告中,推薦了液態鈉的熱物性,如表2所列。

表1 液態鉛鉍和鉛的物性(397.9～2 021 K)[21]Table 1 Property of liquid LBE/Pb (397.9-2 021 K)[21]

表2 液態鈉的物性[22]Table 2 Property of liquid sodium[22]

1.4 求解策略

相對于傳統子通道軟件(如COBRA-EN)采取的牛頓辛普森算法,采用自下而上的方式對整個流體域進行逐層求解并對出口流量進行修正,在求解中低流量時收斂性較好,在求解倒流、回流和極低流量時收斂效果不佳,而SACOS-LMR采用的SIMPLE算法,對整個軸向流體域建立全局矩陣,同時迭代求解,能夠對局部倒流不收斂的情況進行優化,具體求解流程如圖4所示。

圖4 SACOS-LMR程序求解流程Fig.4 Solution process of SACOS-LMR code

2 程序驗證

2.1 37棒束鈉冷瞬態實驗[14]

37棒束鈉冷瞬態實驗是在德國卡爾斯魯厄核研究中心的KNS鈉回路上進行的37棒束失流實驗。圖5為KNS 37棒實驗段及功率分布。參考SNR-300堆芯組件,棒的排列方式為六邊形組件中的三角形排布,采用電加熱棒進行加熱,熱功率分布類似于正弦分布,并由下式計算:

圖5 KNS 37棒實驗段及功率分布[14]Fig.5 KNS 37-pin experiment section and power distribution[14]

(16)

其中:x為單根棒的線功率,W/cm;z為加熱段底部位置的距離,mm;xmax為單根棒的最大線功率,W/cm。

實驗段由200 mm的零功率入口段和900 mm的加熱段以及450 mm的零功率出口段組成。KNS 37棒束的其他幾何參數列于表3。

表3 KNS 37棒實驗段幾何參數Table 3 Geometric parameters of KNS 37

KNS 37棒鈉冷實驗共進行了3種形式的失流實驗(L22、L26和L29),實驗中通過控制泵轉速和閥門來控制實驗段的入口流量。本文選取L22實驗進行程序驗證,分為兩個工況,其具體參數列于表4。對于沸騰前工況的計算,考慮到SACOS-LMR程序目前還未擴展至汽液兩相,因此隨著入口流量不斷減小至鈉沸騰前即停止計算。

表4 KNS 37棒實驗邊界條件Table 4 Boundary conditions of KNS 37

2.2 計算結果

基于上述實驗參數建立KNS 37棒L22實驗穩態工況和沸騰前工況的SACOS-LMR計算模型,計算中考慮到實驗段向周圍環境的熱損失,推薦使用6.7 W/(m2·K)的換熱系數以及40 ℃的環境溫度來計算實驗段的熱損失。

圖6為穩態工況和鈉沸騰前工況中心子通道鈉溫度的軸向變化。可看到,穩態工況和沸騰前工況下程序計算所得中心鈉溫度與實驗值趨勢符合較好,穩態工況最大偏差在20 ℃以內(為進出口溫升的12.2%),沸騰前工況最大偏差在25 ℃以內(為進出口溫升的6.8%)。兩種工況下,在零功率出口段附近的冷卻劑溫度偏差較大,且相較于實驗值都普遍預測偏低。由于金屬堆溫度普遍較高,若采用傳統方法計算溫度誤差,則誤差數值過小難以反映計算的準確性,因此采用相對偏差(RD)表示,如式(17)所示:

圖6 穩態工況和沸騰前工況中心通道鈉溫度軸向變化Fig.6 Variation of axial sodium temperature in central channel of steady case and before boiling

(17)

圖7為穩態工況下壓力沿軸向的變化,可看到程序具備計算格架定位件阻力系數的功能,計算的進出口壓降與實驗值符合較好。圖8為沸騰前工況距離加熱段底部779 mm高度處徑向(90°～270°方向)的溫度分布。可看到,最大偏差在45 ℃以內(為進出口溫升的12.2%),邊通道的偏差最大。對于實驗中90°和270°方向的對稱子通道內的溫度出現不對稱傾斜的情況,可認為是實驗段制造工差引起的。

圖7 穩態工況壓力軸向變化Fig.7 Variation of axial pressure of steady case

圖8 沸騰前工況下1 029 mm截面上子通道溫度分布Fig.8 Temperature distribution at 1 029 mm before boiling

SACOS-LMR對KNS 37棒實驗穩態和瞬態過程的計算結果可表明,程序具備穩態和瞬態的計算功能,且溫度場的相對偏差在15%以內,程序對進出口壓降的預測精度較高,能夠捕捉流量等瞬態過程的變化。

3 程序應用

ALFRED(Advanced Lead Fast Reactor European Demonstrator)是一個歐洲的核反應堆項目。ALFRED反應堆目前處于設計和開發階段,預計未來幾年開始建設。ALFRED反應堆的堆芯由171個燃料組件(FA)組成,根據流量和功率分為4組,如圖9所示,而燃料組件的幾何和物理參數列于表5。

表5 燃料組件的幾何和物理參數Table 5 Geometrical and thermal-hydraulic parameters of ALFRED FA

圖9 ALFRED燃料組件分組[23]Fig.9 Fuel assemble groups in core of ALFRED[23]

基于上述驗證可認為SACOS-LMR具備對液態金屬快堆堆芯進行穩態和瞬態計算與分析的能力。因此采用SACOS-LMR程序對ALFRED堆芯進行熱工水力計算,用于評估其熱工水力特性。

圖10為中央燃料組件活性區出口截面的溫度云圖。可看到,內通道的溫度最高且分布相對均勻,約為490 ℃。邊緣通道的溫度最低,約為473 ℃,而6個角通道的溫度較高,接近483 ℃。最熱的燃料棒為中心燃料棒,這與Petrovich等[23]的計算結果一致。在高度方向上,最熱燃料棒的燃料包殼外表面、包殼內表面、燃料芯塊表面、燃料芯塊中心以及鉛冷卻劑的溫度變化示于圖11。可看到,燃料芯塊的最高溫度出現在0.335 m位置,為1 832.3 ℃,低于2 000 ℃的設計限值。包殼的最高溫度出現在出口處,為544.4 ℃,低于550 ℃的設計限值。

圖10 冷卻劑溫度及最熱棒表面溫度分布Fig.10 Coolant temperature and circumferential temperature of hottest rod

圖11 中心棒軸向溫度變化Fig.11 Axial variation of temperatures in central fuel rod

表6為Petrovich[23]采用ANTEO-LFR子通道軟件與SACOS-LMR軟件計算的最高溫度的對比,二者相對偏差在15%以內。圖12為中心燃料棒在0.335 m處的內部溫度云圖。可看到,氣隙熱阻造成了約210 ℃的溫差,芯塊表面至芯塊中心溫差達到1 000 ℃。

表6 中心組件燃料棒最高溫度對比Table 6 Comparison of maximum temperature of central FA

圖12 中心燃料棒內部溫度分布Fig.12 Temperature distribution inside central fuel rod

圖13為中央燃料組件的質量流量和流速分布。可看到,角通道的流量相對較低,導致局部溫度較高。盡管包殼溫度低于設定的限值,但仍存在一定的安全隱患。可考慮調整幾何參數,即P/W和P/D,從而優化流量分布,在堆芯中實現更均勻的溫度分布。

圖13 出口截面流速和質量流量云圖Fig.13 Distributions of velocity and mass flow rate

為了進一步探究P/W對溫度場和流場的影響,在保證組件對邊距和邊界條件不變的情況下,開展了P/W的敏感性分析。圖14展示了提取參數的子通道編號。圖15為不同P/W下出口處溫度、質量流量和流速的徑向分布。可發現,P/W的變化對內通道影響不大,對邊通道影響較大,對角通道影響最大。隨著P/W的增大,內通道溫度減小,邊通道和角通道溫度增大,流量和流速則呈相反的變化趨勢。此外可發現,當P/W為2.070時,溫度、流量和流速的徑向分布更加平均,對堆芯熱工安全更加有利。圖16為不同P/W下最高溫度的變化。可發現,隨著P/W的增大,芯塊最高溫度下降了6.9 ℃,包殼最高溫度下降了7.8 ℃,冷卻劑最高溫度上升了7.9 ℃,而此時節徑比P/D僅從1.3增加到1.325,足以說明P/W和P/D對堆芯內局部熱點的分布有重要影響。

圖14 子通道編號Fig.14 Measured sub-channel number

圖15 P/W對熱工水力參數的敏感性分析Fig.15 Sensitivity analysis of P/W to thermal-hydraulic parameters

圖16 P/W對燃料棒和冷卻劑最高溫度的敏感性分析Fig.16 Sensitivity analysis of P/W on maximum temperature of fuel rods and coolants

此外,為了分析SACOS-LMR程序的瞬態運行特性,對ALFRED進行了假設的事故場景計算分析,即無保護流失(ULOF)事故。ULOF瞬態是由所有主泵停機引起的。因此,系統逐漸進入一種主冷卻劑質量流量不斷下降的連續過程,直到自然循環模式接管。自然循環堆芯質量流量為39.9 kg/s(額定值的21%)。

圖17為ULOF事故發生后中心燃料組件入口質量流量和出口溫度的變化。芯塊峰值溫度(PFT)在25 s內從1 832 ℃迅速降低至1 395 ℃,300 s內略有上升直至穩定在1 405 ℃,包殼峰值溫度(PCT)和冷卻劑出口最大溫度(PCoT)在25 s內迅速上升至780 ℃和758 ℃,在50 s內降低至652 ℃和638 ℃,然后逐漸穩定不變,該計算結果與ENEA-Relap5計算結果符合良好,但芯塊峰值溫度始終偏低。圖18為ULOF事故后中央燃料組件內最高溫度的變化。

圖17 中心燃料組件入口質量流量和出口溫度的變化Fig.17 Variations of inlet mass flowrate and outlet temperature in central FA

圖18 中心燃料組件的峰值溫度變化Fig.18 Variations of peak temperature in central FA

4 結論

本文基于西安交通大學熱工水力研究室自主開發的適用于壓水堆的子通道程序,通過添加液態金屬快堆特有的模型,如繞絲模型、盒間流模型、液態金屬對流換熱模型等,開發了適用于液態金屬快堆的子通道分析程序SACOS-LMR;同時結合卡爾斯魯厄核研究中心開展的37棒鈉冷瞬態實驗,完成了SACOS-LMR程序的瞬態功能驗證。程序計算的溫度場的相對偏差在15%以內,溫度最大偏差為25 ℃;采用SACOS-LMR程序對ALFRED堆芯開展了熱工安全特性分析,穩態計算結果與ANTEO-LFR的相對偏差在5%左右,預測的包殼最高溫度為544.4 ℃,芯塊最高溫度為1 832.3 ℃,均滿足設計限值要求;發生假想ULOF事故后,包殼峰值溫度在20 s時達到最大,為780 ℃,遠超包殼溫度的設計限值,存在安全隱患,程序計算結果與ENEA-Relap5計算結果整體上符合良好,但芯塊峰值溫度始終偏低。以上計算結果表明,SACOS-LMR子通道程序可用于液態金屬快堆的堆芯設計和熱工水力分析研究,能為中國液態金屬快堆設計和工程應用提供分析工具。