戴 明，朱貴鳳，戴 葉，于世和，鄒 楊，3，余笑寒，3

（1.中國科學(xué)院 上海應(yīng)用物理研究所，上海 201800；2.中國科學(xué)院大學(xué)，北京 100049；3.中國科學(xué)院 核輻射與核能技術(shù)重點(diǎn)實(shí)驗(yàn)室，上海 201800）

氟鹽冷卻球床堆（PB-FHR）經(jīng)10余年的發(fā)展，已逐步由預(yù)概念設(shè)計(jì)走向試驗(yàn)堆基準(zhǔn)設(shè)計(jì)。2012年美國能源部（DOE）啟動合作研究項(xiàng)目（IRP）研究先進(jìn)高溫堆（AHTR）［1］的發(fā)展技術(shù)路線，包括進(jìn)行實(shí)驗(yàn)堆、原型堆和商業(yè)堆的設(shè)計(jì)以及關(guān)鍵技術(shù)、策略需求分析，在該項(xiàng)目中AHTR被命名為氟鹽冷卻高溫堆（FHR），并選用PBFHR（氟鹽冷卻球床堆）作為設(shè)計(jì)基準(zhǔn)［2］。中國科學(xué)院自2011年開始啟動實(shí)施“未來先進(jìn)核裂變能”戰(zhàn)略性先導(dǎo)科技專項(xiàng)，進(jìn)行釷基熔鹽堆核能系統(tǒng)的研究（專項(xiàng)簡稱TMSR），PB-FHR選為該專項(xiàng)試驗(yàn)堆基準(zhǔn)設(shè)計(jì)方案之一。

PB-FHR 中子學(xué)求解與傳統(tǒng)壓水堆方法相比存在一些差異，如柵元能譜計(jì)算要考慮雙重非均勻性及堆芯不同位置泄漏反饋影響［3］，流動球床燃耗計(jì)算時(shí)譜區(qū)內(nèi)不同批次燃料同時(shí)存在且認(rèn)為具有一致的有效微觀截面［4］等。PB-FHR的少群截面計(jì)算是其堆芯擴(kuò)散計(jì)算分析的重要研究前提，它將對擴(kuò)散計(jì)算的準(zhǔn)確性產(chǎn)生重大影響。球床堆復(fù)雜的幾何模型構(gòu)成所謂的雙重非均勻性系統(tǒng)，第1重非均勻性為燃料核和包覆層及石墨基質(zhì)，第2重非均勻性為燃料區(qū)、球石墨殼和氟鹽冷卻劑。雙重非均勻性要求能準(zhǔn)確求解可分辨共振區(qū)的共振吸收［5］。對于可分辨共振能區(qū)一般采用求解基于碰撞概率法的精細(xì)群慢化方程，在ZUT 共振處理軟件中，采用解析的數(shù)值積分精確求解相關(guān)碰撞概率［4］，燃料核逃脫概率與燃料核的總截面相關(guān)，而該截面在共振區(qū)又隨能量改變很大，所以采用精細(xì)群求解。SRAC 的共振處理程序PEACO 能調(diào)用共振區(qū)點(diǎn)連續(xù)截面MCROSS庫直接求解基于碰撞概率法的精細(xì)群慢化方程。對于雙重非均勻性的處理，可采用類似于SCALE 處理雙重非均勻柵元的CSAS［6］的共振計(jì)算兩步法，即先計(jì)算包覆顆粒共振能區(qū)空間均勻化有效微觀共振截面，再利用該有效截面進(jìn)行燃料球柵元計(jì)算。莊坤等［7］使用SRAC 采用該兩步法與MVP-BURN 程序比較了燃料球有效增殖因數(shù)隨燃耗的變化。SRAC也提供直接一步處理雙重非均勻性的功能，相關(guān)碰撞概率的求解沒有ZUT 精確，而是采用了一些假設(shè)。該方法與兩步法無本質(zhì)區(qū)別。

本工作通過與采用連續(xù)能量蒙特卡羅軟件MCNP計(jì)算的參考模型對比，分析確定論軟件SRAC對于準(zhǔn)確加工PB-FHR少群截面的適用性。主要比較參量有無窮增殖因數(shù)（Kinf），少群均勻化總截面、俘獲截面和裂變截面。計(jì)算包括兩個(gè)溫度點(diǎn)：冷態(tài)（300K）和熱態(tài)（1 000K）。

1 計(jì)算程序及方法

1.1 確定論軟件

確定論軟件SRAC［8］是針對各種熱中子堆的中子學(xué)計(jì)算而設(shè)計(jì)的，能加工有效微觀和宏觀群截面，進(jìn)行靜態(tài)柵元和堆芯計(jì)算，能提供堆芯設(shè)計(jì)或?qū)嶒?yàn)分析所需的相關(guān)關(guān)鍵參數(shù)。公用數(shù)據(jù)庫為基于JENDL-3.3庫加工的超過300種核素的多群庫。SRAC 的共振處理程序PEACO 采用基于碰撞概率法直接求解共振能區(qū)超精細(xì)群慢化方程，能很好地處理PB-FHR的球型燃料元件所構(gòu)成的雙重非均勻系統(tǒng)。

本文主要采用SRAC的碰撞概率程序PIJ進(jìn)行通量計(jì)算，采用固定源模式。選擇HOMOSP程序求零幾何曲率下全能群B1方程解，并給出無窮增殖因數(shù)和均勻化能譜。可分辨共振區(qū)采用PEACO 處理，處理能區(qū)為快熱分界能到961.12eV，勒寬間隔為0.000 5，兩區(qū)微觀柵元碰撞率比采用透射截面模型。

1.2 MCNP程序

MCNP程序［9］是由美國洛斯阿拉莫斯國家實(shí)驗(yàn)室開發(fā)的一個(gè)通用的蒙特卡羅程序。MCNP程序采用連續(xù)點(diǎn)截面庫，具有強(qiáng)大的幾何處理能力，在核能新堆型的設(shè)計(jì)上發(fā)揮著很大作用。MCNP 通過使用各種記數(shù)卡統(tǒng)計(jì)大量粒子在介質(zhì)中發(fā)生的輸運(yùn)事件來給出宏觀結(jié)果。本文中，PB-FHR 的燃料與熔鹽等材料溫度非常高，運(yùn)用NJOY 程序?yàn)镸CNP 中涉及的材料做了相應(yīng)的溫度庫以使計(jì)算結(jié)果更加精準(zhǔn)。

MCNP在本文中主要是對球柵元內(nèi)各種反應(yīng)率和通量分別采用FMn卡和Fn 進(jìn)行統(tǒng)計(jì)，兩者相除即為所需的少群截面［10］。為方便與SRAC比較，給出的能譜為107群能譜。由于MCNP未提供統(tǒng)計(jì)出射中子能量及方向信息的計(jì)數(shù)卡［11］，不能直接求得散射矩陣及平均散射余弦，所以散射矩陣和擴(kuò)散系數(shù)未給出相應(yīng)結(jié)果。為統(tǒng)計(jì)所有包覆顆粒（TRISO）的燃料核的總體積，采用VOID 卡、SDEF 卡和F4卡統(tǒng)計(jì)燃料核總體積［9］。

1.3 計(jì)算模型

參考VSOP及MICROX-2軟件對球床模塊高溫氣冷堆（PBMR）球燃料元件建模［12］，PB-FHR 球柵元等效幾何模型如圖1所示，包層冷卻劑的厚度取決于包覆顆粒填充因子。MCNP建模時(shí)彌散的TRISO 采用簡立方柵格排布于燃料區(qū)，燃料球結(jié)構(gòu)也采用圖1模型。少群常數(shù)能群劃分參考文獻(xiàn)［13］給出的石墨慢化堆推薦的8 群劃分方法，取與其相鄰的SRAC的多群庫能群分界能（表1）。相關(guān)模型參數(shù)列于表2。

圖1 PB-FHR球柵元模型Fig.1 Pebble model of PB-FHR

2 結(jié)果與討論

2.1 無窮增殖因數(shù)與少群截面計(jì)算

PB-FHR球柵元無窮增殖因數(shù)與少群截面在冷態(tài)的計(jì)算結(jié)果列于表3。除2群的俘獲截面和1群的裂變截面外，各截面的相對偏差在1.45%以下，特別是包括可分辨共振能區(qū)的3群的俘獲截面相對偏差僅為0.04%。SRAC計(jì)算得到的2群的俘獲截面較MCNP計(jì)算得到的大12.76%，SRAC計(jì)算得到的1群的裂變截面較MCNP計(jì)算得到的小7.08%。SRAC計(jì)算的無窮增殖因數(shù)與MCNP計(jì)算的相比大0.05%。

表1 少群常數(shù)能群劃分Table 1 8-group structure

表2 PB-FHR球柵元主要參數(shù)Table 2 Parameters of lattice of PB-FHR

表3 冷態(tài)時(shí)SRAC計(jì)算的少群截面及無窮增殖因數(shù)與MCNP結(jié)果比較Table 3 Few-group cross sections and infinite multiplication factors comparison at 300K

表4 所列為PB-FHR 球柵元無窮增殖因數(shù)與少群截面在熱態(tài)的計(jì)算結(jié)果。除2群的俘獲截面和1群的裂變截面外，各截面的相對偏差在2.25%以下，3群的俘獲截面相對偏差僅為0.08%。SRAC 計(jì)算得到的2 群的俘獲截面較MCNP計(jì)算得到的大12.74%，SRAC 計(jì)算得到的1群的裂變截面較MCNP 計(jì)算得到的小7.07%。SRAC 計(jì)算的無窮增殖因數(shù)與MCNP計(jì)算的相比大0.04%?？煽闯?，溫度的改變對兩個(gè)計(jì)算的影響是一致的，計(jì)算結(jié)果的差異未擴(kuò)大。SRAC計(jì)算的球柵元的無窮增殖因數(shù)和少群截面與MCNP 結(jié)果基本符合，說明SRAC處理PB-FHR 球柵元雙重非均勻性的方法基本可取。

表4 熱態(tài)時(shí)SRAC計(jì)算的少群截面及無窮增殖因數(shù)與MCNP結(jié)果比較Table 4 Few-group cross sections and infinite multiplication factors comparison at 1 000K

2.2 FLiBe熔鹽

表5列出熱態(tài)下FLiBe熔鹽對1、2群俘獲截面的影響。當(dāng)計(jì)算模型中無FLiBe熔鹽時(shí)，可看到2群的俘獲截面相對偏差由12.74%降低到1.54%，當(dāng)把熔鹽中的F元素由C元素替代后，2群的俘獲截面相對偏差也降為1.08%。說明SRAC對于FLiBe熔鹽中的F 元素處理不當(dāng)是造成2 群俘獲截面相對偏差較大的原因。熔鹽FLiBe中的F 元素為中等質(zhì)量核素，相關(guān)共振能區(qū)位于平滑能量區(qū)（67.4keV～0.82 MeV），對應(yīng)于少群能群結(jié)構(gòu)中的2群，相關(guān)截面在該能區(qū)也較其他輕核的高出1 個(gè)量級，且共振峰更寬。SRAC 帶屏蔽因子的Bondarenko類型截面庫構(gòu)建時(shí)采用窄共振近似，該近似得到的用于并群的能譜能有效適用于重質(zhì)量核素，但對于在高能區(qū)有共振結(jié)構(gòu)的輕中等質(zhì)量核素不太精確。SRAC對F元素采用該近似就造成F 元素共振區(qū)的多群截面上有偏差。同樣，SRAC 對輕核（如C、O、Li和Be）的俘獲處理也不理想，這造成了表5 中無FLiBe熔鹽時(shí)1 群俘獲截面相對偏差達(dá)到23.98%，同樣造成了熔鹽中的F由C代替時(shí)1群俘獲截面相對偏差達(dá)到3.67%。

2.3 TRISO 尺寸

SRAC處理雙重非均勻性時(shí)，對于非PEACO處理能區(qū)（熱區(qū)及快區(qū)E＞961.12eV），燃料區(qū)核素濃度由帶石墨基質(zhì)的TRISO 柵元體積權(quán)重均勻化得到。圖2所示為TRISO幾何按表2中原始尺寸等比例變化后熱態(tài)下的無窮增殖因數(shù)、3群俘獲截面、1群裂變截面及8群裂變截面的計(jì)算結(jié)果。等比例變化保證燃料區(qū)材料體積權(quán)重均勻化后核素原子濃度不變，但相應(yīng)的TRISO 尺寸改變?yōu)楸?計(jì)算模型的0.1～2.5倍。零倍時(shí)即完全均勻分布，無雙重非均勻性效應(yīng)。為考慮MCNP建模中的切球影響，采用MCNP軟件以2.5cm 半徑的60 000 000粒子平面中子源及F4 卡統(tǒng)計(jì)了所有TRISO 的燃料核的總體積。

表5 熱態(tài)下FLiBe對俘獲截面的影響Table 5 Effect of FLiBe salt on cell-homogenized capture cross section at 1 000K

圖2 熱態(tài)下填充比10%時(shí)不同TRISO 大小對結(jié)果的影響Fig.2 Effects of TRISO size on results at 10% TRISO packing factor at 1 000K

TRISO 尺寸由零倍（無雙重非均勻性）增加至2.5倍的過程能很好體現(xiàn)雙重非均勻性對于PB-FHR 柵元計(jì)算的影響。由圖2 可知，TRISO 尺寸越大，含可分辨共振區(qū)的3群的俘獲截面越小，無窮增殖因數(shù)越大，2.5倍尺寸時(shí)無窮增殖因數(shù)較無雙重非均勻性時(shí)偏大約15%，其3群俘獲截面較無雙重非均勻性時(shí)偏小約40%，說明雙重非均勻性效應(yīng)對于PB-FHR球柵元計(jì)算影響很大。同時(shí)，SRAC 計(jì)算的無窮增殖因數(shù)與MCNP結(jié)果相比吻合很好，說明SRAC對雙重非均勻性的處理方法在TRISO尺寸改變時(shí)也能很好適用。圖2中，無窮增殖因數(shù)和3群俘獲截面的偏差隨TRISO 尺寸的波動與MCNP模型中燃料核總體積與理論值的偏差隨TRISO 尺寸的波動規(guī)律一致或相反，說明MCNP 模型在TRISO 尺寸較大時(shí)切球會對結(jié)果產(chǎn)生一定影響。

如圖2所示，SRAC 計(jì)算的1群和8群的裂變截面不會隨TRISO 尺寸變化，與零倍尺寸的結(jié)果相同，這是由于SRAC對于非PEACO 處理能區(qū)（熱區(qū)及快區(qū)E＞961.12eV）的計(jì)算模型為未考慮雙重非均勻性的體積權(quán)重均勻化的模型，即圖2中的零倍尺寸。由MCNP的計(jì)算結(jié)果可看到，當(dāng)TRISO 尺寸增大時(shí)，1 群的裂變截面偏差最大時(shí)約偏小17%，8群的裂變截面最大相對偏差約6.5%，MCNP 計(jì)算結(jié)果隨TRISO 尺寸的波動與其模型中燃料核總體積與理論值的偏差波動規(guī)律一致。1群裂變截面會隨TRISO 尺寸增大而增大，這是由于TRISO 燃料核尺寸增大時(shí)，產(chǎn)生于燃料核內(nèi)的裂變中子逃脫燃料核的概率變小，重核與慢化輕核比例沒有改變，導(dǎo)致裂變中子與慢化輕核碰撞而慢化到其他群的概率降低，裂變中子與重核的碰撞概率增大，1群的裂變截面將增大。8群裂變截面會隨TRISO 尺寸增大而降低是因?yàn)橛奢p核慢化的熱中子進(jìn)入到TRISO 燃料核會有空間自屏，尺寸越大空間自屏效應(yīng)越明顯。所以，1群和8群裂變截面會有較大偏差的原因?yàn)镾RAC在非PEACO 處理能區(qū)未考慮雙重非均勻性效應(yīng)。

3 結(jié)語

本文采用確定論軟件SRAC計(jì)算了PB-FHR球柵元的無窮增殖因數(shù)，少群均勻化總截面、俘獲截面和裂變截面，并使用連續(xù)能量蒙特卡羅軟件MCNP 驗(yàn)證與分析。SRAC 計(jì)算的球柵元的無窮增殖因數(shù)和少群截面與MCNP 結(jié)果基本符合，說明SRAC處理PB-FHR 球柵元雙重非均勻性的方法基本可取。SRAC快群數(shù)據(jù)庫對FLiBe熔鹽中的F 元素共振處理不理想，造成2群俘獲截面偏差較大。TRISO 尺寸變化的計(jì)算結(jié)果說明雙重非均勻性對于PB-FHR球柵元有很大影響，且1群和8群裂變截面會有較大偏差的原因?yàn)镾RAC在非PEACO 處理能區(qū)未考慮雙重非均勻性效應(yīng)。結(jié)果表明：SRAC程序計(jì)算結(jié)果與MCNP吻合，其適用于對PB-FHR 進(jìn)行少群截面加工。

［1］ FORSBERG C W，PETERSON P F，PICKARD P S.Molten-salt-cooled advanced high-temperature reactor for production of hydrogen and electricity［J］.Nuclear Technology，2003，144（1）：289-302.

［2］ FORSBERG C，PETERSON P F，SRIDHARAN K.Fluoride-salt-cooled high-temperature reactors（FHRs）for base-load and peak electricity，grid stabilization，and process heat［R］.Boston： Massachusetts Institute of Technology，2013.

［3］ HUDSON N H，OUGOUAG A M，RAHNEMA F，et al.A pebble bed reactor cross section methodology［J］.Annals of Nuclear Energy，2009，36（1）：1 138-1 150.

［4］ RüTTEN H J，HAAS K A.V S O P（99／05）computer code system for reactor physics and fuel cycle simulation［R］.Jülich：Forschungszentrum Jülich，2005.

［5］ BENDE E E，HOGENBIRK A H.Analytical calculation of the average dancoff factor for a fuel kernel in a pebble bed high-temperature reactor［J］.Nuclear Science and Engineering，1999，133（1）：147-162.

［6］ GOLUOGLU S，LANDERS N F，PETRIE L M，et al.CSAS：Control module for enhanced criticality safety analysis sequences［R］.Oak Ridge：Oak Ridge National Laboratory，2006.

［7］ 莊坤，曹良志，吳宏春，等.球床型固體燃料熔鹽堆計(jì)算軟件開發(fā)［J］.原子能科學(xué)技術(shù)，2013，47（增刊）：463-466.ZHUANG Kun，CAO Liangzhi， WU Hongchun，et al.Software development for pebble bed MSR calculation［J］.Atomic Energy Science and Technology，2013，47（Suppl.）：463-466（in Chinese）.

［8］ OKUMURA K，KUGO T，KANEKO K，et al.SRAC2006：A comprehensive neutronics calculation code system［R］.Ibarakiken：Japan Atomic Energy Agency，2007.

［9］ X-5 Monte Carlo Team.A general Monte Carlo N-particle transport code，Version 5［R］.Los Alamos： Los Alamos National Laboratory，2003.

［10］李滿倉，王侃，姚棟.基于連續(xù)能量蒙特卡羅方法的均勻化群常數(shù)計(jì)算［J］.核科學(xué)與工程，2012，32（4）：306-314.LI Mancang，WANG Kan，YAO Dong.Continuous energy Monte Carlo method based homogenization multi-group constants calculation［J］.Chinese Journal of Nuclear Science and Engineering，2012，32（4）：306-314（in Chinese）.

［11］ILAS G，RAHNEMA F.A Monte Carlo based nodal diffusion model for criticality analysis of spent fuel storage lattices［J］.Annals of Nuclear Energy，2003，30（1）：1 089-1 108.

［12］REITSMA F，STRYDOM G，DE HAAS J B M，et al.The PBMR steady-state and coupled kinetics core thermal-hydraulics benchmark test problems［J］.Nuclear Engineering and Design，2006，236（1）：657-668.

［13］MPHAHLELEA R， OUGOUAGB A M，IVANOV K N，et al.Spectral zone selection methodology for pebble bed reactors［J］.Annals of Nuclear Energy，2011，38（1）：80-87.