【期刊名称】《原子能科学技术》
【年(卷),期】2016(050)007
【总页数】6页(P1245-1250)
【关键词】改进准静态本征函数法;时空中子动力学;因子分解法
【作 者】李明芮;黎浩峰;陈文振;邢晋
新视野大学英语mp3【作者单位】海军工程大学 核能科学与工程系，湖北 武汉 430033;海军工程大学 核能科学与工程系，湖北 武汉 430033; 海军司令部 核安全部，北京 100841;海军工程大学 核能科学与工程系，湖北 武汉 430033;海军工程大学 核能科学与工程系，湖北 武汉 430033; 海军装备部驻重庆地区军事代表局，重庆 400000
【正文语种】中 文
2016年11月19日
【中图分类】TL326
反应堆中子时空动力学方程的求解是一复杂而困难的任务，变量多、维数大、时间变量与空间变量互相耦合，且微分方程具有强烈的刚性。人们无法导出反应堆时空中子动力学方程的精确解析解，但由于解析解对反应堆运行安全具有重要的理论意义和实用价值，因此人们一直在研究其近似解。其基本思路是将空间、能量和时间变量先进行离散，再进行求解[1]。反应堆时空中子动力学方程的解法有数值解法、模项展开法和因子分解法[1-2]。但前两种解法计算量大，无法满足反应堆中子通量快速计算的要求[3-5]。因此，一些学者采用因子分解法求解时空中子动力学方程[6-8]，提出了准静态近似条件下和改进准静态近似条件下不同的计算方法，得到了各种近似条件下的时空中子动力学方程的计算结果[9-11]。在最为复杂和耗时的形状函数计算中，普遍采用归一化方法，为了提高计算精度，就必须进行多次迭代求解，这使得计算复杂化，且增加了计算时间。因此本文研究改进准静态本征函数法，避免使用归一化条件，从而提高计算速度[12]。
不考虑外加中子源时，仅考虑1组缓发中子，多维单速时空中子动力学方程组如下：
其中：N(r,t)和C(r,t)分别为均匀裸堆内r处t时刻的中子通量密度和先驱核浓度；D为扩散常
数；v为中子速率；Σa为堆芯介质的宏观吸收截面；β为缓发中子份额；k∞为无穷增殖因数；λ为先驱核衰变常量。式(1)右端第1项为t时刻每秒因扩散进入r附近单位体积内的中子数；第2项为r处t时刻每秒每单位体积内被介质吸收的中子数；其余两项为考虑缓发中子效应的源项，即t时刻r处每秒单位体积产生的中子数。
改进准静态本征函数法是把中子通量密度函数分解为只与时间有关的幅函数n(t)和形状函数φ(r,t)的乘积，在求解形状函数时，考虑φ(r,t)随时间的变化关系，然后利用本征函数法求解中子通量密度函数的形状函数，而在处理幅函数过程中，将幅函数化简为点堆模型下的形式。
由因子分解法分别将中子通量密度函数和先驱核浓度分解为幅函数和形状函数的乘积：
其中：n(t)、c(t)分别为中子通量密度和先驱核浓度的幅函数；φ(r,t)、g(r,t)分别为中子通量密度和先驱核浓度的形状函数。
将式(3)、(4)代入式(1)，两边同时除以N(r,t)，得：
由文献[1]可假定先驱核浓度的空间分布形状与中子通量密度的空间分布形状相似。根据假
设，这里令：
故式(5)可简化为：
从式(7)可看出，方程左端仅为时间t的函数，方程右端第3、4项为常数项，根据文献[13]得到的形状函数解的形式，是以自然对数为底关于时间变化的幂函数和关于位置变化的余弦式的乘积，可得出方程右端第1、2项只能为常数项，设为-a，则有：
在处理三维形状函数问题时，通常采用柱坐标求解，由式(8)可得到：
利用分离变量法，令：
将式(10)代入式(9)得：
qcc是什么两边同时除以f(r,t)Z(z,t)，整理得：
观察式(12)，r与z是两个彼此无关的独立变量，要使上式对任意的r与z均成立，只有左端为不能含有r的表达式，右端为不能含有z的表达式，即左右两端同时等于常数或关于t的函数。
从文献[13]中平板裸堆的中子通量密度形状函数的结果形式可得出式(12)右端为常数b，即两端同时为b。
由上述分析得出圆柱型堆芯中子通量密度形状函数轴向分布的方程为：
假定在不考虑反射层时，堆芯的上下两端中子通量密度为固定值，即边界条件：
其中：H为堆芯高度；ψ1、ψ2分别为堆芯两端的中子通量密度。
eigendecomposition采用适当的初始条件：倒签提单
观察式(13)的形式、初始条件和边界条件，可采用本征函数法求解式(13)。
先求解对应边界条件的稳态解，此时，形状函数变化与时间无关，代入边界条件可得到稳态解为：
因此非齐次边界问题可化为齐次边界问题：
求对应的齐次方程的本征值问题为：
采用分离变量法，令Z(z,t)=Z(z)T(t)，代入式(18)得：
令式(19)两边同时等于常数，由边界条件可知该常数为负值，设为p，代入边界条件，求得本征函数和本征值为：
利用常数变异法，令：
将式(21)代入式(13)，比较两边的傅里叶系数相等，有：
当a-b=Dv(2n-1)2π2/H2时，可求出T为常数，不符合改进准静态近似要求；当a-b≠Dv(2n-1)2π2/H2时，解得：
其中：
齐次边界条件下，齐次方程的解为：
故采用本征函数法求得的三维中子通量密度的形状函数为：
当反应堆处于稳态时，n取值为1，此时中子通量密度形状函数的变化与时间无关，故可以求得a-b=Dvπ2/H2。
由式(12)可得出圆柱型堆芯中子通量密度径向形状函数方程为：
式(27)的边界条件为：在反应堆堆芯半径R处，中子通量密度为零；中子通量密度有界、连续。
采用适当的初始条件：
采用分离变量法求解式(27)，令：
英语网址livejasmin将式(29)代入式(27)，并令两边同时为常数-μ2，化简得到：
由于R在有限区间上应该有解，因此分离变量后的常数取负值。故式(30)可分解为：
r2R″+rR′+μ2r2R=0
式(31)为零阶贝塞尔方程，其一般解为：
其中，J0、Y0分别为零阶第1类及第2类贝塞尔函数，由边界条件可得到C2必须等于零。这样可得到关于r的本征函数为：
由边界条件求得本征值为：
其中，xn满足零阶第1类贝塞尔函数值为零，即J0(xn)=0。
求解式(32)的结果为：
两个常数的乘积仍为常数，可令两个常数的乘积为Bn，故可以得到径向形状函数的解为：
由贝塞尔展开定理可得出：
当反应堆处于稳态时，n=1又x1=2.405，此时中子通量密度形状函数的变化与时间无关，故可求得b=Dv2.4052/R2。
由a-b=Dvπ2/H2、b=Dv2.4052/R2可求得a=Dvπ2/H2+Dv2.4052/R2=DvB2，将求得的a代入式(7)可得到：
同理，将式(3)、(4)代入式(2)化简得到：
chant是什么意思根据各变量关系，可化简得：
故可得到幅函数方程组：
其中：ρ为反应堆反应性；β为缓发中子总份额；l为瞬发中子平均寿命。
由文献[14-15]推导出中子通量密度的幅函数的近似解为：
其中：n0为初始中子通量密度；ρs为停堆深度，ρs=λl+β-ρ。
将式(26)与(45)、(37)与(45)分别代入式(3)，可得到不考虑外加中子源时，单群时空中子动力学方程轴向和径向的近似解析解分别为：
其中：
migros其中，Z0(z)和f0(r)分别为轴向和径向形状函数的初始条件。
本文通过因子分解法将中子通量密度分解为幅函数和形状函数，利用改进准静态本征函数法求解了改进准静态近似条件下的形状函数，然后将幅函数化简为点堆模型下的形式，最终求解了三维时空中子动力学方程轴向和径向的近似解析解。本文研究的方法克服了文献[13]中在推导平板堆芯形状函数时，不考虑缓发中子份额的局限性，而且克服了求解形状
函数径向随时间变化的困难，因此该方法应用范围更广。另外，式(46)、(47)是显函数表达式，故计算量小，能满足实时甚至超时计算的要求。但本文研究的方法对于临界附近的计算适用性好，对于反应堆物理启动的计算误差较大。在后续研究中如果能找到一种关系式将轴向和径向形状函数结合，而不是简单地相乘，那么就能得到改进准静态下时空中子动力学的解析解，或讨论多群时空中子动力学方程，形状函数的变化与反应性的变化会建立某种关系，这些研究将有着非常重要的意义。
【相关文献】
[1] 蔡章生. 核动力反应堆中子动力学[M]. 北京：国防工业出版社，2005.
[2] 黄祖洽. 核反应堆动力学基础[M]. 北京：北京大学出版社，2007.
[3] ABOANBER A E, NAHLA A A, AL-MUHIAMEED Z I. A novel mathematical model for two-energy groups of the point kinetics reactor dynamics[J]. Progress in Nuclear Energy, 2014, 77: 160-166.
