一种高保真时空中子动力学计算的加速方法与流程
本发明涉及核反应堆堆芯设计和核安全技术领域，具体涉及一种高保真时空中子动力学计算的加速方法。
背景技术：
压水堆核电厂的瞬态和事故工况分析长期基于不考虑空间影响的点堆模型，所谓点堆模型即认为整个核反应堆的功率变化与空间位置无关，显然点堆模型并不适用于功率分布畸变，如弹棒事故、落棒事故、主蒸汽管道破口事故等。近年来随着计算机技术的迅速发展，在点堆模型基础上发展了三维时空中子动力学计算方法，考虑了功率变化在三维空间上的具体分布，在数值模拟的精准性上跨上了一个台阶。但其所采用的计算方法仍是基于传统的两步均匀化方法，无法排除均匀化带来的误差，且不能直接考虑核反应堆内复杂的几何结构，无法直接获得核反应堆内精细的功率分布。
所谓两步均匀化方法，即第一步对构成核反应堆的每一种组件进行单独的计算，第二步对以组件为单位的核反应堆进行计算。在第一步计算中利用全反射边界条件单独地对各种组件进行多能群的中子输运计算，并在组件处于不同的工况下分别进行计算，得到每种组件在不同工况下的多能群的中子通量密度在组件内的分布，进而利用核反应率守恒归并出两能群的、不同工况下、在每种组件内跟空间位置无关的等效均匀化群常数以及不连续因子等物理量；在第二步计算中根据核反应堆所处的工况，对每种组件的等效均匀化群常数利用多项式拟合或多项式插值来获得当前计算工况下的截面参数等，在获得每个组件的截面参数后对核反应堆进行两能群的中子扩散计算，并对各个组件进行功率重构计算，得到三维核反应堆堆芯的功率分布。
在两步均匀化方法计算中只考虑有限种独立的组件计算来计算等效均匀化参数，因此它不能精确考虑在核反应堆中组件之间的相互影响，并且不能精确考虑核反应堆内复杂的几何结构。随着对计算精度要求的不断提高，只有高保真时空中子动力学计算才能根本解决两步均匀化法的问题。
所谓高保真时空中子动力学，就是理论和数值模型不能失真、所得解为真解，计算分辨率和计算精度要高，计算结果精准。要实现高保真时空中子动力学计算，就必须采用一步法对核反应堆进行直接求解。所谓一步法即对核反应堆直接进行几何建模，采用多能群结构，在计算中直接进行中子输运方程的求解，不引入几何近似以及均匀化带来的误差。通过高保真时空中子动力学的模拟可大大提高核反应堆设计的精确性，大幅度减少实验量，提高反应堆的安全性和经济性，对反应堆改进设计和新型反应堆设计均具有重大的理论价值和应用前景。
但是高保真时空中子动力学的应用存在一个非常明显的缺陷，即计算时间过长，计算资源消耗过大。因此对高保真时空中子动力学计算的加速在高保真时空中子动力学计算中起着十分关键的作用，采取高效的加速方法是高保真时空中子动力学计算能够得到实际应用的关键要素之一。
技术实现要素：
为了克服上述现有技术的问题，解决高保真时空中子动力学计算存在的计算时间过长的问题，本发明的目的是提出一种对高保真时空中子动力学计算加速的方法，该方法能在保证计算精度的前提下，能够大幅度减少高保真时空中子动力学计算的时间步数，同时减少在每个时间步上的计算时间，进而大幅缩短整个计算过程的计算耗时。
为了达到上述目的，本发明采用如下技术方案：
一种高保真时空中子动力学计算的加速方法，包括如下步骤：
步骤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群的中子产生截面
φg(r)——r处第g能群的中子标通量密度
φg'(r)——r处第g'能群的中子标通量密度
sf(r)——r处裂变源
keff——输运计算得到的有效增殖因子
由此得到各个平源区的中子通量密度；
5)根据4)中计算所得到的各个平源区的中子通量密度，得到临界状态下的各组缓发中子先驱核密度，具体的公式如下：
式中：
r——空间位置
k——缓发中子先驱核编号
g——能群编号
g——能群总数
keff——输运计算得到的有效增殖因子
ck(r)——r处临界状态第k组缓发中子先驱核密度
βk(r)——r处第k组缓发中子份额
νσf,g——第g群的中子产生截面
λk(r)——r处第k组缓发中子先驱核的衰变常数
φg(r)——r处临界状态第g群的中子通量密度
sf(r)——r处裂变源
6)根据1)、2)、3)得到的信息，对粗网有限差分cmfd方程进行中子共轭计算，得到各个粗网的共轭中子通量密度，具体的计算公式如下：
式中：
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
hu——粗网在u方向的高度
g——当前能群标号
g'——非当前能群标号
g——总能群数
σrg——第g能群的移出截面
σs,g→g'——第g群到第g'群的散射截面
vσf,g——第g能群的中子产生截面
χg'——第g'能群的裂变谱
keff——输运计算得到的有效增殖因子
φg——粗网第g能群的共轭中子通量密度
φg'——粗网第g'能群的共轭中子通量密度
7)将每个平源区的中子通量密度利用因子分解分解成为幅值函数与形状函数的乘积，其中幅值函数为平源区对应粗网的中子通量密度，此时的幅值函数是预估的值，存在较大误差，需要后续进一步校正，利用4)求解得到的临界状态下的中子通量密度求得初始时刻的中子形状函数，具体的计算公式如下：
式中：
g——能群标号
r——空间位置
t——时间
i——粗网编号
ψg(r,t)——t时刻r处第g能群的中子通量密度形状函数
φg(r,0)——临界状态r处第g能群的中子通量密度
ψg(r,0)——临界时刻r处第g能群的中子通量密度形状函数
步骤2：根据输入卡片的描述执行截面扰动，打破核反应堆的临界状态，从而开始时空中子动力学计算；
步骤3：在tn-1,transport至tn,transport时间间隔内求解输运形式的固定源方程，n＝1,n，其中n为整个计算过程划分的输运计算步数，并采用cmfd方法对输运方程进行加速，同时使用全局加速因子对cmfd方程的求解进行加速，求得tn,transport时刻的中子通量密度，具体的计算公式如下：
式中：
ω——角度方向
g——当前能群标号
g'——非当前能群标号
g——能群总数
r——空间位置
n——第n个输运计算步
k——缓发中子标号
keff——输运计算得到的有效增殖因子
ag(r)——r处第g能群的固定源系数
bg(r)——r处第g能群的固定源系数
cg(r)——r处第g能群的固定源系数
vg——第g能群的中子速度
δtn——第n个输运步的步长
χdk,g——第k组缓发中子在第g能群的缓发裂变谱
βk——第k组缓发中子份额
步骤4：利用步骤1和步骤2中求得的平源区的中子通量密度、粗网中子通量密度以及粗网共轭中子通量密度信息，进行校正计算，得到校正后的tn时刻的每个平源区中子标通量密度，具体过程如下：
1)在tn-1,transport至tn,transport时刻的时间间隔内根据设定的cmfd计算时间步长δtcmfd划分cmfd计算所需的时间间隔，使用tn-1,transport时刻的cmfd参数和tn,transport时刻的cmfd参数在计算间隔tn-1,cmfd至tn,cmfd间进行线性插值cmfd参数，其中n＝1,
式中：
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
i——粗网标号
j——平源区标号
g——当前能群标号
g'——非当前能群标号
g——总能群数
ii——i粗网的平源区总数
χg——第g能群的裂变谱
keff——输运计算得到的有效增殖因子
sfi,j——i粗网j平源区的裂变源
vi,j——i粗网内j平源区的体积
采用全局加速因子加速上述cmfd形式的固定源方程，具体如下：
式中：
f——全局加速因子
s——固定源项
lφ(n)——消失项
mφ(n)——产生项
<·>操作表示对全能量、全空间进行积分
φ(n)——n时刻的中子通量密度
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
i——粗网标号
g——当前能群标号
g'——非当前能群标号
g——总能群数
χg——第g能群的裂变谱
φ(n)＝f·φ(n)(公式4‐3)
式中：
φ(n)——第n步粗网的中子通量密度
f——全局加速因子
2)将tn,cmfd时刻粗网的中子通量密度因子分解为幅值函数与形状函数的乘积，进而利用步骤1中求得的粗网的共轭中子通量密度和步骤1中求得的粗网的中子通量密度计算归一化常数c，从而计算tn,cmfd时刻的粗网的中子通量密度形状函数，具体公式如下：
式中：
i——粗网标号
g——能群标号
t——时间变量
np(t)——t时刻预估的幅值函数
n(t)——t时刻的幅值函数
c——归一化常数
v——中子速度
φ*——粗网的共轭中子通量密度
φi,g(0)——临界状态下i粗网的第g能群的中子通量密度
<·>——对全能量、全空间积分
3)分别计算tn-1,cmfd与tn,cmfd时刻的精确点堆参数，具体公式如下：
式中：
t——时间
g——当前能群标号
g'——非当前能群标号
g——总能群数
r——空间位置
k——缓发中子标号
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
hu——粗网在u方向的高度
σrg——第g能群的移出截面
σs,g→g'——第g群到第g'群的散射截面
vσf,g——第g能群的中子产生截面
χg(r)——r处第g能群的裂变谱
keff——输运计算得到的有效增殖因子
φ*(r)——r处粗网的共轭中子通量密度
sf(r,t)——t时刻r处的裂变源
keff——输运计算得到的有效增殖因子
f(t)——t时刻精确点堆参数分母
ρ(t)——t时刻的反应性
βk(r)——r处第k组缓发中子份额
χdk,g(r)——r处第k组缓发中子在第g能群的缓发裂变谱
v(r)——r处中子速度
λ(t)——t时刻等效中子代时间
λk(t)——t时刻第k组缓发中子衰变常数
ck(r,t)——t时刻r处临界状态第k组缓发中子先驱核密度
λk(r)——r处第k组缓发中子先驱核的衰变常数
ck(t)——t时刻第k组缓发中子先驱核浓度
<·>——在全能量、全空间进行积分
4)在tn-1,cmfd与tn,cmfd之间以δtpk为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔tn-1,pk至tn,pk内用tn-1,cmfd与tn,cmfd时间点上的点堆参数进行插值，其中n＝1,
式中：
t——时间变量
i——缓发中子标号
n(t)——t时刻幅值
ci(t)——t时刻第i组缓发中子先驱核浓度
ρ(t)——t时刻反应性
βi(t)——t时刻第i组缓发中子份额
β(t)——t时刻缓发中子份额总和
λ(t)——t时刻等效中子代时间
λi——第i组缓发中子衰变常数
5)当点堆计算的时间点达到tn,cmfd时，校正tn,cmfd时刻的粗网的中子通量密度，具体公式如下：
式中：
i——粗网标号
g——能群标号
t——时间
nc(t)——t时刻校正的幅值函数
np(t)——t时刻预估的幅值函数
6)循环上述1)至5)过程，当tn,cmfd＝tn时，用tn,cmfd时刻的粗网中子通量密度校正tn,transport时刻的每个平源区的中子通量密度，具体公式如下：
式中：
ψg(r,t)——t时刻r处平源区第g能群粗网的形状函数
步骤5：重复执行步骤3与步骤4，直到动力学计算结束为止。
与现有技术相比，本发明具有如下突出优点：
1)采用一步法计算核反应堆中的时空中子动力学方程，不需要进行组件均匀化，可精确处理复杂的非均匀几何问题，大大提高计算精度，达到计算结果高度保真的效果，满足高保真计算的要求。
2)采用了多级的预估校正准静态策略与全局加速因子，能够在保证计算精度的前提下大大减少计算的时间步长，并大大减少每个时间步上的计算时间，显著提高计算效率，对高保真时空中子动力学计算应用到实际的核反应堆设计中去具有重要意义。
附图说明
图1为细网结构图。
图2为堆芯临界计算流程图。
图3为cmfd与全局加速因子加速输运固定源计算流程图。
图4为瞬态多级加速计算流程图。
图5为twigl基准题有无加速计算结果对比图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细的说明：
步骤1：如图2堆芯临界计算流程图所示，进行核反应堆临界计算，采用一步直接输运计算，得到核反应堆处于临界状态时每个能群、每个平源区的中子通量密度以及每个平源区内的每组缓发中子先驱核的密度，具体包括如下步骤：
1)从截面库中读取各个核素的原始多群宏观截面与动力学参数信息；
2)从输入卡片中读取核反应堆的几何信息与计算条件；
3)根据输入卡片中的几何信息建立计算模型：首先根据输入卡片的几何描述得到核反应堆的几何布置；其次根据核反应堆的几何布置建立特征线计算所需的边界条件以及内部特征线的长度信息；为输运计算提供模块化特征信息；
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群的中子产生截面
φg(r)——r处第g能群的中子标通量密度
φg'(r)——r处第g'能群的中子标通量密度
sf(r)——r处裂变源
keff——输运计算得到的有效增殖因子
由此得到各个平源区的中子通量密度；
5)根据4)中计算所得到的各个平源区的中子通量密度，得到临界状态下的各组缓发中子先驱核密度，具体的公式如下：
式中：
r——空间位置
k——缓发中子先驱核编号
g——能群编号
g——能群总数
keff——输运计算得到的有效增殖因子
ck(r)——r处临界状态第k组缓发中子先驱核密度
βk(r)——r处第k组缓发中子份额
νσf,g——第g群的中子产生截面
λk(r)——r处第k组缓发中子先驱核的衰变常数
φg(r)——r处临界状态第g群的中子通量密度
sf(r)——r处裂变源
6)根据1)、2)、3)得到的信息，对粗网有限差分cmfd方程进行中子共轭计算，得到各个粗网的共轭中子通量密度，具体的计算公式如下：
式中：
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
hu——粗网在u方向的高度
g——当前能群标号
g'——非当前能群标号
g——总能群数
σrg——第g能群的移出截面
σs,g→g'——第g群到第g'群的散射截面
vσf,g——第g能群的中子产生截面
χg'——第g'能群的裂变谱
keff——输运计算得到的有效增殖因子
φg——粗网第g能群的共轭中子通量密度
φg'——粗网第g'能群的共轭中子通量密度
7)将每个平源区的中子通量密度利用因子分解分解成为幅值函数与形状函数的乘积，其中幅值函数为平源区对应粗网的中子通量密度，此时的幅值函数是预估的值，存在较大误差，需要后续进一步校正，利用4)求解得到的临界状态下的中子通量密度求得初始时刻的中子形状函数，具体的计算公式如下：
式中：
g——能群标号
r——空间位置
t——时间
i——粗网编号
ψg(r,t)——t时刻r处第g能群的中子通量密度形状函数
φg(r,0)——临界状态r处第g能群的中子通量密度
ψg(r,0)——临界时刻r处第g能群的中子通量密度形状函数
步骤2：根据输入卡片的描述执行截面扰动，打破核反应堆的临界状态，从而开始时空中子动力学计算；
步骤3：在tn-1,transport至tn,transport时间间隔内求解输运形式的固定源方程，n＝1,n，其中n为整个计算过程划分的输运计算步数，并采用cmfd方法对输运方程进行加速，同时使用全局加速因子对cmfd方程的求解进行加速，具体如图3所示，求得tn,transport时刻的中子通量密度，具体的计算公式如下：
式中：
ω——角度方向
g——当前能群标号
g'——非当前能群标号
g——能群总数
r——空间位置
n——第n个输运计算步
k——缓发中子标号
keff——输运计算得到的有效增殖因子
ag(r)——r处第g能群的固定源系数
bg(r)——r处第g能群的固定源系数
cg(r)——r处第g能群的固定源系数
vg——第g能群的中子速度
δtn——第n个输运步的步长
χdk,g——第k组缓发中子在第g能群的缓发裂变谱
βk——第k组缓发中子份额
用来加速上述输运形式的固定源方程的cmfd方程公式如下：
式中：
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
i——粗网标号
j——平源区标号
g——当前能群标号
g'——非当前能群标号
g——总能群数
ii——i粗网的平源区总数
χg——第g能群的裂变谱
keff——输运计算得到的有效增殖因子
sfi,j——i粗网j平源区的裂变源
vi,j——i粗网内j平源区的体积
对于上述cmfd形式的固定源方程，同时采用全局加速因子的技巧进行加速，公式如下：
式中：
f——全局加速因子
s——固定源项
lφ(n)——消失项
mφ(n)——产生项
<·>操作表示对全能量、全空间进行积分
φ(n)——n时刻的中子通量密度
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
i——粗网标号
g——当前能群标号
g'——非当前能群标号
g——总能群数
χg——第g能群的裂变谱
cmfd求解每执行一次，对求解出来的中子通量密度作用全局加速因子来作为下一步迭代使用的中子通量密度，具体公式为：
φ(n)＝f·φ(n)(公式3‐4)
式中：
φ(n)——第n步粗网的中子通量密度
f——全局加速因子
步骤4：如图4所示，利用步骤1和步骤2中求得的平源区的中子通量密度、粗网中子通量密度以及粗网共轭中子通量密度信息，进行校正计算，得到校正后的tn时刻的每个平源区中子标通量密度，具体过程如下：
1)在tn-1,transport至tn,transport时刻的时间间隔内根据设定的cmfd计算时间步长δtcmfd划分cmfd计算所需的时间间隔，使用tn-1,transport时刻的cmfd参数和tn,transport时刻的cmfd参数在计算间隔tn-1,cmfd至tn,cmfd间进行线性插值cmfd参数，其中n＝1,
式中：
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
i——粗网标号
j——平源区标号
g——当前能群标号
g'——非当前能群标号
g——总能群数
ii——i粗网的平源区总数
χg——第g能群的裂变谱
keff——输运计算得到的有效增殖因子
sfi,j——i粗网j平源区的裂变源
vi,j——i粗网内j平源区的体积
采用全局加速因子加速上述cmfd形式的固定源方程，具体如下：
式中：
f——全局加速因子
s——固定源项
lφ(n)——消失项
mφ(n)——产生项
<·>操作表示对全能量、全空间进行积分
φ(n)——n时刻的中子通量密度
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
i——粗网标号
g——当前能群标号
g'——非当前能群标号
g——总能群数
χg——第g能群的裂变谱
φ(n)＝f·φ(n)(公式4‐3)
式中：
φ(n)——第n步粗网的中子通量密度
f——全局加速因子
2)将tn,cmfd时刻粗网的中子通量密度因子分解为幅值函数与形状函数的乘积，进而利用步骤1中求得的粗网的共轭中子通量密度和步骤1中求得的粗网的中子通量密度计算归一化常数c，从而计算tn,cmfd时刻的粗网的中子通量密度形状函数，具体公式如下：
式中：
i——粗网标号
g——能群标号
t——时间变量
np(t)——t时刻预估的幅值函数
n(t)——t时刻的幅值函数
c——归一化常数
v——中子速度
φ*——粗网的共轭中子通量密度
φi,g(0)——临界状态下i粗网的第g能群的中子通量密度
<·>——对全能量、全空间积分
3)分别计算tn-1,cmfd与tn,cmfd时刻的精确点堆参数，具体公式如下：
式中：
t——时间
g——当前能群标号
g'——非当前能群标号
g——总能群数
r——空间位置
k——缓发中子标号
u——坐标轴方向标号
x,y,z——直角坐标系坐标轴方向
hu——粗网在u方向的高度
σrg——第g能群的移出截面
σs,g→g'——第g群到第g'群的散射截面
vσf,g——第g能群的中子产生截面
χg(r)——r处第g能群的裂变谱
keff——输运计算得到的有效增殖因子
φ*(r)——r处粗网的共轭中子通量密度
sf(r,t)——t时刻r处的裂变源
keff——输运计算得到的有效增殖因子
f(t)——t时刻精确点堆参数分母
ρ(t)——t时刻的反应性
βk(r)——r处第k组缓发中子份额
χdk,g(r)——r处第k组缓发中子在第g能群的缓发裂变谱
v(r)——r处中子速度
λ(t)——t时刻等效中子代时间
λk(t)——t时刻第k组缓发中子衰变常数
ck(r,t)——t时刻r处临界状态第k组缓发中子先驱核密度
λk(r)——r处第k组缓发中子先驱核的衰变常数
ck(t)——t时刻第k组缓发中子先驱核浓度
<·>——在全能量、全空间进行积分
4)在tn-1,cmfd与tn,cmfd之间以δtpk为时间步长等距划分点堆计算所需的间隔，在点堆所计算的时间间隔tn-1,pk至tn,pk内用tn-1,cmfd与tn,cmfd时间点上的点堆参数进行插值，其中n＝1,
式中：
t——时间变量
i——缓发中子标号
n(t)——t时刻幅值
ci(t)——t时刻第i组缓发中子先驱核浓度
ρ(t)——t时刻反应性
βi(t)——t时刻第i组缓发中子份额
β(t)——t时刻缓发中子份额总和
λ(t)——t时刻等效中子代时间
λi——第i组缓发中子衰变常数
5)当点堆计算的时间点达到tn,cmfd时，校正tn,cmfd时刻的粗网的中子通量密度，具体公式如下：
式中：
i——粗网标号
g——能群标号
t——时间
nc(t)——t时刻校正的幅值函数
np(t)——t时刻预估的幅值函数
6)循环上述1)至5)过程，当tn,cmfd＝tn时，用tn,cmfd时刻的粗网中子通量密度校正tn,transport时刻的每个平源区的中子通量密度，具体公式如下：
式中：
ψg(r,t)——t时刻r处平源区第g能群粗网的形状函数
步骤5：重复执行步骤3与步骤4，直到动力学计算结束为止。
下面以时空中子动力学基准题twigl的计算结果为例说明本发明的效果：twigl基准题是一个二维两群的时空中子动力学基准题，计算时间从0.0s到0.5s，不采用加速方法的计算时间步长取2.5ms，采用加速方法的计算时间步长取20ms，计算结果如图5所示，从计算结果可以看出，不采用加速方法的计算时间为160分钟，采用加速方法的计算时间为20分钟，采用加速方法后计算时间大大减少，同时计算精度与采用小时间步长的不加速方法相当，此计算结果证明本发明具有相当高的实用价值以及创新性。
- 还没有人留言评论。精彩留言会获得点赞！