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

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

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

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

1 計算方法

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

1.1 隨機球床

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

1) 隨機模型

計算機模擬已成為隨機球床的主要研究方法[14-16]，分為物理和幾何兩類。物理方法如分子動力學和DEM方法重在球體的堆積過程，而幾何方法重在堆積結(jié)果。幾何方法主要分為序列添加和集合重排兩種，具體的實現(xiàn)方法各不相同。本文采用隨機序列添加(RSA)方法，其基本思想為每次迭代過程中生成1個新球到現(xiàn)有系統(tǒng)中。新增球需與已填充球或壁面有3個接觸點。新球的產(chǎn)生采用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)中的判據(jù)L可判斷圓(cc，rc)和圓(cP，rP)是否存在交集：若L>0，存在2個交點；若L=0，存在1個交點；若L<0，沒有交點。5) 若存在交集，交點即為與a、b和c球同時相切的新球球心。球心位置P1,2可通過式(4)、(5)計算得到。cm為兩個新球的中點，α和h為簡化公式引入的系數(shù)。

(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)

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

(6)

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

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

2) 規(guī)則模型

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

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

1.2 燃料球

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

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

(7)

2 TMSR-SF1堆芯

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

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

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

3 結(jié)果和討論

3.1 球床結(jié)構(gòu)

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

3.2 中子物理參數(shù)

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

1)keff

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

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

表2 TMSR-SF1滿裝載冷態(tài)有效增殖因數(shù)計算結(jié)果Table 2 keff for TMSR-SF1 fully-loaded cold state

2) 中子能譜

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

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

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

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

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

4 結(jié)論

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

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

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

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