固态燃料钍基熔盐实验堆(Thorium-based Molten Salt experiment Reactor with Solid Fuel, TMSR-SF1)[1-3]预计是国际上首个将要建造的氟盐冷却球床高温堆。TMSR-SF1与现有反应堆之间的差异对核安全审查提出挑战[4]，为更好地支撑核安全评审工作，其核设计采用通用蒙特卡罗程序。有效缓发中子份额(βeff)和有效中子代时间(Λeff)是反应堆动力学中至关重要的参数，其计算的准确性事关TMSR-SF1的安全分析与动态特性。
βeff为缓发中子对裂变功率的贡献份额。由于缓发中子裂变谱与瞬发中子裂变谱有明显不同，导致单个中子引起的平均裂变几率不同，造成βeff与缓发中子份额(β)有明显差异。Λeff为引起下一代裂变的瞬发中子的平均寿命，而中子代时间(L)则还包含绝大部分只被吸收不裂变的瞬发中子寿命，Λ等于瞬发中子平均寿命(lp)除以keff。由于截面在能量以及空间上的差异，同样导致Λeff与Λ具有很大的不同。可见，中子动态参数的计算都有赖于中子裂变价值的判断，即需要共轭通量的求解。
共轭通量在确定论程序中容易求解，而适合复杂结构的蒙特卡罗方法在求解共轭通量方面具有较大不便，需要考虑粒子的逆向追踪。目前蒙特卡罗方法求解βeff方面有三类：1)基于本征值求解的方法[5-9]，如Prompt方法。Prompt方法通过开关蒙特卡罗多粒子输运程序(Monte Carlo N Particle Transport Code, MCNP)中β项，分别求解两次keff值，其相对差值即为βeff。该方法的近似处理为，在关闭β项时认为瞬发中子裂变谱即为总中子裂变谱，该近似为二阶小量。由于其准确性可靠，计算处理简单，被广泛使用。不足之处在于不能直接分别给出6组βeff数值和有效中子代时间。2)微扰法[10-11]。通过引入微分算符以及虚拟材料可以分别获得βeff和Λeff。其准确性可靠，但涉及的处理工作较大。3)共轭通量法[12]。蒙特卡罗方法的共轭通量求解关键在于使得输运的粒子携带重要的历史信息。通过对大量的粒子历史信息进行统计反过来评价源粒子对某一项的贡献价值。βeff与Λeff可以通过共轭通量与正向通量的积分公式获得。该方法准确性可靠，共轭通量的应用范围较广，但处理过程略复杂。
实际上，通过上述对βeff与Λeff的物理理解以及借鉴共轭通量的处理方法，可以通过粒子携带的历史信息直接统计6组有效缓发中子份额以及有效中子代时间，而不需要先加工共轭通量再求解的过程13]。如：仅对于判定为裂变反应时的粒子历史信息进行统计，求出其平均输运时间，即可得到有效中子代时间；对其源粒子类型分类整理，即可得到6组有效缓发中子份额。
本文先将直接统计方法应用于MCNP中，然后对该方法的准确性进行基准题验证。最后用于TMSR-SF1中子动力学参数计算，并对其变化规律进行分析。
1 直接统计方法直接统计方法在于逆向追踪一个裂变中子的源粒子信息。反过来，如果一个粒子在输运过程中携带了源粒子信息，则可以直接获取信息并加工成中子动态参数。为此，需要先熟悉MCNP中粒子的输运过程以及如何携带信息。
MCNP中粒子的输运主要在判断两类事件：碰撞和穿越栅元的边界。穿越栅元边界时中子信息不发生改变，直接统计法只考虑碰撞事件。以图 1为例，对于初始权重为wt0、初始时刻为t0的瞬发中子类型tp0，经过第一次碰撞后，损失权重到wt1（MCNP默认采用非类比碰撞方法来加大粒子径迹统计数，非类比碰撞认为，粒子按照碰撞中吸收占总碰撞的概率损失其权重，而不是抽样判定为吸收反应或散射反应）。同时判定被吸收的权重部分是否为裂变反应，若判定为是，则进一步判定产生的裂变中子信息并保存，图 1第一次碰撞产生第一组类型缓发中子。第二次的碰撞过程处理相同，并产生第二组类型缓发中子。第三次碰撞被判定为俘获吸收，没有裂变中子产生。第四次发生碰撞后，中子权重wt4已低于权窗下限，需进行轮盘赌注，赌注结果为消失。第四次碰撞的损失权重部分同样需进行吸收类型判定，结果为裂变反应并产生第二组类型缓发中子。在本次中子输运过程中，可以推测到权重为wt0的瞬发中子贡献了(wt0-wt2+wt3-wt4)个裂变中子。通过对每种类型的中子裂变贡献进行统计，便可获得MCNP这一代循环过程中的每组有效缓发中子份额。
|
图 1 MCNP中一个粒子的输运历史 Figure 1 Transportation history of one neutron in MCNP |
由于代与代之间的裂变中子类型具有一定的关联性，在单代粒子数较少的情况下，统计的涨落在代与代间的干扰尤为明显，可能导致结果朝某一方向偏离。一般为了使其充分发展，可追述该裂变中子的前几代源中子类型。本方法中根据测试分析采用追踪前15代源中子。
每一代中有效缓发中子份额的计算公式为：
| $\beta _{{\rm{eff}}}^m = \frac{{\sum\limits_{k = 1}^{15} {\sum\limits_i {\sum\limits_j {(w{t_j} - w{t_{j + 1}}){|_{{\rm{f}},t{p^m}}}} } } }}{{\sum\limits_{m = 1}^7 {\sum\limits_{k = 1}^{15} {\sum\limits_i {\sum\limits_j {(w{t_j} - w{t_{j + 1}}){|_{{\rm{f}},t{p^m}}}} } } } }}$ | (1) |
式中：m为缓发中子组类型，包括瞬发中子在内共7类；j为每一个中子发生的碰撞序号数；i为每一代的中子数；
| ${\Lambda _{{\rm{eff}}}} = \frac{{\sum\limits_i {\sum\limits_j {({t_j} - {t_0}) \times (w{t_j} - w{t_{j + 1}}){|_{\rm{f}}}} } }}{{\sum\limits_i {\sum\limits_j {(w{t_j} - w{t_{j + 1}}){|_{\rm{f}}}} } }}$ | (2) |
多代的累计结果为他们的算术平均值。
在上述计算过程中，需要考虑如何获得中子类型、输运时间等信息。关于输运时间的获取，可以定义一些全局变量，直接统计出发生裂变反应的平均输运时间。而对于中子类型等信息，MCNP本身并没有进行存储，需要考虑中子如何携带。MCNP的裂变中子存储在fso一维向量中，一个中子一般有6个信息，包括三维空间坐标fso(1:3)、能量fso(4)、栅元号fso(5)和粒子序号fso(6)。fso的向量长度为6乘以每代粒子数。为使得裂变中子携带前15代源中子类型信息，可将fso中单个中子信息位从6扩充到21，定义fso(6:20)为其源中子类型信息位。每当判定为发生裂变反应时，将该输运中子的14代源类型信息传递给新的裂变中子，并补充新产生的裂变中子类型，以此完成源中子类型信息的携带与更新。
由于对fso向量长度进行了修改，其同等变量的声明需同样进行修改。考虑新增的全局变量，总共对MCNP进行约30处源代码的修改。
2 基准题验证 2.1 基准题描述为验证直接统计法计算有效缓发中子份额的准确性，选用ICSBEP (International Criticality Safety Benchmark Evaluation Project)数据库基准题14]。其主要为裸燃料球的临界实验，燃料包括235U、233U和239Pu等。模型简介如下：
Godiva (heu-met-fast-001)：94%富集铀裸球，同位素有234U、235U和238U；Jezebel (pu-met-fast- 001)：95%钚裸球，核素有Ga、239Pu、240Pu和241Pu；Skidoo (233U-met-fast-001)：98% 233U裸球，核素有233U、234U、235U、238U；Topsy (heu-met-fast-028)：93%富集铀球外包一层天然铀层；Popsy (pu-met- fast-006)：94%富集度的钚球外包一层厚的天然铀层；Flattop23 (233U-met-fast-006)：98%富集233U球外包一层天然铀层；Big Ten (ieu-met-fast-007)：铀金属圆柱堆芯外包10%富集度235U反射层。
本文计算采用ENDF/B-VII.0库，TMSR-SF1高温下的数据由NJOY核数据处理系统加工而成[15]，MCNP计算每代粒子数100000，共900有效代。
2.2 准确性验证各基准题的计算结果见表 1。第2列数值为实验测得的有效缓发中子份额及测量误差；第3列为本文采用Prompt方法得到的计算值与实验值的相对偏差及统计误差；第4列为直接统计方法计算值与实验值的相对偏差及统计误差。三者的结果比较接近，差异在5%以内，个别例题与实验差异较大，这种结果同样出现在其他文献计算中[5-8]，可能与例题测量误差相关。文献[16]认为，数据库本身提供的缓发中子份额准确性在3%左右。可见，直接统计方法获得的结果基本可靠。
|
表 1 βeff不同计算方法与实验结果比较 Table 1 βeff comparison among calculations and experiment |
直接统计方法需要考虑统计的收敛性。以Godiva为例，如图 2所示，初始循环代βeff的波动较明显，在经历700代之后，趋于稳定，本文结果为1000代时的数值，确保了计算数据的收敛性。
|
图 2 Godiva βeff随中子循环数的收敛情况 Figure 2 Godiva βeff varied with cycle number |
除了总有效缓发中子份额之外，还需要对每组有效缓发中子份额进行比较。Prompt法还不能直接统计每一组缓发中子份额，但通过简单地修改MCNP代码就能实现。主要将总的有效缓发中子份额及裂变能谱抽样直接指定为某一组的有效缓发中子份额和相应裂变能谱。
本文对Godiva模型的6组有效缓发中子份额进行了计算，结果见表 2。Prompt法和直接统计法的各组有效缓发中子份额基本一致。但为了得到较为精确的结果，Prompt法必须投入上亿粒子数，同时还需要进行7次MCNP计算，且6组求和得到的有效缓发中子份额与表 1中Prompt方法获得的具有较大出入。而直接统计法只需不到1000万粒子数即可达到很高的精度，证明直接统计法在计算时间上具有很大的优势。
|
|
表 2 Godiva模型6组有效缓发中子份额 Table 2 6 group of effective delayed neutron fraction of Godiva |
关于有效中子代时间的验证基准题较少，其计算方法也比较直观，本文不做详述，将在TMSR-SF1例题中与中国核动力设计院采用MCNP共轭通量方法的计算结果进行比较。
3 TMSR-SF1动态参数计算 3.1 SF1模型描述TMSR-SF1额定热功率10 MW，堆芯采用球形燃料元件（燃料球），冷却剂为2LiF-BeF2。模型见图 3，设计主要参数见表 3。
|
图 3 固态燃料钍基熔盐实验堆计算模型 Figure 3 Calculation model of TMSR-SF1 |
|
|
表 3 固态燃料钍基熔盐实验堆堆芯总体参数 Table 3 Main parameters of TMSR-SF1 |
球床堆芯采用规则堆积建模（图 4）代替实际的随机堆积。燃料球TRISO (TRIstructural iSOtropic)颗粒采用完整的包覆结构建模，其在燃料区的排列采用简立方结构。球在堆芯的排列采用体心立方结构。
|
图 4 固态燃料熔盐实验堆球床规则堆积建模 (a)球体心立方排列，(b) TRISO简立方排列 Figure 4 Regular packing model of TMSR-SF1 (a) BCC lattice of pebble, (b) Simple cubic lattice of TRISO |
采用直接统计法对TMSR-SF1热态初装载状态下的模型进行了计算，各组有效缓发中子份额与MCNP共轭通量法结果基本一致(±3%)，见表 4。个别较大差异主要体现在第4组中，但仍在统计误差范围内。共轭通量法的误差较大，主要因为其统计的粒子数较少。
|
|
表 4 TMSR-SF1堆芯热态初装载有效缓发中子份额 Table 4 βeff of TMSR-SF1 at the beginning of life |
TMSR-SF1有效缓发中子份额随燃耗时间的变化见表 5。有效缓发中子份额呈现下降的趋势。主要因为239Pu的累积，寿期末239Pu的裂变贡献达到5%。能谱的变化也会导致次级影响，如初期氙毒的影响导致能谱硬化，有效缓发中子份额有所下降。
|
|
表 5 TMSR-SF1堆芯不同燃耗点有效缓发中子份额 Table 5 βeff of TMSR-SF1 at different burnup depths |
TMSR-SF1的有效中子代时间随时间变化见表 6。直接统计法与MCNP共轭通量方法获得的数据基本一致(±3%)，但随燃耗时间具有一定的增加趋势。239Pu在热区的裂变截面较大，或导致其中子代时间的增加。
|
|
表 6 TMSR-SF1堆芯不同燃耗点有效中子代时间 Table 6 Λeff of TMSR-SF1 at different burnup depths (μs) |
此外，MCNP直接输出的TMSR-SF1初态瞬发中子寿命约为890 μs，约是有效中子代时间的1.7倍，主要因为TMSR-SF1堆芯较小，反射层厚度较大，一部分中子在反射层中充分慢化，寿命较长，但其对裂变贡献有限。
TMSR-SF1的能谱较软，其有效中子代时间约是水堆的5倍，但比高温气冷实验堆HTR-10[17]的略低。较长的有效中子代时间，将有利于反应性的控制。
4 结语本文介绍了计算中子动态参数的蒙特卡罗直接统计方法，并将其植入到MCNP程序中。经多个国际基准题的模拟，计算结果与实验偏差在5%以内，验证了直接统计方法在中子动态参数计算上的准确性。运用直接统计方法计算了TMSR-SF1中的6组有效缓发中子份额与有效中子代时间，结果与MCNP共轭通量方法数据基本一致(±3%)，为直接统计方法在TMSR-SF1中进一步应用打下了基础。
| [1] |
冀锐敏, 严睿, 李晓晓, 等. 球形燃料元件中包覆颗粒的分布效应研究[J]. 核技术, 2017, 40(10): 100604. JI Ruimin, YAN Rui, LI Xiaoxiao, et al. Effect of TRISO-particles distributions in pebble fuel[J]. Nuclear Techniques, 2017, 40(10): 100604. DOI:10.11889/j.0253-3219.2017.hjs.40.100604 |
| [2] |
于世和, 严睿, 冀锐敏, 等. PB-FHR的控制棒布局设计及物理效应[J]. 核技术, 2018, 41(1): 010605. YU Shihe, YAN Rui, JI Ruimin, et al. Layout design and physical effect analysis of control rod in PB-FHR[J]. Nuclear Techniques, 2018, 41(1): 010605. DOI:10.11889/j.0253-3219.2018.hjs.41.010605 |
| [3] |
陈家豪, 张海青, 朱智勇. 10 MW固态燃料钍基熔盐堆稳态物理-热工耦合[J]. 核技术, 2017, 40(8): 080603. CHEN Jiahao, ZHANG Haiqing, ZHU Zhiyong. Coupled neutronic and thermal-hydraulic analysis of TMSR-SF1 at steady state[J]. Nuclear Techniques, 2017, 40(8): 080603. DOI:10.11889/j.0253-3219.2017.hjs.40.080603 |
| [4] |
左嘉旭, 高新力, 李朝君, 等. 10 MW固态燃料熔盐实验堆安全分析关键技术初步研究[J]. 核技术, 2017, 40(4): 040601. ZUO Jiaxu, GAO Xinli, LI Chaojun, et al. Preliminary research of safety analysis for 10-MW thorium-based molten salt reactor-solid fuel[J]. Nuclear Techniques, 2017, 40(4): 040601. DOI:10.11889/j.0253-3219.2017.hjs.40.040601 |
| [5] |
Meulekamp R K, Van der Marck S C. Calculating the effective delayed neutron fraction with Monte Carlo[J]. Nuclear Science and Engineering, 2006, 152(2): 142-148. DOI:10.13182/NSE03-107 |
| [6] |
Okajima S, Sakurai T, Lebrat J F, et al. Summary on international benchmark experiments for effective delayed neutron fraction (βeff)[J]. Progress in Nuclear Energy, 2002, 41(1): 285-301. |
| [7] |
Meulekamp R K, Van der Marck S C. Calculating the effective delayed neutron fraction using Monte Carlo techniques[J]. Nuclear Science and Engineering, 2006, 152: 142-148. DOI:10.13182/NSE03-107 |
| [8] |
Carta M, Dulla S, Peluso V, et al. Calculation of the effective delayed neutron fraction by deterministic and Monte Carlo methods[J]. Science and Technology of Nuclear Installations, 2011, 2011(1687-6075): 109-124. |
| [9] |
余纲林, 李成龙, 王侃, 等. 缓发中子有效份额的测量及计算方法[J]. 清华大学学报(自然科学版), 2012, 52(7): 896-900. YU Ganglin, LI Chenglong, WANG Kan, et al. Absolute measurements of the effective delayed neutron fraction[J]. Journal of Tsinghua University (Science and Technology), 2012, 52(7): 896-900. |
| [10] |
Nagaya Y, Mori T. Calculation of effective delayed neutron fraction with Monte Carlo perturbation techniques[J]. Annals of Nuclear Energy, 2011, 38(2): 254-260. |
| [11] |
王冠博, 刘汉刚, 王侃, 等. 反应堆动态参数的蒙特卡罗计算研究[J]. 核动力工程, 2012, 33(2): 116-122. WANG Guanbo, LIU Hangang, WANG Kan, et al. Calculation of reactor kinetic parameters with Monte Carlo method[J]. Nuclear Power Engineering, 2012, 33(2): 116-122. |
| [12] |
汪量子, 姚栋, 王侃. 在蒙特卡罗程序中统计共轭通量并计算中子动力学参数的方法[J]. 原子能科学技术, 2012, 46(5): 579-584. WANG Liangzi, YAO Dong, WANG Kan. Method of tallying adjoint fluence and calculating kinetics parameters in Monte Carlo codes[J]. Atomic Energy Science and Technology, 2012, 46(5): 579-584. |
| [13] |
Nauchi Y, Kameyama T. Proposal of direct calculation of kinetic parameters βeff and based on continuous energy Monte Carlo method[J]. Journal of Nuclear Science and Technology, 2005, 42(6): 503-514. DOI:10.1080/18811248.2004.9726417 |
| [14] |
Briggs J B, Scott L, Nouri A. The international criticality safety benchmark evaluation project[J]. Nuclear Science & Engineering the Journal of the American Nuclear Society, 2002, 145(145): 1-10. |
| [15] |
梅龙伟, 蔡翔舟, 蒋大真, 等. MCNP温度相关热中子散射数据库研制[J]. 核科学与工程, 2013, 33(4): 362-367. MEI Longwei, CAI Xiangzhou, JIANG Dazhen, et al. Development of temperature related thermal neutron scattering database for MCNP[J]. Nuclear Science and Engineering, 2013, 33(4): 362-367. |
| [16] |
Rudstam G, Finck P, Filip A, et al. Delayed neutron data for the major actinides[J]. NEA0WPEC-6 Report, 2002. |
| [17] |
黄晓津, 冯元琨. HTR-10核动力系统反应性扰动下动态特性的仿真研究[J]. 核动力工程, 2002, 23(3): 103-107. HUANG Xiaojin, FENG Yuankun. Simulation and study on reactivity disturbs dynamic character of HTR-10 nuclear power system[J]. Nuclear Power Engineering, 2002, 23(3): 103-107. |