0.   引　言

乏燃料贮存格架用于贮存达到规定燃耗限值的乏燃料组件^[1]，具有组成结构复杂、贮存密度大、工况较多等特点，这也增加了其建模的复杂性。近年来，储存单元采用了中子毒物材料，格架的设计转为密集型，相邻储存单元间的距离大大减小，乏燃料池单位面积的储存量大大增加^[2]，由于其贮存特点，贮存过程中的安全风险依旧存在。而乏燃料贮存水池是贮存乏燃料的主要辅助设施^[3]，因此对乏燃料贮存水池展开临界分析计算对保证其安全性具有重要意义。

MONK是一款可用于反应堆物理及临界安全分析的大型三维粒子输运蒙特卡罗程序。该程序主要通过计算机模拟各种类型中子的产生、输运、死亡来计算中子有效增殖因子（k_eff），同时也可计算中子流、反应率以及边界截面流等参数^[4-5]。

目前，MONK程序广泛应用于乏燃料贮存格架的临界分析计算^[6-9]。RMC为清华大学自主开发的蒙特卡罗程序，在应用到乏燃料贮存格架的临界计算中时，需与MONK程序进行对比验证。其中，乏燃料贮存格架的模型设计算例较多且几何复杂，对于大量乏燃料贮存格架的MONK输入卡，若全部重新手工使用RMC建模，会导致其在建模过程中耗时长，且易于引入错误。因此，为了提高乏燃料贮存格架临界计算的建模效率和建模准确性，便于MONK与RMC程序之间的对比验证，有必要研究更精准且更有效率的转换程序，从而实现自主化蒙特卡罗程序在燃料贮存格架临界分析的应用。

1.   模型转换流程

MONK程序的几何由简单的几何物体构造成不同的子空间（PART或HOLE），再由不同的子空间组装成更大的空间，这个过程会不断重复，直到整个系统被组装完成。MONK程序建模的特点在于：材料部分可定义为原子密度或质量密度，核素直接按元素符号定义，可选择后跟质量数以识别具体同位素；同时，如果定义为原子密度且材料总密度等于各核素密度之和时，材料总密度可简单写为0。几何部分在填充的几何与被填充几何尺寸相同时，则不需在同一位置即可填充，且构建方式支持初始几何绕某一点旋转构成另一几何^[10-12]，仅需定义旋转轴坐标、旋转轴垂直的坐标轴及旋转角度即可。

RMC的空间构成方式与MONK类似，但其几何描述由曲面卡与栅元卡共同定义，系统包括三类基本的几何描述单元，即曲面（SUFACE），栅元（CELL）和空间（UNIVERSE），栅元通过曲面和布尔运算来定义，再由一定数量的栅元构成不同空间^[13]。RMC中材料卡同样支持定义为原子密度或质量密度，但需给定具体总密度，核素按原子序数+质量数定义。几何部分支持对栅元的移动（MOVE），但其空间填充（FILL）需要几何的局部坐标相同，且栅元的几何变换采用变换矩阵的输入方式，仅支持栅元在其所属的局部坐标系中绕xyz轴旋转。

在使用MONK程序构建的乏燃料贮存格架模型中，大量使用了其特定的填充、旋转方式。然而，基于两种程序建模方式的差异，对于该类模型要保证转换前后两种程序的几何模型完全一致，需要对MONK程序输入卡中的几何进行解析，再重新构建成适用于RMC程序的几何，具体步骤为：

第一步：打开MONK输入文件，保留MONK输入文件中的多重变量，并改为适用于RMC多重变量定义的输入方式后将其写入输入卡开头处。

第二步：提取几何部分中几何类型与几何参数，通过基本运算与参数位置来调整构建MACROBODY与SURFACE曲面卡。

第三步：MONK特殊填充中填充几何与被填充几何尺寸相同时局部坐标可不同，根据几何是否被特殊填充分为2种处理方式：①未被特殊填充，则被填充几何的位置不变，直接进行填充；②被特殊填充，首先在填充几何位置构建被填充几何并先进行填充，将填充后的几何平移与旋转至被填充几何的位置，最后，衔接对应曲面组成各个空间生成栅元卡。

第四步：材料部分核素若未定义到具体同位素，则需先利用自然丰度展开，再将核素符号替换为原子序数。MONK程序中材料若按原子密度定义则总密度可简单写为0，因此这种情况下转化后此材料总密度为材料内所有核素密度之和。最后，将转化后的各个材料结合构成材料卡。

2.   程序实现

2.1   材料模块转换

MONK输入卡材料部分中混合物按原子比例或质量比例定义时核素为选择性后跟质量数以确定具体同位素，若未定义到具体同位素，在构成RMC材料卡时可将核素利用自然丰度展开，以此得到各同位素密度（$ \rho_{\mathrm{m}} $），其表达式为：

式中，M_z为核素密度；A为该核素各同位素的对应丰度。

RMC材料卡中，density > 0 表示原子密度，单位为10²⁴原子/cm³；density < 0 表示质量密度，单位为10²⁴ g/cm³。在展开后，按质量密度定义的核素密度前需加负号，单位与MONK相同。此外，MONK程序中在按原子比例定义材料时的材料密度可简写为0，对此转换中将对应材料中所有核素的原子密度相加后得出材料密度写入RMC材料卡对应材料密度处，其表达式为：

式中，$ {\rho }_\mathrm{T} $为材料密度；$ {\rho }_{i} $为材料中各核素密度。

2.2   平移与旋转处理

针对MONK建模中同尺寸几何填充不需在相同位置、且支持几何绕任意轴旋转所形成另一几何的特点，在转化中对此提出解决方式为：先填充栅元，再整体平移或旋转到对应位置。

RMC支持对UNIVERSE和CELL的平移变换和旋转变换。平移变换的表达式为：

式中，$ r=(r_x,r_{y},r_{\textit{z}}) $和$ r'=(r_x',r_y',r_{\textit{z}}') $分别为变换前和变换后的空间任意一点的位置坐标；$ m = ({m_x},{m_y},{m_{\textit{z}}}) $为平移变换向量。

变换过程中首先确定填充几何与被填充几何对应顶点，将两者做差，计算出位置坐标差值，在RMC中的填充几何位置先构建出被填充几何，并进行填充，随后再通过得到的坐标差值将填充后的几何进行栅元平移到被填充几何的位置。

式中，$ p=(p_x,p_y,p_{\textit{z}}) $为被填充几何位置；$ p_{_f}=(p_{fx}, p_{fy},p_{f{\textit{z}}}) $为填充几何位置。

RMC旋转变换可以绕xyz轴之一，其表达式为：

式中，R为旋转变换矩阵。

RMC输入的是旋转变换矩阵的转置矩阵R^T，其参数按照以下方式定义：给定某直角坐标系(x,y,z)，经过旋转变换后得到新坐标(x′,y′,z′),则R^T可表示为：

式中，C_x′x表示x和x′两个坐标轴之间的夹角余弦，依此类推。

针对MONK输入卡中绕某坐标轴旋转所构建的几何，将其对应顶点坐标通过绕点旋转坐标公式进行计算，则可得出旋转后的对应顶点坐标。若是对被填充几何进行绕点旋转，RMC中同样在填充几何位置先构建出被填充几何，并进行填充，随后计算出绕点前被填充几何与填充几何的坐标差值，将得到的坐标差值代入绕点旋转坐标公式进行计算，求解出填充几何位置到绕点旋转后的被填充几何位置需要的移动值（MOVE= x′ y′ z′），可将几何填充后移动至被填充几何旋转后的位置，最后在其栅元中加入对应旋转变换矩阵完成旋转。

以z轴为轴向的绕点旋转坐标公式为：

式中，$ \Delta $x、$ \Delta $y分别为填充几何与被填充几何x、y的坐标差值；x₀、y₀为旋转轴坐标；$ \alpha $为旋转角度；示意图如图1所示。

图  1  移动计算示意图

(x₀,y₀) —旋转轴坐标；(Δx,Δy) —被填充几何坐标与填充几何对应x、y的坐标差值；α—旋转角度

Figure  1.  Schematic Diagram of Mobile Computing

下载: 全尺寸图片 幻灯片

2.3   几何模块转换

MONK程序中定义的栅元分为2类：PART和HOLE，通常采用HOLE填入PART中定义的基本组件内的方式来构建模型。MONK和RMC两程序支持的基本几何体相同，但定义方式有所区别，对应关系如表1所示。

表  1  MONK与RMC宏体定义方式区别表

Table  1.  Difference between MONK and RMC Macro Definition Methods

宏体类型	MONK定义方法	RMC定义方法
球体	SPHERE X0 Y0 Z0 R	SPH X0 Y0 Z0 R
长方体	BOX Xa Ya Za Xl Yl Zl	RPP Xmin Xmax Ymin Ymax Zmin Zmax
圆柱体	(X/Y/Z)ROD Vx Vy Vz R H	RCC Vx Vy Vz hx hy hz R
圆台	(X/Y/Z)CONE Vx Vy Vz RB RT H	TRC Vx Vy Vz hx hy hz RB RT
圆环体	(X/Y/Z)TORUS Vx Vy Vz RH RR	TORUS u v w Vx Vy Vz RH RR
旋转椭球	(X/Y/Z)REL Vx Vy Vz r H	ELL Vx Vy Vz ax ay az r
圆柱扇体	(X/Y/Z)SEC Vx Vy Vz R1 R2 H T1 T2	SEC Vx Vy Vz hx hy hz R1 R2 T1 T2
六棱柱	(X/Y/Z)HEX Vx Vy Vz D H	RHP Vx Vy Vz hx hy hz s1 s2 s3
三棱柱	(X/Y/Z)TRI Xa Ya Za L1 L2 H TH1 TH2	WED Xa Ya Za Ax Ay Az Bx By Bz hx hy hz
注：(X/Y/Z)代表任意轴向；X0 Y0 Z0为宏体中心坐标；R为半径；Xa Ya Za为底面某一顶点坐标；Xl 、Yl 、Zl分别为X、Y、Z方向长度；Xmin、Xmax、Ymin、Ymax、Zmin、Zmax分别为垂直于对应轴的最小最大边界面；Vx Vy Vz 为底面中心坐标；H为轴向高度；hx hy hz为高度向量；RB、RT分别为底面、顶面半径；RH、RR分别为圆环半径、切面圆半径；u v w为圆环方向向量；r为垂直于旋转轴的半轴长度；ax ay az 为椭球体长半轴向量；R1、R2分别为内径、外径；T1、T2分别为坐标轴与圆柱扇体2个侧面的夹角；D为六棱柱中心到第一个面的距离；s1 s2 s3为六棱柱中心到第一个面的向量；L1、L2分别为三棱柱底面顶点到第一、第二条边的距离；TH1 、TH2分别为三棱柱底面第一、第二条边相对于轴的角度；Ax Ay Az、Bx By Bz分别为三棱柱底面顶点到第一、第二条边的向量

下载: 导出CSV 

| 显示表格

2.3.1   PART转换

PART由一组基本几何体构成，按其内部几何体的包含方式分为不同类型，在乏燃料贮存格架的模型中主要使用的类型为NEST、CLUSTER和OVERLAP，区别如图2所示。其区别为：NEST中几何体按定义顺序层层包裹；CLUSTER中除最后一个以外的几何体互相都没有重叠，且最后一个几何体充当容器包裹前面所有的几何体；OVERLAP中每个几何体都与定义列表中前面所有的几何体重叠，且最后一个几何体充当容器包裹前面所有的几何体。

图  2  不同类型PART示意图

Figure  2.  Schematic Diagram of Different Types of PART

下载: 全尺寸图片 幻灯片

模型转换过程为：①提取几何部分的关键词与参数，分别重构组成适用于RMC曲面卡与栅元卡的内容，转化后在曲面卡中使用MACROBODY下的宏体描述，在栅元卡中使用UNIVERSE下的栅元描述，两者皆以1为起始编号；②针对MONK程序中PART的不同类型，在RMC栅元卡中统一使用UNIVERSE描述，识别最后一个NEST定义为UNIVERSE 0，同时配合生成的曲面在栅元卡中通过布尔运算来实现不同类型来的几何表达，如图3所示。

图  3  布尔运算表达实现对应PART的方式

Figure  3.  Boolean Operation Expression Implements the Corresponding PART

下载: 全尺寸图片 幻灯片

在栅元的填充中需有一个最大面充当容器的作用，因此栅元所用曲面为填充栅元的最大曲面，随后提取填充关键词与参数改写入栅元卡，最后在每个栅元末尾添加2.2节中计算出的MOVE与ROTATE值完成栅元卡的生成。

2.3.2   HOLE转换

乏燃料贮存格架的模型中主要使用的HOLE几何类型为PLATE、GLOBE和LATTICE等。PLATE HOLE几何示意图如图4所示，为一维板几何形状，在其他2个方向上是无限的。GLOBE HOLE几何示意图如图5所示，为一组绕z轴的同心圆柱面，圆柱面之间层层包裹，在z方向是无限的。PLATE与GLOBE会被其填充的容器缩减到有限尺寸，一般用于材料划分。LATTICE HOLE几何为含有包层棒的矩形阵列，这些阵列周期性重复形成LATTICE（图6）。

图  4  PLATE示意图（z为轴向）

Figure  4.  Schematic Diagram of PLATE (z is axial)

下载: 全尺寸图片 幻灯片

图  5  GLOBE示意图

Figure  5.  Schematic Diagram of GLOBE

下载: 全尺寸图片 幻灯片

图  6  LATTICE示意图

Figure  6.  Schematic Diagram of LATTICE

下载: 全尺寸图片 幻灯片

同样提取将各个HOLE几何对应参数，分别重构组成适用于RMC曲面卡与栅元卡的内容。转化后在曲面卡中使用SURFACE下的曲面描述，编号衔接宏体部分；在栅元卡中使用UNIVERSE下的栅元描述，编号衔接PART部分。由于HOLE几何在MONK程序中是单独编号的，且最后一个PART为UNIVERSE 0，故HOLE几何转化后的UNIVERSE编号起始为输入卡PART的数量减1，在此UNIVERSE的填充中所填充UNIVERSE的编号为PART数量加此HOLE几何的编号。

PLATE与GLOBE分别转换为RMC曲面卡中的平面与圆柱面，并写入对应参数。随后在其UNIVERSE下的栅元中衔接对应曲面，使用曲面的布尔定义完成栅元的划分并指定对应材料。若此UNIVERSE被填充到重复结构中，需将此UNIVERSE移动至重复结构建立的网格中心，移动值为栅距的一半。

LATTICE HOLE转换为RMC中四边形重复结构，此部分不需生成曲面卡。四边形重复网格建立在xyz坐标系，坐标原点建立在第一个网格的左下角点，因此在建立完成后需整体移至中心位置处，以z轴为轴向的MOVE值计算方法为：

式中，N_x、N_y分别为组件在x、y方向的数量；p_x、p_y分别为x、y方向的栅距。

3.   算例验证与结果分析

为了测试模型转换是否准确，分别使用MONK与RMC展开临界分析计算，将得出的有效增殖因子进行比较来验证结果的正确性。

由于MONK程序控制卡中裂变源指定在裂变材料内，源收敛速度较快，可在20代内收敛；而RMC裂变源需指定具体位置，对于乏燃料贮存格架等较大模型需要更多非活跃代数达到源收敛，因此两程序无法统一代数。MONK设置终止计算条件为标准差小于0.0003，RMC计算条件为总代数500代，非活跃代数200代，计算使用的截面数据为 JEFF3.1 数据库。

结合统计误差(σ)用来确定结果是否在统计上具有一致性，当差异在±3σ区间内（置信度为99.6%）时，可以认为结果在统计上是相同的。σ定义为：

式中，σ₁ 和 σ₂ 分别为MONK和RMC 的k_eff标准偏差。

3.1   乏燃料贮存格架一区模型

乏燃料贮存格架一区又分为4个区域，其中3个区域放置6×7型格架，1个区域放置坠落组件。每个乏燃料组件包含17×17根燃料棒，燃料富集度为1.6%，慢化剂为纯水。该算例的MONK输入卡具有 26 个PART几何、19个HOLE几何、377 个宏体以及8个材料，其中又存在1个含有10个参数的多重变量，代表乏燃料组件坠落至格架拐角不同距离值。

转换后的RMC输入卡具有472个栅元、429个曲面、8个材料以及1个多重变量。图7分别是乏燃料贮存格架一区转换前与转换后的模型示意图，可以看出转换后几何模型并未变化，说明在几何模型方面转换结果可靠。

图  7  乏燃料贮存格架一区转换前后两程序几何示意图

Figure  7.  Geometric Schematic Diagram of the Two Codes before and after the Conversion of the Spent Fuel Storage Rack in Zone 1

下载: 全尺寸图片 幻灯片

由于此算例模型较大，当坠落组件与格架存在一定距离时使用单一点源很难达到收敛，故在RMC中设置初始源为与格架区域同尺寸的长方体源，两程序统一每代30000个粒子。表2 展示了该算例中转换前后两程序每个变量的计算结果。从表2中结果可以看出，转换前后两程序的k_eff偏差最大为90pcm（1pcm=10⁻⁵），最小为7pcm，均在±3σ内，且两程序的k_eff随坠落距离变化趋势符合较好，如图8所示，说明对于k_eff的计算，转换结果可靠。

表  2  乏燃料贮存格架一区转换前后的k_eff

Table  2.  k_eff before and after the Conversion of the Spent Fuel Storage R ack in Zone 1

坠落 距离/cm	MONK keff	MONK 标准差	RMC keff	RMC 标准差	∆keff	3σ
0	0.64060	0.000300	0.64137	0.000202	0.00077	0.00109
1	0.64060	0.000300	0.64110	0.000195	0.00050	0.00107
2	0.64000	0.000300	0.64090	0.000215	0.00090	0.00111
3	0.63990	0.000300	0.63964	0.000212	−0.00026	0.00110
5	0.63890	0.000300	0.63936	0.000204	0.00046	0.00109
7	0.63870	0.000300	0.63877	0.000206	0.00007	0.00109
10	0.63780	0.000300	0.63861	0.000214	0.00081	0.00111
15	0.63790	0.000300	0.63755	0.000214	−0.00035	0.00111
20	0.63730	0.000300	0.63784	0.000219	0.00054	0.00111
30	0.63730	0.000300	0.63790	0.000199	0.00060	0.00108

下载: 导出CSV 

| 显示表格

图  8  乏燃料贮存格架一区k_eff随坠落距离变化趋势图

Figure  8.  Trend Diagram of k_eff in Zone 1 as a Function of Falling Distance

下载: 全尺寸图片 幻灯片

3.2   乏燃料贮存格架二区模型

乏燃料贮存格架二区放置8×10型格架。每个乏燃料组件包含17×17根燃料棒，燃料富集度为1.6%，慢化剂为纯水。该算例的MONK输入卡具有 5 个PART几何、4个HOLE几何、289 个宏体以及13个材料，其中又存在1个含有10个参数的多重变量，代表不同的水密度。

转换后的RMC输入卡具有472个栅元、429个曲面、13个材料以及1个多重变量。图9分别是乏燃料贮存格架二区转换前与转换后的模型示意图，可以看出转换后几何模型并未变化，说明在几何模型方面转换结果可靠。

图  9  乏燃料贮存格架二区转换前后两程序几何示意图

Figure  9.  Geometric Schematic Diagram of the Two Codes before and after the Conversion of the Spent Fuel Storage Rack in Zone 2

下载: 全尺寸图片 幻灯片

RMC中设置初始源为几何中心位置的点源，两程序统一每代20,000个粒子。表3展示了该算例中转换前后两程序每个变量的计算结果。从表3中结果可以看出，转换前后两程序的k_eff偏差最大为110pcm，最小为1pcm，均在±3σ内，且两程序的k_eff随水密度变化趋势符合较好，如图10所示，说明对于k_eff的计算，转换结果可靠。

表  3  乏燃料贮存格架二区转换前后的k_eff

Table  3.  k_eff before and after the Conversion of the Spent Fuel Storage Rack in Zone 2

水密度/(g·cm−3)	MONK keff	MONK 标准差	RMC keff	RMC 标准差	∆keff	3σ
1.0	0.76440	0.000300	0.76439	0.000259	−0.00001	0.00119
0.9	0.75130	0.000300	0.75157	0.000252	0.00027	0.00118
0.8	0.73600	0.000300	0.73710	0.000267	0.00110	0.00120
0.7	0.71840	0.000300	0.71900	0.000249	0.00060	0.00117
0.6	0.69780	0.000300	0.69750	0.000301	−0.00030	0.00127
0.5	0.67060	0.000300	0.67031	0.000263	−0.00029	0.00120
0.4	0.63430	0.000300	0.63536	0.000252	0.00106	0.00118
0.3	0.58620	0.000300	0.58572	0.000256	−0.00048	0.00118
0.2	0.51590	0.000300	0.51529	0.000227	−0.00061	0.00113
0.1	0.41830	0.000300	0.41832	0.000185	0.00002	0.00106

下载: 导出CSV 

| 显示表格

图  10  乏燃料贮存格架二区k_eff随水密度变化趋势图

Figure  10.  Trend Diagram of k_eff in Zone 2 as a Function of Falling Distance

下载: 全尺寸图片 幻灯片

4.   结　论

为了提高乏燃料贮存格架临界计算的建模效率和建模准确性，本文首先研究了MONK和RMC两程序模型构建的特点与差异，针对材料部分利用自然丰度展开未定义到具体同位素的核素，对于MONK中材料密度为0的材料，将材料中所有核素的原子密度相加计算其实际材料密度。其次，针对几何部分提出使用先填充栅元，再整体平移或旋转到对应位置的方式来解决MONK建模中同尺寸几何填充不需在相同位置的问题，具体平移值通过平移与旋转处理得出。采用布尔运算表达衔接对应曲面实现不同类型PART与HOLE的转换。最后，选取乏燃料贮存格架一区、二区为算例验证结果的正确性：一区、二区存在不同数量的几何、材料与变量，运行程序均成功完成转换；将转换前后的几何模型和k_eff进行对比，验证了乏燃料贮存格架临界计算模型转换方法和功能的正确性，可提高乏燃料贮存格架建模效率与准确性，同时验证了自主化蒙特卡罗程序RMC应用到乏燃料贮存格架临界计算中的正确性。