冀銳敏，嚴(yán) 睿，康旭忠，楊 璞，于世和，鄒 楊

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

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

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

1 計(jì)算方法

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

1.1 隨機(jī)球床

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

1) 隨機(jī)模型

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

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

(6)

圖2 球床計(jì)算流程Fig.2 Flowchart for pebble bed generation

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

2) 規(guī)則模型

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

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

1.2 燃料球

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

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

(7)

2 TMSR-SF1堆芯

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

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

3.2 中子物理參數(shù)

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

1)keff

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

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

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

2) 中子能譜

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

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

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

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

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

4 結(jié)論

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

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

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

3) 燃料球中包覆顆粒隨機(jī)分布的影響遠(yuǎn)小于球床隨機(jī)分布的影響。