網(wǎng)格權(quán)窗減方差技術(shù)及其在聚變堆屏蔽分析中應(yīng)用研究

李新梅，鄭華慶，郝麗娟，宋 婧，胡麗琴,3，江 平

(1.合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，安徽合肥230002；2.中國科學(xué)院核能安全技術(shù)研究所，中國科學(xué)院中子輸運理論與輻射安全重點實驗室，安徽合肥230031；3中國科學(xué)技術(shù)大學(xué)，安徽合肥230027)

網(wǎng)格權(quán)窗減方差技術(shù)及其在聚變堆屏蔽分析中應(yīng)用研究

李新梅1,2，鄭華慶2，郝麗娟2，宋 婧2，胡麗琴2,3，江 平1

(1.合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，安徽合肥230002；2.中國科學(xué)院核能安全技術(shù)研究所，中國科學(xué)院中子輸運理論與輻射安全重點實驗室，安徽合肥230031；3中國科學(xué)技術(shù)大學(xué)，安徽合肥230027)

在聚變堆輻射屏蔽計算中，如何有效解決深穿透問題是近年來國際聚變輻射安全領(lǐng)域關(guān)注的焦點之一。針對該問題，本文研究了直角坐標(biāo)系與圓柱坐標(biāo)系下基于網(wǎng)格的權(quán)窗減方差技術(shù)。本文基于超級蒙卡核模擬軟件系統(tǒng)SuperMC實現(xiàn)了該方法，并選取減方差技巧的基準(zhǔn)例題進行測試與分析，初步得出“粗劃真空或密度很小的區(qū)域、細(xì)分密度大的區(qū)域”的網(wǎng)格劃分規(guī)律，能有效提高網(wǎng)格權(quán)窗計算效率?；谠撘?guī)律對聚變屏蔽基準(zhǔn)問題進行對比分析，新的網(wǎng)格劃分與原始網(wǎng)格劃分的計算效率相比，F(xiàn)OM因子提高了1.92倍。減方差技巧的基準(zhǔn)例題和聚變屏蔽基準(zhǔn)問題計算中，SuperMC通量計算結(jié)果與MCNP相比偏差均在0.5%以下，證明了本文中方法的正確性。

基于網(wǎng)格的權(quán)窗；聚變堆；屏蔽計算；減方差技巧；SuperMC

收斂速度慢是蒙特卡羅方法(簡稱蒙卡)在進行聚變反應(yīng)堆[1-5]模擬的瓶頸性限制之一，尤其是對于聚變堆屏蔽計算中深穿透問題。減方差技巧[6]是解決此類問題的主要方法。目前，現(xiàn)有的減方差技巧包括幾何分裂與賭、源偏倚、指數(shù)變換、權(quán)窗等，并在MCNP[6]、FLUKA[7]、Serpent[8]等蒙卡軟件中實現(xiàn)。

多年的實際應(yīng)用表明，目前權(quán)窗[9]是應(yīng)用最廣泛的減方差技巧。權(quán)窗分為基于柵元(cell-based)的權(quán)窗和基于網(wǎng)格(mesh-based)的權(quán)窗兩種模式?；跂旁臋?quán)窗是根據(jù)柵元的重要性設(shè)置權(quán)窗上下限參數(shù)，引導(dǎo)重要粒子朝著感興趣的區(qū)域輸運。為了有效利用柵元權(quán)窗需將尺寸過大柵元劃分成合適的小尺寸柵元。而基于網(wǎng)格的權(quán)窗是根據(jù)網(wǎng)格區(qū)域的重要性設(shè)置網(wǎng)格權(quán)窗上下限，無需劃分真實模型，只需定義劃分網(wǎng)格參數(shù)即可[10]。柵元權(quán)窗一般適用于簡單幾何問題，對于大型、復(fù)雜模型，其適用性較差，尤其針對帶重復(fù)結(jié)構(gòu)的復(fù)雜模型，而基于網(wǎng)格的權(quán)窗恰好可以彌補柵元權(quán)窗的這些缺點。

本文系統(tǒng)地研究了網(wǎng)格權(quán)窗的減方差技巧，并在中國科學(xué)院核能安全技術(shù)研究所·FDS團隊自主研發(fā)的超級蒙卡核模擬軟件系統(tǒng)SuperMC[11-17]的基礎(chǔ)上進行實現(xiàn)，并選取減方差基準(zhǔn)例題對程序進行測試與分析。

1 基于網(wǎng)格的權(quán)窗

目前，國外實現(xiàn)的網(wǎng)格權(quán)窗技巧主要是基于直角坐標(biāo)系和圓柱坐標(biāo)系上建立的[18]。直角坐標(biāo)系下的網(wǎng)格繪制比較方便、直觀。圓柱坐標(biāo)系下的網(wǎng)格可以旋轉(zhuǎn)任意角度進行劃分，比較復(fù)雜但使用廣泛。相比于直角坐標(biāo)系，圓柱坐標(biāo)系的網(wǎng)格權(quán)窗中需將粒子當(dāng)前位置坐標(biāo)轉(zhuǎn)換到直角坐標(biāo)系上，其轉(zhuǎn)換矩陣如下：

通過轉(zhuǎn)換后的坐標(biāo)或者實際坐標(biāo)找到相應(yīng)的網(wǎng)格權(quán)窗下限，繼而進行網(wǎng)格權(quán)窗處理，其基本原理如下：

權(quán)重為wgt的粒子輸運到某網(wǎng)格相應(yīng)的能量區(qū)間時：

1) 當(dāng)wgt

2) 當(dāng)wgt>WU，進行分裂，粒子分裂成int(wgt/WS+ε)個，權(quán)重均為WS；

3) 當(dāng)WL

網(wǎng)格權(quán)窗的使用能將粒子權(quán)重控制在合理的范圍內(nèi)，有效避免極高權(quán)重粒子對計數(shù)的擾動，最終提高計算效率。其流程見圖1。

圖1 網(wǎng)格權(quán)窗的流程圖Fig.1 Flow Chart of Mesh-based Weight Window

由圖1可知，在使用網(wǎng)格權(quán)窗時需考慮粒子當(dāng)前能量和材料下的平均自由程，取抽樣自由程、到下一邊界距離和平均自由程的最小值作為輸運長度。因此，使用網(wǎng)格權(quán)窗技巧修正粒子游動過程時需在下一邊界處、碰撞點以及平均自由程處進行處理。

2 程序測試與結(jié)果分析

本文選取參考文獻[19]中減方差技巧測試的基準(zhǔn)例題，并使用MCNP權(quán)窗產(chǎn)生器[8]生成的網(wǎng)格權(quán)窗下限進行測試。

2.1 長真空圓柱基準(zhǔn)測試?yán)}

2.1.1 模型描述

長真空圓柱測試?yán)}是減方差技巧基準(zhǔn)例題[19-20]，其幾何由半徑100cm、高2010cm的圓柱和半徑210cm、高10cm的圓柱疊加而成(見圖2)。其中，距底面180cm厚度的柵元是密度為2.03g/cm3的重混凝土；180～2000cm為真空柵元；2000～2010cm是密度為0.0203g/cm3的輕混凝土層。能量50% 2MeV與50% 14MeV、方向各向同性的中子點源位于幾何體底層重混凝土中，距底面1.0×10-6cm，計數(shù)區(qū)域為圖中最上層非真空柵元，類型為柵元體通量。由于模型所有柵元均是圓柱型的，因此采用圓柱坐標(biāo)系下的網(wǎng)格劃分。表1中給出了圓柱軸向上不同網(wǎng)格劃分的參數(shù)值。

圖2 長真空圓柱幾何示意圖Fig.2 Figure of Vacuum Cylinder Model

表1 網(wǎng)格劃分方法

續(xù)表

2.1.2 測試結(jié)果與分析

表2中給出了不使用網(wǎng)格權(quán)窗和使用網(wǎng)格權(quán)窗1、2、3的計算結(jié)果，SuperMC與MCNP的相對偏差均小于0.5%。表2中在不使用網(wǎng)格權(quán)窗條件下無計算結(jié)果，這是因為模型中間有很長一段真空材料，屬于典型的深穿透問題，如果不使用減方差技巧，無法得到穩(wěn)定收斂的計算結(jié)果。

表2 長真空圓柱計算結(jié)果

表3中給出了使用網(wǎng)格權(quán)窗1的品質(zhì)因子FOM為7.60，網(wǎng)格權(quán)窗2的FOM為7.54，網(wǎng)格權(quán)窗3的FOM為7.60。比較不使用與使用網(wǎng)格權(quán)窗的計算結(jié)果發(fā)現(xiàn)，有效利用網(wǎng)格權(quán)窗能夠成倍提高計算效率。

表3 長真空圓柱計算效率比較

注：① FOM=1/(σ2T)，σ為標(biāo)準(zhǔn)差，T為計算時間/min。

綜合分析上述三組網(wǎng)格權(quán)窗劃分對計算結(jié)果的影響，可以總結(jié)得到如下網(wǎng)格設(shè)置的規(guī)律：對于真空或者密度很小的材料區(qū)域劃分盡量控制在一個粗網(wǎng)格內(nèi)，無需再細(xì)分；對于材料密度大、密度大的材料區(qū)域劃分則相反。

2.2 聚變屏蔽基準(zhǔn)例題測試

2.2.1 模型描述

聚變屏蔽基準(zhǔn)問題[19]的整個模型在899.16cm×690.85cm×678.18cm的長方體內(nèi)，幾何和材料如圖3所示。中子點源位于中間真空柵元內(nèi)，其坐標(biāo)為(-356.87，232.02，157.40)，能量為14MeV，方向各向同性。在源外部有55.88cm厚的含鐵屏蔽層，計數(shù)區(qū)域為屏蔽層后的柵元，類型為柵元體通量。此模型采用兩種方法劃分：參考文獻[19]中網(wǎng)格劃分(網(wǎng)格1)、根據(jù)上述網(wǎng)格劃分規(guī)律對整個模型進行劃分(網(wǎng)格2)。

圖3 聚變屏蔽模型Fig.3 Fusion Shielding Model

網(wǎng)格2(以y軸方向為例)劃分方法如下：-29.21～0.00cm劃分5份，0.00～178.85cm劃分5份，178.50～196.10cm劃分4份，196.15～223.79cm劃分6份，223.79～259.39cm劃分5份，259.39～279.75cm劃分4份，279.75～407.33cm劃分1份，407.33～661.64cm劃分1份。

2.2.2 測試結(jié)果與分析

表4中給出不使用網(wǎng)格權(quán)窗和使用網(wǎng)格權(quán)窗1、2時的計算結(jié)果，SuperMC與MCNP計算結(jié)果的相對偏差均小于0.5%。可見，無論圓柱坐標(biāo)系下還是直角坐標(biāo)系下使用網(wǎng)格權(quán)窗的計算結(jié)果均與MCNP吻合程度良好。

表4 聚變熱屏蔽計算結(jié)果

續(xù)表

表5中給出不使用網(wǎng)格權(quán)窗的品質(zhì)因子FOM為7.15，使用網(wǎng)格權(quán)窗1后FOM為123.70，計算效率是不使用網(wǎng)格權(quán)窗時的17.30倍；使用網(wǎng)格權(quán)窗2后FOM為237.059，計算效率是不使用網(wǎng)格權(quán)窗時的33.16倍，是網(wǎng)格權(quán)窗1計算效率的1.92倍?？梢姡偨Y(jié)得出的網(wǎng)格劃分規(guī)律在聚變堆屏蔽分析中是有效的。

表5 聚變熱屏蔽計算效率比較

注：① FOM=1/(σ2T)，σ為標(biāo)準(zhǔn)差，T為計算時間/min。計算效率為FOM的比值。

3 總結(jié)

本文對蒙特卡羅粒子輸運模擬中直角坐標(biāo)系與圓柱坐標(biāo)系下的網(wǎng)格權(quán)窗減方差技巧進行了研究。利用減方差技巧的基準(zhǔn)例題進行測試與分析，初步得出“粗劃真空或密度很小的區(qū)域、細(xì)分密度大的區(qū)域”的網(wǎng)格劃分規(guī)律，能有效提高網(wǎng)格權(quán)窗的計算效率?；谠撘?guī)律對聚變屏蔽基準(zhǔn)問題進行對比分析，新的網(wǎng)格劃分與原始網(wǎng)格劃分的計算效率相比FOM因子提高了1.92倍。SuperMC通量計算結(jié)果與MCNP相比偏差均在0.5%以下，證明了本文中方法的正確性。本文網(wǎng)格劃分方法是手工劃分的，可進一步研究網(wǎng)格權(quán)窗的網(wǎng)格自動劃分方法。

致謝

感謝FDS團隊成員對本文工作的幫助與支持。

[1] Y. Wu, FDS Team. CAD-Based Interface Programs for Fusion Neutron Transport Simulation[J]. Fusion Engineering and Design, 2009, 84 (7-11), 1987-1992.

[2] Y. Wu, J. Jiang, M. Wang, M. Jin, et al. A Fusion-Driven Subcritical System Concept Based on Viable Technologies[J].Nuclear Fusion, 2011, 51(10)：103036.

[3] S. Zheng, M. Chen, J. Li, et al. Neutronics Analysis for the Test Blanket Modules Proposed for EAST and ITER[J]. Nuclear Fusion, 47 (2007) 1053-1056.

[4] Y.Chen, Y. Wu. Conceptual Study on High Performance Blanket in a Spherical Tokamak Fusion-DrivenTransmuter[J]. Fusion Engineering and Design, 2000, 49-50：507-512.

[5] Y. Wu, S. Zheng, X. Zhu, et al. Conceptual Design of the Fusion-Driven Subcritical System FDS-I[J]. Fusion Engineering and Design, 2006, 81, Part B： 1305-1311.

[6] X-5 Monte 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.

[7] A. Ferrari, P. R. Sala, A. Fasso, et al. Fluka：A Multi-paricle Transport Code[R]. CERN-2005-010, Oct 12,2005.

[8] J. L. Serpent -A Continuous Energy Monte Carlo Reactor Physics Burn up Calculation Code User’s Manual [M]. March 6, 2013： 131-133.

[9] T. E. Booth. Genesis of the Weight Window and the Weight Window Generator in MCNP-A Personal History[J]. LA-UR-06-5807. 9, July 28, 2006.

[10] 姜宏宇, 王汝贍. 網(wǎng)格權(quán)重窗方法及其在煤漿系統(tǒng)研究中的應(yīng)用[J]. 核電子學(xué)與探測技術(shù), 1998,18(3)： 201-205.

[11] J. Song, G. Sun, Z. Chen, et al. Benchmarking of CAD-based SuperMC with ITER Benchmark Model[J]. Fusion Engineering and Design, DOI： 10.1016/j. fusengdes.2014.5.3.

[12] Y. Wu, J. Song, H. Zheng, et al. CAD-Based Monte Carlo Program for Integrated Simulation of Nuclear System SuperMC[J]. Annals of Nuclear Energy, 82,161-168(2015).

[13] Z.P. Chen, H.Q. Zheng, G.Y. Sun, et al. Preliminary Study on CAD-based Method of Characteristics for Neutron Transport Calculation[J]. Chinese Physics C. (2014) 38 (5), 058201-1 058201-7.

[14] 吳宜燦, 李靜驚, 李瑩 等. 大型集成多功能中子學(xué)計算與分析系統(tǒng)VisualBUS的研究與發(fā)展[J]. 核科學(xué)與工程，2007, 27(4)： 365-373.

[15] 吳宜燦, 胡麗琴, 龍鵬程，等.先進核能軟件發(fā)展與核信息學(xué)實踐[M]. 中國科研信息化藍(lán)皮書, 2013. 北京：科學(xué)出版社, 2013，12： 232-244.

[16] Y. Wu, FDS Team. Design Analysis of the China Dual-Functional Lithium Lead (DFLL) Test Blanket Module in ITER[J]. Fusion Engineering and Design, 2007, 82： 1893-1903.

[17] Y. Wu, X. Zhu, S. Zheng, et al. Neutronics Analysis of Dual-Cooled Waste Transmutation Blanket for the FDS[J]. Fusion Engineering and Design, 2002, 63-64： 133-138.

[18] J. S. Hendricks. Superimposed Mesh Plotting in MCNP [R]. Los Alamos National Laboratory LA-UR-01-1033, September 2001.

[19] C. V. Culbertson, J. S. Hendricks. An Assessment of the MCNP4C Weight Window[R]. Los Alamos National Laboratory LA-13668, December 1999.

[20] T. E. Booth. A Sample Problem for Variance Reduction in MCNP[R]. Los Alamos National Laboratory, LA-10363-MS, October 1985.

Mesh-basedWeightWindowVarianceReductionTechniquesanditsAppliedResearchonFusionReactorShieldingAnalysis

LIXin-mei1,2，ZHENGHua-qing2，HAOLi-juan2，SONGJing2，HuLi-qin2,3，JIANGPing1

(1.Hefei University of Technology, Hefei, Anhui, 230002, China；2.Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei, Anhui, 230031, China；3.University of Science and Technology of China, Hefei, Anhui, 230027, China)

In the field of fusion radiation safety, it is a focus to solve deep penetration problem effectively for the radiation shielding calculation in recent years. In this paper, mesh-based weight window variance reduction techniques under Cartesian and cylindrical coordinate systems were studied and implemented in Super Monte Carlo Program for Nuclear and Radiation Simulation(SuperMC). During the verification with a variance reduction benchmark problem, a regulation of rough mesh division in vacuum or low density regions and fine division in high density regions for improving the computing efficiency was obtained. The fusion shielding benchmark problem was analyzed with a new mesh based on the regulation and the FOM increased by 1.92 times with the new mesh. The difference between the calculation results of SuperMC and MCNP was within 0.5% for the benchmark problem and fusion shielding problem. The accuracy was demonstrated.

Mesh-based weight window；Fusion reactor；Shielding calculation；Variance reduction；SuperMC

2016-12-20

國家ITER 973計劃(2014GB112001)；中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項(XDA03040000)；國家自然科學(xué)基金(11605233、11305203)

李新梅(1988—)，女，安徽人，碩士研究生，現(xiàn)主要從事蒙特卡羅粒子輸運計算研究工作

胡麗琴：liqin.hu@fds.org.cn

TL7

：A

：0258-0918(2017)04-0577-06