代中子时间法求解点堆中子动力学方程

2012-06-26蔡光明阮良成

核科学与工程 2012年4期

蔡光明，阮良成

（中国核电福建福清核电有限公司，福建福清 350318）

点堆中子动力学方程描述反应堆中子密度和反应性之间随时间的变化。求解中子密度随时间的变化，对于反应堆的事故分析、仿真、运行及物理试验都有重要的作用。

几十年来，点堆中子动力学方程求解一直是个研究热点。由于点堆中子动力学方程是个刚性方程[1－3］，即方程中的几个时间常数之间相差几个数量级，因此准确、快速、稳定地求解方程是困难的。得益于现代计算机技术的进步，本文直接采用代中子时间计算法求解点堆中子动力学方程，经过基准例题和动态—逆动态对比计算，验证了模型、程序计算的准确性和稳定性，而计算时间也是可接受的。

1 代中子时间计算法

当反应堆处于动态时，为了简化模型，假设堆内中子注量率各处变化与空间位置无关，即可得到假定“点堆，单群，6组缓发中子”的点堆中子动力学方程：

其中：n（t）为与时间相关的中子密度；βi为第i组缓发中子份额；I为缓发中子平均价值因子；βieff为第i组有效缓发中子份额；βeff为有效缓发中子份额；λi为第i组缓发中子先驱核的衰变常数；ci（t）为第i组先驱核密度；l0为瞬发中子平均寿命；l为瞬发中子平均代时间；S为中子源；keff是反应堆有效增殖系数。

由于计算机技术的发展，本计算模型的中心思想是直接采用瞬发中子平均代时间作为时间增量来计算中子密度的变化，以计算时间为代价来解决方程的刚性问题，来保证点堆中子动力学方程解的准确性和稳定性。因此，令dt＝l，则方程（1）、（2）分别转化为：

则下一代的中子密度、第i组先驱核密度分别为：

其中，反应性ρ定义为

2 程序与验证

根据第1节所述的模型，编制了计算机程序 KNGT（Kinetic Code of Neutron Generation Time method），并用参考文献中的基准例题进行了计算比对。以下是关于程序的介绍与例题比对结果。

2.1 程序编制说明

代中子时间计算法求解点堆中子动力学方程计算机程序KNGT，采用C＋＋语言设计，编译的目标代码执行效率高，可以充分发挥计算机的速度。另外所有的浮点变量全部定义为浮点双精度，尽量避免赋值误差和计算产生的截断误差[3］。

2.2 基准例题验证

例1，计算热堆输入阶跃正反应性后，堆内中子密度n（t）随时间t的变化。其中子动力学参数为：βieff＝0.000 266、0.001 491、0.001 316、0.002 849、0.000 896、0.000 182，λi＝0.012 7、0.031 7、0.115、0.311、1.4、3.87s－1，i＝1～6，l＝2.0E－5s；其他参数为：ρ0＝0，ρ＝0.003，S＝0s－1。计算结果列于表1。

表1 阶跃正反应性输入后的热堆中子密度变化Table1 Neutron density variation in thermal reactor after the input of positive step reactivity

续表

例2，计算热堆输入大阶跃正反应性ρ＝0.007后，堆内中子密度n（t）随时间t的变化，其他参数同例1。计算结果列于表2。

例3，计算热堆输入阶跃负反应性ρ＝－0.007后，堆内中子密度n（t）随时间t的变化，其他参数同例1。计算结果列于表3。

表2 大阶跃正反应性输入后的热堆中子密度变化Table2 Neutron density variation in thermal reactor after the input of great positive step reactivity

表3 阶跃负反应性输入后的热堆中子密度变化Table3 Neutron density variation in thermal reactor after the input of negative step reactivity

例4，计算热堆输入线性反应性后，堆内中子密度n（t）随时间t的变化，计算结果列于表4。线性反应性变化率为0.000 7s－1，其他参数同例1。

例5，计算堆内引入正弦反应性后，堆内中于密度n（t）随时间t的变化。ρ＝ρ0sin（πt／T0），其中，ρ0和T0为系数。计算结果列于表5。βeff＝0.007 9，λ＝0.077s－1，l＝10－8s。

表4 线性反应性输入后的热堆中子密度变化Table4 Neutron density variation in thermal reactor after the input of liner reactivity

表5 复杂反应性输入后的反应堆中子密度峰值Table5 Neutron density peak in thermal reactor after the input of complex reactivity

从以上例1～3几种计算方法的比较，本文方法的计算结果与精确解（埃尔米特法，h＝0.000 1s[4］）比其他方法更吻合；且在计算时间上也并无劣势。例4中本文方法的计算结果也比其他方法更吻合。例5反映了本文方法对复杂反应性计算的稳定性。

3 其他验证

上一节主要用基准例题作了计算比较，也是各位学者一般采用的方法。本节将用动态－逆动态的计算结果互相比较，来做一步验证。

3.1 动态－逆动态计算比较（反应性—中子密度—反应性）

例6，表6是模拟动态刻棒法DRWM的反应性ρdrwm变化，所用的动力学参数同例1，ρ0＝－0.01×10－5，S＝5.0E－10s－1。n是随时间进行的点堆动力学动态计算结果；根据n用ODRM[5］程序进行逆动态计算得到反应性ρodrm，n的计算时间间隔为0.1s。

ρodrm－ρdrwm是反应性通过动态计算得到中子密度再通过逆动态计算得到的反应性的比较，这个偏差最大值为2.4×10－5，相对偏差小于0.3%；±50×10－5范围内，偏差最大值为0.5×10－5，相对偏差最大1%，反应性吻合很好。由此分析可以得出KNGT的计算是正确的，准确的；同时另一方面也验证了逆动态程序ODRM模型与程序的正确性。