崔 勇，陈金根,*，何 帆，李晓晓，蔡翔舟,*

(1.中国科学院 上海应用物理研究所，上海 201800； 2.中国科学院 先进核能创新研究院，上海 201800；3.中国科学院大学，北京 100049)

熔盐堆是6种第4代核能系统候选堆型中唯一使用液态燃料的反应堆，与固态燃料堆不同，熔融的氟化盐混合物(LiF，BeF2，ZrF4，UF4或ThF4)兼作冷却剂与燃料，充满整个一回路(含堆芯)。燃料流动使得缓发中子先驱核(DNP)在堆内的分布受到流场的影响。相应地，缓发中子的价值也会发生改变。部分具有较长半衰期的DNP会随燃料流出堆外并衰变[1]，从而导致有效缓发中子份额(βeff)降低。其次，裂变能直接释放于冷却剂中，因此传统的流动-传热计算方法不再适用。最后，熔盐堆适合采用钍铀燃料循环[2]，因而易裂变材料中可能包含233U和锕系核素(MA)[3]，这些重核的缓发中子数据与采用铀钚燃料循环的固态堆不同。这些特殊性使得熔盐堆βeff的计算较为复杂。

βeff的精确计算对研究熔盐堆的瞬态安全特性具有重要意义。如在失流事故中，缓发中子源项的再分布会向堆内引入正反应性，该反应性的量化需对不同流动工况下的βeff进行评估[3]；在反应性引入事故中，需确定所引入反应性与βeff的关系以得到反应堆的响应特性[4-5]。近年来，国际国内针对DNP的对流[6]、扩散[7]和湍流输运效应[8]对熔盐堆稳态与瞬态特性的影响开展了广泛研究。对于罐式熔盐快堆，由于其堆内无慢化材料，可采用解析方法[3]、数值方法[3,8-9]和蒙特卡罗方法[3]计算其βeff与DNP分布。对于通道式熔盐堆，由于存在石墨、氢化锆等慢化材料，其流动-传热方式具有特殊性[5,10]，对快堆提出的βeff计算方法不再适用。

本文利用解析方法和基于四阶多项式节块展开法与并联多通道模型的三维核热耦合程序TMSR3D[11]，对熔盐实验堆(MSRE)进行稳态计算，得到堆内DNP分布情况，研究入口流量、燃料的堆外流动时间等参数对βeff的影响。

1 计算模型与方法

1.1 计算模型

MSRE堆芯与组件几何结构[12]如图1所示。堆芯由正方形石墨组件构成，组件四周开孔作为燃料盐流道，所有孔道与上下腔室连通。堆芯四周为石墨反射层。表1、2分别列出MSRE的主要设计参数[4,12]和装载235U燃料时的缓发中子数据[4]。详尽的堆芯设计与热工参数参见文献[4]。

表1 MSRE设计参数Table 1 Design parameter of MSRE

a——堆芯纵截面；b——堆芯横截面；c——组件几何图1 MSRE堆芯与组件几何结构Fig.1 Core and assembly configuration of MSRE

1.2 计算方法

熔盐堆稳态下中子学计算模型为：

(1)

(2)

1.2.1解析方法 利用圆柱裸堆几何对MSRE的临界通量进行近似，记MSRE的堆芯半径和高度分别为R和H，可得到圆柱坐标下的通量[13]为：

(3)

其中：J0为零阶贝塞尔函数；φ0为待定系数，由反应堆功率决定。

仅考虑轴向流动，记流速为u，将式(3)代入(2)中得到：

(4)

对上式在(0，z)区间积分得到：

c0exp(-λiz/u)

(5)

其中，c0为堆芯入口z=0处的DNP浓度。

设堆芯中熔盐占堆芯总体积的比值恒为f，且堆芯出口处先驱核浓度均匀混合，其出口浓度ci，out可表示为：

(6)

设熔盐在堆芯外回路中时间为τ，先驱核在堆外回路仅发生衰变，因此入口处先驱核浓度ci，in可表示为：

ci,in=ci,outexp(-λiτ)

(7)

由式(3)可得堆芯内中子产生率N为：

sin(πz/H)2πrfdrdz

(8)

从而得到归一化缓发中子源项Si，d为：

Si,d(r,z)=λici(r,z)/N

(9)

由于熔盐流速u不随半径变化，则得到熔盐在堆芯中的流动时间tc为：

tc=H/u

(10)

将式(5)、(7)和(8)代入式(9)中得到缓发中子源项分布：

Si,d(r,z)=α0,iexp(-λitcz/H)+

(α1,isin(πz/H)-α2,icos(πz/H)+

α2,iexp(-λitcz/H))J0(2.405r/R)

(11)

α0,i=βi·

(12)

(13)

(14)

其中，J1为一阶贝塞尔函数。第i组缓发中子有效份额βi，eff为：

(15)

1.2.2数值方法 熔盐堆稳态G群中子扩散方程和I组DNP守恒方程[11]如下：

g=1,2,…,G;i=1,2,…,I

(16)

(17)

其中，Dg、Σt，g和Σg′g分别为g群扩散系数、总截面和散射矩阵。基于节块展开法和有限差分方法[11]可求解得到稳态下的中子通量分布和DNP分布。

热工水力学计算采用并联多通道模型[5,11,14]，包括燃料盐的质量、动量和能量守恒方程，即：

(18)

(19)

(20)

其中：i和j分别为通道编号和轴向节块编号；M和N分别为轴向总节块数和堆芯总通道数；G、A和W分别为质量流速、等效通道流通面积和质量流量；ρ、f和h为熔盐的密度、摩阻系数和焓；g为重力加速度；D为通道水力直径；Qf和Qg分别为熔盐和石墨体积热源；Δz和Δp分别为节块高度和压降。

对式(18)和(19)耦合求解[14]可得到各通道内的流场和压力场分布。对于单个燃料通道，石墨慢化剂可简化为空心圆柱，只考虑径向导热，其能量守恒方程如下：

(21)

其中，Tg和λg分别为石墨温度和其导热系数。石墨内边界采用熔盐-石墨对流换热条件，外边界为绝热边界条件[5]。根据式(20)得到熔盐温度分布后，可解析求解式(21)得到石墨温度分布[5]。

(22)

其中，V为积分域。

基于TMSR3D的有效缓发中子份额计算流程如图2所示。具体为：1) 采用DRAGON5[15]计算各种组件的少群参数，并进行拟合；2) 物理热工耦合计算；3) 共轭通量计算；4) 根据式(22)计算βeff。

图2 基于数值方法的有效缓发中子份额计算流程Fig.2 Calculation scheme of effective delayed neutron fraction based on numerical method

2 计算结果与分析

2.1 缓发中子先驱核损失份额及其分布

额定工况下，MSRE装载235U时DNP损失份额列于表3。可以看到，本文基于解析方法和数值方法计算所得结果与其他程序[16]均符合较好。对于第2组DNP，由于其具有较大的裂变份额和较小的衰变常量，故损失份额最大。

表3 熔盐实验堆额定工况下缓发中子先驱核损失份额Table 3 Loss of delayed neutron precursor fraction of molten salt reactor experiment under rated condition

基于数值方法计算得到额定工况下DNP浓度的分布，如图3所示，由于熔盐的流动效应，每组DNP浓度的峰值均沿轴向有一定偏移，且流动对半衰期较长的DNP(如第2组)影响较大。这是因为半衰期越长，DNP在衰变前随熔盐行进的距离越长，甚至流出堆芯，造成反应性的损失。

图3 缓发中子先驱核浓度分布Fig.3 Concentration distribution of delayed neutron precursor

图4示出额定工况下最热通道内DNP与缓发中子轴向的分布。其中，DNP浓度由数值程序TMSR3D计算得到，缓发中子分布由解析方法得到。可以看到，第2组DNP归一化浓度高于其他各组DNP。这是由于DNP的稳定分布不仅与裂变份额(β)有关，也与半衰期(1/λ)有关。裂变份额越大，半衰期越长，其稳态浓度越高，即DNP浓度正比于β/λ。由于第2组DNP具有最大的β/λ值，因此其浓度高于其他各组。缓发中子由DNP经β衰变得到，如无燃料流动效应，稳态下各组缓发中子的分布次序应与DNP份额的大小顺序一致。对于液态燃料熔盐堆，还需考虑流动对缓发中子浓度的影响，如对于第5组，尽管其DNP浓度较小(图4a)，然而其具有较短的半衰期，DNP在堆内产生后短期内即衰变为缓发中子，因流动导致流出堆芯损失的缓发中子较少，因而第5组缓发中子具有较大的浓度(图4b)。

图4 额定工况下最热通道内DNP与缓发中子轴向分布Fig.4 Axial distribution of DNP and delayed neutron in the hottest channel at rated condition

2.2 入口流量对βeff的影响

液态燃料熔盐堆中燃料的流动导致引起DNP的空间再分布，部分DNP流出堆芯并衰变，从而对βeff产生影响。图5示出采用解析方法和数值方法计算得到的不同额定流量百分比下堆芯βeff的变化。可以看到，两种方法计算得到的βeff变化趋势基本一致，即入口流量增加，βeff减小。这是由于流量增加，更多的DNP流出堆芯且在堆外衰变，导致了反应性损失，从而引起了βeff的减小。除0%额定流量外，解析方法得到的计算结果始终高于数值程序计算结果，二者最大相对误差σmax为6.26%，额定流量下相对误差σnormal为4.18%。

图5 不同流量下的有效缓发中子份额Fig.5 Effective delayed neutron fraction at different percentages of flow rate

2.3 堆外流动时间对βeff的影响

对于半衰期较长的DNP，其随燃料流出堆芯后可能未发生衰变重新进入堆芯，从而使堆芯反应性及βeff产生波动，这一波动主要与燃料在堆外回路的流动时间有关[17]。图6示出采用解析方法和数值方法计算得到的堆芯βeff随燃料在堆外流动时间的变化，此过程维持入口流量为额定值。可以看到，βeff随燃料在堆外的流动时间的增加而减小，80 s后趋于稳定。这是因为随着燃料在堆外流动时间的增加，更多DNP将在堆外衰变，返回堆芯的DNP减少，从而使βeff减小。由表2可知，第1组DNP的平均寿命最大，约为80.65 s，因此当堆外流动时间超过80.65 s时，所有DNP均将在堆外衰变，因而无任何DNP返回堆芯，βeff降低至最小值，趋于稳定。可以看到，60 s前，两种方法计算结果存在一定偏差，解析方法仍高于数值方法计算结果，二者最大相对误差为10.61%，额定工况下相对误差为4.33%。结合2.2节讨论可发现，解析方法计算得到的βeff相对于数值方法偏大。这是由于，解析方法忽略了堆芯的非均匀性和反射层的作用，采用均匀的流速场，认为堆内所有缓发中子具有相同的价值，即为1，而堆外所有缓发中子价值为0。总的来说，由于燃料的流动性，相对于先驱核衰变产生缓发中子的位置，先驱核产生的位置中子价值是比较高的[3]，这也就是说，解析方法总体上高估了堆内缓发中子的价值，导致βeff计算结果高于数值方法的计算结果。

图6 有效缓发中子份额随燃料在堆外回路流动时间的变化Fig.6 Effective delayed neutron fraction change with fuel residence time out of core

表2 缓发中子先驱核份额及其衰变常量Table 2 Delayed neutron precursor fraction and precursor decay constant

3 结论

有效缓发中子份额βeff的精确计算可为反应堆设计和安全分析提供合理准则和依据。燃料的流动性、特殊的传热方式和钍铀燃料循环的使用导致液态燃料熔盐堆βeff的计算方法与固态燃料反应堆不同。本文针对石墨慢化通道式熔盐堆，分别基于解析方法和数值方法提出了计算βeff的数学模型，计算了MSRE在额定工况下的DNP损失份额和堆内DNP浓度分布，并分析了燃料在堆外流动时间和入口流量对βeff的影响，可得到以下结论。

1) 两种方法均可对DNP行为提供合理描述。解析方法忽略了堆芯的非均匀性和反射层的作用，总体上高估了堆内缓发中子的价值，导致βeff计算结果相对于数值方法偏大。

2) 固定燃料在堆外流动时间，βeff随入口流量的增加而减小。

3) 固定入口流量，βeff随燃料在堆外流动时间的增加而减小，80 s后趋于稳定。