JMCT-S一次屏蔽计算源项生成功能开发

郑 征，黎 辉，丁谦学，上官丹骅

(1.上海核工程研究设计院，上海 200233；2.北京应用物理与计算数学研究所，北京 100190)

JMCT-S一次屏蔽计算源项生成功能开发

郑 征1，黎 辉1，丁谦学1，上官丹骅2

(1.上海核工程研究设计院，上海 200233；2.北京应用物理与计算数学研究所，北京 100190)

一次屏蔽计算对于评估工作人员的辐射剂量、确保反应堆压力容器(RPV)及堆内构件在整个反应堆寿期内的安全性以及防止混凝土屏蔽体及其外部部件和结构被过度活化具有重要的意义。对于一次屏蔽源项的处理，JMCT-S程序自带的源粒子抽样功能无法完全满足其计算需求。本文开发了JMCT-S程序的源项生成程序和源抽样子程序，并在秦山一期和CAP1400一次屏蔽计算模型上进行了验证和应用。数值结果表明，推导的理论模型和开发的程序是正确的，从而为后续提高一次屏蔽设计精度提供了基础。

JMCT-S；源项生成程序；一次屏蔽

一次屏蔽由围绕堆芯的堆内构件(主要包括围板、吊篮等)、水隙[包括围板与吊篮、吊篮与压力容器(RPV)之间的水隙等]、RPV和围绕着RPV的混凝土结构等组成。一次屏蔽具有降低堆芯和结构材料产生的放射性对工作人员的辐射、确保RPV及堆内构件在整个反应堆寿期内安全性、确保混凝土屏蔽体性能的稳定性和完整性，以及防止混凝土一次屏蔽外的部件和结构被过度活化等功能[1]。

国内外通常采用蒙特卡洛方法和离散纵标方法进行一次屏蔽设计。JMCT-S程序[2]是计算复杂三维几何结构中的粒子输运的大型多功能蒙特卡洛程序。JMCT-S自带的源描述功能可以设定多种类型的源粒子位置、能量和方向，常用的源包括点源、面源、体源、各向同性源、各向异性源、单能源和有能量分布的源等等。由于反应堆堆芯几何结构复杂，以及堆芯源项在空间分布(PIN-BY-PIN分布)和能量分布(不同燃耗水平下由于重核素的增加导致不同组件裂变份额产生变化，从而使得不同组件具有不同的裂变谱)上的复杂性，当前JMCT-S程序自带通用源描述无法实现精确模拟，因此有必要开展JMCT-S的源项生成程序和源抽样子程序的研究。

为了进一步完善JMCT-S程序的源项生成功能，提高JMCT-S程序的应用范围，本文研究了JMCT-S程序的源项理论模型，开发了相应的源项生成程序和源抽样子程序，实现了考虑裂变份额的三维全堆芯PIN-BY-PIN源分布，并且在秦山一期计算模型上进行了验证，在CAP1400和CAP1700一次屏蔽设计中进行了应用。

1 理论模型

对于堆芯源粒子抽样，首先根据三维功率分布、裂变核素份额、裂变谱、每次裂变产生粒子数和能量等计算出源粒子的各个参数的累积分布函数(CDF)，然后基于CDF进行抽样，确定源粒子的类型、权重、坐标、方向、能量和所在栅元等参数。其中方向上服从裂变源各向同性分布，因此不需要计算CDF，采用各向同性抽样。

1.1 累积分布函数计算方法

第m个燃料组件的能谱根据该组件的裂变核素份额、每次裂变产生粒子数和裂变谱等参数计算得到[3]：

(1)

式中：g,g′——能群编号；

m——组件编号；

n——核素编号；

p——粒子类型编号；

N——核素数目；

χ′(m,g,p)——第p种粒子、第g群、第m个组件粒子的产生概率；

χ(n,g,p)——第p种粒子、第g群、第n种核素粒子的产生概率。源粒子的能谱采用BUGLE[4]的能群结构；

f(m,n)——第n种核素、第m个组件的粒子裂变份额；

ν(n,p)——第p种粒子、第n种核素每次裂变释放粒子数。

根据各个燃料组件的能谱，可以计算出各个组件的能量CDF：

(2)

式中：E(m,g,p)——第p种粒子、第g群、第m个组件的CDF；

G(p)——第p种粒子的能群数目，p=1表示中子，G(1)=47，p=2表示光子，G(2)=20。

三维粒子源分布计算公式如下：

(3)

式中：S(q,p)——第p种粒子、第q个网格源强，粒子/(cm3·s)；

P(q)——第q个网格真实功率密度，MW/cm3，P(q)=P(i,j,k)，(i,j,k)表示网格编号；

C——单位转换因子，6.24×1012MeV/(s·W)；

K(m)——第m个组件每次裂变释放的能量，MeV。其中f(m,n)和K(m)根据组件燃耗和初始富集度通过六种核素的裂变份额以及每次裂变释放的总能量插值得到。

根据公式(3)计算所得三维PIN-BY-PIN源分布，可以计算得到相应的三维空间各个网格的CDF：

(4)

式中：q,q′,q″——网格编号；

Q——总网格数目；

D(q″,p)——第p种粒子、第q″个网格的CDF；

S(q,p)——第p种粒子、第q个网格的源强，粒子/(cm3·s)；

V(q)——第q个网格的体积，cm3。

堆芯裂变产生第p种粒子数目计算如下：

N′(p)=∑i,j,kS(q,p)·V(q)

(5)

式中：N′(p)——堆芯裂变产生第p种粒子的数目，粒子/s。

经过5年的艰辛努力，西达里亚油田“外练筋骨皮、内练一口气”，各项指标均出现逆转，并实现质的飞跃——原油日产量由2017年1月初的11吨升至目前的65吨，2018年折耗比2017年降低67%，2018年吨油操作成本比2015年降低63%，前三季度实现利润915万元。

根据堆芯裂变产生总的粒子数目，可以计算出粒子类型的CDF供JMCT-S抽样：

(6)

式中：p,p′,p″——粒子类型编号；

P——粒子类型数目，P=1表示中子，P=2表示中子和光子；

F(p)——第p种粒子的CDF。

1.2 源粒子抽样方法

考虑到变量之间的依赖关系，应按照CDF维度从低到高的顺序进行抽样。首先确定粒子类型，再确定粒子空间位置，最后确定粒子能量。方向采用各向同性抽样。

对于粒子类型，选取(0,1]区间内的随机数ξp，若F(p-1)<ξp≤F(p)，则表示源粒子为第p种粒子。

对于空间抽样，选取(0,1]区间内的随机数ξq，若D(q-1,p)<ξq≤D(q,p)，则源粒子在第q个网格内。源粒子的空间位置坐标(x,y,z)计算公式如下：

(7)

式中：x,y,z——源粒子的空间位置坐标，cm；

ξq1,ξq2,ξq3——(0,1]区间内的随机数。

对于能量抽样，选取(0,1]区间内的随机数ξg，若E(g-1,q,p)<ξg≤E(g,q,p)，则源粒子在第g群内。根据下式计算出源粒子的初始能量：

(8)

式中：E′——源粒子的初始能量，MeV；

ξg1——(0,1]区间内的随机数。

2 验证及应用

基于上文介绍的源项理论模型，本文编写了JMCT-S源项生成程序及源抽样子程序，并且在秦山一期计算模型进行了验证，在CAP1400和CAP1700一次屏蔽设计中进行了应用。

2.1 秦山一期辐照监督管中子注量率计算

基于参考文献[5]建立的秦山一期计算模型如图1所示，本文采用新开发的源项程序进行了计算。

图1 秦山一期计算模型示意图Fig.1 Geometry model of Qinshan No.1 reactor

表1给出了在额定功率下第六根辐照监督管中子注量率计算结果。从表1中可以看出，本文计算结果和测量值吻合良好，堆芯中平面结果相对误差在5%以内，上焊缝结果相对误差在20%以内。偏差的主要来源包括数据库、计算模型和程序带来的偏差等。

表1 秦山一期辐照监督管中子注量率Table 1 Capsule neutron fluence rate of Qinshan No.1 reactor

2.2 CAP1400压力容器快中子注量率计算

基于堆芯几何结构等参数，本文建立了JMCT-S的CAP1400计算模型如图2所示。图3给出了压力容器内表面、方位角方向从0°到45°、轴向上从堆芯活性区中平面到堆芯活性区顶部区域的快中子(>1.0MeV)注量率分布。

图2 CAP1400计算模型示意图Fig.2 Geometry model of CAP1400

图3 CAP1400压力容器内表面快中子注量率分布Fig.3 Fast neutron fluence rate distribution on RPV inner surface of CAP1400

基于离散纵标方法的2D/1D注量率综合法在工程上被广泛用于计算RPV中平面高度的快中子注量率[7]，因此选择DORT[8]程序的计算结果作为参考。其中DORT计算采用BUGLE-96数据库，中子划分47群，求积组离散阶数为S8，网格划分数目分别为414×202(二维Rθ计算)，401×231(二维RZ计算)和401(一维R计算)，源项生成模型来自SORCERY程序[9]。

表2给出了压力容器内表面最大快中子注量率。从表中可以看出，JMCT-S和DORT计算结果吻合一致，相对误差在2%以内。

表2 CAP1400压力容器内表面最大快中子注量率Table 2 Maximum fast neutron fluence rate on RPV inner surface of CAP1400

2.3 CAP1700压力容器快中子注量率计算

基于已有的CAP1700计算模型，本文采用新开发的源项程序进行了计算。表3给出了压力容器内表面最大快中子注量率(>1MeV)计算结果。从表中可以看出，当采用相同的235U裂变谱时，本文计算结果和参考结果吻合良好，相对误差在2%左右。当采用混合裂变谱时，本文计算结果比采用235U裂变谱时增大了大约12%。其中TORT计算采用BUGLE-96数据库，中子划分47群，光子划分20群，求积组离散阶数为S8，网格划分数目分别为199×90×99(三维1/8几何RθZ计算)，源项生成模型来自SORCERY程序[9]。

表3 CAP1700压力容器内表面最大快中子注量率Table 3 Maximum fast neutron fluence rate on RPV inner surface of CAP1700

混合裂变谱包括235U、238U、239Pu、240Pu、241Pu和242Pu 6种核素。与235U每次裂变释放的能量相比，238U、239Pu和241Pu每次释放的能量更多；与235U每次裂变释放的中子数目相比，238U、239Pu和240Pu和241Pu每次裂变释放的中子数目更多；与235U裂变谱相比，238U裂变谱更软，239Pu、240Pu和241Pu的裂变谱更硬[3]。这些原因导致采用混合裂变谱时计算结果较采用235U裂变谱时偏大。因此，本文的计算结果是合理的。

3 结论

本文推导了能量、空间和粒子类型等参数的累积分布函数计算公式，给出了源粒子能量、空间和粒子类型等参数的抽样方法，编写了相应的程序模块，并在秦山一期和CAP1400计算模型上进行了初步应用验证。数值结果表明，本文推导的理论模型和开发的程序是正确的，从而为后续提高一次屏蔽设计精度提供了基础。

[1] 傅守信，胡建军.秦山核电二期工程反应堆主屏蔽设计[J]. 核动力工程，2003，24(S1)：64-67.

[2] JMCT-S 用户使用手册[R]. 北京应用物理与计算数学研究所，2015.

[3] A. Vasiliev，H. Ferroukhi，M. A. Zimmermann，et al.，Development of a CASMO-4/SIMULATE-3/MCNPX Calculation Scheme for PWR Fast Neutron Fluence Analysis and Validation against RPV Scraping Test Data[J]. Annals of Nuclear Energy，2007，34：615-627.

[4] BUGLE-96，Coupled 47 Neutron，20 Gamma-Ray Group Cross Section Library Derived from ENDF/B-VI for LWR Shielding and Pressure Vessel Dosimetry Applications[R]. Oak Ridge National Laboratory，1999.

[5] 夏春梅，第六根辐照监督管中子注量计算报告[R]. 上海核工程研究设计院，2014.

[6] 郭玥，秦山320MWe机组反应堆压力容器第六根辐照监督管材料样品中子能谱和快中子注量测定报告[R]. 中国原子能科学研究院反应堆工程研究设计所，2014.

[7] J. A. Kulesza，Comparison of Three-Dimensional Flux Synthesis and Full Three-Dimensional Discrete Ordinates Methods for the Calculation of Reactor Cavity Bioshield Heat Generation Rates[J]. Radiation Transport and Protection，2011，175：228-237.

[8] W. A. Rhoades，R. L. Childs，The DORT Two-Dimensional Discrete Ordinates Transport Code[J]. Nuclear Science and Engineering，1988，99：88-89.

[9] SORCERY User Manual[R]. Westinghouse Electric Company，2006.

[10] W. A. Rhoades，D. B. Simpson，The TORT Three-Di-mensional Discrete Ordinates Neutron/Photon Transport Code[R]. Oak Ridge National Laboratory，1997.

DevelopmentonPrimaryShieldingCalculationSourceGenerationFunctionforJMCT-S

ZHENGZheng1，LIHui1，DINGQian-xue1，SHANGGUANDan-hua2

(1.Shanghai Nuclear Engineering Research and Design Institute，Shanghai，200233，China；2. Institute of Applied Physics and Computational Mathematics，Beijing，100190)

Primary shielding calculation plays an important role in evaluating the radiation dose of personnel，insuring the safety of reactor pressure vessel(RPV)and reactor internals during the whole lifetime，and preventing the concrete shielding as well as the outer components and structures from an excess of activation. The source particle sample function of the JMCT-S code cannot fully satisfy the calculation requirements for generating the source term of primary shielding calculation. In this paper，a source term generation code and a source sample code for JMCT-S are developed. Validation and application on primary shielding calculation model on Qinshan No.1 nuclear power plant and CAP1400 reactor are performed. Numerical results show that the theoretical model and the codes are both correct，and can be used to improve the accuracy of primary shielding calculation.

JMCT-S；source generation code；primary shielding

2017-03-20

郑 征(1985—)，男，陕西人，工程师，博士，现从事核电厂屏蔽设计工作

TL328

A

0258-0918(2017)05-0805-05