賈國斌，伍建輝，陳金根，顧國祥，蔡翔舟，戴 葉,2

(1.中國科學(xué)院 上海應(yīng)用物理研究所，上海 201800；2.中國科學(xué)院 先進核能創(chuàng)新研究院，上海 201800；3.中國科學(xué)院大學(xué)，北京 100049)

液態(tài)熔鹽堆是第4代先進核能系統(tǒng)中唯一使用液體燃料的反應(yīng)堆，與燃料固定于組件的固態(tài)反應(yīng)堆不同，熔鹽既作為冷卻劑又作為燃料，充滿整個一回路(含堆芯)。熔鹽堆有多種設(shè)計堆型，按堆芯中維持鏈式裂變反應(yīng)的平均中子能量不同，可劃分為以法國與俄羅斯分別提出的MSFR[1-4]與MOSART[5]為代表的熔鹽快堆，以及以美國橡樹嶺國家實驗室(ORNL)提出的MSRE[6]與MSBR[7]為代表的熔鹽熱堆。熔鹽快堆通常使用罐式結(jié)構(gòu)，活性區(qū)無慢化劑，堆芯結(jié)構(gòu)簡單。熔鹽熱堆一般采用石墨構(gòu)件作為慢化劑，并采用鎳基合金對其進行固定，控制棒組件根據(jù)需要可布置在活性區(qū)或側(cè)反射層中。相對于熔鹽快堆，熔鹽熱堆具有更復(fù)雜的幾何結(jié)構(gòu)。釷基熔鹽堆核能專項[8]設(shè)計的2 MW液態(tài)熔鹽實驗堆為典型熱堆，出于安全與測試需要，堆芯內(nèi)布置了較多的反應(yīng)性控制設(shè)備與監(jiān)測孔道，從而對中子學(xué)計算工具提出了更高的精確建模要求。蒙特卡羅程序具有精確的幾何建模功能，典型的程序如美國愛達荷國家實驗室(INL)開發(fā)的MCNP[9]、美國ORNL開發(fā)的SCALE[10]以及清華大學(xué)研發(fā)的RMC[11]等。由于以上蒙特卡羅程序可進行精確三維幾何建模，可靠性高，目前釷基熔鹽堆(TMSR, Thorium-based Molten-Salt Reactor)核能系統(tǒng)專項中子學(xué)設(shè)計主要基于蒙特卡羅程序。但此類程序的計算成本較高，相對而言確定論程序可顯著降低計算成本，進而極大縮短堆芯設(shè)計參數(shù)的優(yōu)化時間，是反應(yīng)堆物理分析的重要模擬工具。

目前，針對含有控制棒及實驗孔道等復(fù)雜幾何組件的液態(tài)熔鹽堆開展的確定論程序等效建模與適用性分析等較少，但對包含控制棒組件的輕水堆與重水堆已進行了較多的均勻化分析。加拿大蒙特利爾技術(shù)學(xué)院利用超組件方法[12-13]對重水堆CANDU堆芯中具有特殊結(jié)構(gòu)形式的控制棒組件進行了三維建模與均勻化。西安交通大學(xué)在對輕水堆堆芯圍板和反射層進行宏觀群常數(shù)加工時[14-16]，采用不連續(xù)因子方法確保模型均勻化前后反應(yīng)率與面凈流守恒，計算所得到的堆芯內(nèi)部功率分布相對偏差小于5%，臨界硼質(zhì)量分數(shù)相對偏差小于5×10-5。清華大學(xué)[17-18]通過使用不連續(xù)因子，對側(cè)反射層中包含控制棒的組件進行了均勻化，計算得到的控制棒與吸收球價值和輸運計算所得到的基準解相近，最大絕對誤差約為0.001。上述研究表明，不連續(xù)因子的使用可顯著提高計算精度，但需相應(yīng)輸運程序與擴散程序的支持。

本文基于加拿大蒙特利爾技術(shù)學(xué)院研制的反應(yīng)堆物理數(shù)值計算確定論程序Dragon[19-20]與Donjon[21]，對TMSR專項的液態(tài)熔鹽實驗堆(2 MWth)堆芯進行等效建模與分析。

1 計算模型與方法

1.1 計算模型

液態(tài)熔鹽實驗堆的堆芯如圖1所示。堆芯由活性區(qū)的燃料組件、上下熔鹽腔室、鎳基合金及石墨反射層組成。圖1中紅色為石墨慢化劑，藍色為燃料鹽，綠色為鎳基合金，黑色為控制棒。堆芯活性區(qū)包含1個中子源孔道，8個實驗測試孔道(分成2組，每組含有4個實驗孔道)，6個控制棒孔道?？刂瓢艨椎缽较蛭恢萌鐖D1b所示，其中1號(rod 1)與6號(rod 6)控制棒為調(diào)節(jié)棒組，3號(rod 3)、4號(rod 4)、5號(rod 5)控制棒為補償棒組，2號(rod 2)控制棒為緊急停堆棒。為了展示控制棒在堆芯軸向上的上限位與下限位，圖1中將4號控制棒插入堆芯并處于下限位，其他控制棒處于上限位。可看出，下限位在下熔鹽腔室鎳基合金的正上方，上限位在上熔鹽腔室與燃料組件之間的鎳基合金內(nèi)。堆芯主要參數(shù)列于表1。

a——縱截面；b——橫截面圖1 液態(tài)熔鹽實驗堆堆芯幾何結(jié)構(gòu)Fig.1 Core configuration of liquid molten salt experimental reactor

熔鹽堆活性區(qū)的方形組件邊長為5.1 cm，按照內(nèi)部結(jié)構(gòu)不同分為4種類型，分別為燃料組件、實驗孔道與中子源孔道組件、控制棒孔道組件、控制棒插入孔道組件，其詳細幾何參數(shù)如圖2所示。燃料組件的熔鹽孔道為帶有圓弧的扁長方形。實驗孔道與中子源孔道組件為鎳基合金同心圓環(huán)，四周為燃料鹽，中子源孔道與實驗孔道縱向貫穿整個堆芯活性區(qū)?？刂瓢艨椎赖闹行男A環(huán)為鎳基合金，為控制棒的移動起導(dǎo)向作用?？刂瓢舨牧蠟镚d-Al合金，為同心圓環(huán)結(jié)構(gòu)。

1.2 計算工具與組件等效模型

本文使用傳統(tǒng)的兩步法進行分析，利用加拿大蒙特利爾技術(shù)學(xué)院開發(fā)的反應(yīng)堆物理數(shù)值計算程序包進行計算。首先使用Dragon進行組件計算，得到少群均勻化群常數(shù)(組件均勻化)，然后使用Donjon進行宏觀群常數(shù)的插值與全堆擴散計算。擴散計算時調(diào)用程序包中Trivac程序[22]的本征值求解模塊進行三維擴散計算，最后使用Donjon的功率輸出模塊給出全堆的三維功率分布。

表1 液態(tài)熔鹽實驗堆堆芯的主要參數(shù)Table 1 Key parameter of core of liquid molten salt experimental reactor

為使擴散計算得到的堆芯有效增殖因數(shù)及功率分布與蒙特卡羅計算條件一致，均勻化前后需保證堆芯不同區(qū)域內(nèi)、不同反應(yīng)道的反應(yīng)率守恒[23]。少群常數(shù)計算公式為：

a——燃料組件；b——實驗孔道與中子源孔道組件；c——控制棒孔道組件；d——控制棒插入孔道組件圖2 組件分類Fig.2 Assembly classification

其中：Σx為x類型宏觀群常數(shù)；τ為反應(yīng)率；VⅠ為均勻化Ⅰ區(qū)的體積；ФⅠ為均勻化區(qū)VⅠ的平均中子通量；V為均勻化區(qū)間的體積；r為均勻化區(qū)間的幾何坐標；E為中子能量；Ф為中子通量。

對中子吸收截面較大的部件，如控制棒組件與用于固定石墨組件的鎳基合金，使用超組件方法將周圍的燃料組件包含在一起整體均勻化得到1組群常數(shù)。本文根據(jù)堆芯內(nèi)不同區(qū)域的幾何特征，將其分為6種組件，等效均勻化模型如下。

1) 燃料組件

由于中子輸運程序Dragon不能處理帶有圓弧的熔鹽孔道，將弧形熔鹽孔道等效為等面積的長方形熔鹽孔道。等效方法為保持孔道的寬度(0.417 cm)與面積不變，按照單個熔鹽孔道等面積的方法計算，得到孔道長為2.364 4 cm。使用Dragon程序中GEO模塊下的兩維CAR2D幾何模型，將燃料組件分為5×5的兩維幾何結(jié)構(gòu)，如圖3所示，其中紅色為石墨，藍色為燃料鹽。為了將4個熔鹽孔道劃分為單個區(qū)域的材料，將石墨組件劃分為21個區(qū)域。

圖3 燃料組件均勻化模型Fig.3 Homogenization model of fuel assembly

2) 實驗孔道與中子源孔道組件、控制棒孔道組件

實驗孔道與中子源孔道組件、控制棒孔道組件都具有規(guī)則幾何形狀，因此使用Dragon程序的CARCEL模塊進行精確建模與均勻化，得到相應(yīng)組件群常數(shù)。

3) 控制棒插入孔道組件

對控制棒插入孔道組件進行少群常數(shù)加工時，考慮其對周圍燃料組件能譜影響較大，使用超柵元方法將包圍控制棒的臨近8個燃料組件與控制棒組件同時進行均勻化，均勻化后得到一種材料的少群常數(shù)。計算模型如圖4所示，其中紅色為石墨，藍色為熔鹽管道，黃色為控制棒，綠色為鎳基合金。燃料組件中熔鹽孔道同樣等效為等體積的矩形結(jié)構(gòu)，使用Dragon程序GEO模塊下的兩維CAR2D單獨進行建模，中心的同心圓環(huán)控制棒使用CARCEL單獨精確建模，最終通過CELL模塊將以上二維模型進行組合。

圖4 控制棒插入孔道組件的均勻化模型Fig.4 Homogenization model of control rod assembly with control rod inserted

4) 上下熔鹽腔室

堆芯的上下熔鹽腔室內(nèi)部為燃料鹽，外部使用鎳基合金進行燃料鹽的包容。盡管此區(qū)域的幾何模型復(fù)雜，但中子泄漏較大且不含慢化劑，因此該區(qū)域裂變功率較低，可根據(jù)燃料鹽的等體積方法進行建模。由于用于固定石墨組件的鎳基合金中子吸收截面較大，因此使用三維幾何結(jié)構(gòu)進行建模，計算模型如圖5所示。其中，石墨燃料組件高度為30 cm，用于包容熔鹽腔室的鎳基合金使用實際設(shè)計的厚度參數(shù)3 cm，通過熔鹽等體積的方法，得到熔鹽腔室內(nèi)的熔鹽等效高度為3.86 cm。用于固定石墨組件的鎳基合金使用實際設(shè)計參數(shù)，高度為5 cm。由于Dragon程序在組件均勻化時無法準確處理真空邊界條件(采用真空邊界時，共振自屏計算選擇等效處理SHI模塊時會報錯；如選擇子群共振USS模塊，計算所得到的群常數(shù)與蒙特卡羅統(tǒng)計算出的截面相差較大)，從而本文采用全反射邊界條件進行模擬。為考慮中子泄漏對能譜的影響，使用強中子吸收體10B等效容器外的真空，同時為了將逃逸出堆芯容器的中子全部吸收，保守計算使用厚度為50 cm的10B進行等效。以上模型使用Dragon程序中的CELL模塊進行三維建模。

圖5 上下熔鹽腔室組件均勻化模型Fig.5 Homogenization model of top and bottom molten salt chambers

5) 上熔鹽腔室的控制棒

上熔鹽腔室中控制棒的均勻化模型與上下熔鹽腔室的基本相同，使用三維立方體幾何模型。其中包含石墨燃料組件、用于固定石墨組件的鎳基合金、經(jīng)過等效的熔鹽腔室中的熔鹽(內(nèi)含控制棒)、外側(cè)鎳基合金和使用強吸收體等效的外界真空。使用Dragon程序中的CELL幾何模塊進行三維幾何建模。

6) 側(cè)反射層

為了得到側(cè)反射層石墨的能譜，均勻化模型使用兩維模型，將堆芯徑向方向上靠近石墨的兩個燃料組件與石墨反射層、鎳基合金包殼進行均勻化。同樣為了考慮中子泄漏對鎳基合金的影響，外側(cè)真空使用厚度為50 cm的強吸收材料10B進行等效。以上使用Dragon程序中的CELL幾何進行兩維建模。均勻化后得到反射層石墨組件，用于包容側(cè)反射層石墨的鎳基合金兩種材料的少群常數(shù)。

對以上6種類型組件進行群常數(shù)加工時，使用Dragon程序中的NXT碰撞概率模塊、USS子群共振模塊以及FLU模塊中傳統(tǒng)的基模修正B1模型與輸運修正PNL模型求解擴散系數(shù)。同時，在對補償棒組(3～5號控制棒)與控制棒全部插入計算分析時，首先使用基模修正獲得堆芯keff，然后再將此值代入到輸運方程中進行泄漏修正[13]，程序計算流程與詳細理論公式參考Dragon說明書[20]。

堆芯擴散模型使用等間距劃分的網(wǎng)格模型。堆芯徑向的網(wǎng)格為邊長5.1 cm的正方形網(wǎng)格。對于軸向，活性區(qū)的石墨燃料組件網(wǎng)格高度為10 cm，用于固定燃料組件的鎳基合金、上下熔鹽腔室、包容上下熔鹽腔室熔鹽的鎳基合金的網(wǎng)格高度分別為5、3.86與3 cm。當堆芯單獨插入2號控制棒時，全堆擴散的堆芯模型如圖6所示。由圖6可看出，組件均勻化后的幾何模型與實際堆芯模型(圖1)的不同之處主要在上下熔鹽腔室部分，即將圓弧形上下熔鹽腔室的幾何結(jié)構(gòu)等效為等體積的長方體結(jié)構(gòu)；其次，活性區(qū)與側(cè)反射層交界面處，離散為不連續(xù)的矩形結(jié)構(gòu)。

a——堆芯徑向；b——堆芯軸向圖6 堆芯均勻化模型Fig.6 Homogenization model of core

本文采用Dragon程序進行輸運計算時使用能群數(shù)為361群的精細群結(jié)構(gòu)，選用ENDF/B-Ⅶ.1的微觀數(shù)據(jù)庫，從而保證與本文用于驗證計算的蒙特卡羅程序MCNP5的數(shù)據(jù)庫類型相同?？紤]到液態(tài)熔鹽實驗堆與高溫氣冷堆的石墨慢化能力對快中子的慢化與熱中子的向上散射能力類似，在進行少群能群結(jié)構(gòu)劃分時，重點參考美國賓夕法尼亞大學(xué)對高溫氣冷堆PBMR[24]分析時采用的能群結(jié)構(gòu)[25]。同時考慮到液態(tài)熔鹽中19F與6Li在10～100 keV及100 keV～1 MeV能區(qū)范圍內(nèi)存在共振彈性散射峰，以及初始臨界堆芯中無Pu元素，本文在其能群結(jié)構(gòu)基礎(chǔ)上進行優(yōu)化整理，能群結(jié)構(gòu)列于表2、3。由于PBMR控制棒的組分為B4C，而本文計算的液態(tài)熔鹽實驗堆控制棒組分為Gd-Al合金，為了充分考慮155Gd在2.0 eV與2.6 eV附近的共振吸收峰，分別使用兩種快熱能群分界(2.38 eV與3.93 eV)開展初步敏感性分析。

表2 5種少群常數(shù)能群結(jié)構(gòu)Table 2 Five kinds of broad energy group structure

表3 能群結(jié)構(gòu)劃分的物理依據(jù)Table 3 Physical explanation for recommended group structure

2 計算結(jié)果

2.1 堆芯臨界計算結(jié)果

基于少群能群結(jié)構(gòu)計算了液態(tài)熔鹽實驗堆在運行溫度為920 K時不同控制棒單獨插入與組合插入堆芯時的有效增殖因數(shù)，計算結(jié)果如圖7所示。

圖7中，no rod表示控制棒全部拔出，rod 16表示調(diào)節(jié)棒rod 1和rod 6全部插入，rod 345表示補償棒rod 3、rod 4、rod 5全部插入堆芯，all in表示堆芯內(nèi)6根控制棒全插入堆芯。當快熱能群分界為2.38 eV時，從2群到7群計算得到的最大相對誤差依次為1.2%、0.9%、0.6%、0.5%、0.4%；當使用3.93 eV作為快熱能群分界時，最大相對誤差依次為1.3%、1.0%、0.7%、0.65%、0.48%。由此可知隨能群數(shù)的增加計算精度也不斷增加，當能群結(jié)構(gòu)為7群時，計算精度可滿足小于0.5%的要求。同時由圖7可看出，2.38 eV快熱能群分界具有更高的計算精度，因此采用2.38 eV作為快熱能群分界。此外，最大計算誤差均出現(xiàn)在調(diào)節(jié)棒rod 1與rod 6同時插入堆芯時的工況，其主要原因為控制棒組件均勻化時未考慮側(cè)反射層石墨對能譜的影響，導(dǎo)致計算結(jié)果出現(xiàn)較大偏差。

圖7 堆芯不同運行狀態(tài)下的有效增殖因數(shù)(a)及與MCNP5結(jié)果的相對誤差(b)Fig.7 Effective multiplication factor (a) and relative error compared with MCNP5 result (b) under different operating conditions of core

根據(jù)圖7a不同控制棒插入堆芯的有效增殖因數(shù)，由下式計算控制棒插入堆芯后的積分價值：

計算得到的控制棒價值如圖8所示，由圖8可看出，隨能群數(shù)增加，控制棒價值也不斷接近MCNP5的基準解，當能群數(shù)為7群時，最大絕對誤差為274 pcm。

2.2 功率分布

為驗證Dragon和Donjon程序?qū)θ埯}堆堆芯功率分布模擬的可靠性，使用7群結(jié)構(gòu)，快熱能群分界為2.38 eV，計算得到堆芯軸向與徑向功率分布，并與MCNP5計算結(jié)果進行了比較。控制棒拔出與補償棒組(rod 345)插入堆芯，兩種條件下堆芯活性區(qū)軸向功率分布如圖9所示。由圖9可見：控制棒插入前后軸向功率分布基本沒有變化，即軸向功率峰均在堆芯的中心；靠近上下熔鹽腔室區(qū)域的功率與MCNP5程序計算結(jié)果有較大誤差，最大相對誤差為10%。出現(xiàn)誤差的最主要原因是堆芯幾何問題，即對堆芯進行等效建模時，將上下熔鹽腔室等效為等體積燃料鹽的長方體，從而引起堆芯軸向中子泄漏的增加，降低了此區(qū)域的功率。

圖8 單個控制棒與控制棒組的價值Fig.8 Control rod worth of single rod and rod group

由于確定論計算的全堆模型(圖6)將徑向反射層與燃料石墨組件邊緣等效為不連續(xù)的矩形結(jié)構(gòu)，因此MCNP5程序統(tǒng)計徑向功率時不能與確定論網(wǎng)格完全對應(yīng)。本文使用MCNP5程序統(tǒng)計功率時，使用了更小(邊長為3 cm)的網(wǎng)格進行功率統(tǒng)計。圖10對比了堆芯活性區(qū)軸向在80～90 cm處的徑向功率分布。由圖10可見：控制棒不插入堆芯前，功率峰因子在堆芯中心，且實驗孔道與中子源孔道組件、控制棒孔道組件的功率均低于附近區(qū)域的功率；

圖9 控制棒拔出(a)與補償棒組插入(b)堆芯時的軸向功率分布Fig.9 Axial power distribution with control rod drawn out (a) and compensating rod inserted (b)

當補償棒插入堆芯后，引起徑向功率峰向2號控制棒孔道附近移動。從計算結(jié)果可知，MCNP5程序與確定論計算得到的徑向功率分布基本相同。

a——無控制棒插入，MCNP5；b——無控制棒插入，確定論；c——補償棒插入，MCNP5；d——補償棒插入，確定論圖10 堆芯徑向功率分布對比Fig.10 Comparison of core radial power distribution

2.3 溫度反應(yīng)性系數(shù)

溫度反應(yīng)性系數(shù)是反應(yīng)堆安全分析的重要參數(shù)，為保證反應(yīng)堆具有本征安全性，總溫度反應(yīng)性系數(shù)必須為負。液態(tài)熔鹽實驗堆的總溫度反應(yīng)性系數(shù)可分解為燃料鹽與慢化劑石墨的溫度反應(yīng)性系數(shù)，其中燃料鹽的溫度反應(yīng)性系數(shù)由多普勒效應(yīng)與密度效應(yīng)組成，總的溫度反應(yīng)性系數(shù)計算公式[26-27]為：

由于液態(tài)熔鹽實驗堆燃料鹽的運行溫度范圍為708～973 K，因此本文主要集中對700～1 000 K內(nèi)的溫度反應(yīng)性系數(shù)進行計算分析。計算過程中，控制棒全部處于上限位。圖11示出不同溫度條件下采用7群結(jié)構(gòu)計算所得到的堆芯有效增殖因數(shù)。由圖11可看出，在不同溫度條件下，計算所得到的堆芯整體、燃料鹽、慢化劑的keff結(jié)果與MCNP5基準解的最大相對誤差分別為0.17%、0.14%、0.20%。

圖11 不同溫度條件下的有效增殖因數(shù)(a)及與MCNP5結(jié)果的相對誤差(b)Fig.11 Effective multiplication factor (a) and relative error compared with MCNP5 result (b) at different temperatures

通過計算液態(tài)熔鹽實驗堆燃料鹽與石墨分別在不同溫度下的有效增殖因數(shù)，得到燃料鹽、慢化劑和總的溫度反應(yīng)性系數(shù)，結(jié)果如圖12所示。由圖12可看出，實驗堆總的溫度反應(yīng)性系數(shù)約為-10.5 pcm/K，燃料鹽的溫度反應(yīng)性系數(shù)約為-5.1 pcm/K，慢化劑的溫度反應(yīng)性系數(shù)約為-4.9 pcm/K，與MCNP5計算結(jié)果的最大誤差依次為0.58、0.12和0.58 pcm/K。

3 總結(jié)與展望

本文基于中子輸運程序Dragon與中子擴散程序Donjon，使用傳統(tǒng)的兩步法對液態(tài)熔鹽實驗堆進行了臨界分析，并通過MCNP5程序驗證了Dragon與Donjon程序計算液態(tài)熔鹽實驗堆的適用性。本文使用5種少群能群結(jié)構(gòu)進行了計算分析，并對快熱能群分界進行了敏感性分析，結(jié)果表明使用7群少群結(jié)構(gòu)計算得到的結(jié)果與MCNP5程序計算的基準解最為接近。

圖12 溫度反應(yīng)性系數(shù)與MCNP5程序結(jié)果的比較Fig.12 Comparison of temperature reactivity coefficient with that calculated by MCNP5

本文主要圍繞液態(tài)熔鹽實驗堆初始臨界展開，暫未考慮燃耗計算以及燃料鹽后處理，同時未考慮熔鹽流動對計算結(jié)果的影響。下一步將耦合上述影響因素進一步拓展本工具在熔鹽堆物理分析方面的模擬功能。