CN107145731A

CN107145731A - 一种高保真时空中子动力学计算的加速方法

Info

Publication number: CN107145731A
Application number: CN201710290790.5A
Authority: CN
Inventors: 刘宙宇; 王博; 曹良志; 吴宏春
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2017-04-27
Filing date: 2017-04-27
Publication date: 2017-09-08
Anticipated expiration: 2037-04-27
Also published as: CN107145731B

Abstract

一种高保真时空中子动力学计算的加速方法，包括如下步骤：1、使用一步法计算反应堆临界状态下的中子通量密度与先驱核浓度；2、执行截面扰动，开始时空中子动力学计算；3、在大时间步长内使用一步法进行输运计算，获得预估的中子通量密度；4、在中等时间步内进行CMFD计算，同时采用全局加速因子方法进行加速，在最小时间步内进行点堆计算，使用点堆的幅值去校正粗网通量的幅值，在使用更新后的粗网的幅值去更新细网的中子通量密度；本发明方法能在保证计算精度的前提下，能够大幅度减少高保真时空中子动力学计算的时间步数，同时减少在每个时间步上的计算时间，进而大幅缩短整个计算过程的计算耗时。

Description

一种高保真时空中子动力学计算的加速方法

技术领域

本发明涉及核反应堆堆芯设计和核安全技术领域，具体涉及一种高保真时空中子动力学计算的加速方法。

背景技术

压水堆核电厂的瞬态和事故工况分析长期基于不考虑空间影响的点堆模型，所谓点堆模型即认为整个核反应堆的功率变化与空间位置无关，显然点堆模型并不适用于功率分布畸变，如弹棒事故、落棒事故、主蒸汽管道破口事故等。近年来随着计算机技术的迅速发展，在点堆模型基础上发展了三维时空中子动力学计算方法，考虑了功率变化在三维空间上的具体分布，在数值模拟的精准性上跨上了一个台阶。但其所采用的计算方法仍是基于传统的两步均匀化方法，无法排除均匀化带来的误差，且不能直接考虑核反应堆内复杂的几何结构，无法直接获得核反应堆内精细的功率分布。

所谓两步均匀化方法，即第一步对构成核反应堆的每一种组件进行单独的计算，第二步对以组件为单位的核反应堆进行计算。在第一步计算中利用全反射边界条件单独地对各种组件进行多能群的中子输运计算，并在组件处于不同的工况下分别进行计算，得到每种组件在不同工况下的多能群的中子通量密度在组件内的分布，进而利用核反应率守恒归并出两能群的、不同工况下、在每种组件内跟空间位置无关的等效均匀化群常数以及不连续因子等物理量；在第二步计算中根据核反应堆所处的工况，对每种组件的等效均匀化群常数利用多项式拟合或多项式插值来获得当前计算工况下的截面参数等，在获得每个组件的截面参数后对核反应堆进行两能群的中子扩散计算，并对各个组件进行功率重构计算，得到三维核反应堆堆芯的功率分布。

在两步均匀化方法计算中只考虑有限种独立的组件计算来计算等效均匀化参数，因此它不能精确考虑在核反应堆中组件之间的相互影响，并且不能精确考虑核反应堆内复杂的几何结构。随着对计算精度要求的不断提高，只有高保真时空中子动力学计算才能根本解决两步均匀化法的问题。

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

但是高保真时空中子动力学的应用存在一个非常明显的缺陷，即计算时间过长，计算资源消耗过大。因此对高保真时空中子动力学计算的加速在高保真时空中子动力学计算中起着十分关键的作用，采取高效的加速方法是高保真时空中子动力学计算能够得到实际应用的关键要素之一。

发明内容

为了克服上述现有技术的问题，解决高保真时空中子动力学计算存在的计算时间过长的问题，本发明的目的是提出一种对高保真时空中子动力学计算加速的方法，该方法能在保证计算精度的前提下，能够大幅度减少高保真时空中子动力学计算的时间步数，同时减少在每个时间步上的计算时间，进而大幅缩短整个计算过程的计算耗时。

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

一种高保真时空中子动力学计算的加速方法，包括如下步骤：

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

1)从截面库中读取各个核素的原始多群宏观截面与动力学参数信息；

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

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

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

式中：

Ω——角度方向

——梯度算子

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)中计算所得到的各个平源区的中子通量密度，得到临界状态下的各组缓发中子先驱核密度，具体的公式如下：

式中：

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)将每个平源区的中子通量密度利用因子分解分解成为幅值函数与形状函数的乘积，其中幅值函数为平源区对应粗网的中子通量密度，此时的幅值函数是预估的值，存在较大误差，需要后续进一步校正，利用4)求解得到的临界状态下的中子通量密度求得初始时刻的中子形状函数，具体的计算公式如下：

式中：

g——能群标号

r——空间位置

t——时间

i——粗网编号

——t时刻r处第g能群预估的中子通量密度

——t时刻i粗网处第g能群预估的中子通量密度

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

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

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

步骤2：根据输入卡片的描述执行截面扰动，打破核反应堆的临界状态，从而开始时空中子动力学计算；

步骤3：在t_{n-1,transport}至t_n,transport时间间隔内求解输运形式的固定源方程，n＝1,N，其中N为整个计算过程划分的输运计算步数，并采用CMFD方法对输运方程进行加速，同时使用全局加速因子对CMFD方程的求解进行加速，求得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能群的等效缓发源项

步骤4：利用步骤1和步骤2中求得的平源区的中子通量密度、粗网中子通量密度以及粗网共轭中子通量密度信息，进行校正计算，得到校正后的t_n时刻的每个平源区中子标通量密度，具体过程如下：

1)在t_{n-1,transport}至t_n,transport时刻的时间间隔内根据设定的CMFD计算时间步长Δt_CMFD划分CMFD计算所需的时间间隔，使用t_{n-1,transport}时刻的CMFD参数和t_n,transport时刻的CMFD参数在计算间隔t_n-1,CMFD至t_n,CMFD间进行线性插值CMFD参数，其中n＝1,以Δt_CMFD为时间步长求解CMFD固定源方程，同时对CMFD固定源方程的求解应用全局加速因子，得到t_n,CMFD时刻粗网的中子通量密度具体公式如下：

式中：

u——坐标轴方向标号

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

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

i——粗网标号

j——平源区标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网内第g能群的固定源项

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

Iⁱ——i粗网的平源区总数

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

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

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网j平源区第g能群的中子通量密度

S_F ^i,j——i粗网j平源区的裂变源

V^i,j——i粗网内j平源区的体积

采用全局加速因子加速上述CMFD形式的固定源方程，具体如下：

式中：

f——全局加速因子

S——固定源项

Lφ⁽ⁿ⁾——消失项

Mφ⁽ⁿ⁾——产生项

<·>操作表示对全能量、全空间进行积分

φ⁽ⁿ⁾——n时刻的中子通量密度

u——坐标轴方向标号

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

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

i——粗网标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

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

——i粗网第g能群的固定源系数

φ⁽ⁿ⁾＝f·φ⁽ⁿ⁾ (公式4‐3)

式中：

φ⁽ⁿ⁾——第n步粗网的中子通量密度

f——全局加速因子

2)将t_n,CMFD时刻粗网的中子通量密度因子分解为幅值函数与形状函数的乘积，进而利用步骤1中求得的粗网的共轭中子通量密度和步骤1中求得的粗网的中子通量密度计算归一化常数C，从而计算t_n,CMFD时刻的粗网的中子通量密度形状函数，具体公式如下：

式中：

i——粗网标号

g——能群标号

t——时间变量

——t时刻i粗网的第g能群的预估中子通量密度

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

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

——t时刻i粗网的第g能群的形状函数

C——归一化常数

v——中子速度

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

——t时刻粗网的形状函数

——临界状态下粗网的形状函数

φ_i,g(0)——临界状态下i粗网的第g能群的中子通量密度

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

3)分别计算t_n-1,CMFD与t_n,CMFD时刻的精确点堆参数，具体公式如下：

式中：

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组缓发中子先驱核浓度

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

4)在t_n-1,CMFD与t_n,CMFD之间以Δt_PK为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔t_n-1,PK至t_n,PK内用t_n-1,CMFD与t_n,CMFD时间点上的点堆参数进行插值，其中n＝1,进而求解点堆方程组，得到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组缓发中子衰变常数

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

式中：

i——粗网标号

g——能群标号

t——时间

——t时刻i粗网第g能群校正后的中子通量密度

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

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

——t时刻i粗网第g能群粗网的形状函数

6)循环上述1)至5)过程，当t_n,CMFD＝t_n时，用t_n,CMFD时刻的粗网中子通量密度校正t_n,transport时刻的每个平源区的中子通量密度，具体公式如下：

式中：

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

——t时刻i粗网第g能群校正后的中子通量密度

——t时刻i粗网第g能群预估的中子通量密度

ψ_g(r,t)——t时刻r处平源区第g能群粗网的形状函数

步骤5：重复执行步骤3与步骤4，直到动力学计算结束为止。

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

1)采用一步法计算核反应堆中的时空中子动力学方程，不需要进行组件均匀化，可精确处理复杂的非均匀几何问题，大大提高计算精度，达到计算结果高度保真的效果，满足高保真计算的要求。

2)采用了多级的预估校正准静态策略与全局加速因子，能够在保证计算精度的前提下大大减少计算的时间步长，并大大减少每个时间步上的计算时间，显著提高计算效率，对高保真时空中子动力学计算应用到实际的核反应堆设计中去具有重要意义。

附图说明

图1为细网结构图。

图2为堆芯临界计算流程图。

图3为CMFD与全局加速因子加速输运固定源计算流程图。

图4为瞬态多级加速计算流程图。

图5为TWIGL基准题有无加速计算结果对比图。

具体实施方式

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

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

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

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

式中：

Ω——角度方向

——梯度算子

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——输运计算得到的有效增殖因子

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

式中：

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处裂变源

式中：

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'能群的共轭中子通量密度

式中：

g——能群标号

r——空间位置

t——时间

i——粗网编号

——t时刻r处第g能群预估的中子通量密度

——t时刻i粗网处第g能群预估的中子通量密度

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

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

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

步骤3：在t_{n-1,transport}至t_n,transport时间间隔内求解输运形式的固定源方程，n＝1,N，其中N为整个计算过程划分的输运计算步数，并采用CMFD方法对输运方程进行加速，同时使用全局加速因子对CMFD方程的求解进行加速，具体如图3所示，求得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能群的等效缓发源项

用来加速上述输运形式的固定源方程的CMFD方程公式如下：

式中：

u——坐标轴方向标号

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

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

i——粗网标号

j——平源区标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网内第g能群的固定源项

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

Iⁱ——i粗网的平源区总数

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

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

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网j平源区第g能群的中子通量密度

S_F ^i,j——i粗网j平源区的裂变源

V^i,j——i粗网内j平源区的体积

对于上述CMFD形式的固定源方程，同时采用全局加速因子的技巧进行加速，公式如下：

式中：

f——全局加速因子

S——固定源项

Lφ⁽ⁿ⁾——消失项

Mφ⁽ⁿ⁾——产生项

<·>操作表示对全能量、全空间进行积分

φ⁽ⁿ⁾——n时刻的中子通量密度

u——坐标轴方向标号

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

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

_i——粗网标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

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

——i粗网第g能群的固定源系数

CMFD求解每执行一次，对求解出来的中子通量密度作用全局加速因子来作为下一步迭代使用的中子通量密度，具体公式为：

φ⁽ⁿ⁾＝f·φ⁽ⁿ⁾ (公式3‐4)

式中：

φ⁽ⁿ⁾——第n步粗网的中子通量密度

f——全局加速因子

步骤4：如图4所示，利用步骤1和步骤2中求得的平源区的中子通量密度、粗网中子通量密度以及粗网共轭中子通量密度信息，进行校正计算，得到校正后的t_n时刻的每个平源区中子标通量密度，具体过程如下：

式中：

u——坐标轴方向标号

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

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

i——粗网标号

j——平源区标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网内第g能群的固定源项

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

Iⁱ——i粗网的平源区总数

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

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

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网j平源区第g能群的中子通量密度

S_F ^i,j——i粗网j平源区的裂变源

V^i,j——i粗网内j平源区的体积

式中：

f——全局加速因子

S——固定源项

Lφ⁽ⁿ⁾——消失项

Mφ⁽ⁿ⁾——产生项

<·>操作表示对全能量、全空间进行积分

φ⁽ⁿ⁾——n时刻的中子通量密度

u——坐标轴方向标号

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

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

i——粗网标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

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

——i粗网第g能群的固定源系数

φ⁽ⁿ⁾＝f·φ⁽ⁿ⁾ (公式4‐3)

式中：

φ⁽ⁿ⁾——第n步粗网的中子通量密度

f——全局加速因子

式中：

i——粗网标号

g——能群标号

t——时间变量

——t时刻i粗网的第g能群的预估中子通量密度

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

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

——t时刻i粗网的第g能群的形状函数

C——归一化常数

v——中子速度

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

——t时刻粗网的形状函数

——临界状态下粗网的形状函数

φ_i,g(0)——临界状态下i粗网的第g能群的中子通量密度

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

式中：

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组缓发中子先驱核浓度

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

式中：

t——时间变量

i——缓发中子标号

N(t)——t时刻幅值

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

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

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

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

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

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

式中：

i——粗网标号

g——能群标号

t——时间

——t时刻i粗网第g能群校正后的中子通量密度

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

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

——t时刻i粗网第g能群粗网的形状函数

式中：

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

——t时刻i粗网第g能群校正后的中子通量密度

——t时刻i粗网第g能群预估的中子通量密度

ψ_g(r,t)——t时刻r处平源区第g能群粗网的形状函数

步骤5：重复执行步骤3与步骤4，直到动力学计算结束为止。

下面以时空中子动力学基准题TWIGL的计算结果为例说明本发明的效果：TWIGL基准题是一个二维两群的时空中子动力学基准题，计算时间从0.0s到0.5s，不采用加速方法的计算时间步长取2.5ms，采用加速方法的计算时间步长取20ms，计算结果如图5所示，从计算结果可以看出，不采用加速方法的计算时间为160分钟，采用加速方法的计算时间为20分钟，采用加速方法后计算时间大大减少，同时计算精度与采用小时间步长的不加速方法相当，此计算结果证明本发明具有相当高的实用价值以及创新性。

Claims

1.一种高保真时空中子动力学计算的加速方法，其特征在于：包括如下步骤：

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

式中：

Ω——角度方向

——梯度算子

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——输运计算得到的有效增殖因子

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

式中：

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处裂变源

式中：

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'能群的共轭中子通量密度

式中：

g——能群标号

r——空间位置

t——时间

i——粗网编号

——t时刻r处第g能群预估的中子通量密度

——t时刻i粗网处第g能群预估的中子通量密度

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

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

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

式中：

Ω——角度方向

——梯度算子

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能群的等效缓发源项

用来加速上述输运形式的固定源方程的CMFD方程公式如下：

式中：

u——坐标轴方向标号

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

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

i——粗网标号

j——平源区标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网内第g能群的固定源项

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

Iⁱ——i粗网的平源区总数

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

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

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网j平源区第g能群的中子通量密度

S_F ^i,j——i粗网j平源区的裂变源

V^i,j——i粗网内j平源区的体积

式中：

f——全局加速因子

S——固定源项

Lφ⁽ⁿ⁾——消失项

Mφ⁽ⁿ⁾——产生项

<·>操作表示对全能量、全空间进行积分

φ⁽ⁿ⁾——n时刻的中子通量密度

u——坐标轴方向标号

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

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

i——粗网标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

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

——i粗网第g能群的固定源系数

φ⁽ⁿ⁾＝f·φ⁽ⁿ⁾ (公式3-4)

式中：

φ⁽ⁿ⁾——第n步粗网的中子通量密度

f——全局加速因子

1)在t_{n-1,transport}至t_n,transport时刻的时间间隔内根据设定的CMFD计算时间步长Δt_CMFD划分CMFD计算所需的时间间隔，使用t_{n-1,transport}时刻的CMFD参数和t_n,transport时刻的CMFD参数在计算间隔t_n-1,CMFD至t_n,CMFD间进行线性插值CMFD参数，其中以Δt_CMFD为时间步长求解CMFD固定源方程，同时对CMFD固定源方程的求解应用全局加速因子，得到t_n,CMFD时刻粗网的中子通量密度具体公式如下：

式中：

u——坐标轴方向标号

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

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

i——粗网标号

j——平源区标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网内第g能群的固定源项

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

Iⁱ——i粗网的平源区总数

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

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

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网第g能群的固定源系数

——i粗网j平源区第g能群的固定源系数

——i粗网j平源区第g能群的中子通量密度

S_F ^i,j——i粗网j平源区的裂变源

V^i,j——i粗网内j平源区的体积

式中：

f——全局加速因子

S——固定源项

Lφ⁽ⁿ⁾——消失项

Mφ⁽ⁿ⁾——产生项

<·>操作表示对全能量、全空间进行积分

φ⁽ⁿ⁾——n时刻的中子通量密度

u——坐标轴方向标号

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

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

i——粗网标号

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

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

g——当前能群标号

g'——非当前能群标号

G——总能群数

——i粗网内第g能群的移出截面

——n时刻i粗网的裂变源

——i粗网第g'群到第g群的散射截面

——i粗网第g能群的中子产生截面截面

——n时刻i粗网第g能群的中子通量密度

——n时刻i粗网第g'能群的中子通量密度

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

——i粗网第g能群的固定源系数

φ⁽ⁿ⁾＝f·φ⁽ⁿ⁾ (公式4-3)

式中：

φ⁽ⁿ⁾——第n步粗网的中子通量密度

f——全局加速因子

式中：

i——粗网标号

g——能群标号

t——时间变量

——t时刻i粗网的第g能群的预估中子通量密度

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

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

——t时刻i粗网的第g能群的形状函数

C——归一化常数

v——中子速度

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

——t时刻粗网的形状函数

——临界状态下粗网的形状函数

φ_i,g(0)——临界状态下i粗网的第g能群的中子通量密度

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

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

<mrow> <msup> <mi>&beta;</mi> <mrow> <mi>e</mi> <mi>f</mi> <mi>f</mi> </mrow> </msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mo>&Sigma;</mo> <mi>k</mi> </munder> <msubsup> <mi>&beta;</mi> <mi>k</mi> <mrow> <mi>e</mi> <mi>f</mi> <mi>f</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>&lambda;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo><</mo> <msup> <mi>&phi;</mi> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>&lambda;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>&chi;</mi> <mrow> <mi>d</mi> <mi>k</mi> <mo>,</mo> <mi>g</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>C</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>></mo> </mrow> <mrow> <mo><</mo> <msup> <mi>&phi;</mi> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>&chi;</mi> <mrow> <mi>d</mi> <mi>k</mi> <mo>,</mo> <mi>g</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>C</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>></mo> </mrow> </mfrac> </mrow>

<mrow> <msub> <mi>c</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo><</mo> <msup> <mi>&phi;</mi> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>&chi;</mi> <mrow> <mi>d</mi> <mi>k</mi> <mo>,</mo> <mi>g</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>C</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>></mo> </mrow> <mrow> <mo><</mo> <msup> <mi>&phi;</mi> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>&chi;</mi> <mrow> <mi>d</mi> <mi>k</mi> <mo>,</mo> <mi>g</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>C</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>></mo> <mi>&Lambda;</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>

式中：

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组缓发中子先驱核浓度

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

4)在t_n-1,CMFD与t_n,CMFD之间以Δt_PK为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔t_n-1,PK至t_n,PK内用t_n-1,CMFD与t_n,CMFD时间点上的点堆参数进行插值，其中进而求解点堆方程组，得到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组缓发中子衰变常数

式中：

i——粗网标号

g——能群标号

t——时间

——t时刻i粗网第g能群校正后的中子通量密度

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

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

——t时刻i粗网第g能群粗网的形状函数

式中：

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

——t时刻i粗网第g能群校正后的中子通量密度

——t时刻i粗网第g能群预估的中子通量密度

ψ_g(r,t)——t时刻r处平源区第g能群粗网的形状函数

步骤5：重复执行步骤3与步骤4，直到动力学计算结束为止。