蒙特卡羅粒子輸運權窗減方差方法研究及在ITER屏蔽分析中應用

鄭華慶，宋 婧，郝麗娟，孫光耀，胡麗琴，FDS團隊

（中國科學院核能安全技術研究所，安徽 合肥230031）

蒙特卡羅方法具備處理復雜幾何和模擬精細的物理過程的能力，方法本身和程序結構的簡單性，因而在解決粒子輸運問題方面應用越來越廣泛。屏蔽分析是聚變裝置物理設計與分析的關鍵環節，但蒙特卡羅方法在處理屏蔽問題時，方法本身收斂速度比較慢是一個關鍵的問題，因此往往需要很長的計算時間，在深穿透情況下，這個問題尤其突出。減方差方法可以加快蒙特卡羅方法的收斂速度，因此在應用蒙特卡羅方法解決屏蔽計算問題的過程中，有效地使用減方差方法是必不可少的環節。

現有的蒙特卡羅粒子輸運程序中有四類減方差方法［1］：截斷方法（能量截斷和時間截斷）、總體控制方法（幾何分裂與輪盤賭，能量分裂與輪盤賭，時間分裂與輪盤賭，權截斷和權窗）、修正抽樣方法（指數變換，隱俘獲，強迫碰撞，源偏倚）和部分確定性方法（點探測器，DXTRAN和相關抽樣）。

Booth和 Hendrick提出的權窗方法［2］，在多年的實際應用中逐漸被證實是目前最為有效和最為通用的減方差方法之一。

超級蒙特卡羅計算軟件SuperMC［3］是由FDS團隊自主研發，采用蒙特卡羅方法并耦合確定論方法，定位于基于先進計算機技術實現輸運（穩態粒子輸運、中子動力學）、燃耗與活化（同位素燃耗、材料活化、停機劑量）、多物理耦合（熱工水力學、燃料性能、結構力學）過程的高效能精細模擬，集建模、計算、可視化分析于一體，可廣泛應用于反應堆物理、輻射屏蔽、醫學物理、核探測、高能物理等領域。

本文在中國科學院核能安全技術研究所·FDS團隊前期研究工作的基礎［4－12］上，將權窗減方差方法在SuperMC中實現，并通過基準例題的校驗和國際熱核聚變實驗堆ITER屏蔽分析的應用，證明了方法和程序的正確性和有效性。

1 權窗減方差方法的設計與實現

權窗是相空間（空間－能量、空間－時間）上的分裂與輪盤賭，本文將權窗減方差方法在SuperMC實現的方法如下：

用戶首先給出權窗基本參數CU（權窗上限參數）、CS（賭存活參數）、MXSPIN（最大分裂參數／最大賭存活參數）。一般默認的參數設置是：CU＝5，CS＝3，MXSPIN＝5。

對于每一個相空間柵元有一個權重下限WL（用戶給出）和上限 WU（WU＝WL×CU），當粒子的權重WGT小于下限WL時，發生輪盤賭：賭贏則粒子權重變為 min（CS×WL，WGT×MXSPIN），賭輸則該粒子被殺死；當粒子的權重WGT大于上限WU時，發生分裂：分裂出npa（npa＝min（int（WGT／WU ＋1），MXSPIN）個權重為 WGT／npa的全同粒子；當粒子權重WGT在WU和WL之間時，粒子不發生任何變化。低權重粒子發生輪盤賭，使得不會在權重無意義的粒子上浪費時間；權重高于權窗，粒子發生分裂能有效避免極高權重的粒子對計數的擾動，其中流程如下圖所示：

圖1 權窗減方差流程圖Fig.1 Flow chart of weight windows variance reduction

在SuperMC的用戶圖形界面（GUI）中，如圖2所示，用戶可以通過可視化幾何交互的方式進行權窗參數的設置。在處理復雜的幾何模型時，根據放射源和計數區的位置，用戶可以方便地點擊每個柵元，設置相應柵元的權窗參數。其設置的原則是：粒子由放射源所在柵元開始到計數區為止，考慮其可能經過的所有的柵元，在這些柵元上根據粒子通過該柵元的可能性來設置這些柵元的權窗下限，保證粒子通過合理的引導進入計數區。對于其他對計數區貢獻可以忽略的柵元，可以將權窗下限設為0，讓進入該區域的粒子進行權截斷。

圖2 SuperMC幾何建模功能用戶界面Fig.2 Interface of geometry modeling of SuperMC

2 程序校驗及在ITER屏蔽分析中應用

2.1 混凝土基準校驗例題

2.1.1 模型描述

混凝土是聚變堆典型的屏蔽體材料。計算例題選自參考文獻［13］，其基本的信息如下：200cm厚的混凝土屏蔽層，材料密度為ρ＝2.03g／cm3，被分成20層，每層是10cm，如圖3所示。放射源是位于（X＝0cm；Y＝0cm；Z＝1.0×10－6cm），單能（14MeV）、單方向（0，0，1）的中子點源。計數柵元是混凝土圓柱最上層的柵元，計數類型為體通量計算。三組權窗的柵元權重下限設置參見表1。

圖3 混凝土基準測試例題幾何示意圖Fig.3 Geometry of concrete benchmark case

表1 不同權窗的柵元權重下限設置Table 1 Weight low limit in each cell of different weight windows

2.1.2 測試結果

SuperMC與MCNP使用相同的數據庫HENDL［14］，計算結果如表2所示，在沒有添加權窗、添加權窗1、添加權窗2和添加權窗3，這三種計算條件下，SuperMC通量計算結果與參考程序MCNP的計算結果基本吻合，偏差均小于0.5%。

計算效率的比較如表3所示，在沒有添加權窗的情況下，計算品質因子FOM（Figure of Merit）僅為0.593，添加權窗1的情況FOM為4.109，計算效率提升了近6倍；添加權窗2的情況FOM為66.077，計算效率提升了110倍；添加權窗3的情況FOM為90.240，計算效率提升了150倍。

表2 混凝土測試例題通量計算結果Table 2 Flux results of concrete benchmark case

表3 混凝土測試例題計算效率比較Table 3 Efficiency comparison of FOM of concrete benchmark case

2.2 ITER－T426屏蔽分析應用例題

2.2.1 模型描述

國際熱核聚變實驗堆ITER是目前全球規模最大的國際科研合作項目之一，ITER裝置是一個能產生大規模核聚變反應的超導托卡馬克。本文選用ITER－T426例題［15］作為應用例題。該例題基于 FNG（Frascati neutron generator）裝置，設計一個由不銹鋼和有機玻璃組成的屏蔽結構，模擬ITER真空室。本文選用ITER－T426例題，如圖4所示。源是位于底層 柵 元 （X＝0cm；Y＝ 0.001cm；Z＝0cm）處，單能（14MeV）、單方向（0，1，0）的中子點源。計數區域是中心點位于（X＝0cm；Y＝74.93cm；Z＝0cm），半徑為1.5cm，高度為5cm的不銹鋼圓柱，計數類型為體通量計算。

圖4 ITER－T426例題在SuperMC的顯示圖Fig.4 Geometry view of ITER－T426in SuperMC

2.2.2 測試結果

SuperMC與MCNP使用相同的核數據庫HENDL，計算結果如表4所示，在沒有添加權窗、添加權窗兩種計算條件下，SuperMC通量計算結果與參考程序MCNP的結果基本吻合，偏差均小于0.2%。

計算效率的比較如表5所示，在沒有添加權窗的情況下，計算品質因子FOM僅為59.707，添加權窗的情況FOM 為905.199，計算效率提升了14倍。

表4 ITER－T426例題計算結果Table 4 Flux results of ITER－T426case

表5 ITER－T426例題計算效率比較Table 5 Efficiency comparison of FOM of ITER－T426case

3 總結

本文對蒙特卡羅粒子輸運模擬中權窗減方差方法進行了研究：給出了關鍵參數設置的原則和程序實現的思路，并基于超級蒙特卡羅計算軟件SuperMC進行實現。通過混凝土基準例題的校驗和ITER屏蔽分析的應用，證明了方法和程序的正確性和有效性，通過合理設置權窗參數可以提高計算效率。

［1］ X－5Monte Carlo Team，MCNP－A General Monte Carlo N－Particle Transport Code，Version 5［R］.Los Alamos National Laboratory Report LA－UR－03－1987，April 24，2003.

［2］ Booth T E，Hendricks J S.Importance Estimation in Forward Monte Carlo Calculations ［J］. Nuclear Technology Fusion，1984，5：90－100.

［3］ 孫光耀，宋婧，鄭華慶，等 .超級蒙特卡羅軟件SuperMC2.0中子輸運計算校驗［J］.原子能科學技術，2013，47（s2）：520－525.

［4］ 吳宜燦，胡麗琴，龍鵬程，等 .先進核能軟件發展與核信息學實踐［J］.中國科研信息化藍皮書2013，2013.12：232－244.

［5］ Li Y，Lu L，Ding A P，et al.Benchmarking of MCAM 4.0with the ITER 3DModel［J］.Fusion Engineering and Design，2007，82（15－24）：2861－2866.

［6］ 吳宜燦，李靜驚，李瑩，等 .大型集成多功能中子學計算與分析系統VisualBUS的研究與發展［J］.核科學與工程，2007，27（5）：72－85.

［7］ Wu Y C，FDS Team.CAD－based Interface Programs for Fusion Neutron Transport Simulation ［J］.Fusion Engineering and Design，2009，84（7－11）：1987－1992.

［8］ 吳宜燦，李瑩，盧磊，等 .蒙特卡羅粒子輸運計算自動建模程序系統的研究與發展［J］.核科學與工程，2006，26（1）：20－27.

［9］ Wu Y C，Xie Z S，Fischer U.A Discrete Ordinates Nodal Method for One－dimensional Neutron Transport Calculation in Curvilinear Geometries［J］.Nuclear Science and Engineering，1999，133（3）：350－357.

［10］ Wu Y C，FDS Team.Conceptual Design Activities of FDS Series Fusion Power Plants in China［J］.Fusion Engineering and Design，2006，81（23－24）：2713－2718.

［11］ Wu Y C，FDS Team.Design Analysis of the China Dual－functional Lithium Lead （DFLL）Test Blanket Module in ITER［J］.Fusion Engineering and Design，2007，82（15－24）：1893－1903.

［12］ Wu Y C，FDS Team.Overview of Liquid Lithium Lead Breeder Blanket Program in China ［J］. Fusion Engineering and Design，2011，86（9－11）：2343－2346.

［13］ Booth T E，MCNP Variance Reduction Examples［R］.Los Alamos National Laboratory，Diagnostic Applications Group X－5，Mail Stop F663，December 23，2004.

［14］ Zou J，He Z Z，Zeng Q，et al.Development and Testing of Multigroup Library with Correction of Self－shielding Effects in Fusion－Fission Hybrid Reactor［J］.Fusion Engineering and Design，2010，85（7－9）：1587－1590.

［15］ Batistoni P.Experimental Validation of Shutdown Dose Rates Calculations Inside ITER Cryostat［J］.Fusion Engineering and Design，2001，58－59：613－616.