CN103049610B - 一种反应堆设计方法 - Google Patents

Abstract

本发明一种反应堆设计方法，包括首先测得缓发中子产额并写入数据库中，然后计算出缓发中子有效份额和中子代时间，再计算反应堆的动态特性，直到符合设计标准等步骤。本发明不受反应堆堆型、堆芯材料等的限制，可以广泛应用在各种类型的反应堆上，不需要修改蒙卡源程序，具有良好的可实现性。

Description

一种反应堆设计方法

技术领域

本发明属于反应堆设计的方法，具体涉及采用蒙特卡罗方法计算反应堆的两个重要动态参数——缓发中子有效份额（β_eff）和中子代时间（Λ），并将其应用于反应堆设计的方法。

背景技术

缓发中子有效份额（β_eff）和中子代时间（Λ）是反应堆设计任务所必须的两个动态参数，它们对于反应堆的瞬态特性分析十分重要。近年来，β_eff和Λ的计算越来越多的采用蒙特卡罗（蒙卡）方法，这是因为蒙卡方法具有物理图像直观、可描述任意几何等传统计算方法不具备的优点。但是，连续点截面的蒙特卡罗方法不能直接求解β_eff和Λ定义式中的通量伴随函数，因此在反应堆设计任务中，使用蒙卡方法计算动态参数仍存在着困难，同时也是反应堆设计工作一直想攻克的难题。

计算缓发中子有效份额β_eff和Λ的蒙卡方法主要有以下几种：（1）构造通量伴随函数的替代函数。例如Nauchi提出用裂变产生的下一代中子数来代替中子伴随函数，从而计算β_eff和Λ；（2）微扰法，例如Nagaya和B.Verboomen根据相应的微扰理论分别计算出β_eff和Λ。（3）其它方法，包括计算β_eff的k本征值法和计算Λ的瞬发中子密度衰减法等。这些方法或者需要增加复杂的抽样，不易实现，或者存在着较大误差，难以满足反应堆设计的动态参数计算要求。

基于上述原因，在反应堆设计任务中，采用蒙卡方法计算β_eff和Λ方面还存在问题，为了建立具有自主知识产权、实现方便且计算结果可靠的蒙卡方法。本发明从微扰理论和瞬发动力学理论出发，成功研制了计算β_eff和Λ的可靠稳定的蒙卡方法，并将其用于反应堆的设计工作。

发明内容

本发明目的在于提供一种实现方便、结果可靠的计算β_eff和Λ的蒙卡方法用于反应堆设计，解决了目前采用蒙卡方法计算β_eff和Λ的难以实现或者计算误差大的问题。

一种反应堆设计的方法，其特殊之处在于：

1】根据已有基准反应堆测量裂变核素的缓发中子产额，并写入数据库中，具体如下：

1.1】在已有基准反应堆中，使用探测器测量缓发中子先驱核的半衰期、中子谱和缓发中子的产额；

1.2】将测得的缓发中子产额、缓发中子谱、半衰期按一定格式记入数据库中；

2】计算缓发中子有效份额β_eff，具体如下：

2.1】将数据库中的裂变核素的缓发中子产额修改为测量值的（1+a）倍形成新的数据库，使用新的数据库计算出对应的反应堆有效增殖因子k(a)，并计算出采用原始数据库的反应堆有效增殖因子k(0)；

2.2】根据泰勒展开，采用数值微分方法，由k(a)的值计算出有效增殖因子k(a)在0点处的导数k'(0)：

k′(0)=(k(a)-k(-a))/2a；

2.3】根据k(0)和k'(0)，利用微扰公式计算出缓发中子有效份额β_eff：

${\mathrm{\β}}_{\mathrm{eff}}=\frac{1}{k\left(0\right)}{k}^{\′}\left(0\right);$

3】计算中子代时间Λ，具体如下；

3.1】使用蒙卡程序，计算出仅考虑瞬发中子的反应堆有效增殖因子k_p，由公式ρ_p=（1-1/k_p）计算出反应性ρ_p；

3.2】采用标准源，使用蒙卡程序记录反应堆中子通量密度随时间变化曲线，并以指数公式拟合曲线，得到α值，所述指数公式为φ(t)=φ_oexp(αt)，式中的φ_o和φ(t)分别是开始拟合时刻和t时刻的中子通量密度；

3.3】由计算公式Λ=ρ_P/α计算出中子代时间Λ；

4】由计算出的参数β_eff和Λ，结合动态分析模型，计算反应堆重要参数包括功率和最高温度随时间的变化，并由计算结果判断反应堆是否符合设计标准；如果反应堆不符合设计标准，则通过调整堆芯水铀比、堆芯体积、反射层材料和厚度来调整β_eff和Λ，重复步骤2和3，直至根据β_eff和Λ值计算出的重要参数符合反应堆的设计标准。

本发明的优点：

1、本发明不受反应堆堆型、堆芯材料等的限制，可以广泛应用在各种类型的反应堆上。

2、本发明不需要修改蒙卡源程序，具有良好的可实现性。

3、本发明具有计算误差小、稳定性好的特点。

附图说明

图1是本发明技术方案图

图2是蒙卡程序记录的归一化中子通量密度衰减曲线；

图3是归一化中子通量密度的指数拟合曲线；

图4是功率随时间变化曲线；

图5是温度随时间变化曲线。

具体实施方式

本发明的技术路线图如图1所示。首先测得缓发中子产额并写入数据库中，然后计算出缓发中子有效份额和中子代时间，再计算反应堆的动态特性，直到符合设计标准，具体环节介绍如下：

1、缓发中子的产额测量和数据库制作

使用探测器测量各缓发中子先驱核的半衰期和中子谱等参数，拟合得到每组缓发中子的产额和每组缓发中子先驱核的衰变常数，目前缓发中子多分为六组；将缓发中子产额、中子谱、衰变常数等参数按一定格式制成数据库，如ENDF数据库。

2、计算缓发中子有效份额β_eff

2.1有效增殖因子计算

使用原始数据库计算出反应堆的有效增殖因子k(0)，再将数据库中的缓发中子产额的测量值修改成原先的（1+a）倍和（1-a）倍，使用新数据库求出反应堆的有效增殖因子k(a)和k(-a)。。

2.2由数值微分的方法计算k(a)在0点处的导数k'(0)

采用数值微分的方法求出k'(0)，如采用二阶泰勒展开，可以得到k'(0)计算公式：

$k^{\′}\left(0\right)=\frac{k\left(a\right)-k(-a)}{2a}---\left(1\right)$

2.3由微扰公式计算出缓发中子有效份额β_eff

将反应堆的缓发中子产额增加a倍，可得微扰公式：

${\mathrm{\β}}_{\mathrm{eff}}=\frac{1}{k\left(0\right)}\frac{\mathrm{dk}\left(a\right)}{\mathrm{da}}{|}_{a=0}---\left(2\right)$

也可以写成为：

${\mathrm{\β}}_{\mathrm{eff}}=\frac{1}{k\left(0\right)}{k}^{\′}\left(0\right)---\left(3\right)$

根据已经求出的k'(0)和k(0)即可求出β_eff

2.4微扰公式的证明

中子输运方程为：

$\mathrm{L\φ}(r,E,\mathrm{\Ω})=\frac{1}{k\left(0\right)}{S}_f(r,E,\mathrm{\Ω})---\left(4\right)$

式中，r,E,Ω分别代表中子的位置、能量和运动方向，k(0)是反应堆的有效增殖因子，φ(r,E,Ω)代表中子在相空间点(r,E,Ω)处的中子通量密度，Lφ(r,E,Ω)为中子消失项，S_f(r,E,Ω)为中子产生项。

中子消失项表达式：

$\mathrm{L\φ}(r,E,\mathrm{\Ω})=\mathrm{\Ω}\·\▿\mathrm{\φ}(r,E,\mathrm{\Ω})+{\mathrm{\Σ}}_t(r,E)\mathrm{\φ}(r,E,\mathrm{\Ω})$

$-{\∫\mathrm{d\Ω}}^{\′}{\∫\mathrm{dE}}^{\′}{\mathrm{\Σ}}_s(r,E^{\′},{\mathrm{\Ω}}^{\′}\→E,\mathrm{\Ω})\mathrm{\φ}(r,E^{\′},{\mathrm{\Ω}}^{\′})---\left(5\right)$

Σ_t(r,E)为中子在相空间点（r，E）处的裂变截面，∑_s(r,E，Ω→E′,Ω′)为从（E,Ω）到（E′,Ω′）的散射截面。式（5）右边第一项代表泄漏出相空间点（r,E，Ω）的中子，第二项代表移出相空间点（r,E，Ω）的中子，第三项代表散射到相空间点（r,E，Ω）的中子。

源项：

$S_f(r,E,\mathrm{\Ω})=S_f^p(r,E,\mathrm{\Ω})+\underset{i}{\mathrm{\Σ}}{S}_{f,i}^d(r,E,\mathrm{\Ω})---\left(6\right)$

（6）式右边第一项为瞬发中子源项，第二项为缓发中子总源项。

瞬发中子源项：

$S_f^p(r,E,\mathrm{\Ω})=\underset{m}{\mathrm{\Σ}}\∫{\mathrm{dE}}^{\′}{\∫\mathrm{d\Ω}}^{\′}{\mathrm{\χ}}_m^p(E\←E^{\′},\mathrm{\Ω})v_m^p\left(E^{\′}\right){\mathrm{\Σ}}_{f,m}\left({r,E}^{\′}\right)\mathrm{\φ}({r,E}^{\′},{\mathrm{\Ω}}^{\′})---\left(7\right)$

式（7）中，为第m种裂变核素的瞬发中子裂变谱，为第m种核素裂变产生的瞬发中子数，∑_f，m(r,E′)为第m种裂变核素的裂变截面。

第i组缓发中子源项：

$S_{f,i}^d(r,E,\mathrm{\Ω})=\underset{m}{\mathrm{\Σ}}{\∫\mathrm{dE}}^{\′}{\∫\mathrm{d\Ω}}^{\′}{\mathrm{\χ}}_{m,i}^d(E\←E^{\′},\mathrm{\Ω})v_{m,i}^d\left(E^{\′}\right){\mathrm{\Σ}}_{f,m}\left({r,E}^{\′}\right)\mathrm{\φ}(r,E^{\′},{\mathrm{\Ω}}^{\′})---\left(8\right)$

为第m种核素裂变产生的第i组缓发中子的裂变谱，为第m种核素裂变产生的第i组缓发中子数。

输运方程的伴随方程，即价值方程为：

$L^{+}{\mathrm{\φ}}^{+}(r,E,\mathrm{\Ω})=\frac{1}{k\left(0\right)}{\mathrm{\Σ}}_m\∫{\mathrm{dE}}^{\′}\∫{\mathrm{d\Ω}}^{\′}[{\mathrm{\χ}}_m^p(E^{\′}\←E,{\mathrm{\Ω}}^{\′})v_m^p\left(E\right)+\underset{i}{\mathrm{\Σ}}{\mathrm{\χ}}_{m,i}^d(E^{\′}\←E,{\mathrm{\Ω}}^{\′})v_{m,i}^d\left(E\right)]---\left(9\right)$

${\×\mathrm{\Σ}}_{f,m}(r,E){\mathrm{\φ}}^{+}(r,E^{\′},{\mathrm{\Ω}}^{\′})$

式中，φ⁺(r,E，Ω)是相空间点（r,E，Ω）的中子价值，L⁺φ⁺(r,E,Ω)代表中子价值消失项，式（9）右边代表中子价值产生项

中子价值消失项：

$L^{+}{\mathrm{\φ}}^{+}(r,E,\mathrm{\Ω})=-{\mathrm{\Ω}\·\▿\mathrm{\φ}}^{+}(r,E,\mathrm{\Ω})+E_t(r,E){\mathrm{\φ}}^{+}(r,E,\mathrm{\Ω})$

（10）

$-{\∫\mathrm{d\Ω}}^{\′}{\∫\mathrm{dE}}^{\′}{\mathrm{\Σ}}_s(r,E^{\′},{\mathrm{\Ω}}^{\′}\←E,\mathrm{\Ω}){\mathrm{\φ}}^{+}(r,E^{\′},{\mathrm{\Ω}}^{\′})$

上式中，右边第一项代表中子沿方向无碰撞飞行时，中子价值的增加量，第二项代表发生碰撞的中子价值移出量，第三项代表散射进来的中子价值增加量。

如果在反应堆中加一个微扰，将缓发中子产额增加a倍，那么可以得到扰动后中子输运方程：

$\mathrm{L\φ}(r,E,\mathrm{\Ω},a)=\frac{1}{k\left(a\right)}{S}_f(r,E,\mathrm{\Ω},a)---\left(11\right)$

$S_f(r,E,\mathrm{\Ω},a)=S_f^p(r,E,\mathrm{\Ω},a)+(1+a)\underset{i}{\mathrm{\Σ}}{S}_{f,i}^d(r,E,\mathrm{\Ω},a)---\left(12\right)$

$S_f^p(r,E,\mathrm{\Ω},a)=\underset{m}{\mathrm{\Σ}}{\∫\mathrm{dE}}^{\′}{\∫\mathrm{d\Ω}}^{\′}{\mathrm{\χ}}_m^p(E\←E^{\′},\mathrm{\Ω})v_m^p\left(E^{\′}\right){\mathrm{\Σ}}_{f,m}(r,E^{\′})\×\mathrm{\φ}(r,E^{\′},{\mathrm{\Ω}}^{\′},a)---\left(13\right)$

$S_{f,i}^d(r,E,\mathrm{\Ω},a)=\underset{m}{\mathrm{\Σ}}{\∫\mathrm{dE}}^{\′}\∫d{\mathrm{\Ω}}^{\′}{\mathrm{\χ}}_{m,i}^d(E\←E^{\′},\mathrm{\Ω})v_{m,i}^d\left(E^{\′}\right){\mathrm{\Σ}}_{f,m}(r,E^{\′})\×\mathrm{\φ}(r,E^{\′},{\mathrm{\Ω}}^{\′},a)---\left(14\right)$

$\mathrm{L\φ}(r,E,\mathrm{\Ω},a)=\mathrm{\Ω}\·\▿\mathrm{\φ}(r,E,\mathrm{\Ω},a)+{\mathrm{\Σ}}_t(r,E)\mathrm{\φ}(r,E,\mathrm{\Ω},a)$

（15）

$-\∫{\mathrm{d\Ω}}^{\′}{\∫\mathrm{dE}}^{\′}{\mathrm{\Σ}}_s(r,E,\mathrm{\Ω}\←E^{\′},{\mathrm{\Ω}}^{\′})\mathrm{\φ}(r,E^{\′},{\mathrm{\Ω}}^{\′},a)$

方程（11）两边都乘上φ⁺(r,E，Ω)，方程（9）两边同乘φ(r,E,Ω,a)，分别对相空间(r,E,Ω)积分，再相减，可得：

$<{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})\&CenterDot;\mathrm{L\&phi;}(r,E,\mathrm{\&Omega;},a)>-<\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a)\&CenterDot;L^{+}{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>$

$=\frac{1}{k\left(a\right)}<S_f(r,E,\mathrm{\&Omega;},a){\&CenterDot;\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>---\left(16\right)$

将式（10）和式（15）代入（16）式左边，得到：

（16）式左边=

$<\mathrm{\&Omega;}\&CenterDot;\&dtri;\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a){\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})+\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a)\&CenterDot;\mathrm{\&Omega;}\&CenterDot;\&dtri;{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>$

$+<{\mathrm{\&Sigma;}}_t(r,E)\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a){\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})-{\mathrm{\&Sigma;}}_t(r,E){\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a)>$

$=0$

（16）式左边=0，那么（16）式右边=0，这样就有：

$\frac{1}{k\left(a\right)}<S_f(r,E,\mathrm{\&Omega;},a)\&CenterDot;{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>$

$=\&Integral;\mathrm{dr}\&Integral;\mathrm{dE}\&Integral;\mathrm{d\&Omega;}\frac{1}{k\left(0\right)}{\mathrm{\&Sigma;}}_m{\&Integral;\mathrm{dE}}^{\&prime;}{\&Integral;\mathrm{d\&Omega;}}^{\&prime;}[{\mathrm{\&chi;}}_m^p(E^{\&prime;}\&LeftArrow;E,{\mathrm{\&Omega;}}^{\&prime;})v_m^p\left(E\right)$

$+\underset{i}{\mathrm{\&Sigma;}}{\mathrm{\&chi;}}_{m,i}^d(E^{\&prime;}\&LeftArrow;E,{\mathrm{\&Omega;}}^{\&prime;})v_{m,i}^d\left(E\right)]{\&times;\mathrm{\&Sigma;}}_{f,m}(r,E){\&times;\mathrm{\&phi;}}^{+}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;})\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a)---\left(18\right)$

$=\frac{1}{k\left(0\right)}{\mathrm{\&Sigma;}}_m\&Integral;d{E}^{\&prime;}\&Integral;d{\mathrm{\&Omega;}}^{\&prime;}{\mathrm{\&phi;}}^{+}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;})\&Integral;\mathrm{dr}\&Integral;\mathrm{dE}\&Integral;\mathrm{d\&Omega;}[{\mathrm{\&chi;}}_m^p(E^{\&prime;}\&LeftArrow;E,{\mathrm{\&Omega;}}^{\&prime;})v_m^p\left(E\right)+$

$+\underset{i}{\mathrm{\&Sigma;}}{\mathrm{\&chi;}}_{m,i}^d(E^{\&prime;}\&LeftArrow;E,{\mathrm{\&Omega;}}^{\&prime;})v_{m,i}^d\left(E\right)]\&times;{\mathrm{\&Sigma;}}_{f,m}(r,E)\&times;\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},a)$

$=\frac{1}{k\left(0\right)}<S_f^p(r,E,\mathrm{\&Omega;},a){\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})+\underset{i}{\mathrm{\&Sigma;}}{S}_{f,i}^d(r,E,\mathrm{\&Omega;},a){\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>$

由（18）式，得：

$\frac{1}{k\left(0\right)}\frac{k\left(a\right)-k\left(0\right)}{a}=\frac{<\underset{i}{\mathrm{\&Sigma;}}{S}_{f,i}^d(r,E,\mathrm{\&Omega;},a)\&CenterDot;{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>}{<[S_f^p(r,E,\mathrm{\&Omega;},a)+\underset{i}{\mathrm{\&Sigma;}}{S}_f^d(r,E,\mathrm{\&Omega;},a)]{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})>}---\left(19\right)$

再由β_eff的定义式：

${\mathrm{\&beta;}}_{\mathrm{eff}}=\frac{<{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})\underset{i}{\mathrm{\&Sigma;}}{S}_{f,i}^d(r,E,\mathrm{\&Omega;})>}{<{\mathrm{\&phi;}}^{+}(r,E,\mathrm{\&Omega;})S_f(r,E,\mathrm{\&Omega;})>}---\left(20\right)$

在a→0时，就证明了计算β_eff的微扰理论：

${\mathrm{\&beta;}}_{\mathrm{eff}}=\frac{1}{k\left(0\right)}\frac{\mathrm{dk}\left(a\right)}{\mathrm{da}}{|}_{a=0}---\left(21\right)$

3、计算中子代时间

3.1仅考虑瞬发中子的反应堆有效增殖因子k_p和反应性ρ_p的计算

使用蒙卡程序，计算出仅考虑瞬发中子的反应堆有效增殖因子k_p，并由公式ρ_p=（1-1/k_p）计算出反应性ρ_p。

3.2使用蒙卡程序记录反应堆中子通量密度随时间的变化曲线

采用标准源，使用蒙卡程序记录反应堆中子通量密度随时间的变化曲线，之后采用指数公式φ(t)=φ_o exp(αt)拟合曲线得到α值，所述指数公式中的φ_o和φ(t)分别是开始拟合时刻和t时刻的中子通量密度。

3.3采用本发明的公式计算中子代时间

本发明提出的计算公式为：

$\mathrm{\&Lambda;}=\frac{{\mathrm{\&rho;}}_P}{\mathrm{\&alpha;}}---\left(22\right)$

3.4、本发明公式（22）的理论推导

中子代时间的定义式为：

$\mathrm{\&Lambda;}=\frac{1}{F\left(t\right)}<\frac{\mathrm{\&phi;}(r,E,\mathrm{\&Omega;},t)}{\mathrm{\&upsi;}}{\mathrm{\&phi;}}_o^{+}(r,E,\mathrm{\&Omega;})>---\left(23\right)$

式中，υ是中子速度，φ(r,E,Ω,t)是t时刻的中子通量密度，φ_o ⁺(r,E,Ω)为临界反应堆的中子价值， $F\left(t\right)=\&Integral;\&Integral;\&Integral;\tilde{x}\left(E\right){\mathrm{v\&Sigma;}}_f(r,E^{\&prime;})\mathrm{\&phi;}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;},t){\mathrm{\&phi;}}_o^{+}(r,E,\mathrm{\&Omega;})\mathrm{drd}{\mathrm{\&Omega;}}^{\&prime;}d{E}^{\&prime;}\mathrm{d\&Omega;dE},$ 代表的是所有中子的总价值。

式（22）的推导分为以下两部分。

3.4.1仅考虑瞬发中子的反应堆的反应性ρ_p的定义

仅考虑瞬发中子的临界反应堆中子输运方程为：

$\mathrm{\&Omega;}\&CenterDot;\&dtri;{\mathrm{\&phi;}}_o(r,E,\mathrm{\&Omega;})+{\mathrm{\&Sigma;}}_t(r,E){\mathrm{\&phi;}}_o(r,E,\mathrm{\&Omega;})=$

$\&Integral;\&Integral;{\mathrm{\&Sigma;}}_s(r,E^{\&prime;})f_s(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;}\&RightArrow;E,\mathrm{\&Omega;}){\mathrm{\&phi;}}_o(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;})d{\mathrm{\&Omega;}}^{\&prime;}d{E}^{\&prime;}---\left(24\right)$

$+\&Integral;\&Integral;{\widetilde{\mathrm{\&chi;}}}_p\left(E^{\&prime;}\right)v{\mathrm{\&Sigma;}}_f(r,E^{\&prime;}){\mathrm{\&phi;}}_o(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;})d{\mathrm{\&Omega;}}^{\&prime;}d{E}^{\&prime;}$

式中，φ_o(r,E，Ω)为临界状态的中子通量密度，为瞬发中子裂变谱。

仅考虑瞬发中子的临界反应堆中子价值方程：

$-\mathrm{\&Omega;}\&CenterDot;\&dtri;{{\mathrm{\&phi;}}_o}^{+}(r,E,\mathrm{\&Omega;})+{\mathrm{\&Sigma;}}_t(r,E){{\mathrm{\&phi;}}_o}^{+}(r,E,\mathrm{\&Omega;})$

$={\mathrm{\&Sigma;}}_s(r,E)\&Integral;\&Integral;f_s(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;}\&LeftArrow;E,\mathrm{\&Omega;}){{\mathrm{\&phi;}}_o}^{+}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;}){\mathrm{d\&Omega;}}^{\&prime;}d{E}^{\&prime;}---\left(25\right)$

$+v{\mathrm{\&Sigma;}}_f(r,E)\&Integral;\&Integral;{\widetilde{\mathrm{\&chi;}}}_p\left(E^{\&prime;}\right){{\mathrm{\&phi;}}_o}^{+}(r,E^{\&prime;}{\mathrm{\&Omega;}}^{\&prime;})d{\mathrm{\&Omega;}}^{\&prime;}d{E}^{\&prime;}$

在反应堆引入微扰，Σ_t，∑_s f_s，v∑_f，φ分别变为∑_t ^*， φ^*。

扰动后的反应堆中子输运方程为：

$\mathrm{\&Omega;}\&CenterDot;\&dtri;{\mathrm{\&phi;}}^{*}(r,E,\mathrm{\&Omega;},t)+{\mathrm{\&Sigma;}}_t^{*}(r,E){\mathrm{\&phi;}}^{*}(r,E,\mathrm{\&Omega;},t)=$

$\&Integral;\&Integral;{\mathrm{\&Sigma;}}_s^{*}(r,E^{\&prime;})f_s^{*}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;}\&RightArrow;E,\mathrm{\&Omega;}){\mathrm{\&phi;}}^{*}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;},t)d{\mathrm{\&Omega;}}^{\&prime;}d{E}^{\&prime;}---\left(26\right)$

$+\frac{1}{k^{*}}\&Integral;\&Integral;{\widetilde{\mathrm{\&chi;}}}_p\left(E^{\&prime;}\right)v{\mathrm{\&Sigma;}}_f^{*}(r,E^{\&prime;}){\mathrm{\&phi;}}^{*}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;},t){\mathrm{d\&Omega;}}^{\&prime;}{\mathrm{dE}}^{\&prime;}$

（26）式两边同乘φ_o ⁺(r,E,Ω)，（25）式两边同乘φ^*(r,E,Ω,t)，并对整个相空间(r,E,Ω)积分，再相减，用<>代表积分符∫∫∫dVdEdΩ，得到：

$<{{\mathrm{\&phi;}}_o}^{+}{\mathrm{\&Omega;}\&CenterDot;\&dtri;\mathrm{\&phi;}}^{*}+{\mathrm{\&phi;}}^{*}{\mathrm{\&Omega;}\&CenterDot;\&dtri;\mathrm{\&phi;}}_o^{+}>+<({\mathrm{\&Sigma;}}_t^{*}-{\mathrm{\&Sigma;}}_t){{\mathrm{\&phi;}}_o}^{+}{\mathrm{\&phi;}}^{*}>$

$=\&Integral;\&Integral;<({\mathrm{\&Sigma;}}_s^{*}{f}_s^{*}-{\mathrm{\&Sigma;}}_s{f}_s){\mathrm{\&phi;}}^{*}>{{\mathrm{\&phi;}}_o}^{+}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;}){\mathrm{dE}}^{\&prime;}d{\mathrm{\&Omega;}}^{\&prime;}+---\left(27\right)$

$\&Integral;\&Integral;<{\widetilde{\mathrm{\&chi;}}}_p\left(E\right)(\frac{v^{*}{{\mathrm{\&Sigma;}}_f}^{*}}{k^{*}}-v{\mathrm{\&Sigma;}}_f){\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{+}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;})>d{E}^{\&prime;}\mathrm{d\&Omega;}$

令ΔΣ_t＝Σ_t ^*-Σ_t，Δ（vΣ_f）＝v^*Σ_f ^*v∑_f，Δk=k^*-1，并使用代替经变换可得到（27）式左边第一项为0，左边第二项为∫∫<ΔΣ_tφ^*φ_o′⁺>dE'dΩ'，右边第一项为∫∫<Δ(∑_sf_s)φ^*φ_o′⁺>dE'dΩ'。

经整理，可得仅考虑瞬发中子的反应堆的反应性ρ_p的严格定义：

${\mathrm{\&rho;}}_P=\frac{\mathrm{\&Delta;k}}{k^{*}}=\frac{1}{\&Integral;\&Integral;<{\widetilde{\mathrm{\&chi;}}}_p\left(E\right)v{\mathrm{\&Sigma;}}_f^{*}{\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}}---\left(28\right)$

$\&times;(\&Integral;\&Integral;<{\widetilde{\mathrm{\&chi;}}}_p\left(E\right)\mathrm{\&Delta;}\left({\mathrm{v\&Sigma;}}_f\right){\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}d{\mathrm{\&Omega;}}^{\&prime;}+\&Integral;\&Integral;<\mathrm{\&Delta;}\left({\mathrm{\&Sigma;}}_s{f}_s\right){\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}-<{\mathrm{\&Delta;\&Sigma;}}_t{\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{+}>)$

令φ^*＝φ_o+Δφ,φ_o为临界状态下的中子通量密度。代入式（28），注意到式（28）的每个含φ^*的子式中，都有另一个无穷小量，略去所有二阶小量，可以得到反应性ρ_p近似定义：

${\mathrm{\&rho;}}_P\&ap;\frac{1}{\&Integral;\&Integral;<{\widetilde{\mathrm{\&chi;}}}_p\left(E\right)\mathrm{\&Delta;}\left({\mathrm{v\&Sigma;}}_f\right){\mathrm{\&phi;}}_o{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}}---\left(29\right)$

$\&times;(\&Integral;\&Integral;<{\widetilde{\mathrm{\&chi;}}}_p\left(E\right)\mathrm{\&Delta;}\left({\mathrm{v\&Sigma;}}_f\right){\mathrm{\&phi;}}_o{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}+\&Integral;\&Integral;<\mathrm{\&Delta;}\left({\mathrm{\&Sigma;}}_s{f}_s\right){\mathrm{\&phi;}}_o{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}-<\mathrm{\&Delta;}{\mathrm{\&Sigma;}}_t{\mathrm{\&phi;}}_o{{\mathrm{\&phi;}}_o}^{+}>)$

3.4.2瞬发动力学方程推导

首先可以将通量φ^*(r,E，Ω,t)写成为幅因子n(t)和形状因子ψ(r,E,Ω,t)的乘积：

φ^*(r,E,Ω,t)＝n(t)·ψ(r,E,Ω,t)    （30）

当系统从临界状态受到轻微的扰动时，形状因子随时间变化很小，有： $\frac{\&PartialD;}{\&PartialD;t}<\frac{1}{\mathrm{\&upsi;}}{\mathrm{\&phi;}}_o^{+}(r,E,\mathrm{\&Omega;})\mathrm{\&psi;}(r,E,\mathrm{\&Omega;},t)>=0,$ υ为中子速度。

这样就有：

$<\frac{1}{\mathrm{\&upsi;}}{\mathrm{\&phi;}}_o^{+}(r,E,\mathrm{\&Omega;})\frac{{\&PartialD;\mathrm{\&phi;}}^{*}(r,E,\mathrm{\&Omega;},t)}{\&PartialD;t}>=\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}<\frac{1}{\mathrm{\&upsi;}}{\mathrm{\&phi;}}_o^{+}(r,E,\mathrm{\&Omega;})\mathrm{\&psi;}(r,E,\mathrm{\&Omega;},t)>---\left(31\right)$

中子输运方程：

$\frac{1}{\mathrm{\&upsi;}}\frac{{\&PartialD;\mathrm{\&phi;}}^{*}(r,E,\mathrm{\&Omega;},t)}{\&PartialD;t}+\mathrm{\&Omega;}\&CenterDot;\&dtri;{\mathrm{\&phi;}}^{*}(r,E,\mathrm{\&Omega;},t)+{\mathrm{\&Sigma;}}_t^{*}(r,E){\mathrm{\&phi;}}^{*}(r,E,\mathrm{\&Omega;},t)$

$=\&Integral;\&Integral;{\mathrm{\&Sigma;}}_s^{*}(r,E^{\&prime;})f_s^{*}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;}\&RightArrow;E,\mathrm{\&Omega;}){\mathrm{\&phi;}}^{*}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;},t){\mathrm{d\&Omega;}}^{\&prime;}{\mathrm{dE}}^{\&prime;}---\left(32\right)$

$+\&Integral;\&Integral;\widetilde{\mathrm{\&chi;}}\left(E^{\&prime;}\right)v^{*}{\mathrm{\&Sigma;}}_f^{*}(r,E^{\&prime;}){\mathrm{\&phi;}}^{*}(r,E^{\&prime;},{\mathrm{\&Omega;}}^{\&prime;},t){\mathrm{d\&Omega;}}^{\&prime;}{\mathrm{dE}}^{\&prime;}$

用φ_o ⁺(r,E,Ω)乘（32）式，用φ^*(r,E，Ω,t)乘（26）式，并对整个相空间(r,E,Ω)积分，相减后可得：

$<\frac{1}{\mathrm{\&upsi;}}\frac{{\&PartialD;\mathrm{\&phi;}}^{*}}{\&PartialD;t}{{\mathrm{\&phi;}}_o}^{+}>+<{{\mathrm{\&phi;}}_o}^{+}{\mathrm{\&Omega;}\&CenterDot;\&dtri;\mathrm{\&phi;}}^{*}+{\mathrm{\&phi;}}^{*}\mathrm{\&Omega;}\&CenterDot;\&dtri;\mathrm{\&phi;}>+<({\mathrm{\&Sigma;}}_t^{*}-{\mathrm{\&Sigma;}}_t){{\mathrm{\&phi;}}_o}^{+}{\mathrm{\&phi;}}^{*}>=---\left(33\right)$

$\&Integral;\&Integral;<({\mathrm{\&Sigma;}}_s^{*}{f}_s^{*}-{\mathrm{\&Sigma;}}_s{f}_s){\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}+\&Integral;\&Integral;<\widetilde{\mathrm{\&chi;}}\left(E\right)(v^{*}{{\mathrm{\&Sigma;}}_f}^{*}-{\mathrm{v\&Sigma;}}_f){\mathrm{\&phi;}}^{*}{{\mathrm{\&phi;}}_o}^{\&prime;+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}$

经整理，得到：

$\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}=\frac{(\&Integral;\&Integral;<\mathrm{\&Delta;}\left({\mathrm{\&Sigma;}}_s{f}_s\right){{\mathrm{\&psi;\&phi;}}_o}^{+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}+\&Integral;\&Integral;<\widetilde{\mathrm{\&chi;}}\left(E\right)\mathrm{\&Delta;}\left({\mathrm{v\&Sigma;}}_f\right){{\mathrm{\&psi;\&phi;}}_o}^{+}>{\mathrm{dE}}^{\&prime;}{\mathrm{d\&Omega;}}^{\&prime;}-<{\mathrm{\&Delta;\&Sigma;}}_t{{\mathrm{\&psi;\&phi;}}_o}^{+}>)}{<\frac{1}{\mathrm{\&upsi;}}{{\mathrm{\&psi;\&phi;}}_o}^{+}>}n\left(t\right)---\left(34\right)$

由中子代时间定义（23）和反应性定义式（29），可得瞬发动力学方程：

$\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}=\frac{{\mathrm{\&rho;}}_P}{\mathrm{\&Lambda;}}n\left(t\right)---\left(35\right)$

在与临界状态差别不大的情况下，形状因子变化很小，有：

ψ(r,E,Ω,t)≈ψ_o(r,E,Ω)    （36）

那么由式（32）和（37）就有：

$\frac{{\mathrm{d\&phi;}}^{*}\left(t\right)}{\mathrm{dt}}=\frac{{\mathrm{\&rho;}}_P}{\mathrm{\&Lambda;}}{\mathrm{\&phi;}}^{*}\left(t\right)---\left(37\right)$

可得：

${\mathrm{\&phi;}}^{*}\left(t\right)={\mathrm{\&phi;}}_o\exp \left(\frac{{\mathrm{\&rho;}}_P}{\mathrm{\&Lambda;}}t\right)---\left(38\right)$

令：

$\mathrm{\&alpha;}=\frac{{\mathrm{\&rho;}}_P}{\mathrm{\&Lambda;}}---\left(39\right)$

可得：

$\mathrm{\&Lambda;}=\frac{{\mathrm{\&rho;}}_P}{\mathrm{\&alpha;}}---\left(40\right)$

4、计算反应堆的重要参数，判断是否满足设计标准

由计算出的参数β_eff和Λ，选择一定的动态分析模型，如在西安脉冲堆的计算中选择了Nordheim-Fuchs模型（该模型将在实施例中详述，值得一提的是，无论采用何种动态分析模型，β_eff和Λ都是所必须的参数，本发明也因此具有很强的实用性），计算反应堆重要参数包括功率、最高温度等随时间的变化曲线，并由计算结果判断反应堆是否符合设计标准；如果反应堆不符合设计标准，则通过调整堆芯水铀比、堆芯体积、反射层材料和厚度来调整β_eff和Λ，重复步骤2和3，直至根据β_eff和Λ值计算出的重要参数符合反应堆的设计标准。

实施例：

以西安脉冲堆为例，介绍本发明的应用。

第一步，使用探测器测量各缓发中子先驱核的半衰期和中子谱等参数，拟合得到六组缓发中子的产额和相应缓发中子先驱核的衰变常数，并将缓发中子产额、中子谱、衰变常数等参数按一定格式写入数据库。

第二步，采用未经修改的原始数据库，计算出反应堆的有效增殖因子k(0)。再将数据库中的裂变核素²³⁵U、²³⁸U的缓发中子产额提高和降低0.25倍，使用蒙卡程序调用相应数据库，计算出反应堆有效增殖因子k(0.25)和k(-0.25)，再采用公式k′(0)＝2(k(0.25)-k(-0.25))计算出k′(0)，然后根据微扰公式计算出β_eff。β_eff计算结果稳定在7.267‰。

第三步，使用MCNP程序，计算出仅考虑瞬发中子的反应堆的有效增殖因子k_p=0.99500±0.00005,计算出对应的反应性ρ_P=-（5.025±0.05）‰。使用蒙卡程序记录堆芯中子通量密度随着时间的变化，如附图2所示，在曲线的开始阶段有一段不稳定的区域，曲线的拟合需要略过不稳定区域，本发明在2ms-8ms的时间区间内以函数φ(t)=φ_o exp(αt)拟合中子通量密度随时间的变化曲线，如附图3所示，拟合函数曲线与中子通量密度衰减曲线基本上是重合的，拟合得到α=-（137.5±0.25）s^-1。根据公式计算出中子代时间Λ=（36.55±0.36）μs。

第四步，在西安脉冲堆引入反应性ρ₀=0.02482以发射脉冲，分析反应堆在引入最大脉冲反应性时是否符合设计标准。脉冲棒在95.8ms内匀速弹出堆芯，动态分析模型采用Nordheim-Fuchs模型。Nordheim-Fuchs模型的核心方程：

$\frac{\mathrm{dP}\left(t\right)}{\mathrm{dt}}=\frac{\mathrm{\&rho;}\left(t\right)-{\mathrm{\&beta;}}_{\mathrm{eff}}}{\mathrm{\&Lambda;}}P\left(t\right)---\left(41\right)$

式中，P(t)是功率，ρ(t)是反应性。

结合基本的传热原理，即可得到反应堆设计中几个非常重要的参数的变化曲线，附图4、5分别是功率和平均最高温度随着时间的变化曲线。考虑燃料芯块内部的功率分布和反应堆的功率峰因子，计算得到的燃料最热点的最高温度为828.4℃。

第五步，判断重要参数是否符合设计标准。西安脉冲堆的设计安全限值为：在包壳温度大于500℃时，燃料芯块最大温度限值为970℃，在包壳温度小于500℃时，燃料芯块最大温度限值为1150℃。因此，在引入最大反应性ρ₀=0.02482时，西安脉冲反应堆是符合设计安全标准的。如果经计算，反应堆不符合设计标准，则需调整堆芯结构、材料等。具体可以通过调整堆芯水铀比、堆芯体积、反射层材料和厚度等来调整，然后再根据本发明的方法得到反应堆重要参数值，直到反应堆重要参数的值满足设计标准。

Claims (1)

1.一种反应堆设计方法，其特征在于：

1】根据已有基准反应堆测量裂变核素的缓发中子产额，并写入数据库中，具体如下：

1.1】在已有基准反应堆中，使用探测器测量缓发中子先驱核的半衰期、中子谱和缓发中子的产额；

1.2】将测得的缓发中子产额、缓发中子谱、半衰期按一定格式记入数据库中；

2】计算缓发中子有效份额β_eff，具体如下：

2.1】将数据库中的裂变核素的缓发中子产额修改为测量值的(1+a)倍形成新的数据库，使用新的数据库计算出对应的反应堆有效增殖因子k(a)，将数据库中的裂变核素的缓发中子产额修改为测量值的(1-a)倍形成新的数据库，使用新的数据库计算出对应的反应堆有效增殖因子k(-a)，并计算出采用原始数据库的反应堆有效增殖因子k(0)；

2.2】根据泰勒展开，采用数值微分方法，由k(a)的值计算出有效增殖因子k(a)在0点处的导数k'(0)：

k′(0)＝(k(a)-k(-a))/2a；

2.3】根据k(0)和k'(0)，利用微扰公式计算出缓发中子有效份额β_eff：

${\mathrm{\&beta;}}_{\mathrm{eff}}=\frac{1}{K\left(0\right)}{k}^{\&prime;}\left(0\right);$

3】计算中子代时间Λ，具体如下；

3.1】使用蒙卡程序，计算出仅考虑瞬发中子的反应堆有效增殖因子k_p，由公式ρ_p＝(1-1/k_p)计算出反应性ρ_p；

3.2】采用标准源，使用蒙卡程序记录反应堆中子通量密度随时间变化曲线，并以指数公式拟合曲线，得到α值，所述指数公式为φ(t)＝φ_oexp(αt)，式中的φ_o和φ(t)分别是开始拟合时刻和t时刻的中子通量密度；

3.3】由计算公式Λ＝ρ_P/α计算出中子代时间Λ；

4】由计算出的参数β_eff和Λ，结合动态分析模型，计算反应堆重要参数包括功率和最高温度随时间的变化，并由计算结果判断反应堆是否符合设计标准；如果反应堆不符合设计标准，则通过调整堆芯水铀比、堆芯体积、反射层材料和厚度来调整β_eff和Λ，重复步骤2和3，直至根据β_eff和Λ值计算出的重要参数符合反应堆的设计标准。