蒙特卡羅自動建模中的空腔簡化算法研究

俞盛朋，程夢云，龍鵬程，胡麗琴

（中國科學院 核能安全技術(shù)研究所，中國科學院 中子輸運理論與輻射安全重點實驗室，安徽 合肥 230031）

蒙特卡羅粒子輸運計算方法通過在幾何空間中追蹤粒子與材料之間發(fā)生的物理反應來近乎真實地模擬粒子的運動規(guī)律。為能精確高效進行蒙特卡羅粒子輸運計算，需準確簡潔描述蒙特卡羅計算的幾何模型。傳統(tǒng)的手工蒙特卡羅建模方法費時費力，并可能引入幾何錯誤，限制了蒙特卡羅程序的使用。為解決蒙特卡羅的建模困難，國內(nèi)外已研發(fā)了一系列基于CAD的蒙特卡羅自動建模工具，如德國KIT研發(fā)的McCad［1］、日本JAEA研發(fā)的GEOMIT［2］和FDS團隊研發(fā)的SuperMC／MCAM［3-4］。其中SuperMC／MCAM（multi-physics coupling analysis modeling program）是FDS團隊研發(fā)的超級蒙特卡羅核計算仿真軟件系統(tǒng)SuperMC［5-6］的建模模塊，是一款基于現(xiàn)代軟件工程理念［7］設計開發(fā)的多物理耦合分析自動建模軟件系統(tǒng)，其目的是解決核分析中的多物理耦合建模挑戰(zhàn)。目前SuperMC／MCAM 已廣泛應用于國內(nèi)外核相關(guān)領(lǐng)域中［8-17］。

出于計算效率和正確性的考慮，在典型蒙特卡羅程序，如MCNP［18］中，所有的空間包括實體間空隙（空腔）均被描述為柵元。由于一般的工程CAD 模型僅建立代表實體的幾何，對于空腔幾何，需要由建模軟件自動生成，否則需要大量的人力手工創(chuàng)建代表空腔的CAD 模型［3］，并且在手工建立空腔模型的過程中仍有導入錯誤的風險［19］。由于CAD 模型的轉(zhuǎn)換表達式相對固定，自動建模軟件中自動空腔生成方法的優(yōu)劣決定了整體模型轉(zhuǎn)換效率以及最終蒙特卡羅輸運計算的效率。

目前在自動建模軟件中已實現(xiàn)了基本的自動空腔填充方法，如在SuperMC／MCAM 中已實現(xiàn)了基于CAD實體的自動空腔填充算法［3-4，19］，該方法能自動使用子空間表達式減去與其相交的實體表達式描述子空間內(nèi)的空腔，但其對空腔的描述復雜度受限于所相交實體的描述復雜度，在復雜模型中無法生成高效簡潔的空腔描述，限制了其在復雜模型上的應用。為進一步優(yōu)化空腔生成的效率，本文發(fā)展一種基于凸實體的自動空腔生成算法并基于SuperMC／MCAM 最新版本SuperMC／MCAM5.2進行實現(xiàn)，基于國際熱核聚變實驗堆ITER 中子學基準模型，對新的空腔生成算法進行測試。

1 基于凸實體的空腔生成算法

基于凸實體的空腔簡化算法是使用較實體幾何體更小單位的幾何體進行空腔相交計算，通過描述更少的相交幾何體部分使空腔描述更簡潔，以提高空腔生成的速度及輸運計算中粒子徑跡跟蹤的效率。該算法的基本流程如圖1所示，對其描述如下：

1）創(chuàng)建包裹所有實體的長方體a，形成初始的空腔長方體集合A＝｛a｝，假設所有的實體s組成實體集合S。

2）將所有實體切分為凸實體c，形成整體凸實體集合，若凸實體由原實體外殼面組成，則凸實體為正凸實體，記其集合為Cout，否則凸實體屬于內(nèi)凸實體集合Cin。凸實體可描述實體，如對 于s∈S，s＝（∪c∩s＝candc∈CoutDesc（c））-（∪c∩s＝sandc∈CinDesc（c）），其中，Desc（c）指實 體c的蒙特卡羅幾何描述。

圖1 基于凸實體的空腔生成算法流程Fig.1 Flow chart of void generation algorithm based on convex solids

3）執(zhí)行如下步驟，直到A＝?。

（1）從空腔長方體集合中取出空腔長方體a∈A，檢查與a 相交的兩類凸實體，形成兩類凸實體集合：與a相交的外殼集合Cint，out＝｛c｜c∈Coutand c∩a＝?｝以及與a 相交的內(nèi)殼集合Cint，in＝｛c｜c∈Cinand c∩a≠?｝，a 中的空腔部分設為空腔va，其表達式為：Desc（va）＝（Desc（a）- （∪c∈CoutDesc（c）））∪（Desc（a）∩（∪c∈Cint，inDesc（c）））。

（2）判斷a 內(nèi)空腔的描述Desc（va），若達到預定的要求（空腔表達式整體描述長度小于閾值），則生成該空腔的表達式，并將a 從A 中刪除并轉(zhuǎn)到步驟1，否則轉(zhuǎn)到步驟3。

（3）取a最長一邊作為切分邊，以通過該邊中點且垂直于該邊的平面切分空腔長方體形成空腔長方體anew1和anew2，加入A，轉(zhuǎn)到步驟1。

相比SuperMC／MCAM 中原空腔算法，基于凸實體的空腔生成算法（新算法）提升了模型轉(zhuǎn)化的效率，主要原因如下：1）新算法只需判斷凸實體與空腔長方體是否相交，而原算法需判斷空腔長方體與整個實體相交與否，相交的判斷和參與判斷的實體復雜度成正比，所以相比原算法，新算法減少了空腔生成過程的幾何計算量，提升了空腔生成效率；2）新算法采用空腔長方體減凸實體描述方式，該描述方式較原描述方式減少了實體中不與空腔長方體相交的部分的描述，故新算法的空腔描述復雜度降低，因此采用較少的分割即可達到與原算法同樣的空腔復雜度要求，這也提升了空腔生成效率。同時，新空腔生成算法相比原算法提升了蒙特卡羅計算效率，主要原因如下：1）新空腔算法相比原空腔算法，建立的空腔的極限復雜度更低，降低了粒子與空腔柵元單次計算的復雜度；2）新空腔算法相比原空腔算法，在相同復雜度下，減少了粒子所需計算的空腔數(shù)量，減少了計算次數(shù)。

2 測試驗證

本文在SuperMC／MCAM 中實現(xiàn)新的空腔算法。為證實新空腔算法的正確性以及效率提升，采用ITER Alite模型，通過MCNP建模與MCNP計算兩方面對比新算法與基本算法，進行正確性和效率的測試。圖2為ITER Alite模型，該模型實體部分包含4 194個實體柵元，32 683 個 面，是ITER 于2009 年 發(fā) 布 的 用 于ITER 中子學分析的基準模型。該模型包括了托卡馬克本體和部分生物屏蔽層，考慮到模型的計算效率和托卡馬克裝置的對稱性，該模型僅一個40°的扇形區(qū)域，通過反射面模擬全反應堆裝置計算。

圖2 SuperMC／MCAM 中的ITER Alite三維視圖Fig.2 3Dview of ITER Alite model in SuperMC／MCAM

2.1 正確性測試

原空腔生成算法的正確性已得到證明且已應 用于ITER 建 模 分 析 中［10，17，19］，并 已 成 功 應用于ITER Alite的自動建模中［20］。為驗證新算法的正確性，本文基于ITER Alite的CAD模型，分別使用基于凸實體的空腔生成算法和原空腔生成算法生成MCNP模型，在兩個模型上分別進行MCNP計算。在計算模型中，于坐標（600，0，0）處設置各向同性能量為14 MeV的點源，進行輸運計算并統(tǒng)計模型各柵元中的粒子徑跡長度。使用原空腔建模算法和新空腔建模算法的計算結(jié)果的相對偏差如圖3所示。

在兩種空腔建模算法得到的模型的計算結(jié)果中，柵元徑跡長度的計算結(jié)果的最大相對偏差為2.80×10-4，在允許范圍內(nèi)。實驗結(jié)果證實了新空腔算法建模的正確性。

圖3 兩種空腔算法建模后計算結(jié)果的相對偏差Fig.3 Relative deviation of calculation result of models generated by two void generation algorithms

2.2 效率測試

為對比新、舊空腔建模算法的建模效率，分別使用不同的復雜度參數(shù)生成一系列模型。其中空腔復雜度指描述空腔使用的面的數(shù)量，此參數(shù)限制了生成的空腔的最大復雜度。兩種空腔生成算法的空腔生成時間隨空腔復雜度變化的趨勢如圖4所示。

圖4 兩種算法空腔建模時間與空腔復雜度的關(guān)系Fig.4 Relationship between void modeling time and void complexity of two algorithms

空腔復雜度為500以下時，新空腔生成算法的建模時間均小于原空腔算法，在二者均能達到的空腔復雜度區(qū)域，新空腔生成算法的效率是原算法的1.1～25倍；原空腔建模算法在平均空腔復雜度為251時已達極限，而新空腔算法的建模極限是空腔復雜度為49，這證實了新空腔算法對空腔復雜度的可擴展性。

為驗證新、舊空腔算法所建立模型的計算效率，對比了同樣空腔復雜度情況下新、舊算法生成的空腔數(shù)量，從模型整體復雜度的角度來側(cè)面對比兩種算法生成模型的計算效率。圖5為新、舊算法的空腔復雜度與空腔數(shù)量的關(guān)系。

圖5 兩種算法空腔數(shù)量與空腔復雜度的關(guān)系Fig.5 Relationship between void number and void complexity of two algorithms

在同樣的空腔復雜度下，新算法的空腔數(shù)量均低于舊算法。在二者有交集的空腔復雜度范圍內(nèi)，新空腔生成算法的空腔數(shù)量是原空腔算法的1／3～1／60。

為正面驗證新空腔建模算法得到的模型在計算性能上的提升，對不同復雜度的模型的計算速度進行測試。測試均使用了真空材料，點源為設置在（600，0，0）處能量為14 MeV 的中子源，模擬1千萬粒子歷史的測試結(jié)果如圖6所示。

圖6 兩種算法空腔復雜度與計算時間的關(guān)系Fig.6 Relationship between void complexity andcalculation time of two algorithms

從圖6可看出，在同樣的空腔復雜度下，由于新空腔算法可生成更少的空腔數(shù)量，其計算效率優(yōu)于原空腔算法；而由于新空腔算法可在很好地控制空腔數(shù)量的前提下得到更低的空腔復雜度，所以其全局最優(yōu)計算效率要高于原空腔算法。試驗得到，新空腔算法生成的模型中計算速度最快的時間為26.9min，而原空腔算法生成的模型最快計算時間為64.91 min，新空腔算法提高了1.4倍計算效率。

3 總結(jié)與展望

本文提出了一種新的基于凸實體的自動空腔生成算法，將參與空腔描述的單位從柵元縮小為組成柵元的凸實體，這使得空腔自動建模過程可在控制空腔數(shù)量的前提下得到較優(yōu)的空腔復雜度。該方法實現(xiàn)后已集成進SuperMC／MCAM 中。通過測試對比發(fā)現(xiàn)，對應不同的空腔復雜度，新算法的建模效率提升最多達25倍，使生成的空腔復雜度簡化至1／60。對比兩種算法各自生成的最高效的計算模型，新算法較原算法提升了1.4 倍計算效率。

感謝FDS團隊其他成員提供的幫助和支持。

［1］ TSIGE-TAMIRAT H，F(xiàn)ISCHER U，SERIKOV A，et al.Automatic generation and validation of an ITER neutronics model from CAD data［J］.Fusion Engineering and Design，2007，82：1 956-1 959.

［2］ SATO S，IIDA H，NISHITANI T.Development of CAD／MCNP interface program prototype for fusion reactor nuclear analysis［J］.Fusion Engineering and Design，2006，81：2 767-2 772.

［3］ WU Y.CAD-based interface programs for fusion neutron transport simulation［J］.Fusion Engineering and Design，2009，84（7-11）：1 987-1 992.

［4］ 吳宜燦，李瑩，盧磊，等.蒙特卡羅粒子輸運計算自動建模程序系統(tǒng)的研究與發(fā)展［J］.核科學與工程，2006，26（1）：20-27.WU Yican，LI Ying，LU Lei，et al.Research and development of the automatic modeling system for Monte Carlo particle transport simulation［J］.Chinese Journal of Nuclear Science and Engineering，2006，26（1）：20-27（in Chinese）.

［5］ WU Y C，SONG J，ZHENG H Q，et al.CADbased Monte Carlo program for integrated simulation of nuclear system SuperMC［J］.Annals of Nuclear Energy，doi：10.1016／j.anucene.2014.08.058.

［6］ SONG Jing，SUN Guangyao，CHEN Zhenping，et al.Benchmarking of CAD-based SuperMC with ITER benchmark model［J］.Fusion Engineering and Design，2014，89（11）：2 499-2 503.

［7］ 吳宜燦，胡麗琴，龍鵬程，等.先進核能軟件發(fā)展與核信息學實踐［M］∥中國科研信息化藍皮書.北京：科學出版社，2013：232-244.

［8］ HU H，WU Y，CHEN M，et al.Benchmarking of SNAM with the ITER 3D model［J］.Fusion Engineering and Design，2007，82：2 867-2 871.

［9］ BOTH J P，MAZZOLO A，PENELIAU Y，et al.User manual for version 4.3of the TRIPOLI-4 Monte-Carlo method particle transport computer code［M］.France：CEA，2003.

［10］LI Y，LU L，DING A，et al.Benchmarking of MCAM 4.0with the ITER 3Dmodel［J］.Fusion Engineering and Design，2007，82（15-24）：2 861-2 866.

［11］吳宜燦，李靜驚，李瑩，等.大型集成多功能中子學計算與分析系統(tǒng)VisualBUS 的研究與發(fā)展［J］.核科學與工程，2007，27（4）：365-373.WU Yican，LI Jingjing，LI Ying，et al.An integrated multi-functional neutronics calculation and analysis code system：VisualBUS［J］.Chinese Journal of Nuclear Science and Engineering，2007，27（4）：365-373（in Chinese）.

［12］WU Y，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）：2 713-2 718.

［13］WU Y，QIAN J，YU J.The fusion-driven hybrid system and its material selection［J］.Journal of Nuclear Materials，2002，307-311：1 629-1 636.

［14］WU Y，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：1 893-1 903.

［15］CHEN Y，WU Y.Conceptual study on high performance blanket in a spherical Tokamak fusiondriven transmitter［J］.Fusion Engineering and Design，2000，49-50：507-512.

［16］吳宜燦，柏云清，宋勇，等.中國鉛基研究反應堆概念設計研究［J］.核科學與工程，2014，34（2）：56-63.WU Yican，BAI Yunqing，SONG Yong，et al.Conceptual design of China Lead-based Research Reactor CLEAR-Ⅰ［J］.Nuclear Science and Engineering，2014，34（2）：56-63（in Chinese）.

［17］王國忠，黨同強，熊健，等.MCAM4.8 在ITER建筑大廳中子學建模中的應用［J］.核科學與工程，2011，31（4）：351-355.WANG Guozhong，DANG Tongqiang，XIONG Jian，et al.Application of MCAM4.8in creating neutronics model for ITER building［J］.Chinese Journal of Nuclear Science and Engineering，2011，31（4）：351-355（in Chinese）.

［18］BRIESMEISTER J F.MCNP：A general Monte Carlo N-particle transport code，version 4C，LA-13709-M［M］.US：LANL，2000.

［19］李瑩.蒙特卡羅粒子輸運計算自動建模系統(tǒng)預處理關(guān)鍵技術(shù)研究［D］.合肥：中國科學院合肥物質(zhì)科學研究院，2009.

［20］LU L，LEE Y K，ZHANG J，et al.Development of Monte Carlo automatic modeling functions of MCAM for TRIPOLI-ITER application［J］.Nuclear Instruments and Methods in Physics Research Section A，2009，605（3）：384-387.