基于CAD技术的特征线中子输运计算程序开发

陈珍平，王电喜，何 桃，王国忠，郑华庆，FDS团队

（1.中国科学技术大学，安徽 合肥 230027；2.中国科学院核能安全技术研究所，安徽 合肥230031）

中子输运问题一直是核反应堆物理分析与辐射屏蔽设计研究的重点问题之一。随着反应堆技术的快速发展，在未来进行先进反应堆设计与分析时对中子学模拟计算程序的要求也越来越高，主要体现在计算精度、计算效率以及处理几何复杂度等方面，因此开发具有自主知识产权的、快速精确的中子学计算程序将具有非常重要的意义。特征线方法[1］（Method of Characteristics，MOC）是近年来在中子输运计算与数值模拟研究中热点研究的方法之一。该方法是从中子输运方程的全微分形式出发，沿着中子飞行轨迹（称为特征线）进行积分求解的一种确定论方法。由于该方法同时具有碰撞概率法和离散纵标法的优点，且理论上不受计算模型几何形状的限制，近十多年来已逐渐成为国际上求解中子输运方程的主要确定论方法之一。

长期以来，阻碍特征线方法能否广泛应用的主要问题在于能否将特征线理论与有效的几何处理方法结合起来。本文在FDS团队自主研发的大型集成多功能中子学计算与分析系统VisualBUS[2－3］框架下，利用核与辐射输运计算自动建模软件 MCAM[4－6］的几何处理能力进行特征线方法几何预处理，可高效、快速地完成几何建模、网格剖分和特征线生成等几何处理功能。最后，利用基准例题对特征线数值计算方法和程序进行了充分地检验，其结果验证了本研究几何处理方法和计算程序的正确性与可靠性。

1 特征线理论

特征线方法是将中子输运方程沿着中子飞行轨迹和特定的离散方向进行输运求解的一种方法。在本文中，对方位角离散采用Filippone[7］提出的离散方法；对极角离散耦合了LO离散和TY离散[8］方法；对几何模型区域进行非结构网格剖分，划分为许多平源近似区，在每个平源区里的中子源项被认为是常数。

图1 特征线示意图Fig.1 The characteristic lines of MOC

根据前面的基本假设，单群特征线中子输运方程可表示为：

其中：

s是特征线长度，Σt，i是区域i的宏观总截面，Ψi，k（s，Ωm）是区域i沿着特征线方向长度为s处的中子角通量密度，Qf（Ωm）是区域i的中子源项。

方程（1）为一阶线性微分方程，可以解析求解得到：

式（2）中Ψini，k（Ωm）是特征线k入射到区域i里入射点处的入射角通量密度。

根据方程（2）可以得到，特征线k穿出区域i的出射角通量密度为：

式（3）中si，k是特征线k穿过区域i的总长度；将特征线k上所有点的角通量密度进行积分平均，则可以获得特征线k上任一点的平均角通量密度：

同理，将穿过区域i的所有特征线上的平均通量密度在i区进行体积加权平均，可获得区域i的平均角通量密度：

式中δAk为特征线间距。

将区域i各个离散方向的角通量密度进行加权求和，得到区域i的平均中子标通量密度：

式中ωm是方向Ωm对应的离散权重，M为空间离散方向总数目。

2 几何预处理方法

2.1 几何建模

FDS团队自主研发的核与辐射输运计算自动建模软件MCAM实现了多种商用工程CAD软件（CATIA、AutoCAD、UG等）与多种辐射输运计算程序（MCNP、TRIPOLI[9］等）之间的接口。一方面，可以直接由计算机辅助设计（CAD）模型生成完整的辐射输运计算程序的输入模型，包括空腔模型、材料、源和计数信息等；另一方面，可以解析辐射输运计算程序的输入模型，生成CAD模型并可视化，供分析和修正。MCAM具有非常强大的几何建模能力，能够方便地建立各种复杂几何的CAD模型。目前MCAM已经在一些复杂核装置（如FDS系列聚变反应堆[10］、ITER 模型[11－13］等）的中子学建模和计算分析中得到了广泛的应用。

本文基于MCAM提供的造型功能和内嵌的布尔运算功能，在MCAM平台上进行二次开发了MOC几何建模模块，通过调用该模块的接口，可以快速、精确地进行几何建模，且能方便地构造出如反应堆燃料棒、燃料组件和反应堆堆芯等各种常见的几何模型。同时，利用MOC建模模块所构造的几何模型可以通过调用MCAM底层的渲染接口在MCAM的GUI中渲染出来（如图2），可非常清晰、直观地表达出模型的几何结构，便于用户进行模型正确性检查。同时，还为几何模型添加了材料属性编辑功能，即用户可以对图2中的几何模型进行材料属性添加，则MOC程序在进行输运计算时，可以直接根据几何模型的材料属性来获取材料的各种截面信息。最终构造好的几何模型（已添加材料属性）可保存为标准的CAD文件格式（sat，step，igs等），MOC 程序做中子输运计算时可以直接将上述文件导入进行输运计算，即本文开发的MOC程序支持CAD文件直接导入进行输运计算的功能。

图2 MCAM几何建模Fig.2 Geometry modeling of MCAM

2.2 特征线跟踪

由前面的特征线公式（1）～（5）可知，在进行特征线中子输运计算时，需要不断地进行特征线跟踪来获取模型的几何和材料等相关信息。所以，MOC程序在进行中子输运计算时所需要的几何和材料等信息主要以特征线为载体，主要包括：特征线间距、特征线段长度、特征线所在网格编号、对应的材料号以及边界信息等。本文主要基于MCAM提供的射线跟踪功能，二次开发了MOC程序特征线自动处理模块，图3表示了整个MOC程序中几何建模模块、射线跟踪模块以及物理计算模块等各模块间的相互关系。

图3 MOC程序模块关系图Fig.3 The relationship of different modules in MOC

在开始特征线跟踪之前，程序根据用户在外部配置文件中提供的方位角数目，自动进行方位角离散并计算每个方位角对应的特征线间距；然后根据网格划分参数对模型进行网格划分；最后根据已离散的方位角数值，通过调用MCAM的射线跟踪接口来产生所有方位角下对应的特征线信息，并将之保存用于MOC程序的输运计算。由于MCAM提供了射线跟踪功能，并且计算精度和效率能保持在合理范围内，在二次开发时可直接继承这些功能进行应用程序的开发，大大提高了应用程序的开发效率。

3 数值验证

本文提出了基于MCAM进行MOC几何预处理的方法，并开发了特征线中子输运计算MOC程序，为了验证程序的正确性与可靠性，本文给出了一维7群例题和二维C5G7－UO2组件基准例题的数值验证结果。

3.1 一维7群例题

本文采用文献[14］中定义的一维7群问题，其几何尺寸和材料布置如图4所示，各种材料的7群截面来自C5G7－MOX基准题[15］。由于该问题材料间能谱差异大、泄漏较强，因此可以很好的验证MOC程序的正确性。

图4 一维7群例题几何布置（cm）Fig.4 The layout of 1D7Gproblem（cm）

该问题的参考解由文献[14］中的程序PEACH给出。表2给出了与PEACH程序计算结果的比较，其中本文采用平源菱形差分特征线法给出，平源区长度为0.1cm，方位角数目为4个，特征线间距约为0.05cm，从表1中可以看出本文计算的结果与参考值相比，keff误差约为－3pcm，结果符合良好。

表1 一维7群例题MOC与PEACH结果比较Table 1 Comparision of kefffor 1D7G problem between MOC and PEACH

表2 二维C5G7－UO2组件例题kinf计算结果比较Table 2 Comparision of kinffor C5G7－UO2assembly results

在特征线方法中，影响计算结果精度的主要因素包括划分网格尺寸大小、方位角数目、极角数目、特征线间距等。针对本例题，本文就前两个因素做了敏感性分析。图5给出了在不同网格划分、相同角度离散（TY－2极角和4个方位角）、相同特征线间距（0.05cm）情况下keff误差敏感性分析。从图5可看出，网格大小对keff精度会有较大影响；当网格的尺寸小于1.0cm时，MOC程序能够获得非常满意的计算精度（小于70pcm）。图6给出了在相同网格尺寸（0.1cm）、相同极角离散（TY－2极角）、相同特征线间距（0.05cm）和不同方位角数目下keff误差敏感性分析。从图6中可以看出，当方位角数目N≥4时，继续增加方位角数目对最终结果影响不大，误差绝对值小于3pcm，故采用4个方位角已经可以获得比较满意的计算精度。

图5 不同网格大小keff误差分析Fig.5 kefferror analysis with different mesh dimensions

图6 不同方位角数目keff误差分析Fig.6 kefferror analysis with different angle dimensions

3.2 二维C5G7－UO2组件基准例题

二维C5G7－UO2基准题是国际经济合作与发展组织（OECD）发布的二维C5G7基准例题[15］的二氧化铀燃料组件，组件燃料布置与几何结构如图7所示。上图为17×17个栅元组成的UO2组件，包括UO2燃料栅元、导向管栅元和裂变室栅元。栅元为两区非均匀结构，各栅元的间距为1.26cm，栅元中的燃料芯块和导向管均匀化后的半径为0.54cm，组件四周采用全反射边界条件。同时，采用图7右下角的四分之一UO2燃料组件模型[16］进行进一步计算，通过和单组件计算结果对比，以进一步验证程序的正确性。

图7 2DC5G7－UO2燃料组件几何布置（cm）Fig.7 The layout of the 2DC5G7benchmark with UO2assembly （cm）

本文MOC程序计算时各材料区的截面信息采用文献[15］中二维C5G7模型提供的宏观截面数据。本文对1／4UO2组件和1／1UO2组件分别进行了kinf计算。在计算时，将模型进行非结构网格剖分；对空间角离散采用LO－2优化两极角和6方位角离散；特征线间距为0.02cm，计算结果如表2所示。其计算结果与 程 序 GALAXY[16］、GMVP[17］的 误 差 约 为－0.15%，说明与其他两程序计算结果吻合良好；同时，MOC程序计算的1／4UO2组件和1／1UO2组 件 的kinf分 别 为 1.331 345 和1.331 782，相对误差约为0.03%，两模型计算结果符合良好，证明了本文MOC程序的正确性与可靠性。

总结

几何预处理方法，即基于MCAM的几何处理引擎，通过C＋＋语言面向对象二次开发了MOC几何预处理模块，可以非常方便地进行特征线几何建模以及特征线信息快速获取，这样便将特征线理论与MCAM有效地结合起来。本文基于MCAM二次开发了特征线中子输运计算MOC程序，并利用相关基准例题对程序进行了验证，计算结果与国内外其他相关计算程序的结果吻合良好，表明本文方法和程序的可行性、正确性与可靠性。

致谢

本文在开展研究过程中，得到了西安交通大学刘宙宇博士在特征线理论方面的帮助与指导，在此深表感谢。

本文提出了一种基于CAD技术的特征线

[1]Askew R.A Characteristics for Mulation of the Neutron Transport Equation in Complicated Geometries[R］.Report AEEW－R－1108，UK Atomic Energy Authority，1972.

[2]吴宜灿，李静惊，李莹，等.大型集成多功能中子学计算与分析系统VisualBUS的研究与发展[J］.核科学与工程，2007，27（5）：72－85.

[3]吴宜灿，胡丽琴，龙鹏程，等.先进核能系统设计分析软件与数据库研发进展[J］.核科学与工程，2010，30（1）：55－64.

[4]Wu Y C，FDS Team.CAD－based Interface Prog－rams for Fusion Neutron Transport Simulation [J］.Fusion Engineering and Design，2009，84：1987－1992.

[5]吴宜灿，李莹，卢磊，等.蒙特卡罗粒子输运计算自动建模程序系统的研究与发展[J］.核科学与工程，2006，26（01）：20－27.

[6]李莹，曾勤，卢磊.利用ITER基准模型对 MCAM4.2进行检验[J］.核科学与工程，2008，28（1）：47－50.

[7]Filippone W L， Woolf S，Lavigne R J.Particle Transport Calculation with the Method of Streaming Rays[J］.Nuclear Science and Engineering，1981，77：119.

[8]Yamamoto A，Tabuchi M，Sugimura N，et al.Derivation of Optimum Polar Angle Quadrature Set for the Method of Characteristics Based on Approximation Error for the Bickley Function [J］. Journal of Nuclear and Technology，2007，44（2）：129－136.

[9]张俊军，曾勤，王国忠，等.蒙特卡罗程序TRIPOLI自动建模方法研究[J］.核科学与工程，2010，20（3）：84－88.

[10]Wu Y，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]曾勤，卢磊，李莹，等.蒙特卡罗粒子输运计算自动建模程序MCAM在ITER核分析建模中的应用[J］.原子核物理评论，2006，23（2）：138－141.

[12]王国忠，党同强，熊健，等.MCAM4.8在ITER建筑大厅中子学建模中的应用[J］.核科学与工程，2011，31（4）：351－355.

[13]熊健，王国忠，王电喜，等.MCAM 在ITER装置TRIPOLI三维中子学建模中的应用[J］.核科学与工程，2011，31（2）：162－167.

[14]汤春桃.中子输运方程特征线解法及嵌入式组件均匀化方法的研究 [D］.上海：上海交通大学，2009.

[15]Lewis E E，Simith M A，Tsoulfanidis N，et al.Benchmark Specification for Deterministic 2－D／3－D MOX Fuel Assembly Transport Calculations without Spatial Homogenization （C5G7MOX）[R］.NEA／NSC／DOC，2001.

[16]Yamaji K，Matsumoto H，Sato D，et al.Simple and Efficient Parallelization Method for MOC Calculation[J］.Journal of Nuclear Science and Technology，2010，47（1）：90－102.

[17]Mori T，Nakagawa T.MVP／GMVP：General Purpose Monte Carlo codes for Neutron and Photon Transport Calculations Based on Continuous Energy and Multigroup Methods [R］.Japan Atomic Energy Research Institute，1994.