RELAP5缓发中子先驱核输运模型扩展及验证

李 锐 程懋松 戴志敏,

1（上海科技大学 上海201210）

2（中国科学院上海应用物理研究所 上海201800）

熔盐堆作为第四代核反应堆国际论坛Generation IV International Forum（GIF）提出的6种反应堆概念之一，因为其固有安全性、防核扩散能力以及较高的燃料利用率，在国际上受到越来越多研究者的关注［1-2］。在熔盐堆的设计和分析过程中，反应堆系统分析程序具有重要的作用，其开发和验证是国内外重点关注内容之一。但由于熔盐堆独特的中子动力学和热工水力学特性，传统的固态堆系统分析程序不能适用于熔盐堆，所以需要为熔盐堆开发新的系统分析程序，或者在已有的固态堆系统分析程序基础上进行改进［3］。美国爱达荷国家实验室开发的RELAP5是适用于轻水堆的系统分析程序，具有极强的建模灵活性，涵盖了绝大部分事故瞬态模型［4-6］。因此，在RELAP5基础上进行适应性修改是获得液态燃料熔盐堆安全分析工具的有效办法。

为了开展熔盐堆瞬态和安全分析，中国科学院上海应用物理研究所对RELAP5/Mod4.0程序进行了大量改进。焦小伟［7］在RELAP5/Mod4.0程序中添加了熔盐物性以及适用于球床堆的Wakao换热关系式，对固态熔盐堆进行了事故工况下的瞬态分析。姜淑颖等［8］针对RELAP5/Mod4.0程序无法计算流动熔盐与不凝气体接触的相关问题进行了改进和验证。施承斌等［9-10］采用体积比法和衰变法的缓发中子先驱核守恒方程扩展了RELAP5/Mod4.0固态点堆模型，使用时滞型Runge-Kutta数值方法求解，并进行了基准题验证。

因为在液态燃料熔盐堆中，堆芯产生的缓发中子先驱核（Delayed Neutron Precursors，DNP）会随着燃料盐的流动在整个一回路中循环并衰变，导致反应堆的中子动力学受到燃料盐流速的影响［11］，而RELAP5/Mod4.0程序原本的点堆中子动力学模型无法计算这种效应，所以需要针对这一点对RELAP5/Mod4.0程序进行适用性修改。虽然文献［9］已经扩展了RELAP5/Mod4.0程序的点堆模型，使程序能够准确计算液态燃料熔盐堆稳态时DNP流动对中子动力学带来的影响，但是当燃料盐流速变化时程序对DNP的计算仍不够准确，同时使用时滞型Runge-Kutta法会使RELAP5/Mod4.0程序无法自适应计算步长，导致计算效率降低。因此，在上述工作的基础上，依据液态燃料熔盐堆安全分析需求，为提高RELAP5/Mod4.0在液态燃料熔盐堆安全分析中的适用性和精确度，基于中国科学院上海应用物理研究所现有的RELAP5/Mod4.0程序，采用一维缓发中子先驱核输运模型和二阶Godunov数值计算方法进一步改进其点堆动力学模型，并以熔盐实验堆（Molten Salt Reactor Experiment，MSRE）实验数据作为基准，进行程序修改后的正确性验证。

1 数学模型

在RELAP5/Mod4.0的点堆模型中，堆芯中子被视作均匀分布。在之前的工作中，分别使用体积比法和衰变法对原有的点堆模型进行了修改，考虑了DNP在一回路中的分布对中子动力学的影响［9］。为了进一步提高缓发中子先驱核输运计算的精度和计算效率，在原来的零维体积比和衰变模型基础上，新增一维DNP输运方程，并和零维堆芯中子守恒方程共同组成的新点堆动力学模型。

堆芯中子守恒方程保持不变，根据堆芯中子守恒原理，有：

式中：t代表时间；p代表中子通量；r代表反应性；β代表缓发中子份额；l代表平均中子代时间；Nd代表缓发中子组数；λi代表第i组DNP的衰变常数；cc，i代表堆芯内第i组DNP的密度。

一维缓发中子先驱核守恒一般方程表示为：

式中：ci代表控制体中第i组DNP的密度；A代表通流面积；vf代表燃料盐流速；x代表流动方向；Dm，i代表扩散系数；∇·(D m，i∇c i(t))代表扩散项；βi代表第i组缓发中子的份额。

由于液体中的扩散系数极小，所以式（2）中的扩散项可以忽略。因此，式（2）忽略扩散项，分别表示成堆芯和堆芯外回路两个一维缓发中子先驱核输运方程：

式（1）、（3）和（4）共同构成了采用一维缓发中子先驱核输运模型的液态燃料熔盐堆点堆动力学模型。

2 数值求解方法

在包括一维缓发中子先驱核输运模型的新点堆中子动力学模型中，采用Rider提出的二阶Godunov数值方法求解DNP输运方程［12］。在以下推导中，下标或上标中的K、L和M代表控制体，j和j+1代表连接点，它们的关系如图1所示。

图1 控制体K、L、M以及连接点j和j+1的关系Fig.1 Relationships between control volume K,L,M and junction j and j+1

基于有限体积法，根据高斯定理，式（3）可被写作：

式（5）的数值解可由式（6）得到：

式中：VL代表控制体L的体积。表示连接点j第n时间步DNP的流量，可由式（7）得到：

式中：Φ(γ，1)代表限幅器；S j+1代表连接点j+1处DNP在流动方向上的梯度；θL代表不连续检测器；ωL代表Courant数项；(1+θLωL)代表压缩项，由以下方程求得：

3 程序修改及验证

3.1 RELAP5程序修改

在RELAP5/Mod4.0程序中，每一个时间步内都会先进行热构件模型的计算，然后进行热工水力模型计算，最后在计算点堆中子动力学时同时计算中子和DNP数量。修改点堆动力学模型之后，将零维中子守恒方程和一维DNP守恒方程分为两步计算，先使用本时间步的热构件、热工水力参数、以及前一时间步的DNP参数求解零维中子守恒方程，然后使用本时间步的中子数和热工水力参数计算一维DNP分布，计算流程如图2所示。

3.2 程序验证

美国ORNL进行的MSRE作为目前世界上唯一成功的并且有详细实验记录的熔盐堆，其运行过程进行了大量实验，积累了大量实验数据［13］。因此，改进的RELAP5/Mod4.0程序选择MSRE的启泵、停泵、自然循环和反应性引入的功率响应实验作为基准题对程序进行评估验证。

3.2.1 MSRE简介

图2 RELAP5计算流程图Fig.2 The calculation flowchart of RELAP5

MSRE是ORNL为发展民用熔盐堆技术而搭建的小型熔盐堆，满载功率为8 MW［14］。MSRE堆芯由石墨棒组件排列而成，石墨棒表面被加工有凹槽，形成熔盐流道。熔盐在堆芯被加热后被一回路泵送往双熔盐换热器，将热量传给二回路熔盐，然后流回堆芯。二回路熔盐从双熔盐换热器流出后进入空气散热器，将热量传给外部环境，然后再被二回路泵送回双熔盐换热器［15］。

3.2.2 RELAP5模型

使用改进的RELAP5/Mod4.0程序建模MSRE系统，采用的具体设计与运行参数如表1所示［15］，对应的堆芯区域划分图和系统节块图见图4和图5，建模所用的MSRE熔盐、石墨和结构材料参数如表2所示［17］，中子动力学参数如表3所示［18］。由于在4个基准题验证过程中熔盐的流量均为已知参数，所以，为了方便建模和最小化热工水力误差造成的影响，忽略了泵模型，将泵内熔盐的体积合并入管道体积，并采用RELAP5/Mod4.0中的时间控制连接点直接以流量随时间的变化值作为输入参数。另外，在4个基准题中，启泵和停泵是以235U为燃料，采用235U的中子动力学参数；自然对流循环和反应性引入实验以233U为燃料，使用233U的中子动力学参数。

图3 MSRE系统简要示意图[16]Fig.3 Simplified schematic of the MSRE system[16]

3.2.3 启泵和停泵实验验证

在MSRE启泵瞬态实验过程中，反应堆在零功率下运行，因此温度反馈可以忽略。初始状态下，燃料盐处于静止状态，堆芯处于临界状态。一回路泵从t=0时刻开始从零加速，燃料盐流速变化如图6所示［19］。随着熔盐流速增加，DNP流出堆芯，造成缓发中子损失。为保证堆芯临界，控制系统根据反馈开始提升控制棒引入反应性，并在第4～12 s之间达到2.941×10-4s-1的引入速度上限。随后当一回路中的DNP开始再次进入堆芯时，引入的反应性随之降低，在几次波动后趋于稳定。改进后的RELAP5/Mod4.0采用拟合的燃料熔盐流速作为输入（图6），分别计算了有超临界限制和无超临界限制条件下MSRE启泵50 s内控制棒补偿的反应性随时间的变化，并与MSRE实验值和文献［8］体积比模型（volume ratio model）和衰变模型（decay model）的模拟结果进行了对比，如图7所示。相比于体积比模型和衰变模型，采用一维DNP输运方程的新点堆模型更贴近实验值。在有超临界限制情况下，数值结果与实验结果吻合良好；在无限制条件下，数值结果出现一个尖峰，高于实验值。

在停泵实验中，同样由于堆芯处于零功率，温度反馈被忽略。初始状态下，堆芯处于临界稳态，燃料盐处于额定流速下。一回路泵在t=0时刻开始降低转速，燃料盐流量在20 s内逐渐降低为零，如图8所示［19］。由于燃料盐流量的减小导致DNP损失减小，因此控制棒引入的反应性也减小。RELAP5/Mod4.0模拟结果与MSRE实验值对比如图9所示，数值结果与实验结果吻合良好，并且比体积比模型和衰变模型更加准确。

表1 MSRE设计与运行参数Table 1 MSER design and operation parameters

图4 MSER堆芯区域划分Fig.4 Core dividing scheme of MSRE

图5 MSRE的RELAP5模型节块Fig.5 RELAP5 nodalization of MSRE

表2 MSRE热物理性参数Table 2 MSRE thermo-physical properties

表3 233U和235U的中子动力学参数表Table 3 Neutron kinetic parameters of 233U and 235U

3.2.4 自然循环实验验证

在自然循环实验中，控制棒棒位不变。一回路泵处于关闭状态，整个系统由手动控制的冷却系统驱动［20］。实验开始前系统处于稳态，堆芯功率为4.1 kW，燃料盐流量极低。实验开始后冷却系统开始提升散热量，并在每一次散热量提升后等待堆芯再次达到临界稳态后再进行下一次散热量提升。由于冷却系统散热量提升，导致流回堆芯的燃料盐温度降低，引起负温度反馈增加，导致堆芯功率增加和堆芯出口温度提升。堆芯进出口温差增大，引起燃料盐流量增加。改进后的RELAP5/Mod4.0程序计算结果与MSRE实验结果对比如图10所示，数值结果与实验结果吻合良好。

图6 启泵过程燃料盐流量随时间变化Fig.6 Variation of fuel salt flow during fuel pump startup transient

图7 启泵过程控制棒引入的反应性变化Fig.7 Variation of reactivity inserted during fuel pump startup transient

图8 停泵过程燃料盐流量变化Fig.8 Variation of fuel salt flow during fuel pump coast-down transient

图9 停泵过程引入的反应性变化Fig.9 Variation of reactivity inserted during fuel pump coastdown transient

图10 自然循环实验堆芯功率随时间变化Fig.10 Variation of core power during the natural circulation experiment

3.2.5 反应性引入的功率响应实验验证

ORNL执行了MSRE装载233U燃料情况下反应堆功率对反应性引入的响应实验［21］。分别采用MSRE在1 MW、5 MW和8 MW工况下反应性引入的反应堆功率响应实验来验证改进后的RELAP5/Mod4.0程序。

MSRE运行在8 MW功率工况下，阶跃引入1.39×10-4反应性，堆芯功率随时间的变化如图11所示。MSRE运行在5 MW功率工况下，阶跃引入1.9×10-4反应性，堆芯功率随时间的变化如图12所示。MSRE运行在1 MW功率工况下，阶跃引入2.48×10-4反应性，堆芯功率随时间的变化如图13所示。

如图11～13结果所示，在不同稳态功率运行工况下，引入相应的反应性，堆芯功率响应的RELAP5/Mod4.0计算值与MSRE的实验值趋势上基本一致，且RELAP5/Mod4.0计算值具有更高的保守度，满足安全分析的要求。

图11 8 MW下堆功率对1.39×10-4反应性阶跃变化的时间响应Fig.11 Power response to a 1.39×10-4 step insertion of reactivity at 8 MW

图12 5 MW下堆功率对1.9×10-4反应性阶跃变化的时间响应Fig.12 Power response to a 1.9×10-4 step insertion of reactivity at 5 MW

图13 1 MW下堆功率对2.48×10-4反应性阶跃变化的时间响应Fig.13 Power response to a 2.48×10-4 step insertion of reactivity at 1 MW

4 结语

为了进一步提高RELAP5/Mod4.0程序在液态燃料熔盐堆瞬态分析和安全评估中的适用性和准确性，新增了一维缓发中子先驱核输运模型和相应的带有梯度限制器的二阶Godunov数值计算方法，并使用MSRE启停泵、自然循环和反应性引入下堆芯功率响应实验基准题对改进的RELAP5/Mod4.0程序进行了对比验证。结果表明：RELAP5/Mod4.0程序计算值与MSRE实验值吻合良好，验证了改进后程序的正确性和适用性。改进的RELAP5/Mod4.0程序将为液态燃料熔盐堆瞬态和安全特性研究、安全分析和安全许可申请提供更加可靠的支撑工具，对液态燃料熔盐堆的设计和分析具有重要意义。