本发明涉及核反应堆设计和反应堆物理计算领域，具体涉及一种pin-by-pin分析压水堆瞬态的方法。
背景技术：
在反应堆运行过程中，需要实时掌控堆内的真实情况，以便于操作员可以及时进行反应性调节，避免事故的发生，或在事故发生后控制事故的范围和等级，保障工作人员与周围群众的生命财产安全。这需要研究一种实时且准确分析压水堆堆芯瞬态过程的方法，即中子动态学计算。堆芯与事件相关的现象有很多种，根据其过程时间常数的不同大致可以分成三类：短时间瞬态、中等时间瞬态和长时间瞬态。短时间瞬态持续时间通常为毫秒到秒量级，主要是由于人为地或者事故地引入反应性而导致的中子通量密度迅速改变，比如瞬态实验，启动，停堆过程等。中等时间瞬态通常为持续数小时或者数天的过程，主要是热中子堆中裂变产物(xe135，sm149等)的积累、消耗以及β衰变所导致的瞬态过程。长时间瞬态的时间跨度达到数个月甚至数年，主要是由于易裂变产物的积累和燃耗，和可裂变核素的积累、燃耗以及衰变。
用确定论方法直接求解中子输运方程是很耗时的，无法满足实时性要求，因此传统方法都是求解角度p1近似的中子扩散方程，空间上通常采用组件均匀化思想，先求得组件均匀化的中子通量密度，再通过功率重构得到组件内部的功率分布。但这样求得的中子通量密度和功率分布往往与实际的真实结果误差较大，特别是对于局部的非均匀效应。sp3近似方程是相对于中子扩散方程更加精确的中子输运方程形式，通过空间上pin-by-pin(以栅元为典型网格大小)求解sp3近似方程可以获得每个栅元内部的中子通量密度和功率分布。西安交通大学necp实验室已经开发堆芯程序efen(指数函数展开节块方法)可以求解sp3近似方程。
因此，关键在于如何pin-by-pin求解瞬态的sp3近似方程，得到每个栅元内部的中子学信息，实时反映和分析反应堆运行瞬态过程中堆芯内部的真实情况。
技术实现要素：
为了克服上述现有技术存在的问题，本发明的目的是为了给出一种pin-by-pin分析压水堆瞬态的方法，采用差分方法和积分方法将瞬态中子输运方程离散成稳态中子输运方程，再将它转化成稳态sp3近似方程，再在每个时间步下pin-by-pin迭代求解稳态方程，得到各个时间步的栅元内部的中子通量密度和功率分布，同时计算时间依旧满足实际运用需求。
为了达到上述目的，本发明采取了以下技术方案予以实施：
一种pin-by-pin分析压水堆瞬态的方法，步骤如下：
步骤1：计算堆芯初始状态，作为瞬态分析的初始点；得到每个栅元内部的零阶和二阶中子通量密度和修正后的多群微观总截面，微观散射截面和微观裂变截面；具体包括如下步骤：
1)读取输入卡片中堆芯几何参数进行几何建模，读取核素的原始宏观总截面，宏观散射截面，缓发中子份额，衰变常数，能群信息，中子能谱，单次裂变产生中子数，读取瞬态时间步长信息和扰动信息；
2)基于输入参数和堆芯几何计算初始的缓发中子先驱核浓度、sp3稳态方程的系数和源项；计算公式如下所示：
式中:
g'——碰撞前中子能群；
t0——初始时间；
x——中子位置；
sf(t0)——初始时刻的裂变源；
ν(x)——单次裂变产生中子数；
σf,g'(x,t0)——初始时刻在x处的第g'群中子宏观裂变截面；
ν(x)σf,g'(x,t0)——初始时刻在x处的第g'群中子宏观产生截面；
式中：
i——缓发中子组数；
nd——缓发中子总组数；
ci(x,t0)——初始时刻在x处的第i组缓发中子先驱核浓度；
λi——第i组缓发中子衰变常数；
βi(x)——在x处的第i组缓发中子份额；
式中：
sd——缓发中子源项；
χd,g,i(x)——在x处的第i组缓发中子能谱；
式中：
d1——等效扩散系数；
σt,g(x)——在x处的第g群中子宏观总截面；
σr,g(x)＝σt,g(x)-σs,g→g(x)公式(5)
式中：
σr,g(x)——在x处的第g群中子宏观移出截面；
σs,g→g(x)——在x处的第g群中子宏观自散射截面；
式中：
sg——等效中子源；
μ0——初始角度；
g——碰撞后中子能群；
σs,g'→g(x,μ0)——在x处角度为μ0的中子从第g'群散射到第g群的宏观截面；
keff——有效增殖系数；
χg(x)——第g群中子能谱；
3)利用指数函数展开节块方法进行pin-by-pin稳态计算，得到初始时间步每个栅元内部的零阶和二阶中子通量密度和有效增值系数；sp3稳态方程如下所示：
式中：
4)利用有效增值系数修正原始截面，得到更新后的截面参数,如公式所示：
步骤2：对下一个时间步，考虑外部扰动信息，同时利用上一个时间步的中子通量密度更新截面；得到sp3稳态方程的系数和源项；
式中：
式中：
vg(x)——在x处的第g群中子速度；
δtn＝tn+1-tn——第n个到第n+1个时间步的步长；
sf(tn+1)——第n+1个时间步在x处的裂变源；
ci(x,tn)——tn时刻在x处的第i组缓发中子先驱核浓度；
式中：
γ——第n个时间步长与第n-1个时间步长的比值；
δtn-1＝tn-tn-1——第n-1个到第n个时间步的步长；
κ0(λi)——自定义函数；
κ1(λi)——自定义函数；
κ2(λi)——自定义函数；
步骤3：利用指数函数展开节块方法空间上pin-by-pin数值求解时间离散的第tn+1步的稳态sp3方程，迭代求解收敛后得到新时间步下栅元的零阶和二阶中子通量密度；稳态sp3方程如下所示：
式中：
堆芯的有效增殖系数和通量在相邻两次时间步的结果误差小于
用户给定收敛限即认为收敛；收敛准则公式如下所示：
式中：
errork——有效增值系数误差；
keff(tn+1)——在tn+1时刻计算得到的有效增值系数；
keff(tn)——在tn时刻计算得到的有效增值系数；
errorf0——0阶中子通量密度矩误差；
errorf2——2阶中子通量密度矩误差；
步骤4：更新缓发中子先驱核浓度，公式如下所示：
公式(15)
式中：
ci(x,tn+1)——tn+1时刻在x处的第i组缓发中子先驱核浓度；
sf(tn)——第n个时间步在x处的裂变源；
sf(tn-1)——第n个时间步在x处的裂变源；
步骤5：判断当前时间步是否为最后时间步，若不是，则重复2-4；若是，则停止计算，输出瞬态参数；为反应堆操作员的及时操作提供及时准确的反应堆运行瞬态信息。
本发明与现有技术相比，具有如下优点：
1.空间上采用pin-by-pin求解，可以得到每个栅元的中子通量密度和功率分布，相对于组件均匀化方法，不需要考虑功率重构引入误差。
2.角度上采用sp3近似，可以得到零阶和二阶中子通量密度，相对于扩散近似只能得到零阶中子通量密度，保留了更多的堆芯信息。同时求解时间并没有显著增加。
附图说明
图1为pin-by-pin分析压水堆瞬态的操作流程图。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步详细说明：
本发明采用差分方法和积分方法将瞬态中子输运方程离散成稳态中子输运方程，再将它转化成稳态sp3近似方程，再在每个时间步下pin-by-pin迭代求解稳态方程，得到各个时间步的栅元内部的中子通量密度和功率分布，同时计算时间依旧满足实际运用需求。
为了达到上述目的，本发明采取了以下技术方案予以实施：
一种pin-by-pin分析压水堆瞬态的方法，步骤如下：
步骤1：计算堆芯初始状态，作为瞬态分析的初始点；得到每个栅元内部的零阶和二阶中子通量密度和修正后的多群微观总截面，微观散射截面和微观裂变截面；具体包括如下步骤：
1)读取输入卡片中堆芯几何参数进行几何建模，读取核素的原始宏观总截面，宏观散射截面，缓发中子份额，衰变常数，能群信息，中子能谱，单次裂变产生中子数，读取瞬态时间步长信息和扰动信息；
2)基于输入参数和堆芯几何计算初始的缓发中子先驱核浓度、sp3稳态方程的系数和源项；计算公式如下所示：
式中:
g'——碰撞前中子能群；
t0——初始时间；
x——中子位置；
sf(t0)——初始时刻的裂变源；
ν(x)——单次裂变产生中子数；
σf,g'(x,t0)——初始时刻在x处的第g'群中子宏观裂变截面；
ν(x)σf,g'(x,t0)——初始时刻在x处的第g'群中子宏观产生截面；
式中：
i——缓发中子组数；
nd——缓发中子总组数；
ci(x,t0)——初始时刻在x处的第i组缓发中子先驱核浓度；
λi——第i组缓发中子衰变常数；
βi(x)——在x处的第i组缓发中子份额；
式中：
sd——缓发中子源项；
χd,g,i(x)——在x处的第i组缓发中子能谱；
式中：
d1——等效扩散系数；
σt,g(x)——在x处的第g群中子宏观总截面；
σr,g(x)＝σt,g(x)-σs,g→g(x)公式(5)
式中：
σr,g(x)——在x处的第g群中子宏观移出截面；
σs,g→g(x)——在x处的第g群中子宏观自散射截面；
式中：
sg——等效中子源；
μ0——初始角度；
g——碰撞后中子能群；
σs,g'→g(x,μ0)——在x处角度为μ0的中子从第g'群散射到第g群的宏观截面；
keff——有效增殖系数；
χg(x)——第g群中子能谱；
3)利用efen程序(指数函数展开节块方法)进行pin-by-pin稳态计算，得到初始时间步每个栅元内部的零阶和二阶中子通量密度和有效增值系数；sp3稳态方程如下所示：
式中：
4)利用有效增值系数修正原始截面，得到更新后的截面参数,如公式
所示：
步骤2：对下一个时间步，考虑外部扰动信息，同时利用上一个时间步的中子通量密度更新截面；得到sp3稳态方程的系数和源项；
式中：
式中：
vg(x)——在x处的第g群中子速度；
δtn＝tn+1-tn——第n个到第n+1个时间步的步长；
sf(tn+1)——第n+1个时间步在x处的裂变源；
ci(x,tn)——tn时刻在x处的第i组缓发中子先驱核浓度；
式中：
γ——第n个时间步长与第n-1个时间步长的比值；
δtn-1＝tn-tn-1——第n-1个到第n个时间步的步长；
κ0(λi)——自定义函数；
κ1(λi)——自定义函数；
κ2(λi)——自定义函数；
步骤3：利用指数函数展开节块方法(efen程序)空间上pin-by-pin数值求解时间离散的第tn+1步的稳态sp3方程，迭代求解收敛后得到新时间步下栅元的零阶和二阶中子通量密度；稳态sp3方程如下所示：
式中：
堆芯的有效增殖系数和通量在相邻两次时间步的结果误差小于用户给定收敛限即可认为收敛。收敛准则公式如下所示：
式中：
errork——有效增值系数误差；
keff(tn+1)——在tn+1时刻计算得到的有效增值系数；
keff(tn)——在tn时刻计算得到的有效增值系数；
errorf0——0阶中子通量密度矩误差；
errorf2——2阶中子通量密度矩误差；
步骤4：更新缓发中子先驱核浓度，公式如下所示：
式中：
ci(x,tn+1)——tn+1时刻在x处的第i组缓发中子先驱核浓度；
sf(tn)——第n个时间步在x处的裂变源；
sf(tn-1)——第n个时间步在x处的裂变源；
步骤5：判断当前时间步是否为最后时间步，若不是，则重复2-4；若是，则停止计算，输出瞬态参数；为反应堆操作员的及时操作提供及时准确的反应堆运行瞬态信息。