Analytical Solution Study on Transient Neutron Transport in Homogeneous Bare Reactor Based on Neutron Telegraph Equation
-
摘要: 斐克定律作为扩散理论的基础，其本质是在中子输运方程的P1近似方程（电报方程）的基础上进一步忽略了中子流密度关于时间的导数项，因此在瞬态过程中无法准确描述实际的中子动力学过程。本文基于单能中子电报动力学方程，针对一维无限大平板裸堆中子动力学问题，利用分离变量法推导了其解析解，并与扩散中子动力学方程的解析解进行了对比分析。研究发现：在瞬态变化过程中，电报方程的解与扩散方程的解相比，空间项仍保持余弦函数的形式，但是时间项的变化则更为复杂，一是时间项的级数组合形式受到问题的几何和材料的影响，二是高阶谐波表现出随时间振荡变化的规律。该研究结果可以为后续基于中子电报方程的数值理论研究提供参照和依据。Abstract: The Fick's law, as the foundation of the diffusion theory, further neglects the time derivative term of the neutron current density in the P1 approximation equation (telegraph equation) of the neutron transport equation. Therefore, it is difficult to accurately describe the actual neutron kinetic behavior in the transient process. In this paper, based on the monoenergetic neutron telegraph equation, an analytical solution is derived for the one-dimensional infinite slab bare reactor neutron kinetics problem using the method of separation of variables, and it is compared and analyzed with the analytical solution of the neutron diffusion equation. The study reveals that during transient changes, the spatial term of the telegraph equation still maintains the form of a cosine function compared to the diffusion equation's solution, but the variation of the temporal term is more complex. Firstly, the combination form of the temporal term's orders is influenced by the geometry and materials of the problem. Secondly, higher-order harmonics exhibit oscillatory changes with time. These research findings can provide references and foundations for subsequent numerical theoretical studies based on the neutron telegraph equation.
-
Key words:
- Neutron kinetics /
- P1 approximation /
- Telegraph equation /
- Diffusion equation /
- Analytical solution
-
0. 引 言
为了保障核反应堆的安全稳定运行，需要对堆芯各种工况下的状态进行准确预测和模拟，尤其是瞬态安全分析。对反应堆中子动力学过程的准确预测，不仅能为堆芯核设计提供更加可靠的理论依据，而且能为实际反应堆的安全运行提供参照和指导，预防堆芯瞬态事故的发生。
球谐函数方法是中子输运理论中发展最早的一种计算方法，其中最简单且应用最广泛的是P1近似及其产生的扩散近似。传统的瞬态计算理论模型主要是基于扩散方程的点堆动力学模型或者组件均匀化的节块扩散中子动力学计算方法。在单速且中子源是各向同性的稳态问题中，P1近似和扩散近似是等价的[1]。因此，经常将P1近似直接称为扩散近似。然而事实上对于中子动力学问题，扩散近似需要在P1近似的基础上进一步忽略中子流密度关于时间的导数项，从而得到斐克定律；对于考虑多能群中子的情形还需额外引入输运修正近似[2]。这种近似在理论上势必会影响堆芯瞬态计算的可靠性。
在外文文献中，P1近似方程又被称之为中子电报方程。电报方程这一概念被应用在各个领域[3]，但在核反应堆物理分析中却鲜有直接应用。相关文献[4-5]均指出，扩散方程由于在电报方程的基础上忽略了中子流密度关于时间的导数项，使得抛物型的扩散方程所表征的中子传播速度为无限大，尽管扩散方程中显含的速度v（v 为实际的中子速率）为有限值，这一点在脉冲中子实验中有所印证。这显然与实际不符，因为已知典型的热中子速度为2200 m/s (0.0253 eV)，为有限值。而电报方程为双曲型偏微分方程，在数学上表征了有限的中子传播速度，尽管由此给出的中子传播速度的值是有误的($v/\sqrt 3 $)，但是其相较于扩散方程中的无限中子传播速度，已经有了很好的改进。因此，有必要基于直接的P1方程（电报方程）研究和分析中子动力学过程，并对比分析扩散中子动力学方程的区别。
近年来，Altahhan等人开始将电报方程模型应用在核反应堆动力学分析中，以克服扩散方程的一些缺陷。基于中子电报方程的点堆动力学模型（TNRK）被建立，并与扩散点堆动力学模型（DPRK）关于反应性引入实验做了较为充分的比较[5-8]。基于Gen-Foam程序，时空多群中子电报方程也被应用于熔盐堆系统的动力学分析中[9]。在上述电报方程点堆动力学模型的基础上，考虑两群中子的点堆模型求解[10]以及伴随敏感性分析也随之被提出[11]。总之，针对点堆问题，电报方程与扩散方程的对比分析已经比较成熟，在电报方程的点堆中子动力学分析过程中定义了中子弛豫时间这一新的物理量，并认为其为有限值，其值受中子传输介质的影响；此时基于扩散方程的点堆中子动力学即为电报方程在中子弛豫时间等于0时的特殊情形。相比之下，针对同时考虑空间变量的中子动力学问题，尚未发现相应的理论分析，其中的数学物理机制尚不明确。
因此，本文针对一维均匀平板裸堆动力学问题，推导了单能中子电报方程的解析解，并与已知扩散方程给出的解析解进行了对比分析，旨在从理论上揭示电报方程和扩散方程在时空中子动力学过程中刻画的中子通量密度与中子流密度的变化机理。
1. 电报方程
针对单能时空中子输运方程，应用P1近似后，可得电报方程：
$$ \begin{split} &\frac{1}{v}\frac{{\partial \phi ({\boldsymbol{r}},t)}}{{\partial t}} + \nabla \cdot {\boldsymbol{J}}({\boldsymbol{r}},t) + {\varSigma _{\mathrm{a}}}({\boldsymbol{r}},t)\phi ({\boldsymbol{r}},t) \\ & =\nu {\varSigma _{\mathrm{f}}}({\boldsymbol{r}},t)\phi ({\boldsymbol{r}},t) \end{split}$$ (1) $$ \frac{1}{v}\frac{{\partial {\boldsymbol{J}}({\boldsymbol{r}},t)}}{{\partial t}} + \frac{1}{3}\nabla \phi ({\boldsymbol{r}},t) + {\varSigma _{{\mathrm{tr}}}}({\boldsymbol{r}},t){\boldsymbol{J}}({\boldsymbol{r}},t) = 0 $$ (2) 式中，v为中子飞行速率，cm/s；$\phi ({\boldsymbol{r}},t)$和$ \boldsymbol{J}(\boldsymbol{r},t) $分别为在t时刻空间位置r处的中子（标）通量密度和中子流密度，cm−2 · s−1；${\varSigma _{\mathrm{a}}}({\boldsymbol{r}},t)$为中子宏观吸收截面，cm−1；$ \nu \varSigma_{{\mathrm{f}}}({\boldsymbol{r}}, t) $为裂变中子宏观产生截面，cm−1；${\varSigma _{{\mathrm{tr}}}}({\boldsymbol{r}},t)$称为宏观输运截面，cm−1。
2. 扩散方程
若忽略式（2）中的中子流密度关于时间的导数项，即假设：
$$ \frac{1}{{v{\varSigma _{{\mathrm{tr}}}}({\boldsymbol{r}},t)}}\frac{{\partial {\boldsymbol{J}}({\boldsymbol{r}},t)}}{{\partial t}} = 0 $$ (3) 令扩散系数$D({\boldsymbol{r}}, t) $为：
$$ D({\boldsymbol{r}}, t)=\frac{1}{3 \varSigma_{{\mathrm{t r}}}({\boldsymbol{r}}, t)} $$ (4) 得到斐克定律：
$$ {\boldsymbol{J}}({\boldsymbol{r}},t) = - D({\boldsymbol{r}},t)\nabla \phi ({\boldsymbol{r}},t) $$ (5) 将式（5）代入式（1），即可得到扩散方程：
$$ \frac{1}{v}\frac{{\partial \phi ({\boldsymbol{r}},t)}}{{\partial t}} = \nabla D\nabla \phi ({\boldsymbol{r}},t) + (\nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})\phi ({\boldsymbol{r}},t) $$ (6) 此处省略了截面参数中的空间变量和时间变量。
3. 平板裸堆问题的扩散方程解析解
针对厚度为h（包含外推距离，边界中子通量密度为0）的平板均匀裸堆问题（几何描述见图1），在截面参数不随时间变化和不考虑缓发中子的情况下，文献[12]给出了单能中子扩散方程[式（6）]的解析解：
$$ \phi(x,t)=\sum\limits_{n=1}^{\infty}A_n\mathrm{e}^{(k_n-1)t/l_n}\cos (B_nx),n=1,2,\cdots $$ (7) $$ \left\{ \begin{gathered} {l_n} = \frac{{{L^2}}}{{Dv(1 + {L^2}B_n^2)}} \\ {L^2} = \frac{D}{{{\varSigma _{\mathrm{a}}}}} \\ {B_n} = \frac{{(2n - 1)\pi }}{h} \\ {k_n} = \frac{{{k_\infty }}}{{1 + {L^2}B_n^2}} \\ \end{gathered} \right. $$ 式中，${k_\infty }$为无限介质增殖系数；An为待定系数，可由初始条件确定。
4. 平板裸堆问题的电报方程解析解
针对该平板均匀裸堆问题，单能中子电报方程为：
$$ \frac{1}{v}\frac{{\partial \phi (x,t)}}{{\partial t}} = (\nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})\phi (x,t) - \frac{{\partial {\boldsymbol{J}}(x,t)}}{{\partial t}} $$ (8) $$ - \frac{1}{v}\frac{{\partial {\boldsymbol{J}}(x,t)}}{{\partial t}} = {\varSigma _{{\mathrm{tr}}}}({\boldsymbol{r}},t){\boldsymbol{J}}(x,t) + \frac{1}{3}\frac{{\partial \phi (x,t)}}{{\partial t}} $$ (9) 对应的初始条件为：
$$ \phi (x,0) = {\varphi _0}(x),{\text{ }}{\boldsymbol{J}}(x,t) = {{\boldsymbol{J}}_0}(x) $$ (10) 边界条件为：
$$ \phi \left( - \frac{h}{2},t\right) = 0,{\text{ }}\phi \left(\frac{h}{2},t\right) = 0 $$ (11) 对式（8）两边关于时间变量t求偏导，对式（9）两边关于空间变量x求偏导（忽略时间和空间之间的依赖性），并结合原式（8）可以消去J(x,t)项，得到关于中子通量密度的二阶偏微分方程：
$$ \frac{\tau}{Dv}\frac{\partial^2\phi}{\partial t^2}+\left(\frac{1}{Dv}-\tau B_{\mathrm{m}}^2\right)\frac{\partial\phi}{\partial t}=\frac{\partial^2\phi}{\partial x^2}+B_{\mathrm{m}}^2\phi $$ (12) $$ \tau = \frac{{3D}}{v} = \frac{1}{{v{\varSigma _{{\mathrm{tr}}}}}} $$ $$ B_\mathrm{m}^2 = \frac{{\nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}}}}{D} $$ 式中，$ \tau $为中子弛豫时间，s：$ B_{\mathrm{m}}^2 $为核反应堆的材料曲率，cm−2。
该方程也是中子电报方程的一种形式[5]，可以看出，当τ=0时，式（12）会退化成中子扩散方程。
接下来利用分离变量法对式（12）进行解析求解，设中子通量密度函数可分离为形状函数X(x)和幅函数T(t)的乘积：
$$ \phi (x,t) = X (x)T(t) $$ (13) 并将之代入式（12）可得：
$$\begin{split} &\frac{\tau}{Dv}\frac{1}{T(t)}\frac{\mathrm{d}^2T}{\mathrm{d}t^2}+\left(\frac{1}{Dv}-\tau B_{\mathrm{m}}^2\right)\frac{1}{T(t)}\frac{\mathrm{d}T}{\mathrm{d}t}-B_{\mathrm{m}}^2 \\ &=\frac{1}{X(x)}\frac{\mathrm{d}^2X}{\mathrm{d}x^2}=-B^2 \end{split} $$ (14) 式中，B2为任意常数。
根据空间函数与时间函数的无关性，可得关于空间变量和时间变量的方程分别如下：
$$ X''(x) + {B^2}X(x) = 0 $$ (15) $$ T''(t)+\left(\frac{1}{\tau}-DvB_{\mathrm{m}}^2\right)T'(t)+\frac{Dv}{\tau}(B^2-B_{\mathrm{m}}^2)T(t)=0 $$ (16) 由式（10）可得式（15）的边界条件：
$$ X( - h/2) = X(h/2) = 0 $$ (17) 式（15）及其边界条件构成一个特征值问题，当且仅当：
$$ {B_n} = \frac{{(2n - 1)\pi }}{h},{\text{ }}n = 1,2,3 \cdots $$ 该特征值问题才有特征解：
$$ X_{n}(x)=\cos \left(B_{n} x\right) $$ (18) 所以，式（12）的通解可表示为：
$$ \phi (x,t) = \sum\limits_{n = 1}^\infty {{T_n}(t)\cos ({B_n}x)} $$ (19) 式（16）为二阶常系数齐次线性微分方程，其对应的特征方程为：
$$ q^2+\left(\frac{1}{\tau}-DvB_{\mathrm{m}}^2\right)q+\frac{Dv}{\tau}(B^2-B_{\mathrm{m}}^2)=0 $$ (20) 其特征根为：
$$ \begin{split} q_1^{(n)} =\;& - \frac{1}{{2\tau }}(1 - \tau DvB_{\rm{m}}^2) + \\ \;&\frac{1}{{2\tau }}(1 + \tau DvB_{\rm{m}}^2){\left[1 - \frac{{4\tau DvB_n^2}}{{(1 + \tau DvB_{\rm{m}}^2)}}\right]^{1/2}} \end{split} $$ (21) $$ \begin{split} q_2^{(n)} =\;& - \frac{1}{{2\tau }}(1 - \tau DvB_{\rm{m}}^2) - \\ \;&\frac{1}{{2\tau }}(1 + \tau DvB_{\rm{m}}^2){\left[1 - \frac{{4\tau DvB_n^2}}{{(1 + \tau DvB_{\rm{m}}^2)}}\right]^{1/2}} \end{split} $$ (22) 根据特征根的取值情况，可以确定原方程通解的形式。
（1）当：
$$ 1 - \frac{{4\tau DvB_n^2}}{{(1 + \tau DvB_{\rm{m}}^2)}} = 1 - \frac{{4B_n^2}}{{3({\varSigma _{{\mathrm{tr}}}} + \nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})}} > 0 $$ 即$ {B_n} < \sqrt 3 ({\varSigma _{{\mathrm{tr}}}} + \nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})/2 $时，式（16）的解可表示为：
$$ T_n(t)=A_n\mathrm{e}^{q_1^{(n)}t}+B_n\mathrm{e}^{q_2^{(n)}t} $$ (23) （2）当$ {B_n} = \sqrt 3 ({\varSigma _{{\mathrm{tr}}}} + \nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})/2 $时，式（16）的解可表示为：
$$ {T_n}(t) = ({A_n}^\prime + {B_n}^\prime t){{\mathrm{e}}^{\alpha t}} $$ (24) $$ \alpha = - \frac{1}{{2\tau }}(1 - \tau DvB_{\rm{m}}^2) = \frac{v}{2}(\nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}} - {\varSigma _{{\mathrm{tr}}}}) $$ （3）当$ {B_n} > \sqrt 3 ({\varSigma _{{\mathrm{tr}}}} + \nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})/2 $时，式（16）的解可表示为：
$$ T_n(t)=(A_n^{''}\cos\beta_nt+B_n^{''}\sin\beta_nt)\mathrm{e}^{\alpha t} $$ (25) $$ \begin{gathered} {\beta _n} = \frac{1}{{2\tau }}(1 + \tau DvB_{\rm{m}}^2){\left[1 - \frac{{4\tau DvB_n^2}}{{(1 + \tau DvB_{\rm{m}}^2)}}\right]^{1/2}} \\ = \frac{v}{2}\left[\frac{4}{3}B_n^2 - {({\varSigma _{{\mathrm{tr}}}} + \nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})^2}\right] \\ \end{gathered} $$ 由此存在常数N，使得方程的通解具有如下形式：
$$ \begin{split} & \phi (x,t) = \sum\limits_{n = 1}^N {({A_n}{{\mathrm{e}}^{q_1^{(n)}t}} + {B_n}{{\mathrm{e}}^{q_2^{(n)}t}})\cos ({B_n}x)} + \\ &\qquad \sum\limits_{n = N + 1}^\infty {({A_n}^{\prime \prime }\cos {\beta _n}t + {B_n}^{\prime \prime }\sin {\beta _n}t){{\mathrm{e}}^{\alpha t}}\cos ({B_n}x)} \end{split} $$ (26) 需要指出的是，由于式（24）所示的情况较为特殊，仅在该式两边恰好严格相等时才出现，所以通解式（26）中并未包含式（24）这种特殊形式。另外，可以看到，N的大小取决于Bn与$ \sqrt 3 ({\varSigma _{{\mathrm{tr}}}} + \nu {\varSigma _{\mathrm{f}}} - {\varSigma _{\mathrm{a}}})/2 $的相对大小，前者只与问题的几何尺寸相关，而后者只与问题的材料（截面参数）相关。
5. 解析解的对比分析
对比扩散方程解析解[式（7）]和电报方程解析解[式（26）]，可以得到以下结论：
（1）在空间分布的函数组分上，二者是一致的，都是余弦函数系。
（2）在时间项上，电报方程解析解的构成更为复杂。一方面，时变函数的形式会在某个临界阶数N前后出现从指数函数到三角函数与指数函数组合的变化；另一方面，各阶谐波指数衰减的速率也与扩散方程的解析解不同。
（3）易知$q_1^{(1)} > \alpha > q_2^{(1)}$，若$q_1^{(1)} = 0$，则$\phi (x,t)$会随着时间的增加趋于稳定，即此时的核反应堆处于临界状态。将$q_1^{(1)} = 0$代入式（21）可得核反应堆的临界条件[式（27）]，等式右边即为核反应堆的几何曲率，该结论与扩散模型下的单群理论的临界方程相同[12]。
$$ B_\mathrm{m}^2 = B_1^2 = {\left(\frac{\pi }{h}\right)^2} $$ (27) （4）在阶跃反应性问题中，基波和各阶谐波的幅值随时间的变化规律不同。在扩散解析解中，每一阶波函数的幅值都是按照单个指数函数的规律发生变化；且阶数越高，变化越快；变化最慢的基波在临界情况下可以形成稳定分布。相比之下，在电报解析解中，n>N的高阶谐波幅值整体也是按照指数函数进行衰减的，且不同阶的谐波衰减指数相同，但每一阶谐波还会在指数衰减范围内按照两个三角函数的组合进行振荡变化，且谐波阶数越高，振荡频率越大；低阶波（n<N）的幅值按照两个指数函数的组合发生变化，临界条件下，基波中的两个指数函数中，有一个会成为时间上的常数，另一个会随着时间逐渐衰减，使基波形成稳定分布的形状。
6. 结 论
本文基于单能中子电报方程，即中子输运方程的P1近似方程，针对一维无限大均匀平板裸堆中子动力学问题，推导了电报方程的解析解，并与扩散方程给出的解析解作了定性的对比分析。结果表明，电报方程由于在P1方程中保留了中子流密度关于时间的导数项，即消去了作为扩散理论基础的斐克定律的假设，其中子动力学行为表现得更为复杂，在解析解中不仅时间项的构成受到问题的材料和几何的影响，而且在其高阶谐波中还表现出随时间振荡变化的规律。
未来将基于本文研究结果进一步开展关于电报方程数值计算模型的研究，并探究其在实际反应堆堆芯瞬态问题中相较于扩散方程结果所带来的改进。
-
[1] 谢仲生,邓力. 中子输运理论数值计算方法[M]. 西安: 西北工业大学出版社,2005: 71-72. [2] 吴宏春,郑友琦,曹良志,等. 中子输运方程确定论数值方法[M]. 北京: 中国原子能出版社,2018: 100-102. [3] HEIZLER S I. Asymptotic telegrapher’s equation ( P1) approximation for the transport equation[J]. Nuclear Science and Engineering, 2010, 166(1): 17-35. doi: 10.13182/NSE09-77 [4] WEISS G H. Some applications of persistent random walks and the telegrapher's equation[J]. Physica A:Statistical Mechanics and its Applications, 2002, 311(3-4): 381-410. doi: 10.1016/S0378-4371(02)00805-1 [5] ALTAHHAN M R, NAGY M S, ABOU-GABAL H H, et al. Formulation of a point reactor kinetics model based on the neutron telegraph equation[J]. Annals of Nuclear Energy, 2016, 91: 176-188. doi: 10.1016/j.anucene.2016.01.011 [6] ALTAHHAN M R, ABOANBER A E, ABOU-GABAL H H. Analytical solution of the Telegraph Point Reactor Kinetics model during the cold start-up of a nuclear reactor[J]. Annals of Nuclear Energy, 2017, 109: 574-582. doi: 10.1016/j.anucene.2017.06.001 [7] ALTAHHAN M R, ABOANBER A E, ABOU-GABAL H H, et al. Response of the point-reactor telegraph kinetics to time varying reactivities[J]. Progress in Nuclear Energy, 2017, 98: 109-122. doi: 10.1016/j.pnucene.2017.03.008 [8] ALTAHHAN M R, AVRAMOVA M N, IVANOV K N. Reactor dynamics based on the neutron telegraph equation[J]. Transactions of the American Nuclear Society, 2017, 117: 753-756. [9] ALTAHHAN M R, BALESTRA P, HOU J, et al. Implementation of the multigroup telegraph based P1 approximation and comparison to the multigroup diffusion based P1 approximation in Gen-Foam[C]//Proceedings of the Fourth International Conference on Physics and Technology of Reactors and Applications. Marrakech, Morocco: GMTR, 2018. [10] ABOANBER A E, NAHLA A A, EL-MHLAWY A M, et al. An efficient exponential representation for solving the two-energy group point telegraph kinetics model[J]. Annals of Nuclear Energy, 2022, 166: 108698. doi: 10.1016/j.anucene.2021.108698 [11] ZAREI M. An adjoint sensitivity analysis of the Telegrapher's neutron kinetic equations[J]. Annals of Nuclear Energy, 2022, 167: 108839. doi: 10.1016/j.anucene.2021.108839 [12] 吴宏春,曹良志,郑友琦,等. 核反应堆物理[M]. 北京: 中国原子能出版社,2014: 51-54. -
下载:
下载:
计量
- 文章访问数: 523
- HTML全文浏览量: 48
- PDF下载量: 27
- 被引次数: 0
