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

蔡光明，阮良成

（中国核电 福建福清核电有限公司，福建 福清 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模型与程序的正确性。

表6 动态－逆动态计算反应性比较Table6 Comparison between kinetic－inverse kinetic contrast calculation on reactivity

续表

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

例7，图1中I是某反应堆稀释刻棒的中子密度变化实测曲线，ρ是用ODRM程序进行逆动态计算得到反应性，n是用KNGT进行的点堆动力学动态计算得到的。I、ρ、n计算的时间间隔为1s。

n－I是中子密度通过逆动态计算再通过动态计算的比较，这个偏差最大值为0.082，最大相对偏差小于0.11%；中子密度吻合很好，图1中子密度n、I曲线重叠无法区分。由此分析可以再次得出KNGT的计算是正确的、准确的；同时另一方面也再次验证了逆动态程序ODRM模型与程序的正确性。

图1 刻棒中的中子密度与反应性的变化Fig.1 Variation of neutron density and reactivity in rod measurement test

4 小结

对于用代中子时间法编制的程序求解点堆中子动力学方程，经过例题验证及动态－逆动态的验证，有如下结论：

（1）经过基准例题比较计算验证，模型是正确的，程序KNGT是正确的，计算结果是准确的。

（2）由于模型使用了最小的时间参数，未采用特别的假设，因此模型的计算结果是稳定的，可用于各种复杂动态计算。

（3）在当前计算机硬件条件下，程序运行是足够快的，可用于热堆的实时仿真计算。

（4）使用动态－逆动态相互计算比较，更进一步验证了动态程序KNGT的计算正确性与稳定性；另一方面也验证了逆动态程序ODRM计算的正确性。

需要另外指出的是，例6中中子源S≠0，这在参考文献中鲜有类似的基准例题验证。除了本文的方法验证之外，希望得到其他学者的验证。

[1]廖忠岳，李铁军，刘长欣 .反应堆中子动力学方程仿真算法及其应用[J］.核动力工程，1995，16（4）：306－311.

[2]黎浩峰，等 .用高阶泰勒多项式积分方法求解点堆中子动力学方程[J］.原子能科学技术，2008，42：162－168.

[3]经荥清，等 .点动态学方程几种解法的比较[J］.清华大学学报（自然科学版），1995，35（3）：7－12.

[4]陈东，刘振华 .求解点堆中子动力学方程的三角插值法[J］.原子能科学技术，2010，44（10）：1195－1120.

[5]蔡光明 .落棒法测量所有控制棒全插时的棒价值[J］.核科学与工程，2005，25（2）：154－158.