说明书一种反应堆设计方法
技术领域
本发明属于反应堆设计的方法，具体涉及采用蒙特卡罗方法计算反应堆的两个重要动态参数——缓发中子有效份额（β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：
3】计算中子代时间Λ，具体如下；
3.1】使用蒙卡程序，计算出仅考虑瞬发中子的反应堆有效增殖因子kp，由公式ρp=（1‑1/kp）计算出反应性ρ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)计算公式：
2.3由微扰公式计算出缓发中子有效份额βeff
将反应堆的缓发中子产额增加a倍，可得微扰公式：
也可以写成为：
根据已经求出的k'(0)和k(0)即可求出βeff
2.4微扰公式的证明
中子输运方程为：
式中，r,E,Ω分别代表中子的位置、能量和运动方向，k(0)是反应堆的有效增殖因子，φ(r,E,Ω)代表中子在相空间点(r,E,Ω)处的中子通量密度，Lφ(r,E,Ω)为中子消失项，Sf(r,E,Ω)为中子产生项。
中子消失项表达式：
Σt(r,E)为中子在相空间点（r，E）处的裂变截面，∑s(r,E，Ω→E′,Ω′)为从（E,Ω）到（E′,Ω′）的散射截面。式（5）右边第一项代表泄漏出相空间点（r,E，Ω）的中子，第二项代表移出相空间点（r,E，Ω）的中子，第三项代表散射到相空间点（r,E，Ω）的中子。
源项：
（6）式右边第一项为瞬发中子源项，第二项为缓发中子总源项。
瞬发中子源项：
式（7）中，为第m种裂变核素的瞬发中子裂变谱，为第m种核素裂变产生的瞬发中子数，∑f，m(r,E′)为第m种裂变核素的裂变截面。
第i组缓发中子源项：
为第m种核素裂变产生的第i组缓发中子的裂变谱，为第m种核素裂变产生的第i组缓发中子数。
输运方程的伴随方程，即价值方程为：
式中，φ+(r,E，Ω)是相空间点（r,E，Ω）的中子价值，L+φ+(r,E,Ω)代表中子价值消失项，式（9）右边代表中子价值产生项
中子价值消失项：
（10）
上式中，右边第一项代表中子沿方向无碰撞飞行时，中子价值的增加量，第二项代表发生碰撞的中子价值移出量，第三项代表散射进来的中子价值增加量。
如果在反应堆中加一个微扰，将缓发中子产额增加a倍，那么可以得到扰动后中子输运方程：
（15）
方程（11）两边都乘上φ+(r,E，Ω)，方程（9）两边同乘φ(r,E,Ω,a)，分别对相空间(r,E,Ω)积分，再相减，可得：
将式（10）和式（15）代入（16）式左边，得到：
（16）式左边=
（16）式左边=0，那么（16）式右边=0，这样就有：
由（18）式，得：
再由βeff的定义式：
在a→0时，就证明了计算βeff的微扰理论：
3、计算中子代时间
3.1仅考虑瞬发中子的反应堆有效增殖因子kp和反应性ρp的计算
使用蒙卡程序，计算出仅考虑瞬发中子的反应堆有效增殖因子kp，并由公式ρp=（1‑1/kp）计算出反应性ρp。
3.2使用蒙卡程序记录反应堆中子通量密度随时间的变化曲线
采用标准源，使用蒙卡程序记录反应堆中子通量密度随时间的变化曲线，之后采用指数公式φ(t)=φo exp(αt)拟合曲线得到α值，所述指数公式中的φo和φ(t)分别是开始拟合时刻和t时刻的中子通量密度。
3.3采用本发明的公式计算中子代时间
本发明提出的计算公式为：
3.4、本发明公式（22）的理论推导
中子代时间的定义式为：
式中，υ是中子速度，φ(r,E,Ω,t)是t时刻的中子通量密度，φo+(r,E,Ω)为临界反应堆的中子价值，代表的是所有中子的总价值。
式（22）的推导分为以下两部分。
3.4.1仅考虑瞬发中子的反应堆的反应性ρp的定义
仅考虑瞬发中子的临界反应堆中子输运方程为：
式中，φo(r,E，Ω)为临界状态的中子通量密度，为瞬发中子裂变谱。
仅考虑瞬发中子的临界反应堆中子价值方程：
在反应堆引入微扰，Σt，∑s fs，v∑f，φ分别变为∑t*，φ*。
扰动后的反应堆中子输运方程为：
（26）式两边同乘φo+(r,E,Ω)，（25）式两边同乘φ*(r,E,Ω,t)，并对整个相空间(r,E,Ω)积分，再相减，用<>代表积分符∫∫∫dVdEdΩ，得到：
令ΔΣt＝Σt*‑Σt，Δ（vΣf）＝v*Σf*v∑f，Δk=k*‑1，并使用代替经变换可得到（27）式左边第一项为0，左边第二项为∫∫<ΔΣtφ*φo′+>dE'dΩ'，右边第一项为∫∫<Δ(∑sfs)φ*φo′+>dE'dΩ'。
经整理，可得仅考虑瞬发中子的反应堆的反应性ρp的严格定义：
令φ*＝φo+Δφ,φo为临界状态下的中子通量密度。代入式（28），注意到式（28）的每个含φ*的子式中，都有另一个无穷小量，略去所有二阶小量，可以得到反应性ρp近似定义：
3.4.2瞬发动力学方程推导
首先可以将通量φ*(r,E，Ω,t)写成为幅因子n(t)和形状因子ψ(r,E,Ω,t)的乘积：
φ*(r,E,Ω,t)＝n(t)·ψ(r,E,Ω,t) （30）
当系统从临界状态受到轻微的扰动时，形状因子随时间变化很小，有：υ为中子速度。
这样就有：
中子输运方程：
用φo+(r,E,Ω)乘（32）式，用φ*(r,E，Ω,t)乘（26）式，并对整个相空间(r,E,Ω)积分，相减后可得：
经整理，得到：
由中子代时间定义（23）和反应性定义式（29），可得瞬发动力学方程：
在与临界状态差别不大的情况下，形状因子变化很小，有：
ψ(r,E,Ω,t)≈ψo(r,E,Ω) （36）
那么由式（32）和（37）就有：
可得：
令：
可得：
4、计算反应堆的重要参数，判断是否满足设计标准
由计算出的参数βeff和Λ，选择一定的动态分析模型，如在西安脉冲堆的计算中选择了Nordheim‑Fuchs模型（该模型将在实施例中详述，值得一提的是，无论采用何种动态分析模型，βeff和Λ都是所必须的参数，本发明也因此具有很强的实用性），计算反应堆重要参数包括功率、最高温度等随时间的变化曲线，并由计算结果判断反应堆是否符合设计标准；如果反应堆不符合设计标准，则通过调整堆芯水铀比、堆芯体积、反射层材料和厚度来调整βeff和Λ，重复步骤2和3，直至根据βeff和Λ值计算出的重要参数符合反应堆的设计标准。
实施例：
以西安脉冲堆为例，介绍本发明的应用。
第一步，使用探测器测量各缓发中子先驱核的半衰期和中子谱等参数，拟合得到六组缓发中子的产额和相应缓发中子先驱核的衰变常数，并将缓发中子产额、中子谱、衰变常数等参数按一定格式写入数据库。
第二步，采用未经修改的原始数据库，计算出反应堆的有效增殖因子k(0)。再将数据库中的裂变核素235U、238U的缓发中子产额提高和降低0.25倍，使用蒙卡程序调用相应数据库，计算出反应堆有效增殖因子k(0.25)和k(‑0.25)，再采用公式k′(0)＝2(k(0.25)‑k(‑0.25))计算出k′(0)，然后根据微扰公式计算出βeff。βeff计算结果稳定在7.267‰。
第三步，使用MCNP程序，计算出仅考虑瞬发中子的反应堆的有效增殖因子kp=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=0.02482以发射脉冲，分析反应堆在引入最大脉冲反应性时是否符合设计标准。脉冲棒在95.8ms内匀速弹出堆芯，动态分析模型采用Nordheim‑Fuchs模型。Nordheim‑Fuchs模型的核心方程：
式中，P(t)是功率，ρ(t)是反应性。
结合基本的传热原理，即可得到反应堆设计中几个非常重要的参数的变化曲线，附图4、5分别是功率和平均最高温度随着时间的变化曲线。考虑燃料芯块内部的功率分布和反应堆的功率峰因子，计算得到的燃料最热点的最高温度为828.4℃。
第五步，判断重要参数是否符合设计标准。西安脉冲堆的设计安全限值为：在包壳温度大于500℃时，燃料芯块最大温度限值为970℃，在包壳温度小于500℃时，燃料芯块最大温度限值为1150℃。因此，在引入最大反应性ρ0=0.02482时，西安脉冲反应堆是符合设计安全标准的。如果经计算，反应堆不符合设计标准，则需调整堆芯结构、材料等。具体可以通过调整堆芯水铀比、堆芯体积、反射层材料和厚度等来调整，然后再根据本发明的方法得到反应堆重要参数值，直到反应堆重要参数的值满足设计标准。