蒙特卡羅粒子輸運(yùn)權(quán)窗減方差方法研究及在ITER屏蔽分析中應(yīng)用

鄭華慶，宋 婧，郝麗娟，孫光耀，胡麗琴，F(xiàn)DS團(tuán)隊(duì)

（中國(guó)科學(xué)院核能安全技術(shù)研究所，安徽 合肥230031）

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

現(xiàn)有的蒙特卡羅粒子輸運(yùn)程序中有四類(lèi)減方差方法［1］：截?cái)喾椒ǎ芰拷財(cái)嗪蜁r(shí)間截?cái)啵⒖傮w控制方法（幾何分裂與輪盤(pán)賭，能量分裂與輪盤(pán)賭，時(shí)間分裂與輪盤(pán)賭，權(quán)截?cái)嗪蜋?quán)窗）、修正抽樣方法（指數(shù)變換，隱俘獲，強(qiáng)迫碰撞，源偏倚）和部分確定性方法（點(diǎn)探測(cè)器，DXTRAN和相關(guān)抽樣）。

Booth和 Hendrick提出的權(quán)窗方法［2］，在多年的實(shí)際應(yīng)用中逐漸被證實(shí)是目前最為有效和最為通用的減方差方法之一。

超級(jí)蒙特卡羅計(jì)算軟件SuperMC［3］是由FDS團(tuán)隊(duì)自主研發(fā)，采用蒙特卡羅方法并耦合確定論方法，定位于基于先進(jìn)計(jì)算機(jī)技術(shù)實(shí)現(xiàn)輸運(yùn)（穩(wěn)態(tài)粒子輸運(yùn)、中子動(dòng)力學(xué)）、燃耗與活化（同位素燃耗、材料活化、停機(jī)劑量）、多物理耦合（熱工水力學(xué)、燃料性能、結(jié)構(gòu)力學(xué)）過(guò)程的高效能精細(xì)模擬，集建模、計(jì)算、可視化分析于一體，可廣泛應(yīng)用于反應(yīng)堆物理、輻射屏蔽、醫(yī)學(xué)物理、核探測(cè)、高能物理等領(lǐng)域。

本文在中國(guó)科學(xué)院核能安全技術(shù)研究所·FDS團(tuán)隊(duì)前期研究工作的基礎(chǔ)［4－12］上，將權(quán)窗減方差方法在SuperMC中實(shí)現(xiàn)，并通過(guò)基準(zhǔn)例題的校驗(yàn)和國(guó)際熱核聚變實(shí)驗(yàn)堆ITER屏蔽分析的應(yīng)用，證明了方法和程序的正確性和有效性。

1 權(quán)窗減方差方法的設(shè)計(jì)與實(shí)現(xiàn)

權(quán)窗是相空間（空間－能量、空間－時(shí)間）上的分裂與輪盤(pán)賭，本文將權(quán)窗減方差方法在SuperMC實(shí)現(xiàn)的方法如下：

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

對(duì)于每一個(gè)相空間柵元有一個(gè)權(quán)重下限WL（用戶(hù)給出）和上限 WU（WU＝WL×CU），當(dāng)粒子的權(quán)重WGT小于下限WL時(shí)，發(fā)生輪盤(pán)賭：賭贏則粒子權(quán)重變?yōu)?min（CS×WL，WGT×MXSPIN），賭輸則該粒子被殺死；當(dāng)粒子的權(quán)重WGT大于上限WU時(shí)，發(fā)生分裂：分裂出npa（npa＝min（int（WGT／WU ＋1），MXSPIN）個(gè)權(quán)重為 WGT／npa的全同粒子；當(dāng)粒子權(quán)重WGT在WU和WL之間時(shí)，粒子不發(fā)生任何變化。低權(quán)重粒子發(fā)生輪盤(pán)賭，使得不會(huì)在權(quán)重?zé)o意義的粒子上浪費(fèi)時(shí)間；權(quán)重高于權(quán)窗，粒子發(fā)生分裂能有效避免極高權(quán)重的粒子對(duì)計(jì)數(shù)的擾動(dòng)，其中流程如下圖所示：

圖1 權(quán)窗減方差流程圖Fig.1 Flow chart of weight windows variance reduction

在SuperMC的用戶(hù)圖形界面（GUI）中，如圖2所示，用戶(hù)可以通過(guò)可視化幾何交互的方式進(jìn)行權(quán)窗參數(shù)的設(shè)置。在處理復(fù)雜的幾何模型時(shí)，根據(jù)放射源和計(jì)數(shù)區(qū)的位置，用戶(hù)可以方便地點(diǎn)擊每個(gè)柵元，設(shè)置相應(yīng)柵元的權(quán)窗參數(shù)。其設(shè)置的原則是：粒子由放射源所在柵元開(kāi)始到計(jì)數(shù)區(qū)為止，考慮其可能經(jīng)過(guò)的所有的柵元，在這些柵元上根據(jù)粒子通過(guò)該柵元的可能性來(lái)設(shè)置這些柵元的權(quán)窗下限，保證粒子通過(guò)合理的引導(dǎo)進(jìn)入計(jì)數(shù)區(qū)。對(duì)于其他對(duì)計(jì)數(shù)區(qū)貢獻(xiàn)可以忽略的柵元，可以將權(quán)窗下限設(shè)為0，讓進(jìn)入該區(qū)域的粒子進(jìn)行權(quán)截?cái)唷?/p>

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

2 程序校驗(yàn)及在ITER屏蔽分析中應(yīng)用

2.1 混凝土基準(zhǔn)校驗(yàn)例題

2.1.1 模型描述

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

圖3 混凝土基準(zhǔn)測(cè)試?yán)}幾何示意圖Fig.3 Geometry of concrete benchmark case

表1 不同權(quán)窗的柵元權(quán)重下限設(shè)置Table 1 Weight low limit in each cell of different weight windows

2.1.2 測(cè)試結(jié)果

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

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

表2 混凝土測(cè)試?yán)}通量計(jì)算結(jié)果Table 2 Flux results of concrete benchmark case

表3 混凝土測(cè)試?yán)}計(jì)算效率比較Table 3 Efficiency comparison of FOM of concrete benchmark case

2.2 ITER－T426屏蔽分析應(yīng)用例題

2.2.1 模型描述

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

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

2.2.2 測(cè)試結(jié)果

SuperMC與MCNP使用相同的核數(shù)據(jù)庫(kù)HENDL，計(jì)算結(jié)果如表4所示，在沒(méi)有添加權(quán)窗、添加權(quán)窗兩種計(jì)算條件下，SuperMC通量計(jì)算結(jié)果與參考程序MCNP的結(jié)果基本吻合，偏差均小于0.2%。

計(jì)算效率的比較如表5所示，在沒(méi)有添加權(quán)窗的情況下，計(jì)算品質(zhì)因子FOM僅為59.707，添加權(quán)窗的情況FOM 為905.199，計(jì)算效率提升了14倍。

表4 ITER－T426例題計(jì)算結(jié)果Table 4 Flux results of ITER－T426case

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

3 總結(jié)

本文對(duì)蒙特卡羅粒子輸運(yùn)模擬中權(quán)窗減方差方法進(jìn)行了研究：給出了關(guān)鍵參數(shù)設(shè)置的原則和程序?qū)崿F(xiàn)的思路，并基于超級(jí)蒙特卡羅計(jì)算軟件SuperMC進(jìn)行實(shí)現(xiàn)。通過(guò)混凝土基準(zhǔn)例題的校驗(yàn)和ITER屏蔽分析的應(yīng)用，證明了方法和程序的正確性和有效性，通過(guò)合理設(shè)置權(quán)窗參數(shù)可以提高計(jì)算效率。

［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］ 孫光耀，宋婧，鄭華慶，等 .超級(jí)蒙特卡羅軟件SuperMC2.0中子輸運(yùn)計(jì)算校驗(yàn)［J］.原子能科學(xué)技術(shù)，2013，47（s2）：520－525.

［4］ 吳宜燦，胡麗琴，龍鵬程，等 .先進(jìn)核能軟件發(fā)展與核信息學(xué)實(shí)踐［J］.中國(guó)科研信息化藍(lán)皮書(shū)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］ 吳宜燦，李靜驚，李瑩，等 .大型集成多功能中子學(xué)計(jì)算與分析系統(tǒng)VisualBUS的研究與發(fā)展［J］.核科學(xué)與工程，2007，27（5）：72－85.

［7］ Wu Y C，F(xiàn)DS Team.CAD－based Interface Programs for Fusion Neutron Transport Simulation ［J］.Fusion Engineering and Design，2009，84（7－11）：1987－1992.

［8］ 吳宜燦，李瑩，盧磊，等 .蒙特卡羅粒子輸運(yùn)計(jì)算自動(dòng)建模程序系統(tǒng)的研究與發(fā)展［J］.核科學(xué)與工程，2006，26（1）：20－27.

［9］ Wu Y C，Xie Z S，F(xiàn)ischer 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，F(xiàn)DS 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，F(xiàn)DS 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，F(xiàn)DS 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.