模块化粒子输运程序包PHEN的开发与应用

朱金辉，陶应龙，卓 俊，谢红刚，左应红，商 鹏，韦 源，牛胜利

（西北核技术研究所，西安710024；强脉冲辐射环境模拟与效应国家重点实验室，西安710024）

经过半个多世纪的发展，特别是近年来随着高性能计算技术的飞速提升，蒙特卡罗（MC）方法已成为解决粒子输运问题的首选工具［1］。目前，比较成熟的粒子输运大型MC计算程序主要有MCNP、EGS、PENELOPE 和 GEANT4等［2－4］。这些程序功能强大、通用性好，被广泛应用于核技术相关领域的模拟计算中，为核工程技术及辐射环境的研究与发展发挥了重要的支撑作用［5－6］。但是针对一些特定模拟计算问题，如粒子次级反应过程中的特定参数获取、物理过程之间的耦合计算等，这些程序皆无法满足需求，且由于这些程序结构复杂和庞大，在其基础上进行二次开发实现特定功能比较困难。

为此，本文采用模块化的程序设计方法，在前期开发的中子－光子输运程序TOPAN［7］的基础上，开发了能够模拟光子、质子、电子和中子与物质相互作用的主要过程，粒子在几何中输运及输运结果统计记录的模块化程序，形成了模块化程序包PHEN（the transportation program of photon，hadron，electron，and neutron）。该程序包中的各个程序模块能够较方便地被调用，可根据用户需要开发特定功能程序，也能够较方便地实现与其他物理过程的耦合计算。

1 PHEN程序包基本功能介绍

PHEN程序包采用FORTRAN语言编制，基于MC方法模拟各个要素［8］，能够模拟光子、质子、电子和中子与物质相互作用及其输运问题，包含粒子输运所需的各基本功能模块。与常见的MCNP及GEANT4程序相比，PHEN程序包具有3个特点：

1）清晰的树结构几何建模功能

在几何建模方面，PHEN程序包采用树结构描述几何模型［4］。树结构几何模型的基本思想是利用子体和母体的关系建立具有代（母子体）关系：几何结构的树状关系图，在几何相关计算中只需处理本代和具有直接子母体关系的部分，从而减小复杂几何结构的几何算法计算量，理论上可从O（N）量级降低到O（logkN）量级，其中，N为几何体总数，k为树结构总代数。

树结构几何的基本思想如图1所示。根据图1左侧的结构可以建立图1右侧的树状结构，其中，⑨是③和⑧的母体，⑧是⑦和④的母体……，当粒子位置在体⑧中时，求粒子前向与结构界面交点的程序，只需处理⑧的子体（⑦和④）及其与母体（⑨）的交界面，不需要遍历模型内所有的结构交界面，理论上用于计算粒子输运前方界面交点的时间将减少约2／3。

图1 树结构几何示例Fig.1Sketch of tree structure

PHEN程序包的粒子输运几何建模模块能够建立各种标准体模型（包含各种二次曲面围成的体）及其各种组合模型。采用PHEN程序的树结构几何建模方法，建模过程中从最小的子体开始，逐级建模，几何结构清晰，可以较容易地建立较为复杂的几何模型。图2为利用PHEN几何模块建立的计算模型示意图，图中不同颜色代表不同的材料。

图2 利用PHEN程序几何模块建立的计算模型Fig.2Acalculation model established by PHEN geographic module

PHEN程序包的几何模块，还可以用于射线追踪模拟计算，比如可见光的输运模拟等［9］。

2）带电粒子输运的功能拓展

在带电粒子输运的MC模拟中，低能带电粒子的输运一直是模拟的难点之一。在PHEN程序包的研发过程中，研究了低能电子直接MC模拟方法［10］，并将电子输运的直接模拟方法和压缩历史方法相结合，建立了适合于能量在50eV～1GeV范围的电子输运MC方法，即混合模拟方法。把低能电子与物质的相互作用分为能量损失较少的“软碰撞”和能量损失较大的“硬碰撞”。其中，“软碰撞”发生次数特别多，采用压缩历史方法模拟；“硬碰撞”发生次数较少，采用直接模拟方法模拟。

在低能质子与物质的相互作用模拟中，非弹性散射物理过程的模拟采用预复合模型［5］，能够模拟能量为10keV的质子在常用材料中的输运过程。另外，针对中子与含氢材料产生反冲质子的物理现象，在中子弹性散射模块中添加了产生次级质子的模型，能够模拟由弹性碰撞产生的反冲质子，结合质子输运模块，能够实现中子输运－反冲质子－质子输运过程的一体化模拟。

3）可以定制的数据统计功能［3］

根据MC模拟的特点，并借鉴相关研究工作［11］的思路，采用了粒子标识方法记录粒子及其次级粒子的输运历史，并根据标识变量的值进行统计，以达到求解特定物理量的目的。除了记录射线种类、位置、方向、能量、权重等射线基本参数外，专门设置了粒子标识数组用来标识射线产生的位置、作用核素、反应类型等感兴趣参数。利用粒子标识数组，可以对具体特定属性的粒子进行统计分析，得到某种材料或某种反应产生的粒子数对总计数的贡献，如特定材料（n，2n）反应对中子通量的贡献，光子康普顿散射对光子点通量的贡献等。通过粒子标识和用户定义响应函数，还可以扩充数据统计模块功能。

2 中子与典型核素产生的次级γ射线特征参数计算

获得不同能量中子与特定核素相互作用过程中产生的次级γ射线特征参数，对于分析中子－γ联合屏蔽问题具有重要意义。在模拟过程中，希望直接记录中子与原子核相互作用产生的次级γ信息，而通用的MC模拟程序，如MCNP，在粒子输运模拟的统计记录过程中难以排除介质在后续输运过程中对次级γ射线参数的影响。为此，在PHEN程序包的中子及γ输运过程模块的基础上进行二次开发，能够计算不同能量的中子与典型核素相互作用后产生的次级γ的产额及平均能量。

计算中不模拟粒子的几何输运过程，直接模拟特定能量中子与特定原子核的相互作用过程，统计产生的次级γ射线的产额及能量。设入射中子总数为Nn，所产生的次级γ的总数为Nγ，次级γ总能量为Eγ，tot。设中子能量为En，统计得到能量为En的中子产生的次级γ的平均能量Eγ，ave为

定义次级γ产额Y为

其中，Y表示中子与原子核相互作用过程产生次级γ的概率。

针对典型核素H、C、N、Fe及U，计算得到了不同能量中子作用下产生的次级γ产额及平均能量，如图3所示。

图3 不同能量中子在典型核素中产生的次级γ产额及平均能量Fig.3Yield and average energy of secondaryγ－rays generated by neutron interacting with typical nuclei

从图3可以看到，不同核素与中子相互作用的次级γ产额和平均能量存在很大差异，通常高Z核素的次级γ产额很高，但平均能量相对较小。

3 基于指向概率法的点通量角度谱计算

在常见粒子输运模拟程序中，MCNP程序具备中子及光子点通量计数功能，采用指向概率方法记录点通量，通过计算粒子与介质原子核反应时指向记录点的虚次级粒子对记录点的通量贡献，得到总的点通量计数。在具体应用中，用户有时还希望得到点通量的角度分布信息，但MCNP程序目前不具备该功能。

利用模块化程序包PHEN，基于MCNP程序中的粒子点通量计数方法，开发得到了包含角度分布信息的点通量计数功能。

为了便于记录分析，需预先设定参考方向。进行点通量计算时，在记录虚次级粒子的点通量贡献的同时，计算虚次级粒子运动方向与参考方向的夹角，并对其进行统计记录。

为了验证开发的点通量角度谱分析功能的正确性，利用MCNP程序球面上面通量角度谱的对称性进行对比计算，计算模型如图4所示。模型中材料为铁。首先利用MCNP程序记录3个球面S1，S2，S3上的面通量及其角度谱；其次利用PHEN程序包记录3个点P1，P2，P3的点通量及其角度谱，参考方向为x轴正方向；最后利用PHEN程序包记录S1，S2，S3上的面通量及其角度谱。球心能量为1MeV中子产生的点通量和面通量计算结果，如图5所示。

图4 计算模型Fig.4Calculation model

图5 中子通量角度谱结果Fig.5Angle spectra of neutron fluence

根据对称性分析，上述两种方法计算的中子通量角度谱应该一致。从图5可以看出，PHEN程序包计算给出的面通量和点通量及其角度分布与MCNP程序给出的结果基本一致，验证了本文开发的点通量角度谱分析功能的正确性。需要说明的是，本文开发的点通量角度谱分析功能也适用于非对称情况下的角度谱分析，用户只需根据研究需求，预先指定参考方向，即可得到所需的角度分布信息。

4 中子输运时间特性模拟计算

在超临界系统分析中，除了系统裂变中子的有效增殖因子keff，系统的时间增殖特性也是研究人员关心的重要方面。MCNP5程序具有成熟的裂变驿站方法，可以方便地给出系统的keff，但难以直接给出时间增殖特性。

本文利用模块化程序包PHEN，开发了粒子输运的时间驿站功能，能够计算给出系统的时间增殖特性。同时参照MCNP5程序的裂变驿站方法，也能够计算超临界系统裂变中子的有效增殖因子［12］。

针对图6所示的超临界系统模型，利用PHEN程序包和MCNP5程序计算得到的有效增殖因子，如表1所列。可见，2种程序的计算结果基本一致。

图6 超临界系统计算模型Fig.6Calculation model of a supercritical system

采用时间驿站功能，对图6所示模型计算了不同参数条件下，系统内中子总数随时间的变化关系，如图7所示。

图7 中子总数随时间的变化Fig.7Sum number of neutron vs.time

从图7可以看出，PHEN程序包能够有效处理与时间相关的粒子输运问题，可用于超临界系统中子输运模拟，能够给出系统裂变中子的有效增殖因子及时间增殖特性。

5 总结与展望

针对特定问题需要，开发了能够模拟粒子输运过程的模块化程序包PHEN，具备基本的中子、光子、电子和质子耦合输运蒙特卡罗模拟能力。PHEN程序包中的模块化程序能够较方便地被调用，既可以根据用户需要开发特定功能程序，也能够较方便地实现与其他物理过程的耦合计算。PHEN程序包的基本功能模块可用于粒子输运模拟技术研究和应用研究。

在后续工作中，将不断改进完善PHEN程序包，如增加光电子耦合输运中的厚靶韧致辐射模型、增加切伦科夫辐射模型、改进截面数据处理算法等；同时，要不断提高程序包的接口适应性，为粒子输运相关数值模拟工作提供更加完备好用的粒子输运模块化程序包。