隨機球床模擬方法對氟鹽冷卻高溫堆中子學參數的影響分析

冀銳敏，嚴 睿，康旭忠，楊 璞，于世和，鄒 楊

(中國科學院 上海應用物理研究所，上海 201800)

中國科學院于2011年啟動了戰略性先導專項“未來先進核能系統——釷基熔鹽堆(TMSR)核能系統”[1]，目標包括建造1座固態燃料釷基熔鹽堆(TMSR-SF1)。TMSR-SF1是熱功率為10 MW的氟鹽冷卻高溫堆(FHR)，是21世紀初美國新提出的第4代先進反應堆，結合了高溫氣冷堆和熔鹽堆的優點[2]。TMSR-SF1采用直徑6 cm的球形燃料和2LiF-BeF2冷卻劑。由于密度小于熔鹽，燃料球從堆芯底部投入，靠浮力上升后逐步構筑堆芯。已開展的相關堆積實驗[3-6]表明，TMSR-SF1堆積密度低于HTR-10，且不確定性較大。

TMSR-SF1堆芯結構是物理和熱工設計的重要基礎，也是其不確定性的主要來源。為解決該問題，本文采用蒙特卡羅程序(MCNP)[7]基于TMSR-SF1設計完成高保真物理分析。TMSR-SF1中子物理計算面臨很大挑戰：上萬顆燃料球在堆芯內隨機排布以及每顆燃料球的燃料區域內均有上萬個包覆顆粒隨機排布。蒙特卡羅程序通常采用3種方法完成球床堆[8-11]隨機結構的模擬：1) 建立簡立方(SC)、體心立方(BCC)等規則結構完成隨機結構的近似(方法1)。HTR-10首次臨界計算中曾采用了該方法，并獲得了與實驗較為吻合的結果[12-13]。2) 不提前精確定位，在計算中對其隨機定位[14-17]，MCNP尚不具備該功能(方法2)。3) 用數值方法模擬隨機結構(方法3)。通常采用離散單元法(DEM)和軟球模型進行受力分析，獲取各燃料球的精確位置[11,17]。此方法有3點不足：1) 計算模型導致球發生形變；2) 輸入參數中有些依據經驗輸入，如恢復系數、摩擦系數等；3) 所需的計算資源多。綜合以上因素，本文使用方法1和方法3模擬TMSR-SF1堆芯結構。方法3選擇工業領域內礦團堆積中使用的RSA方法[18-19]，其基本思想是逐個加入新球，直至完成給定區域的填充。新球位置的計算參考Halo方法[20]，根據現有球床計算出所有可能的新球的位置。隨后對上述位置進行判斷，排除超出邊界、發生重疊和不能獲得穩定支撐的球體位置[21]，以獲取穩定的堆積。為加速計算，對計算區域進行網格化，對各球體位置進行像素化。基于上述兩種方法建模，分析不同燃料球床以及燃料球描述模型對TMSR-SF1主要中子物理參數的影響，如keff、堆芯能譜、控制棒價值和溫度系數等。

1 計算方法

TMSR-SF1堆芯結構非常復雜，故臨界計算選用MCNP(版本6-6.1.1)，截面數據庫采用ENDF/B-Ⅵ.0。

1.1 隨機球床

本文分別采用隨機模型和規則模型來模擬隨機球床。

1) 隨機模型

計算機模擬已成為隨機球床的主要研究方法[14-16]，分為物理和幾何兩類。物理方法如分子動力學和DEM方法重在球體的堆積過程，而幾何方法重在堆積結果。幾何方法主要分為序列添加和集合重排兩種，具體的實現方法各不相同。本文采用隨機序列添加(RSA)方法，其基本思想為每次迭代過程中生成1個新球到現有系統中。新增球需與已填充球或壁面有3個接觸點。新球的產生采用Halo方法[20]：1) 與a球相切的所有球心可形成1個球面，以a球的球心為原點，2r(r為球半徑)為半徑，簡稱a球Halo(圖1a)。2) a球和b球的Halo交集，即與a球和b球均相切的球心集合，可形成1個圓(圖1b)，標識為圓(cc，rc)。向量cc表示圓心位置，rc為圓半徑，可由式(1)計算得到。向量a和b為a球和b球的球心位置。3) 與球a、b

a——a球Halo；b——a球和b球Halo的交集；c——c球Halo與圓(cc,rc)的交集；d——球心集合{P1,P2}圖1 Halo方法示意圖Fig.1 Halo algorithm scheme

和c均相切的球心可通過投影計算。c球Halo在圓(cc，rc)所屬平面的截面同樣為1個圓(cP，rP)，如圖1c所示。圓心位置cP和圓半徑rP可由式(2)計算得到。向量c為c球球心位置，向量n為a球到b球方向的單位向量，λ為向量ccc在向量n上的投影。4) 式(3)中的判據L可判斷圓(cc，rc)和圓(cP，rP)是否存在交集：若L>0，存在2個交點；若L=0，存在1個交點；若L<0，沒有交點。5) 若存在交集，交點即為與a、b和c球同時相切的新球球心。球心位置P1,2可通過式(4)、(5)計算得到。cm為兩個新球的中點，α和h為簡化公式引入的系數。

(1)

λ=(c-cc)n,n=(a-b)/|a-b|

(2)

L=rP+rc-|cc-cP|

(3)

cm=(1-α)cc+αcP

(4)

P1,2=cm±h(n×cccm)/|n×cccm|

(5)

計算出的新球還需進行有效性判斷：不得超出堆芯活性區；不能與已有燃料球發生重疊；落位應穩定[20-21]。在上述過程中需多次檢索給定球的鄰球。為加速檢索，初始將整個堆芯區域進行網格化，預分為若干立方體區域，立方體的邊長取燃料球的半徑r。根據燃料球位置產生與之對應的像素。假設堆芯原點位置為O=[Ox，Oy，Oz]，位置為P=[Px，Py，Pz]的燃料球可定位為(i，j，k)，如式(6)所示。

(6)

圖2 球床計算流程Fig.2 Flowchart for pebble bed generation

隨機球床的計算流程如圖2所示。簡述如下：輸入燃料球半徑和堆芯活性區信息，將堆芯網格化，建立已填充燃料球序列和待填充燃料球序列，用于存儲燃料球編號和位置信息；隨機產生放置于堆芯活性區頂部的1個球；利用Halo方法結合已有燃料球計算所有可能的新球位置；對其進行有效性判斷；將有效的燃料球放置于待填充球序列；依據填充的密集程度，對待填充燃料球進行排序；取出排名第一的球至已填充序列；根據新填充的燃料球重復上述過程，直至待填充序列為空集，即堆芯活性區已填滿燃料球。

2) 規則模型

HTR-10在首次臨界計算中曾采用一些規則的模型來模擬隨機堆積的燃料球床，如BCC模型、FCC模型、CHPOP模型等[12-13]，并獲得了與實驗較為一致的結果。本文使用BCC模型進行等效處理，球床模型細節示于圖3。

a——周期性規則填充的柵格；b——燃料球堆積水平方向；c——燃料球堆積豎直方向圖3 燃料球堆積的規則等效模型Fig.3 Equivalent regular model of fuel pebble packing

1.2 燃料球

圖4 燃料球模型與URAN擾動原理示意圖Fig.4 Schematic of fuel pebble and URAN vibration principle

受限于MCNP程序中對柵元、曲面等編號數量的限制，并未對燃料球中的包覆顆粒進行隨機位置的建模，而將其進行規則的建模，然后通過MCNP中的URAN卡對其位置進行隨機擾動[7]。URAN卡實現擾動的基本原理圖示于圖4。詳細步驟如下：1) 對燃料球中包覆顆粒詳細建模；2) 將包覆顆粒填充至燃料填充區域中劃分好的簡單立方格中，簡單立方格的尺寸由燃料球中包覆顆粒的堆積密度確定；3) 設定URAN卡，指定擾動的柵元編號以及擾動范圍；4) 在計算過程中，待抽樣到該柵元時，抽樣執行隨機擾動，將包覆顆粒隨機移動(式(7))。其中，[x，y，z]為擾動后的位置；[x0，y0，z0]為擾動前的位置；ξ1、ξ2和ξ3為0～1的隨機數；δx、δy和δz為指定的擾動范圍。值得指出的是，該擾動受擾動范圍的限定，通常限制于劃分好的簡單立方網格內。即實現了燃料球內包覆顆粒有限的隨機分布，與真實的分布仍有一定的差異。單燃料球中包覆顆粒分析表明，相比于完全的隨機，這種簡單的處理方式會導致邊界處包覆顆粒的切割，引起無限增殖因數低估[10]。

(7)

2 TMSR-SF1堆芯

TMSR-SF1使用全陶瓷型包覆顆粒制成的球形燃料元件、氟鋰鈹熔鹽(FLiBe)作為冷卻劑、石墨作為慢化劑和堆芯的結構材料，熱功率為10 MW，其堆芯結構示意圖如圖5a、b所示，主要由堆芯活性區和石墨反射層組成。堆芯活性區由反射層的石墨構件圍成，形成上下圓臺和中間圓柱區域。圓柱直徑為135 cm，高為180 cm；上下圓臺最小直徑為30 cm；圓臺斜面與水平面的夾角為30°；反射層外圍形狀為圓柱體，直徑為285 cm，高為306 cm，反射層外的堆芯圍筒厚度為2 cm。堆芯活性區內燃料球隨機堆積，燃料球之間的空隙形成不規則的流道，供冷卻劑自下向上流動，帶走裂變產生的熱量。靠近活性區的石墨側反射層中布置16個孔道，其中間隔布置的12根控制棒作為反應性控制系統，排空熔鹽作為第2套停堆系統。燃料球由中心燃料填充區域和外部石墨殼組成，如圖5c所示。TMSR-SF1設計中所用燃料球的基本參數列于表1。燃料填充區域內由包覆燃料顆粒(TRISO)和石墨基體混合而成。TRISO由UO2核心和4層包覆層組成。

a——堆芯三維圖；b——堆芯俯視圖；c——燃料球示意圖圖5 TMSR-SF1堆芯結構示意圖Fig.5 Schematic diagram of TMSR-SF1

表1 TMSR-SF1球形燃料元件參數Table 1 TMSR-SF1 fuel element characteristic parameter

3 結果和討論

3.1 球床結構

球床堆積模擬中的一個關鍵問題是球體堆積結構的描述。除堆芯平均堆積密度外，還分別使用投影以及軸向和徑向堆積密度分布[22]兩種方法給出兩種不同建模方法所獲球床的結構細節，如圖6所示。圖6a、b為RSA方法獲得的隨機球床在底面和側面的投影，圖6c、d為規則堆積方法獲得的規則球床在底面和側面的投影，圖6e、f為與圖6a、b對應的徑向和軸向堆積密度分布，圖6g、h為與圖6c、d對應的徑向和軸向堆積密度分布。可見，隨機球床在靠近活性區邊界處的區域以及球床頂部區域的堆積密度較高，軸向堆積密度隨球床高度的降低而有所降低；規則球床在靠近活性區邊界處的區域堆積密度較低，軸向堆積密度則為規律的鋸齒狀。對比可見，盡管兩組球床模型的堆積密度均接近60%，但球床的細節結構有很大不同。

3.2 中子物理參數

使用3種建模方法計算滿裝載時的冷態反應性：1) RSA方法模擬隨機堆積球床，依據各燃料球位置利用MCNP實現詳細建模(方法1)；2) 等效規則球床，MCNP直接建模，通過不同的堆積密度擬合出給定堆積密度的情況(方法2)；3) 規則球床等效隨機球床，在相同堆積密度的前提下，BCC結構填充原點有差異，以統計活性區邊界對燃料球的不同切割效應(方法3)。為計算包覆顆粒隨機分布的影響，在上述3種輸入的前提下，通過MCNP程序URAN卡實現所有燃料球中包覆顆粒的隨機移動，分別計算球隨機分布以及包覆顆粒隨機分布對計算結果的影響。

1)keff

MCNP計算中，每代粒子數100 000，循環總數1 000代，有效循環500代，計算得到有效增殖因數keff的標準偏差小于0.000 10。上述3種方法的計算結果列于表2。由表2可見，方法1的計算結果最大，方法2、3的計算結果非常接近。平均堆積密度為59.58%時，隨機堆積模型結果較規則堆積模型結果約大500 pcm，接近0.5%。3種計算方法中由包覆顆粒隨機分布導致的有效增殖因數差異分別為+31、+26、-21 pcm。由此可見，包覆顆粒的隨機分布對滿裝載冷態有效增殖因數的影響較小，相對于燃料球隨機分布其影響完全可忽略。

圖6 球床堆積密度分布對比Fig.6 Comparison of packing density distribution of pebble bed

表2 TMSR-SF1滿裝載冷態有效增殖因數計算結果Table 2 keff for TMSR-SF1 fully-loaded cold state

2) 中子能譜

對于隨機球床，距活性區邊界約3個燃料球直徑范圍內的堆積密度出現較大振蕩，因此將堆芯活性區(R<67.5 cm，R為模型徑向距離)、活性區中心區域(R<50 cm)、活性區邊界區域(50 cm

3) 控制棒價值和溫度系數

為進一步量化不同模型的影響，基于第10組算例和前2種方法計算控制棒總價值和溫度系數，結果列于表3。由表3可見，相對于隨機模型，規則模型會高估控制棒價值約5.78%，而溫度系數的結果較為一致。控制棒價值高估與堆芯有效增殖因數低估的原因一致，燃料球排布的簡化增加了堆芯中子的泄漏。

圖7 堆芯各區域中子通量密度分布對比Fig.7 Comparison of neutron flux density in different regions

表3 TMSR-SF1滿裝載冷態控制棒價值和溫度系數Table 3 Control rod worth and temperature coefficient for TMSR-SF1 fully-loaded cold state

4 結論

基于TMSR-SF1設計模型，利用RSA方法(方法1)實現了隨機球床的模擬，與BCC結構的規則等效球床進行了球床堆積(方法2)情況對比；基于MCNP程序中URAN卡完成了燃料球內包覆顆粒隨機分布的模擬；針對上述計算模型進行了中子物理關鍵參數分析，得到如下主要結論。

1) 盡管隨機球床和規則球床的平均堆積密度均接近60%，但球床結構的細節有很大不同。隨機球床在靠近活性區邊界約3個燃料球直徑區域內形成較為規則的分布，徑向堆積密度在該區域內有一定的振蕩，在堆芯區域則較為平穩。隨機球床的軸向堆積密度從頂部到底部逐步降低。規則球床的徑向和軸向堆積密度均呈現鋸齒狀分布。

2) 由于球床結構的不同，規則球床模型會造成更多的中子泄漏，導致堆芯有效增殖因數的低估和布置于反射層的控制棒價值的高估，但對溫度系數的影響可忽略。

3) 燃料球中包覆顆粒隨機分布的影響遠小于球床隨機分布的影響。