CN107122545A - 一种精确计算核反应堆内时空中子分布的方法 - Google Patents

Abstract

一种精确计算核反应堆时空中子分布的方法，包括如下步骤：1、使用一步法直接求解稳态中子输运方程获得稳态状态下的中子通量密度分布以及先驱核浓度；2、在0.25ms内求解时空中子输运方程，通过0.25时刻的中子通量密度与中子通量密度形状函数获得0.25ms的幅值，此幅值作为点堆计算的初始值；3、在大步长内使用一步法计算时空中子输运方程，获得预估解，在极小步长内求解点堆方程，用点堆的幅值去校正中子通量密度；本发明通过改进现有的预估校正准静态方法，在不显著增加计算时间的前提下，能够精确捕捉核反应堆在极短时间内中子的迅速变化，彻底克服传统预估校正准静态方法的缺陷，使其在高保真时空中子动力学计算中真正发挥作用。

Description

一种精确计算核反应堆内时空中子分布的方法

技术领域

本发明涉及核反应堆堆芯设计和核安全技术领域，具体涉及一种精确计算核反应堆内时空中子分布的方法。

背景技术

在核电厂的设计中，核反应堆的设计一直处于最核心的地位，而在核反应堆的设计中，精确模拟核反应堆在瞬态过程中时在不同时间点以及不同位置处的中子分布一直是一个巨大的挑战。目前数值模拟是核反应堆设计的最主要的手段，因此为了精确模拟核反应堆在瞬态过程中时在不同时间点以及不同位置处的中子分布，数值模拟计算的方法便显得尤为重要。

近年来随着计算机计算能力的大幅度提高，在数值模拟中采用高保真时空中子动力学模型已经成为国际的主流。高保真时空中子动力学能够保证理论模型和计算模型不失真，计算所得解为真解，计算分辨率和计算精度高，计算结果精确，能够达到高度保真的要求。要实现高保真时空中子动力学的计算，就必须采用一步法对核反应堆直接进行几何建模，采用多能群结构，通过精确求解中子输运方程，在计算过程中不引入几何近似以及传统两步法均匀化带来的误差。通过高保真时空中子动力学的模拟可大大提高核反应堆设计的精确性，大幅度减少实验量，提高反应堆的安全性和经济性，对核反应堆改进设计和新型反应堆设计均具有重大的理论价值和应用前景。

但巨额的计算时间成为了限制高保真时空中子动力学计算发展的最主要瓶颈。国际上目前主要流行通过预估校正准静态方法处理时空中子动力学中与时间相关的因素。在预估校正准静态方法中，通过在较大时间步长内以高计算代价直接求解时空中子输运方程获得预估的中子分布，此为预估，再在此求解输运方程的步长内以较小步长和较小计算代价求解点堆方程获得校正的中子分布的幅值，最后在输运计算的步长末通过校正的幅值去修正预估的中子分布的幅值，此为校正。但是预估校正准静态方法在处理极短时间内核反应堆内中子迅速变化的情况时会出现较大的误差，导致最终计算结果存在较大误差，因此极大地限制了预估校正准静态方法在高保真时空中子动力学中的应用。

因此，充分利用预估校正准静态极大节省计算时间的优势，并通过一定的手段克服预估校正准静态在处理极短时间内核反应堆内中子迅速变化时失准的缺陷，才能使预估校正准静态方法在高保真时空中子动力学的模拟中发挥真正的作用。

发明内容

为了克服上述现有技术的问题，解决在核反应堆功率急剧变化的情况下，预估校正准静态方法计算失准的问题，本发明的目的是提出一种精确计算核反应堆内时空中子分布的方法，该方法在现有预估校正准静态方法的基础上，通过改进现有的预估校正准静态方法，在不显著增加计算时间的前提下，能够精确捕捉核反应堆在极短时间内中子的迅速变化，彻底克服传统预估校正准静态方法的缺陷，使其在高保真时空中子动力学计算中真正发挥作用。

为了达到上述目的，本发明采用如下技术方案：

一种精确计算核反应堆时空中子分布的方法，包括如下步骤：

步骤1：进行核反应堆稳态计算，采用一步法直接进行输运计算，得到核反应堆处于稳态时每个能群、每个平源区的中子通量密度以及每个平源区内的每组缓发中子先驱核初始的浓度，具体包括如下步骤：

1)从截面文件中读取各个材料的原始多群宏观截面信息与动力学参数信息；

2)从输入卡片中读取核反应堆的几何信息与计算条件信息；

3)根据输入卡片中读取的几何信息进行直接的几何建模：首先根据输入卡片中的几何描述得到核反应堆的几何布置；其次根据核反应堆的几何布置建立特征线方法计算所需的边界条件以及内部特征线的长度信息，为中子输运计算提供模块化特征线信息；

4)根据1)、2)、3)中得到的信息采用特征线方法进行中子输运计算，得到各个平源区的中子通量密度分布，具体的计算公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群编号

g'——非当前能群编号

G——能群总数

r——空间位置

Σ_t,g(r)——r处第g群的宏观总截面

Σ_s,g'→g(r)——r处g'能群到g能群的散射截面

χ_g(r)——r处第g能群的裂变谱

νΣ_f,g——第g群的中子产生截面

——r,Ω处第g能群中子角通量密度

φ_g(r)——r处第g能群的中子标通量密度

φ_g'(r)——r处第g'能群的中子标通量密度

S_F(r)——r处裂变源

k_eff——输运计算得到的有效增殖因子

由此得到各个平源区的中子通量密度；

5)根据4)中计算所得到的各个平源区的中子通量密度以及1)中读取的动力学参数信息，得到稳态下的各组缓发中子先驱核密度，具体的公式如下：

式中：

r——空间位置

k——缓发中子先驱核编号

g——能群编号

G——能群总数

k_eff——输运计算得到的有效增殖因子

C_k(r)——r处临界状态第k组缓发中子先驱核密度

β_k(r)——r处第k组缓发中子份额

νΣ_f,g——第g群的中子产生截面

λ_k(r)——r处第k组缓发中子先驱核的衰变常数

φ_g(r)——r处临界状态第g群的中子通量密度

S_F(r)——r处裂变源

6)根据1)、2)、3)得到的信息，对粗网有限差分CMFD方程进行中子共轭计算，得到各个粗网的共轭中子通量密度，具体的计算公式如下：

式中：

u——坐标轴方向标号

x,y,z——直角坐标系坐标轴方向

h_u——粗网在u方向的高度

——粗网在u方向右边界的净中子流

——粗网在u方向左边界的净中子流

g——当前能群标号

g'——非当前能群标号

G——总能群数

Σ_rg——第g能群的移出截面

Σ_s,g→g'——第g群到第g'群的散射截面

vΣ_f,g——第g能群的中子产生截面

χ_g'——第g'能群的裂变谱

k_eff——输运计算得到的有效增殖因子

φ_g——粗网第g能群的共轭中子通量密度

φ_g'——粗网第g'能群的共轭中子通量密度

7)将每个平源区的中子通量密度利用因子分解分解成为幅值函数与形状函数的乘积，此时幅值函数为1，利用4)求解得到的临界状态下的中子通量密度求得初始时刻的中子形状函数，具体的计算公式如下：

式中：

g——能群标号

r——空间位置

φ_g(r,0)——稳态时r处第g能群的中子通量密度

ψ_g(r,0)——稳态时r处第g能群的中子通量密度形状函数

N(0)——稳态时中子通量密度的幅值函数

步骤2：读取输入卡片扰动信息，根据得到的扰动信息执行截面扰动，打破核反应堆的稳态状态，从而开始核反应堆时空中子动力学计算；

步骤3：执行时空中子动力学计算的第一个时间步，以0.25ms为输运计算时间步长，进行中子输运形式的固定源计算，得到0.25ms时刻的中子通量密度以及中子通量密度形状函数，并由中子通量密度以及中子通量密度形状函数得到中子通量密度幅值，具体包括如下步骤：

1)以0.25ms为输运计算时间步长进行中子输运形式的固定源计算，得到0.25ms时刻的中子通量密度，具体公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群标号

g'——非当前能群标号

G——能群总数

r——空间位置

n——第n个输运计算步

k——缓发中子标号

——n时刻r处第g群的宏观总截面

——n时刻r处g'能群到g能群的散射截面

——n时刻r处第g能群的裂变谱

——n时刻第g群的中子产生截面

——n时刻r,Ω处第g能群中子角通量密度

——n时刻r处第g能群的中子标通量密度

——n时刻r处第g'能群的中子标通量密度

——n-1时刻r处第g能群的中子标通量密度

k_eff——输运计算得到的有效增殖因子

——n时刻r处第g能群的固定源

——n时刻r处的裂变源

A_g(r)——r处第g能群的固定源系数

B_g(r)——r处第g能群的固定源系数

C_g(r)——r处第g能群的固定源系数

v_g——第g能群的中子速度

Δt_n——第n个输运步的步长

——n时刻第g能群的等效缓发份额

χ_dk,g——第k组缓发中子在第g能群的缓发裂变谱

β_k——第k组缓发中子份额

——n-1时刻r处第g能群的等效缓发源项

2)计算归一化常数，归一化常数是求解得到0.25ms时刻处中子通量密度形状函数的必须要素，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

φ^*——粗网的共轭中子通量密度

v——中子速度

——t时刻的中子通量密度形状函数

——稳态下的中子通量密度形状函数

φ_i,g(0)——稳态下i平源区的第g能群的中子通量密度

C——归一化常数

3)利用1)和2)中计算得到的中子通量密度与归一化常数，得到0.25ms时刻的中子通量密度形状函数，具体计算公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

——t时刻i平源区的第g能群中子通量密度形状函数

φ_i,g(t)——t时刻i平源区的第g能群的中子通量密度

N(t)——t时刻的幅值函数

φ(t)——t时刻的中子通量密度

C——归一化常数

v——中子速度

φ^*——粗网的共轭中子通量密度

——t时刻中子通量密度的形状函数

<·>——对全能量、全空间积分

4)由1)计算得到的中子通量密度与3)中计算得到的中子通量密度形状函数得到0.25ms时刻的中子通量密度幅值，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

N(t)——t时刻的中子通量密度幅值函数

φ_i,g(t)——t时刻i平源区的第g能群的中子通量密度

——t时刻i平源区的第g能群的中子通量密度形状函数

步骤4：在t_{n-1,transport}至t_n,transport时间间隔内以远大于步骤3中0.25ms的时间步长进行求解输运形式的固定源方程，n＝1,N，其中N为整个计算过程划分的输运计算步数，求得t_n,transport时刻的中子通量密度，此步骤为预估校正准静态中的预估计算步，具体的计算公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群标号

g'——非当前能群标号

G——能群总数

r——空间位置

n——第n个输运计算步

k——缓发中子标号

——n时刻r处第g群的宏观总截面

——n时刻r处g'能群到g能群的散射截面

——n时刻r处第g能群的裂变谱

——n时刻第g群的中子产生截面

——n时刻r,Ω处第g能群中子角通量密度

——n时刻r处第g能群的中子标通量密度

——n时刻r处第g'能群的中子标通量密度

——n-1时刻r处第g能群的中子标通量密度

k_eff——输运计算得到的有效增殖因子

——n时刻r处第g能群的固定源

——n时刻r处的裂变源

A_g(r)——r处第g能群的固定源系数

B_g(r)——r处第g能群的固定源系数

C_g(r)——r处第g能群的固定源系数

v_g——第g能群的中子速度

Δt_n——第n个输运步的步长

——n时刻第g能群的等效缓发份额

χ_dk,g——第k组缓发中子在第g能群的缓发裂变谱

β_k——第k组缓发中子份额

——n-1时刻r处第g能群的等效缓发源项

步骤5：将步骤4中得到的中子通量密度通过因子分解分解为中子通量密度幅值函数与中子通量密度形状函数的乘积，在t_{n-1,transport}与t_n,transport两个时间点计算点堆参数，在t_{n-1,transport}与t_n,transport之间对点堆参数进行线性插值，在t_{n-1,transport}与t_n,transport之间求解点堆方程组，得到t_n,transport时刻的点堆幅值，利用t_n,transport时刻的点堆幅值对t_n,transport时刻的中子通量密度进行幅值的校正，此为预估校正准静态中的校正步，具体包含如下步骤：

1)将步骤4中得到的中子通量密度通过因子分解分解为中子通量密度幅值函数与中子通量密度形状函数的乘积，进而利用步骤1中求得的粗网的共轭中子通量密度，求解出t_n,transport的中子通量密度形状函数，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

——t时刻i平源区的第g能群的预估中子通量密度

N^p(t)——t时刻预估的中子通量密度幅值函数

N(t)——t时刻的幅值函数

——t时刻i平源区的第g能群的中子通量密度形状函数

C——归一化常数

v——中子速度

φ^*——粗网的共轭中子通量密度

——t时刻的中子通量密度形状函数

——稳态下的中子通量密度形状函数

φ_i,g(0)——稳态下i平源区的第g能群的中子通量密度

<·>——对全能量、全空间积分

2)分别计算t_{n-1,transport}与t_n,transport时刻的精确点堆参数，具体公式如下:

F(t)＝<φ^*(r)χ_g(r)S_F(r,t)>

式中：

t——时间

g——当前能群标号

g'——非当前能群标号

G——总能群数

r——空间位置

k——缓发中子标号

u——坐标轴方向标号

x,y,z——直角坐标系坐标轴方向

h_u——粗网在u方向的高度

——粗网第g群在u方向右边界的净中子流

——粗网第g群在u方向左边界的净中子流

Σ_rg——第g能群的移出截面

Σ_s,g→g'——第g群到第g'群的散射截面

vΣ_f,g——第g能群的中子产生截面

χ_g(r)——r处第g能群的裂变谱

k_eff——输运计算得到的有效增殖因子

φ^*(r)——r处粗网的共轭中子通量密度

——t时刻r处的中子通量密度形状函数

S_F(r,t)——t时刻r处的裂变源

k_eff——输运计算得到的有效增殖因子

F(t)——t时刻精确点堆参数分母

ρ(t)——t时刻的反应性

β_k(r)——r处第k组缓发中子份额

χ_dk,g(r)——r处第k组缓发中子在第g能群的缓发裂变谱

——t时刻第k组缓发中子份额

——t时刻缓发中子份额总和

v(r)——r处中子速度

Λ(t)——t时刻等效中子代时间

λ_k(t)——t时刻第k组缓发中子衰变常数

C_k(r,t)——t时刻r处临界状态第k组缓发中子先驱核密度

λ_k(r)——r处第k组缓发中子先驱核的衰变常数

c_k(t)——t时刻第k组缓发中子先驱核浓度

<·>——在全能量、全空间进行积分

3)在t_{n-1,transport}与t_n,transport之间以Δt_PK为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔t_n-1,PK至t_n,PK内用t_{n-1,transport}与t_n,transport时间点上的点堆参数进行插值，其中进而求解点堆方程组，得到t_n,PK时刻的幅值N^c，具体公式如下：

式中：

t——时间变量

i——缓发中子标号

N(t)——t时刻幅值

c_i(t)——t时刻第i组缓发中子先驱核浓度

ρ(t)——t时刻反应性

β_i(t)——t时刻第i组缓发中子份额

β(t)——t时刻缓发中子份额总和

Λ(t)——t时刻等效中子代时间

λ_i——第i组缓发中子衰变常数

4)当点堆计算的时间点达到t_n,transport时，校正t_n,transport时刻的中子通量密度，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间

——t时刻i平源区第g能群校正后的中子通量密度

N^c(t)——t时刻校正的幅值函数

N^p(t)——t时刻预估的幅值函数

——t时刻i平源区第g能群的中子通量密度形状函数

步骤6：重复执行步骤4与步骤5，直到时空中子动力学计算结束。

与现有技术相比，本发明具有如下突出优点：

1)彻底消除了传统预估校正准静态方法在处理核反应堆功率急剧变化时计算结果偏差较大的缺陷，在不显著增加计算时间的前提下大幅度提高计算精度，对预估校正准静态方法在高保真时空中子动力学计算中得到真正应用具有极重要作用；

2)能够实现变步长的计算策略，能够实现在核反应堆功率急剧变化的时间点采取小时间步长进行计算，在核反应堆功率变化平缓的时间点采取大时间步长进行计算，能够灵活充分应用计算机资源，提高计算的经济性。

附图说明

图1为堆芯稳态计算流程。

图2为输运修正点堆幅值流程图。

图3为点堆修正输运幅值流程图。

图4为C5G7‐TD基准题功率变化对比图。

图5为C5G7‐TD基准题误差对比图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细的说明：

步骤1：如图1所示，进行核反应堆稳态计算，采用一步法直接进行输运计算，得到核反应堆处于稳态时每个能群、每个平源区的中子通量密度以及每个平源区内的每组缓发中子先驱核初始的浓度，具体包括如下步骤：

1)从截面文件中读取各个材料的原始多群宏观截面信息与动力学参数信息；

2)从输入卡片中读取核反应堆的几何信息与计算条件信息；

3)根据输入卡片中读取的几何信息进行直接的几何建模：首先根据输入卡片中的几何描述得到核反应堆的几何布置；其次根据核反应堆的几何布置建立特征线方法计算所需的边界条件以及内部特征线的长度信息，为中子输运计算提供模块化特征线信息；

4)根据1)、2)、3)中得到的信息采用特征线方法进行中子输运计算，得到各个平源区的中子通量密度分布，具体的计算公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群编号

g'——非当前能群编号

G——能群总数

r——空间位置

Σ_t,g(r)——r处第g群的宏观总截面

Σ_s,g'→g(r)——r处g'能群到g能群的散射截面

χ_g(r)——r处第g能群的裂变谱

νΣ_f,g——第g群的中子产生截面

——r,Ω处第g能群中子角通量密度

φ_g(r)——r处第g能群的中子标通量密度

φ_g'(r)——r处第g'能群的中子标通量密度

S_F(r)——r处裂变源

k_eff——输运计算得到的有效增殖因子

由此得到各个平源区的中子通量密度；

5)根据4)中计算所得到的各个平源区的中子通量密度以及1)中读取的动力学参数信息，得到稳态下的各组缓发中子先驱核密度，具体的公式如下：

式中：

r——空间位置

k——缓发中子先驱核编号

g——能群编号

G——能群总数

k_eff——输运计算得到的有效增殖因子

C_k(r)——r处临界状态第k组缓发中子先驱核密度

β_k(r)——r处第k组缓发中子份额

νΣ_f,g——第g群的中子产生截面

λ_k(r)——r处第k组缓发中子先驱核的衰变常数

φ_g(r)——r处临界状态第g群的中子通量密度

S_F(r)——r处裂变源

6)根据1)、2)、3)得到的信息，对粗网有限差分CMFD方程进行中子共轭计算，得到各个粗网的共轭中子通量密度，具体的计算公式如下：

式中：

u——坐标轴方向标号

x,y,z——直角坐标系坐标轴方向

h_u——粗网在u方向的高度

——粗网在u方向右边界的净中子流

——粗网在u方向左边界的净中子流

g——当前能群标号

g'——非当前能群标号

G——总能群数

Σ_rg——第g能群的移出截面

Σ_s,g→g'——第g群到第g'群的散射截面

vΣ_f,g——第g能群的中子产生截面

χ_g'——第g'能群的裂变谱

k_eff——输运计算得到的有效增殖因子

φ_g——粗网第g能群的共轭中子通量密度

φ_g'——粗网第g'能群的共轭中子通量密度

7)将每个平源区的中子通量密度利用因子分解分解成为幅值函数与形状函数的乘积，此时幅值函数为1，利用4)求解得到的临界状态下的中子通量密度求得初始时刻的中子形状函数，具体的计算公式如下：

式中：

g——能群标号

r——空间位置

φ_g(r,0)——稳态时r处第g能群的中子通量密度

ψ_g(r,0)——稳态时r处第g能群的中子通量密度形状函数

N(0)——稳态时中子通量密度的幅值函数

步骤2：读取输入卡片扰动信息，根据得到的扰动信息执行截面扰动，打破核反应堆的稳态状态，从而开始核反应堆时空中子动力学计算；

步骤3：如图2所示，执行时空中子动力学计算的第一个时间步，以0.25ms为输运计算时间步长，进行中子输运形式的固定源计算，得到0.25ms时刻的中子通量密度以及中子通量密度形状函数，并由中子通量密度以及中子通量密度形状函数得到中子通量密度幅值，具体包括如下步骤：

1)以0.25ms为输运计算时间步长进行中子输运形式的固定源计算，得到0.25ms时刻的中子通量密度，具体公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群标号

g'——非当前能群标号

G——能群总数

r——空间位置

n——第n个输运计算步

k——缓发中子标号

——n时刻r处第g群的宏观总截面

——n时刻r处g'能群到g能群的散射截面

——n时刻r处第g能群的裂变谱

——n时刻第g群的中子产生截面

——n时刻r,Ω处第g能群中子角通量密度

——n时刻r处第g能群的中子标通量密度

——n时刻r处第g'能群的中子标通量密度

——n-1时刻r处第g能群的中子标通量密度

k_eff——输运计算得到的有效增殖因子

——n时刻r处第g能群的固定源

——n时刻r处的裂变源

A_g(r)——r处第g能群的固定源系数

B_g(r)——r处第g能群的固定源系数

C_g(r)——r处第g能群的固定源系数

v_g——第g能群的中子速度

Δt_n——第n个输运步的步长

——n时刻第g能群的等效缓发份额

χ_dk,g——第k组缓发中子在第g能群的缓发裂变谱

β_k——第k组缓发中子份额

——n-1时刻r处第g能群的等效缓发源项

2)计算归一化常数，归一化常数是求解得到0.25ms时刻处中子通量密度形状函数的必须要素，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

φ^*——粗网的共轭中子通量密度

v——中子速度

——t时刻的中子通量密度形状函数

——稳态下的中子通量密度形状函数

φ_i,g(0)——稳态下i平源区的第g能群的中子通量密度

C——归一化常数

3)利用1)和2)中计算得到的中子通量密度与归一化常数，得到0.25ms时刻的中子通量密度形状函数，具体计算公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

——t时刻i平源区的第g能群中子通量密度形状函数

φ_i,g(t)——t时刻i平源区的第g能群的中子通量密度

N(t)——t时刻的幅值函数

φ(t)——t时刻的中子通量密度

C——归一化常数

v——中子速度

φ^*——粗网的共轭中子通量密度

——t时刻中子通量密度的形状函数

<·>——对全能量、全空间积分

4)由1)计算得到的中子通量密度与3)中计算得到的中子通量密度形状函数得到0.25ms时刻的中子通量密度幅值，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

N(t)——t时刻的中子通量密度幅值函数

φ_i,g(t)——t时刻i平源区的第g能群的中子通量密度

——t时刻i平源区的第g能群的中子通量密度形状函数

步骤4：如图3所示，在t_{n-1,transport}至t_n,transport时间间隔内以远大于步骤3中0.25ms的时间步长进行求解输运形式的固定源方程，n＝1,N，其中N为整个计算过程划分的输运计算步数，求得t_n,transport时刻的中子通量密度，此步骤为预估校正准静态中的预估计算步，具体的计算公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群标号

g'——非当前能群标号

G——能群总数

r——空间位置

n——第n个输运计算步

k——缓发中子标号

——n时刻r处第g群的宏观总截面

——n时刻r处g'能群到g能群的散射截面

——n时刻r处第g能群的裂变谱

——n时刻第g群的中子产生截面

——n时刻r,Ω处第g能群中子角通量密度

——n时刻r处第g能群的中子标通量密度

——n时刻r处第g'能群的中子标通量密度

——n-1时刻r处第g能群的中子标通量密度

k_eff——输运计算得到的有效增殖因子

——n时刻r处第g能群的固定源

——n时刻r处的裂变源

A_g(r)——r处第g能群的固定源系数

B_g(r)——r处第g能群的固定源系数

C_g(r)——r处第g能群的固定源系数

v_g——第g能群的中子速度

Δt_n——第n个输运步的步长

——n时刻第g能群的等效缓发份额

χ_dk,g——第k组缓发中子在第g能群的缓发裂变谱

β_k——第k组缓发中子份额

——n-1时刻r处第g能群的等效缓发源项

步骤5：将步骤4中得到的中子通量密度通过因子分解分解为中子通量密度幅值函数与中子通量密度形状函数的乘积，在t_{n-1,transport}与t_n,transport两个时间点计算点堆参数，在t_{n-1,transport}与t_n,transport之间对点堆参数进行线性插值，在t_{n-1,transport}与t_n,transport之间求解点堆方程组，得到t_n,transport时刻的点堆幅值，利用t_n,transport时刻的点堆幅值对t_n,transport时刻的中子通量密度进行幅值的校正，此为预估校正准静态中的校正步，具体包含如下步骤：

1)将步骤4中得到的中子通量密度通过因子分解分解为中子通量密度幅值函数与中子通量密度形状函数的乘积，进而利用步骤1中求得的粗网的共轭中子通量密度，求解出t_n,transport的中子通量密度形状函数，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

——t时刻i平源区的第g能群的预估中子通量密度

N^p(t)——t时刻预估的中子通量密度幅值函数

N(t)——t时刻的幅值函数

——t时刻i平源区的第g能群的中子通量密度形状函数

C——归一化常数

v——中子速度

φ^*——粗网的共轭中子通量密度

——t时刻的中子通量密度形状函数

——稳态下的中子通量密度形状函数

φ_i,g(0)——稳态下i平源区的第g能群的中子通量密度

<·>——对全能量、全空间积分

2)分别计算t_{n-1,transport}与t_n,transport时刻的精确点堆参数，具体公式如下:

F(t)＝<φ^*(r)χ_g(r)S_F(r,t)>

式中：

t——时间

g——当前能群标号

g'——非当前能群标号

G——总能群数

r——空间位置

k——缓发中子标号

u——坐标轴方向标号

x,y,z——直角坐标系坐标轴方向

h_u——粗网在u方向的高度

——粗网第g群在u方向右边界的净中子流

——粗网第g群在u方向左边界的净中子流

Σ_rg——第g能群的移出截面

Σ_s,g→g'——第g群到第g'群的散射截面

vΣ_f,g——第g能群的中子产生截面

χ_g(r)——r处第g能群的裂变谱

k_eff——输运计算得到的有效增殖因子

φ^*(r)——r处粗网的共轭中子通量密度

——t时刻r处的中子通量密度形状函数

S_F(r,t)——t时刻r处的裂变源

k_eff——输运计算得到的有效增殖因子

F(t)——t时刻精确点堆参数分母

ρ(t)——t时刻的反应性

β_k(r)——r处第k组缓发中子份额

χ_dk,g(r)——r处第k组缓发中子在第g能群的缓发裂变谱

——t时刻第k组缓发中子份额

——t时刻缓发中子份额总和

v(r)——r处中子速度

Λ(t)——t时刻等效中子代时间

λ_k(t)——t时刻第k组缓发中子衰变常数

C_k(r,t)——t时刻r处临界状态第k组缓发中子先驱核密度

λ_k(r)——r处第k组缓发中子先驱核的衰变常数

c_k(t)——t时刻第k组缓发中子先驱核浓度

<·>——在全能量、全空间进行积分

3)在t_{n-1,transport}与t_n,transport之间以Δt_PK为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔t_n-1,PK至t_n,PK内用t_{n-1,transport}与t_n,transport时间点上的点堆参数进行插值，其中进而求解点堆方程组，得到t_n,PK时刻的幅值N^c，具体公式如下：

式中：

t——时间变量

i——缓发中子标号

N(t)——t时刻幅值

c_i(t)——t时刻第i组缓发中子先驱核浓度

ρ(t)——t时刻反应性

β_i(t)——t时刻第i组缓发中子份额

β(t)——t时刻缓发中子份额总和

Λ(t)——t时刻等效中子代时间

λ_i——第i组缓发中子衰变常数

4)当点堆计算的时间点达到t_n,transport时，校正t_n,transport时刻的中子通量密度，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间

——t时刻i平源区第g能群校正后的中子通量密度

N^c(t)——t时刻校正的幅值函数

N^p(t)——t时刻预估的幅值函数

——t时刻i平源区第g能群的中子通量密度形状函数

步骤6：重复执行步骤4与步骤5，直到时空中子动力学计算结束。

下面以时空中子动力学基准题C5G7‐TD系列的TD0‐5的计算结果为例说明本发明的效果：C5G7‐TD基准题系列的TD0‐5是一个二维七群非均匀的时空中子动力学基准题，基准解为以0.25ms为步长的全隐式向后差分计算所得，以2ms为时间步长进行传统预估校正准静态计算与本发明修正的预估校正准静态计算，总功率变化对比如图4所示，传统预估校正准静态与本发明的修正预估校正准静态和基准解的误差如图5所示。从计算结果可以看出，本发明在基本不增加计算时间的前提下极大地改善了传统预估校正准静态的精度，此计算结果证明本发明具有相当高的实用价值以及创新性。

Claims (1)

1.一种精确计算核反应堆时空中子分布的方法，其特征在于：包括如下步骤：

步骤1：进行核反应堆稳态计算，采用一步法直接进行输运计算，得到核反应堆处于稳态时每个能群、每个平源区的中子通量密度以及每个平源区内的每组缓发中子先驱核初始的浓度，具体包括如下步骤：

1)从截面文件中读取各个材料的原始多群宏观截面信息与动力学参数信息；

2)从输入卡片中读取核反应堆的几何信息与计算条件信息；

3)根据输入卡片中读取的几何信息进行直接的几何建模：首先根据输入卡片中的几何描述得到核反应堆的几何布置；其次根据核反应堆的几何布置建立特征线方法计算所需的边界条件以及内部特征线的长度信息，为中子输运计算提供模块化特征线信息；

4)根据1)、2)、3)中得到的信息采用特征线方法进行中子输运计算，得到各个平源区的中子通量密度分布，具体的计算公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群编号

g'——非当前能群编号

G——能群总数

r——空间位置

Σ_t,g(r)——r处第g群的宏观总截面

Σ_s,g'→g(r)——r处g'能群到g能群的散射截面

χ_g(r)——r处第g能群的裂变谱

νΣ_f,g——第g群的中子产生截面

——r,Ω处第g能群中子角通量密度

φ_g(r)——r处第g能群的中子标通量密度

φ_g'(r)——r处第g'能群的中子标通量密度

S_F(r)——r处裂变源

k_eff——输运计算得到的有效增殖因子

由此得到各个平源区的中子通量密度；

5)根据4)中计算所得到的各个平源区的中子通量密度以及1)中读取的动力学参数信息，得到稳态下的各组缓发中子先驱核密度，具体的公式如下：

式中：

r——空间位置

k——缓发中子先驱核编号

g——能群编号

G——能群总数

k_eff——输运计算得到的有效增殖因子

C_k(r)——r处临界状态第k组缓发中子先驱核密度

β_k(r)——r处第k组缓发中子份额

νΣ_f,g——第g群的中子产生截面

λ_k(r)——r处第k组缓发中子先驱核的衰变常数

φ_g(r)——r处临界状态第g群的中子通量密度

S_F(r)——r处裂变源

6)根据1)、2)、3)得到的信息，对粗网有限差分CMFD方程进行中子共轭计算，得到各个粗网的共轭中子通量密度，具体的计算公式如下：

式中：

u——坐标轴方向标号

x,y,z——直角坐标系坐标轴方向

h_u——粗网在u方向的高度

——粗网在u方向右边界的净中子流

——粗网在u方向左边界的净中子流

g——当前能群标号

g'——非当前能群标号

G——总能群数

Σ_rg——第g能群的移出截面

Σ_s,g→g'——第g群到第g'群的散射截面

vΣ_f,g——第g能群的中子产生截面

χ_g'——第g'能群的裂变谱

k_eff——输运计算得到的有效增殖因子

φ_g——粗网第g能群的共轭中子通量密度

φ_g'——粗网第g'能群的共轭中子通量密度

7)将每个平源区的中子通量密度利用因子分解分解成为幅值函数与形状函数的乘积，此时幅值函数为1，利用4)求解得到的临界状态下的中子通量密度求得初始时刻的中子形状函数，具体的计算公式如下：

式中：

g——能群标号

r——空间位置

φ_g(r,0)——稳态时r处第g能群的中子通量密度

ψ_g(r,0)——稳态时r处第g能群的中子通量密度形状函数

N(0)——稳态时中子通量密度的幅值函数

步骤2：读取输入卡片扰动信息，根据得到的扰动信息执行截面扰动，打破核反应堆的稳态状态，从而开始核反应堆时空中子动力学计算；

步骤3：执行时空中子动力学计算的第一个时间步，以0.25ms为输运计算时间步长，进行中子输运形式的固定源计算，得到0.25ms时刻的中子通量密度以及中子通量密度形状函数，并由中子通量密度以及中子通量密度形状函数得到中子通量密度幅值，具体包括如下步骤：

1)以0.25ms为输运计算时间步长进行中子输运形式的固定源计算，得到0.25ms时刻的中子通量密度，具体公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群标号

g'——非当前能群标号

G——能群总数

r——空间位置

n——第n个输运计算步

k——缓发中子标号

——n时刻r处第g群的宏观总截面

——n时刻r处g'能群到g能群的散射截面

——n时刻r处第g能群的裂变谱

——n时刻第g群的中子产生截面

——n时刻r,Ω处第g能群中子角通量密度

——n时刻r处第g能群的中子标通量密度

——n时刻r处第g'能群的中子标通量密度

——n-1时刻r处第g能群的中子标通量密度

k_eff——输运计算得到的有效增殖因子

——n时刻r处第g能群的固定源

——n时刻r处的裂变源

A_g(r)——r处第g能群的固定源系数

B_g(r)——r处第g能群的固定源系数

C_g(r)——r处第g能群的固定源系数

v_g——第g能群的中子速度

Δt_n——第n个输运步的步长

——n时刻第g能群的等效缓发份额

χ_dk,g——第k组缓发中子在第g能群的缓发裂变谱

β_k——第k组缓发中子份额

——n-1时刻r处第g能群的等效缓发源项

2)计算归一化常数，归一化常数是求解得到0.25ms时刻处中子通量密度形状函数的必须要素，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

φ^*——粗网的共轭中子通量密度

v——中子速度

——t时刻的中子通量密度形状函数

——稳态下的中子通量密度形状函数

φ_i,g(0)——稳态下i平源区的第g能群的中子通量密度

C——归一化常数

3)利用1)和2)中计算得到的中子通量密度与归一化常数，得到0.25ms时刻的中子通量密度形状函数，具体计算公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

——t时刻i平源区的第g能群中子通量密度形状函数

φ_i,g(t)——t时刻i平源区的第g能群的中子通量密度

N(t)——t时刻的幅值函数

φ(t)——t时刻的中子通量密度

C——归一化常数

v——中子速度

φ*——粗网的共轭中子通量密度

——t时刻中子通量密度的形状函数

<·>——对全能量、全空间积分

4)由1)计算得到的中子通量密度与3)中计算得到的中子通量密度形状函数得到0.25ms时刻的中子通量密度幅值，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

N(t)——t时刻的中子通量密度幅值函数

φ_i,g(t)——t时刻i平源区的第g能群的中子通量密度

——t时刻i平源区的第g能群的中子通量密度形状函数

步骤4：在t_{n-1,transport}至t_n,transport时间间隔内以远大于步骤3中0.25ms的时间步长进行求解输运形式的固定源方程，n＝1,N，其中N为整个计算过程划分的输运计算步数，求得t_n,transport时刻的中子通量密度，此步骤为预估校正准静态中的预估计算步，具体的计算公式如下：

式中：

Ω——角度方向

——梯度算子

g——当前能群标号

g'——非当前能群标号

G——能群总数

r——空间位置

n——第n个输运计算步

k——缓发中子标号

——n时刻r处第g群的宏观总截面

——n时刻r处g'能群到g能群的散射截面

——n时刻r处第g能群的裂变谱

——n时刻第g群的中子产生截面

——n时刻r,Ω处第g能群中子角通量密度

——n时刻r处第g能群的中子标通量密度

——n时刻r处第g'能群的中子标通量密度

——n-1时刻r处第g能群的中子标通量密度

k_eff——输运计算得到的有效增殖因子

——n时刻r处第g能群的固定源

——n时刻r处的裂变源

A_g(r)——r处第g能群的固定源系数

B_g(r)——r处第g能群的固定源系数

C_g(r)——r处第g能群的固定源系数

v_g——第g能群的中子速度

Δt_n——第n个输运步的步长

——n时刻第g能群的等效缓发份额

χ_dk,g——第k组缓发中子在第g能群的缓发裂变谱

β_k——第k组缓发中子份额

——n-1时刻r处第g能群的等效缓发源项

步骤5：将步骤4中得到的中子通量密度通过因子分解分解为中子通量密度幅值函数与中子通量密度形状函数的乘积，在t_{n-1,transport}与t_n,transport两个时间点计算点堆参数，在t_{n-1,transport}与t_n,transport之间对点堆参数进行线性插值，在t_{n-1,transport}与t_n,transport之间求解点堆方程组，得到t_n,transport时刻的点堆幅值，利用t_n,transport时刻的点堆幅值对t_n,transport时刻的中子通量密度进行幅值的校正，此为预估校正准静态中的校正步，具体包含如下步骤：

1)将步骤4中得到的中子通量密度通过因子分解分解为中子通量密度幅值函数与中子通量密度形状函数的乘积，进而利用步骤1中求得的粗网的共轭中子通量密度，求解出t_n,transport的中子通量密度形状函数，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间变量

——t时刻i平源区的第g能群的预估中子通量密度

N^p(t)——t时刻预估的中子通量密度幅值函数

N(t)——t时刻的幅值函数

——t时刻i平源区的第g能群的中子通量密度形状函数

C——归一化常数

v——中子速度

φ*——粗网的共轭中子通量密度

——t时刻的中子通量密度形状函数

——稳态下的中子通量密度形状函数

φ_i,g(0)——稳态下i平源区的第g能群的中子通量密度

<·>——对全能量、全空间积分

2)分别计算t_{n-1,transport}与t_n,transport时刻的精确点堆参数，具体公式如下:

式中：

t——时间

g——当前能群标号

g'——非当前能群标号

G——总能群数

r——空间位置

k——缓发中子标号

u——坐标轴方向标号

x,y,z——直角坐标系坐标轴方向

h_u——粗网在u方向的高度

——粗网第g群在u方向右边界的净中子流

——粗网第g群在u方向左边界的净中子流

Σ_rg——第g能群的移出截面

Σ_s,g→g'——第g群到第g'群的散射截面

vΣ_f,g——第g能群的中子产生截面

χ_g(r)——r处第g能群的裂变谱

k_eff——输运计算得到的有效增殖因子

φ^*(r)——r处粗网的共轭中子通量密度

——t时刻r处的中子通量密度形状函数

S_F(r,t)——t时刻r处的裂变源

k_eff——输运计算得到的有效增殖因子

F(t)——t时刻精确点堆参数分母

ρ(t)——t时刻的反应性

β_k(r)——r处第k组缓发中子份额

χ_dk,g(r)——r处第k组缓发中子在第g能群的缓发裂变谱

——t时刻第k组缓发中子份额

——t时刻缓发中子份额总和

v(r)——r处中子速度

Λ(t)——t时刻等效中子代时间

λ_k(t)——t时刻第k组缓发中子衰变常数

C_k(r,t)——t时刻r处临界状态第k组缓发中子先驱核密度

λ_k(r)——r处第k组缓发中子先驱核的衰变常数

c_k(t)——t时刻第k组缓发中子先驱核浓度

<·>——在全能量、全空间进行积分

3)在t_{n-1,transport}与t_n,transport之间以Δt_PK为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔t_n-1,PK至t_n,PK内用t_{n-1,transport}与t_n,transport时间点上的点堆参数进行插值，其中进而求解点堆方程组，得到t_n,PK时刻的幅值N^c，具体公式如下：

式中：

t——时间变量

i——缓发中子标号

N(t)——t时刻幅值

c_i(t)——t时刻第i组缓发中子先驱核浓度

ρ(t)——t时刻反应性

β_i(t)——t时刻第i组缓发中子份额

β(t)——t时刻缓发中子份额总和

Λ(t)——t时刻等效中子代时间

λ_i——第i组缓发中子衰变常数

4)当点堆计算的时间点达到t_n,transport时，校正t_n,transport时刻的中子通量密度，具体公式如下：

式中：

i——平源区标号

g——能群标号

t——时间

——t时刻i平源区第g能群校正后的中子通量密度

N^c(t)——t时刻校正的幅值函数

N^p(t)——t时刻预估的幅值函数

——t时刻i平源区第g能群的中子通量密度形状函数

步骤6：重复执行步骤4与步骤5，直到时空中子动力学计算结束。