蒙特卡罗自动建模中的空腔简化算法研究

俞盛朋，程梦云，龙鹏程，胡丽琴

（中国科学院 核能安全技术研究所，中国科学院 中子输运理论与辐射安全重点实验室，安徽 合肥 230031）

蒙特卡罗粒子输运计算方法通过在几何空间中追踪粒子与材料之间发生的物理反应来近乎真实地模拟粒子的运动规律。为能精确高效进行蒙特卡罗粒子输运计算，需准确简洁描述蒙特卡罗计算的几何模型。传统的手工蒙特卡罗建模方法费时费力，并可能引入几何错误，限制了蒙特卡罗程序的使用。为解决蒙特卡罗的建模困难，国内外已研发了一系列基于CAD的蒙特卡罗自动建模工具，如德国KIT研发的McCad［1］、日本JAEA研发的GEOMIT［2］和FDS团队研发的SuperMC／MCAM［3-4］。其中SuperMC／MCAM（multi-physics coupling analysis modeling program）是FDS团队研发的超级蒙特卡罗核计算仿真软件系统SuperMC［5-6］的建模模块，是一款基于现代软件工程理念［7］设计开发的多物理耦合分析自动建模软件系统，其目的是解决核分析中的多物理耦合建模挑战。目前SuperMC／MCAM 已广泛应用于国内外核相关领域中［8-17］。

出于计算效率和正确性的考虑，在典型蒙特卡罗程序，如MCNP［18］中，所有的空间包括实体间空隙（空腔）均被描述为栅元。由于一般的工程CAD 模型仅建立代表实体的几何，对于空腔几何，需要由建模软件自动生成，否则需要大量的人力手工创建代表空腔的CAD 模型［3］，并且在手工建立空腔模型的过程中仍有导入错误的风险［19］。由于CAD 模型的转换表达式相对固定，自动建模软件中自动空腔生成方法的优劣决定了整体模型转换效率以及最终蒙特卡罗输运计算的效率。

目前在自动建模软件中已实现了基本的自动空腔填充方法，如在SuperMC／MCAM 中已实现了基于CAD实体的自动空腔填充算法［3-4，19］，该方法能自动使用子空间表达式减去与其相交的实体表达式描述子空间内的空腔，但其对空腔的描述复杂度受限于所相交实体的描述复杂度，在复杂模型中无法生成高效简洁的空腔描述，限制了其在复杂模型上的应用。为进一步优化空腔生成的效率，本文发展一种基于凸实体的自动空腔生成算法并基于SuperMC／MCAM 最新版本SuperMC／MCAM5.2进行实现，基于国际热核聚变实验堆ITER 中子学基准模型，对新的空腔生成算法进行测试。

1 基于凸实体的空腔生成算法

基于凸实体的空腔简化算法是使用较实体几何体更小单位的几何体进行空腔相交计算，通过描述更少的相交几何体部分使空腔描述更简洁，以提高空腔生成的速度及输运计算中粒子径迹跟踪的效率。该算法的基本流程如图1所示，对其描述如下：

1）创建包裹所有实体的长方体a，形成初始的空腔长方体集合A＝｛a｝，假设所有的实体s组成实体集合S。

2）将所有实体切分为凸实体c，形成整体凸实体集合，若凸实体由原实体外壳面组成，则凸实体为正凸实体，记其集合为Cout，否则凸实体属于内凸实体集合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）执行如下步骤，直到A＝∅。

（1）从空腔长方体集合中取出空腔长方体a∈A，检查与a 相交的两类凸实体，形成两类凸实体集合：与a相交的外壳集合Cint，out＝｛c｜c∈Coutand c∩a＝∅｝以及与a 相交的内壳集合Cint，in＝｛c｜c∈Cinand c∩a≠∅｝，a 中的空腔部分设为空腔va，其表达式为：Desc（va）＝（Desc（a）- （∪c∈CoutDesc（c）））∪（Desc（a）∩（∪c∈Cint，inDesc（c）））。

（2）判断a 内空腔的描述Desc（va），若达到预定的要求（空腔表达式整体描述长度小于阈值），则生成该空腔的表达式，并将a 从A 中删除并转到步骤1，否则转到步骤3。

（3）取a最长一边作为切分边，以通过该边中点且垂直于该边的平面切分空腔长方体形成空腔长方体anew1和anew2，加入A，转到步骤1。

相比SuperMC／MCAM 中原空腔算法，基于凸实体的空腔生成算法（新算法）提升了模型转化的效率，主要原因如下：1）新算法只需判断凸实体与空腔长方体是否相交，而原算法需判断空腔长方体与整个实体相交与否，相交的判断和参与判断的实体复杂度成正比，所以相比原算法，新算法减少了空腔生成过程的几何计算量，提升了空腔生成效率；2）新算法采用空腔长方体减凸实体描述方式，该描述方式较原描述方式减少了实体中不与空腔长方体相交的部分的描述，故新算法的空腔描述复杂度降低，因此采用较少的分割即可达到与原算法同样的空腔复杂度要求，这也提升了空腔生成效率。同时，新空腔生成算法相比原算法提升了蒙特卡罗计算效率，主要原因如下：1）新空腔算法相比原空腔算法，建立的空腔的极限复杂度更低，降低了粒子与空腔栅元单次计算的复杂度；2）新空腔算法相比原空腔算法，在相同复杂度下，减少了粒子所需计算的空腔数量，减少了计算次数。

2 测试验证

本文在SuperMC／MCAM 中实现新的空腔算法。为证实新空腔算法的正确性以及效率提升，采用ITER Alite模型，通过MCNP建模与MCNP计算两方面对比新算法与基本算法，进行正确性和效率的测试。图2为ITER Alite模型，该模型实体部分包含4 194个实体栅元，32 683 个 面，是ITER 于2009 年 发 布 的 用 于ITER 中子学分析的基准模型。该模型包括了托卡马克本体和部分生物屏蔽层，考虑到模型的计算效率和托卡马克装置的对称性，该模型仅一个40°的扇形区域，通过反射面模拟全反应堆装置计算。

图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的点源，进行输运计算并统计模型各栅元中的粒子径迹长度。使用原空腔建模算法和新空腔建模算法的计算结果的相对偏差如图3所示。

在两种空腔建模算法得到的模型的计算结果中，栅元径迹长度的计算结果的最大相对偏差为2.80×10-4，在允许范围内。实验结果证实了新空腔算法建模的正确性。

图3 两种空腔算法建模后计算结果的相对偏差Fig.3 Relative deviation of calculation result of models generated by two void generation algorithms

2.2 效率测试

为对比新、旧空腔建模算法的建模效率，分别使用不同的复杂度参数生成一系列模型。其中空腔复杂度指描述空腔使用的面的数量，此参数限制了生成的空腔的最大复杂度。两种空腔生成算法的空腔生成时间随空腔复杂度变化的趋势如图4所示。

图4 两种算法空腔建模时间与空腔复杂度的关系Fig.4 Relationship between void modeling time and void complexity of two algorithms

空腔复杂度为500以下时，新空腔生成算法的建模时间均小于原空腔算法，在二者均能达到的空腔复杂度区域，新空腔生成算法的效率是原算法的1.1～25倍；原空腔建模算法在平均空腔复杂度为251时已达极限，而新空腔算法的建模极限是空腔复杂度为49，这证实了新空腔算法对空腔复杂度的可扩展性。

为验证新、旧空腔算法所建立模型的计算效率，对比了同样空腔复杂度情况下新、旧算法生成的空腔数量，从模型整体复杂度的角度来侧面对比两种算法生成模型的计算效率。图5为新、旧算法的空腔复杂度与空腔数量的关系。

图5 两种算法空腔数量与空腔复杂度的关系Fig.5 Relationship between void number and void complexity of two algorithms

在同样的空腔复杂度下，新算法的空腔数量均低于旧算法。在二者有交集的空腔复杂度范围内，新空腔生成算法的空腔数量是原空腔算法的1／3～1／60。

为正面验证新空腔建模算法得到的模型在计算性能上的提升，对不同复杂度的模型的计算速度进行测试。测试均使用了真空材料，点源为设置在（600，0，0）处能量为14 MeV 的中子源，模拟1千万粒子历史的测试结果如图6所示。

图6 两种算法空腔复杂度与计算时间的关系Fig.6 Relationship between void complexity andcalculation time of two algorithms

从图6可看出，在同样的空腔复杂度下，由于新空腔算法可生成更少的空腔数量，其计算效率优于原空腔算法；而由于新空腔算法可在很好地控制空腔数量的前提下得到更低的空腔复杂度，所以其全局最优计算效率要高于原空腔算法。试验得到，新空腔算法生成的模型中计算速度最快的时间为26.9min，而原空腔算法生成的模型最快计算时间为64.91 min，新空腔算法提高了1.4倍计算效率。

3 总结与展望

本文提出了一种新的基于凸实体的自动空腔生成算法，将参与空腔描述的单位从栅元缩小为组成栅元的凸实体，这使得空腔自动建模过程可在控制空腔数量的前提下得到较优的空腔复杂度。该方法实现后已集成进SuperMC／MCAM 中。通过测试对比发现，对应不同的空腔复杂度，新算法的建模效率提升最多达25倍，使生成的空腔复杂度简化至1／60。对比两种算法各自生成的最高效的计算模型，新算法较原算法提升了1.4 倍计算效率。

感谢FDS团队其他成员提供的帮助和支持。

［1］ TSIGE-TAMIRAT H，FISCHER 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］ 吴宜灿，李莹，卢磊，等.蒙特卡罗粒子输运计算自动建模程序系统的研究与发展［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］ 吴宜灿，胡丽琴，龙鹏程，等.先进核能软件发展与核信息学实践［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］吴宜灿，李静惊，李莹，等.大型集成多功能中子学计算与分析系统VisualBUS 的研究与发展［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，FDS 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，FDS 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］李莹.蒙特卡罗粒子输运计算自动建模系统预处理关键技术研究［D］.合肥：中国科学院合肥物质科学研究院，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.