三维并行程序JSNT对HBR-2装置的屏蔽计算与分析

杨 超，程汤培，邓 力，付元光，温丽丽

(1.北京应用物理与计算数学研究所，北京 100094；2.中物院高性能数值模拟软件中心，北京 100088)

基于并行自适应结构网格应用支撑软件框架JASMIN[1]，中物院高性能数值模拟软件中心研制了三维并行离散纵标(SN)屏蔽计算程序JSNT[2]，主要应用于核电厂的屏蔽计算。JSNT程序基于区域分解实现了大规模并行计算和多种迭代加速算法[3]，具有可视化建模和后处理的功能，拥有对诸如反应堆、厂房等大型装置进行精确建模和高分辨率计算的能力[4]。H. B. Robinson 2号机组(HBR-2)是热功率为2 300 MW的压水堆装置，对外公开发布了辐照监督管处和中子剂量测量仪处的比活度测量值，是屏蔽计算程序校验的经典算例。本文利用JSNT对HBR-2装置进行屏蔽计算与分析。首先利用可视化建模程序JLAMT建立HBR-2的精细模型，然后基于HBR-2基准报告[5]中pin-by-pin的三维功率分布，利用源项转换程序SourceTrans计算得到屏蔽计算所需的三维裂变中子源项，最后利用JSNT程序计算辐照监督管处和中子剂量测量仪处的中子通量密度分布以及基准报告中给出的6个核素的放射性比活度，并与实验测量值进行对比。

1 离散纵标方法

稳态中子输运方程[6]：

E′)ψ(r,E′,Ω′)dE′dΩ′+S(r,E,Ω)

(1)

式中：ψ为中子角通量密度，cm-2·s-1；r为位置向量；E为能量变量；Ω为方向向量；Σt为宏观总截面，cm-1；Σs为从能量E′、角度Ω′散射到E和Ω的宏观散射截面，cm-1；Σf为宏观裂变截面，cm-1；χ为中子裂变谱；S为外源源强，cm-3·s-1。

对式(1)角度方向上采用SN方法离散，空间上采用差分，能量上采用分群近似，则有：

ωmμm(Ai+1/2,j,kψi+1/2,j,k,m,g-Ai-1/2,j,kψi-1/2,j,k,m,g)+

ωmξm(Bi,j+1/2,kψi,j+1/2,k,m,g-Bi,j-1/2,kψi,j-1/2,k,m,g)+

ωmηm(Ci,j,k+1/2ψi,j,k+1/2,m,g-Ci,j,k-1/2ψi,j,k-1/2,m,g)+

(Ai+1/2,j,k-Ai-1/2,j,k)(αm+1/2ψi,j,k,m+1/2,g-

αm-1/2ψi,j,k,m-1/2,g)+ωmΣt,i,j,k,gVi,j,kψi,j,k,m,g=

ωmVi,j,kQi,j,k,m,g

(2)

式中：ωm为方向m上的权重；g为离散后的能群号；i、j、k为离散后的空间网格索引；A、B、C分别为网格单元与R、θ、Z方向正交的网格面的面积，cm2；V为计算网格体积，cm3；Q为总源强(裂变源、散射源和外源之和)，cm-3·s-1；α为待定系数。

JSNT采用源迭代方法对式(2)进行求解，外迭代用于更新裂变源及特征值，内迭代用于更新散射源。基于JASMIN框架、计算单元以及边界上反射方向间数据依赖关系创建有向非循环图(DAG)，根据DAG实现区域分解的大规模并行计算。

2 源分布计算

对于裂变源项的计算，首先根据堆芯的三维棒功率分布、每次裂变产生中子数和能量、裂变份额以及裂变谱等参数计算棒上的源强，然后将棒与计算网格进行映射。对于直角坐标，采用体积权重方法进行映射，如图1a所示；对于柱坐标系，将棒离散为带权的点，它们在棒内均匀分布，然后根据柱计算网格所覆盖的离散点的数目，计算得到该计算网格内的源强，如图1b所示。

三维空间棒网格源分布计算公式如下：

S(i,j,k,g)=P(i,j,k)·C′·B′(m)/K′(m)·

(3)

式中：P为棒网格功率密度，W/cm3；C′为单位转换因子，即将MeV转换为W，6.24×1012MeV·s-1·W-1；B′为组件偏倚因子；m为组件号；K为每次裂变释放的能量，MeV；n为核素编号；f′为裂变份额；ν为每次裂变释放粒子数。

a——直角坐标网格；b—— 柱坐标网格图1 不同坐标系下的映射Fig.1 Remapping process of different mesh ordinates

3 计算模型与计算结果分析

3.1 计算模型

HBR-2快中子注量率实验基准题由ORNL在1997年发布，其目的是验证辐射屏蔽输运计算程序和数据库的正确性。HBR-2基准题报告[5]给出了详细的反应堆几何、材料、堆芯功率分布和反应堆运行历史等数据，并提供了第9循环末辐照监督管和辐照剂量测量仪处比活度的测量数据。HBR-2热功率为2 300 MW，堆芯装有157盒燃料组件，每盒燃料组件内有204个燃料元件，燃料元件布置方式为15×15。利用JLAMT软件对HBR-2装置进行建模，JLAMT是北京应用物理与计算数学研究所开发的一款面向粒子输运程序的自动、可视建模软件[7]，可将三维模型自动输出转换为GDML文件与JSNT对接。HBR-2整体几何结构如图2所示，堆芯中平面处的剖面如图3所示，堆芯外围结构包括围板、吊篮、热屏、反应堆压力容器和生物屏蔽等部件。

3.2 中子源项计算

基于HBR-2基准报告中提供的第9循环的pin-by-pin平均功率(图4a)，利用源转换程序SourceTrans将pin-by-pin功率转换为pin-by-pin裂变中子源，然后映射到RθZ圆柱几何的计算网格上。为便于与基准报告中的计算结果比较，采用基准报告中提供的47群中子裂变能谱，即采用235U和239Pu的均匀混合裂变谱[5]，其能群结构与47群的BUGLE-96截面库一致。图4b示出了堆芯中平面处第18群的裂变源的分布，可看出其分布与pin-by-pin平均功率分布是一致的。

图2 HBR-2几何结构图Fig.2 Configuration of HBR-2

图3 HBR-2剖面图Fig.3 Schematic section of HBR-2

a——功率分布；b——源分布图4 pin-by-pin功率分布和第18群裂变源分布Fig.4 Distributions of pin-by-pin power and fission neutron source for the 18th group

3.3 计算条件

基于47群的BUGLE-96截面库[8]计算HBR-2装置的辐照监督管处和中子剂量测量仪处的中子通量密度以及放射性比活度，并分析两种角度离散(S8和S16)和2种空间离散对计算结果的影响。2种空间离散的方式如下。

1) CASE1。201×70×96离散[9]，R方向堆芯网格由内到外依次为5.0、2.0、1.0、0.8 cm，结构材料为0.8～1.0 cm，探测位置为0.5～0.8 cm；θ方向在探测位置网格由0.5°逐渐递减到0.1°，其他位置为0.5°，R-θ剖面如图5a所示；Z方向测点位置为1.0 cm，其他位置为5.0 cm，R-Z剖面如图6a所示。

图5 中平面处R-θ方向网格Fig.5 Mesh of R-θ direction at core mid-plane

图6 R-Z方向网格Fig.6 Mesh of R-Z direction

2) CASE2。274×195×104离散，R方向堆芯网格由内到外依次为2.0、0.8 cm，结构材料为0.8～1.0 cm，探测位置附近为0.5 cm；θ方向探测位置为0.125°，其他位置为0.5°，R-θ剖面如图5b所示；Z方向测点位置为1.0 cm，其他位置为5.0 cm，R-Z剖面如图6b所示。

3.4 计算结果与分析

1) JSNT的计算结果与分析

针对201×70×96计算网格，截面采用P3和方向离散采用S8的方案，JSNT在某浪潮集群上利用1个节点完成并行计算，该集群1个计算节点上共有24个处理器，其计算时间为0.8 h。可发现，JSNT通过并行较大地提高了计算效率，降低了计算时间。图7示出了中子通量密度的三维空间分布，可看出中子通量密度沿径向逐渐衰减。

图7 中子通量密度分布Fig.7 Distribution of neutron flux density

2) JSNT与DORT计算结果的比较

图8、9分别示出了两计算点处的47群中子通量密度JSNT与DORT的计算结果的比较，其中DORT是基于BUGLE-96截面计算的结果[5]。在辐照监督管处，最大相对偏差小于20%，大于1 MeV的能群相对偏差均小于5%；在中子剂量测量仪处，大于1 MeV的能群相对偏差小于10%，在低能区存在较大偏差。

图8 辐照监督管处的中子通量密度分布Fig.8 Distribution of neutron flux density at location of surveillance capsule

图9 中子剂量测量仪处的中子通量密度分布Fig.9 Distribution of neutron flux density at location of neutron dosimeter

由于计算放射性比活度的反应均是阈能反应，低能区偏差大对结果影响较小。

表1、2列出了辐照监督管处和中子剂量测量仪处的比活度的计算值与测量值的比值。随着角度离散数和计算网格数的增加，计算结果趋于测量值。在辐照监督管处，JSNT与DORT计算的比活度较一致，与测量的相对偏差均小于20%。在中子剂量测量仪处，除237Np和238U的裂变反应外，计算的比活度与测量值之间的相对偏差也小于20%。与DORT计算结果相比，JSNT计算的237Np和238U的裂变反应的比活度偏低，这是由于DORT的计算中子通量在高能区相对较大引起的。在中子剂量测量仪处，237Np和238U的裂变反应比活度与测量值偏差较大，其原因可能是与堆芯距离较远，中子通量密度相对较小，使实验测量难度增加，且裂变反应在工程上是通过测量它的裂变产物来得到反应率的，其测量的不确定度本身比较大；另外屏蔽截面数据库和响应截面的精度也是导致偏差较大的原因之一。

表1 辐照监督管处比活度的计算值与测量值的比值Table 1 Ratios of calculated-to-measured specific activities at location of surveillance capsule

表2 中子剂量测量仪处比活度的计算值与测量值的比值Table 2 Ratios of calculated-to-measured specific activities at location of neutron dosimeter

4 结论

基于HBR-2基准，利用自主研发的并行程序JSNT完成了该装置的屏蔽计算。计算结果表明，当采用较多的网格数时，在辐照监督管处各核素的放射性比活度的计算值和测量值的相对偏差均小于20%，满足屏蔽计算的工程需求。但在中子剂量测量仪处，237Np和238U的裂变反应比活度的计算值和测量值存在较大偏差，有待进一步研究。