氟鹽冷卻球床堆球柵元少群截面計算

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

（1.中國科學院 上海應用物理研究所，上海 201800；2.中國科學院大學，北京 100049；3.中國科學院 核輻射與核能技術重點實驗室，上海 201800）

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

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

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

1 計算程序及方法

1.1 確定論軟件

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

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

1.2 MCNP程序

MCNP程序［9］是由美國洛斯阿拉莫斯國家實驗室開發的一個通用的蒙特卡羅程序。MCNP程序采用連續點截面庫，具有強大的幾何處理能力，在核能新堆型的設計上發揮著很大作用。MCNP 通過使用各種記數卡統計大量粒子在介質中發生的輸運事件來給出宏觀結果。本文中，PB-FHR 的燃料與熔鹽等材料溫度非常高，運用NJOY 程序為MCNP 中涉及的材料做了相應的溫度庫以使計算結果更加精準。

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

1.3 計算模型

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

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

2 結果與討論

2.1 無窮增殖因數與少群截面計算

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

表1 少群常數能群劃分Table 1 8-group structure

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

表3 冷態時SRAC計算的少群截面及無窮增殖因數與MCNP結果比較Table 3 Few-group cross sections and infinite multiplication factors comparison at 300K

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

表4 熱態時SRAC計算的少群截面及無窮增殖因數與MCNP結果比較Table 4 Few-group cross sections and infinite multiplication factors comparison at 1 000K

2.2 FLiBe熔鹽

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

2.3 TRISO 尺寸

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

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

圖2 熱態下填充比10%時不同TRISO 大小對結果的影響Fig.2 Effects of TRISO size on results at 10% TRISO packing factor at 1 000K

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

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

3 結語

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

［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］.原子能科學技術，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］李滿倉，王侃，姚棟.基于連續能量蒙特卡羅方法的均勻化群常數計算［J］.核科學與工程，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.