Burnup lib compression method based on pseudo decay nuclides definement
-
摘要: 燃耗数据库的构建决定了燃耗和衰变热计算的准确性,评价核数据库中的燃耗信息过于复杂,导致燃耗矩阵规模大,刚性强,计算效率低。从燃耗数据库的基本组成出发,考虑燃耗数据库中,各核素及其转化关系对中子学计算精度和重要核素核子密度计算精度的影响,并作为燃耗库压缩的依据。对于因裂变产物压缩而损失的衰变热计算精度,通过非线性最小二乘优化算法拟合衰变释热函数,构造伪衰变核代替裂变产物衰变热计算,以保持衰变热的计算精度。验证结果表明,原精细燃耗库中有超过1 500种核素,经压缩后保留不足200种核素。压缩后的燃耗数据库在有效增殖因子计算和核子密度计算中并未引入明显偏差。在衰变热计算方面,伪衰变核对于衰变热计算精度有显著的复原效果,对总功率贡献的计算偏差小于0.5%,满足衰变热计算精度的需求。Abstract: The construction of the burnup lib determines the accuracy of burnup and decay heat calculations. The evaluation of burnup information in the nuclear lib is complex, leading to a large, rigid, and inefficient burnup matrix. This paper begins with the basic composition of the burnup lib, considering the impact of each nuclide and its transformation relationships on the accuracy of neutronics calculations and target nuclide nuclear density calculations, which serves as the basis for the compression of the burnup lib. To address the decay heat calculation accuracy loss caused by the compression of fission products, a nonlinear least squares optimization algorithm is used to fit the decay heat release function, and pseudo-decay nuclides are constructed to replace the fission product decay heat calculation, thereby maintaining the accuracy of decay heat calculations. Verification results show that the original detailed burnup lib contains more than 1 500 nuclides, which are reduced to fewer than 200 nuclides after compression. The compressed burnup lib does not introduce significant deviations in the calculation of the effective multiplication factor and nuclear density. In terms of decay heat calculations, the pseudo-decay nuclides significantly restore the decay heat calculation accuracy, with the contribution of decay heat to total power having a calculation deviation of less than 0.5%, meeting the required accuracy for decay heat calculations.
-
燃耗和衰变热计算是核反应堆安全分析的重要环节,在计算中需要考虑裂变、俘获、衰变等诸多反应对核素的相互转化以及能量释放的影响。燃耗数据库描述核素的性质以及燃耗过程中核素之间的转化关系,它的构建决定了计算精度。由评价核数据库[1]定义的精细燃耗数据库中,包含了数千种核素以及核素间的转化关系,对计算和储存性能提出较高需求。若能将燃耗数据库压缩,则可提高计算与储存效率。燃耗库的压缩算法可分为定性分析和定量分析,早期的压缩算法多为经验或半经验性质的定性分析方法[2-3]。近年来,国内外学者提出定量分析方法,如贡献矩阵法[4]、奇异值分解法[5]、贡献函数法[6]和双约束压缩方法[7]等。这类方法较早期的定性分析方法相比,能更好保留计算精度和提高计算效率,但仍存在不足。前两种方法在分析中,局限于特定的燃耗矩阵,也未考虑绝对反应率影响。贡献函数方法仅限于裂变产物的压缩,并对所有非保留核素进行衰变平衡态的近似处理。双约束压缩方法在压缩时缺少对核素中子学性质的分析,采用伪核素提升压缩燃耗数据库计算精度。同时,以上方法在压缩了大量裂变产物,难以保持衰变热计算精度。
为考虑衰变热计算精度,在压缩过程中会遇到以下难点:(1)对衰变热有贡献的裂变产物多达数百种;(2)裂变产物的半衰期分布在毫秒至数十年;(3)裂变产物受中子辐照影响,衰变特性会发生变化。基于上述问题,本文从定量化分析方法[8]出发,构造伪衰变核加入燃耗数据库,替代裂变产物衰变热计算以保持计算精度。
1. 理论模型
1.1 燃耗数据库压缩方法框架
在燃耗计算时,基于评价核数据库求解燃耗方程(1)[9],需考虑数千种核素的衰变、裂变和吸收反应等信息,因而对计算资源和效率提出较高要求。然而在真实堆芯中,众多核素及转化关系对堆芯反应性以及原子核密度场影响较小,如短寿期核素在产生后便迅速发生衰变。事实上,核素及其转化关系也是燃耗库的基本组成。因此,可从燃耗库基本组成出发,评价各核素及其转化关系对堆芯反应性和核子密度的相对重要性,保留重要的核素及其转化关系,从而实现燃耗库压缩。
{dNi(r,t)dt=βi−1Ni−1(r,t)−(λidecay+G∑g=1σa,g,iφg(r,t))Ni(r,t)+Fiβi−1={λdecayi−1①G∑g=1σa,g,i−1φg(r,t)②Fi=G∑g′=1∑i′γi,i′σf,g′,i′φg′(r,t)Ni′(r,t) (1) 式中:
{N_i}({{\boldsymbol{r}}},t) 为在位置r时刻t处第i核素的核子密度(cm−3);{\lambda _i}^{{\text{decay}}} 为第i核素衰变常数(s−1);{\sigma _{{\text{a}},g,i}} 为第i核素在第g能群下的微观吸收截面(cm2);{\varphi _g}(r,t) 为在能群g中,时刻t与空间r处的中子通量密度(cm2·s1);{\sigma _{{\text{f}},g,i}} 为第i核素在第g能群下的微观裂变截面(cm2);{\gamma _{i,{i^{'}}}} 为第{i'}{\text{ }} 易裂变核裂变时对第i核素的产额(%)。\;{ \beta _{i - 1}} 存在两种形式:1)当第i核素由第i−1核素发生衰变产生时,选择式①代入计算;2)当第i核素由第i−1核素发生中子反应产生时,选择式②代入计算。为准确评估燃耗库基本组成对堆芯反应性的影响,本文从自主开发的中子学-燃耗输运耦合程序出发,其中,中子学功能由变分节块法中子输运程序VITAS[10]提供,采用CRAM[11]方法实现燃耗求解功能。燃耗库的压缩框架如图1所示,大体上分为两部分:(1)选取特征问题进行燃耗计算,基于燃耗结果进行衰变热分析,筛选出对衰变热计算有重要贡献的核素与临界安全和源项计算中的重要核素共同作为目标核素;对燃耗库基本组成逐一进行基于目标核素的重要性分析评价,保留相对重要的核素及其转化关系形成适用于中子学计算的压缩燃耗库。(2)将堆芯中的裂变过程拆解为理想化的裂变系统,用伪衰变核替代各裂变系统中的裂变产物,并推导衰变释热函数;根据真实衰变释热数据拟合出衰变释热函数中的参数,进而构造伪衰变核并加入压缩燃耗库,达到保留衰变热计算精度的目的。
1.2 基本压缩操作的重要性二元组
从燃耗库基本组成出发,分别执行以下四类基本压缩操作:(1)中子反应道的删除,删除特定中子反应道的转化关系;(2)衰变反应道的删除,删除特定衰变反应道的转化关系;(3)核素的删除,删除特定核素以及与该核素相关的中子反应道和衰变反应道转化关系;(4)核素的衰变平衡态处理,删除核素及附属于该核素的中子反应道转换关系,并令该核素的产生的转化关系指向衰变子核,同时更新产生转化关系的裂变产额、反应道分支比等信息。在特定问题中,执行基本压缩操作后从三个方面对作展开评价:中子产生精度损失,中子吸收精度损失和目标核素的核子密度计算精度损失。在辐照时间T内,中子产生精度损失可定义为
{L_{{\text{fission}}}} = \quad \dfrac{{\left| {{{\rho '}_{{\text{fission}}}} - {\rho _{{\text{fission}}}}} \right|}}{{{\rho _{{\text{fission}}}}}} (2) 式中:
{\rho _{{\text{fission}}}} 和{\rho '_{{\text{fission}}}} 分别表示执行基本压缩操作前后的中子产生密度,中子产生密度定义如下{\rho }_{\text{fission}}=\displaystyle\sum _{l\in {\bf{AN}}}{\overline{N}}_{l}T{\displaystyle\int }_{0}^{\infty }{\nu }_{l}(E){\sigma }_{l,\text{f}}(E)\phi (E)\text{d}E (3) 式中:AN代表锕系核素所组成的集合;
{\bar N_l} 为第l号核素在辐照时间T内,积分平均的核子密度(cm−3),采用切比雪夫有利近似法[9](CRAM)计算;νl(E)为每次裂变平均释放中子数;σl,f (E)为第l号核素微观裂变截面(cm2);ϕ(E)为中子通量密度(cm2·s)。同理,中子吸收精度损失定义为
\left\{ \begin{array}{l} {\rho }_{\text{absorption}}=\displaystyle\sum _{l=1}^{n}{\overline{N}}_{l}T{\displaystyle\int }_{0}^{\infty }{\sigma }_{l,\text{a}}(E)\phi (E)\text{d}E\\ {L}_{\text{absorption}}=\dfrac{\left|{{\rho }^{\prime }}_{\text{absorption}}-{\rho }_{\text{absorption}}\right|}{{\rho }_{\text{absorption}}}\end{array} \right. (4) 式中:n是核素总数;σl,a(E)为第l号核素微观吸收截面(cm2);
{\rho _{{\text{absorption}}}} 和{\rho '_{{\text{absorption}}}} 分别表示执行基本压缩操作前后的中子消失密度。目标核素的核子密度计算精度损失定义为
{L}_{\text{nuclide}}=\displaystyle\sum _{l=1}^{m}\text{ }\displaystyle\sum _{d\in {{\bf{TN}}}_{l}}\dfrac{\left|{{N}^{\prime }}_{l,d}-{N}_{l,d}\right|}{{N}_{l,d}} (5) 式中:m是基于目标核素衰变特性设定衰变时间的数目,TNl代表在第l个衰变时间点下,定义的目标核素的集合,目标核素根据用户需求而变化;
{N_{l,d}} 和{N'_{l,d}} 分别表示在执行基本压缩操作前后、第d号核素在第l衰变时间下的核子密度。最终,用于评价基本压缩操作的重要性二元组可定义如下
\left\{ \begin{array}{l} {s}_{1}=\underset{i\in {\bf{I}},j\in {\bf{J}},k\in {\bf{K}}}{\mathrm{max}}\left({L}_{\text{fission},i,j,k}\right)+\underset{i\in {\bf{I}},j\in {{\bf{J}}}_{i},k\in {\bf{K}}}{\mathrm{max}}\left({L}_{\text{absorption},i,j,k}\right)\\ {s}_{2}=\underset{i\in {\bf{I}},j\in {\bf{J}},k\in {\bf{K}}}{\mathrm{max}}\left({L}_{\text{nuclide},i,j,k}\right)\end{array} \right. (6) 式中:I、J、K分别代表特征问题、燃耗区、燃耗步所组成的集合。
{s_1},{s_2} 分别表示基本压缩操作中的中子学计算精度损失指标和对目标核素核子密度计算精度损失指标。设定合适的截断值{\varepsilon _1},{\varepsilon _2} ,当{s_1} > {\varepsilon _1},{s_2} > {\varepsilon _2} 时,执行基本压缩操作。1.3 衰变热先驱核拟合函数
参照缓发中子先驱核的思路,将裂变产物的衰变热计算用少量的伪衰变核代替(伪衰变核或称为衰变热先驱核),以达到简化衰变热计算的目的。在理想情况下,裂变产物不发生中子反应,衰变特性不发生变化(称之为衰变假设),可通过构造数量有限的衰变热先驱核,代替裂变产物的衰变热计算。事实上,在中子辐照过程中,受中子辐照的裂变产物(NIFP)的衰变特性会因发生中子反应而改变,可称为中子辐照效应。筛选出受中子辐照效应显著影响的NIFP核素集合,并假设剩余裂变产物不受中子辐照效应的影响,仅发生衰变反应,拟合剩余裂变产物的衰变释热曲线,进而构造衰变热先驱核。
衰变热先驱核参数的拟合从裂变系统出发,裂变系统指代锕系元素发生一次裂变反应而形成的裂变产物集合(如U-235热中子裂变,U-238快中子裂变等众多裂变系统)。裂变产物的衰变释热函数可表示为
{P}_{\text{fp}}(t)\approx {\displaystyle\int }_{0}^{T}\displaystyle\sum _{i\in \text{FS}}{R}_{i}(\tau ){P}_{i}(t-\tau )\text{d}\tau +\displaystyle\sum _{j\in \text{NIFP}}{Q}_{j}^{\text{decay}}{\lambda }_{j}{N}_{j}(t) (7) 式中:FS指代裂变系统组成的集合;T是积分时间;
{R_i}\left( \tau \right) 是第i裂变系统在\tau 时刻的裂变率;{P_i}(\tau ) 是第i裂变系统扣除NIFP核素贡献后,忽略中子辐照效应的衰变释热率。Q_j^{{\text{decay}}} 为第j裂变产物的衰变释热(MeV)。{\lambda _j} 为第j裂变产物的衰变常数(s−1);{N_j} 为第j裂变产物的核子密度(cm−3)。为了便于衰变热先驱核参数的推导和拟合,约定衰变热先驱核的衰变热为1 MeV,那么裂变产物的衰变释热函数可表示为
\begin{split} {P}_{\text{fp}}(t)\approx & {\displaystyle\int }_{0}^{T}\displaystyle\sum _{i\in \text{FS}}{R}_{i}(\tau )\displaystyle\sum _{j\in \text{DHP}}{\lambda }_{j}{\beta }_{i,j}{\text{e}}^{-{\lambda }_{j}(t-\tau )}\text{d}\tau +\displaystyle\sum _{k\in \text{NIFP}}{Q}_{k}^{\text{decay}}{\lambda }_{k}{N}_{k}(t)=\\ & {\displaystyle\int }_{0}^{T}\displaystyle\sum _{j\in \text{DHP}}{\lambda }_{j}\displaystyle\sum _{i\in \text{FS}}{R}_{i}(\tau ){\beta }_{i,j}{\text{e}}^{-{\lambda }_{j}(t-\tau )}\text{d}\tau +\displaystyle\sum _{k\in \text{NIFP}}{Q}_{k}^{\text{decay}}{\lambda }_{k}{N}_{k}(t)=\\ & \displaystyle\sum _{j\in \text{DHP}}{\lambda }_{j}{N}_{j}(t)+\displaystyle\sum _{k\in \text{NIFP}}{Q}_{k}^{\text{decay}}{\lambda }_{k}{N}_{k}(t) \end{split} (8) 式中:DHP指代衰变热先驱核构成的集合;
\;{ \beta _{i,j}} 指代第j衰变热先驱核对第i裂变系统的裂变产额;裂变产物的衰变热由NIFP核素与衰变热先驱核共同贡献,但首先需确定各衰变热先驱核的衰变常数与裂变份额。接下来讨论衰变热先驱核的参数拟合方法。首先选取典型裂变系统集合(Typical Fission Systems,TFS),考虑到各个系统之间的差异性,定义TFS中裂变产物衰变释热函数为
P_{\operatorname{mix}}(t)=\displaystyle\sum_{i \in \mathrm{TFS}} \omega_{i} P_{i}(t) (9) 式中:ωi是第i裂变系统的权重系数;Pi(t)是第i裂变系统中,扣除NIFP核素贡献的衰变释热率。
由式(8)可知,扣除NIFP核素贡献后,TFS中的衰变释热函数也可表示为
{\tilde{P}}_{\text{mix}}(t,{\boldsymbol{\lambda}} ,{\boldsymbol{\beta}} )=\displaystyle\sum _{i\in \text{DHP}}{\tilde{{\boldsymbol{\lambda}} }}_{i}{\tilde{{\boldsymbol{\beta}} }}_{i}{\text{e}}^{-{\tilde{{\boldsymbol{\lambda}} }}_{i}t} (10) 式中:
{\tilde{{\boldsymbol{\lambda}} }}_{i} 和{\tilde{\;{\boldsymbol{\beta}} }}_{i} 分别表示衰变热先驱核的衰变常数与裂变份额。1.4 参数拟合方法
式(10)为扣除NIFP核素贡献的TFS衰变释热拟合函数形式,只需拟合出函数中的参数即可完成衰变热先驱核的构造。本类搜索拟合问题可归类为非线性回归[12]问题,本质上是需要求解一组向量使得残差平方和函数值最小,即
\left\{ \begin{array}{l} \boldsymbol{x}=\operatorname{argmin}_{\boldsymbol{x}}\{\boldsymbol{R}(\boldsymbol{x})\} \\ \boldsymbol{R}(\boldsymbol{x})=\dfrac{1}{2} \displaystyle\sum_{i=1}^P\left(\mathrm{r}_i(\boldsymbol{x})\right)^2=\dfrac{1}{2}\|\boldsymbol{r}(\boldsymbol{x})\|^2=\dfrac{1}{2} \boldsymbol{r}(\boldsymbol{x})^{\top} \boldsymbol{r}(\boldsymbol{x}) .\end{array} \right. (11) 式中:P为考虑TFS衰变时间点的数目;x为待拟合的向量参数,在本文中x代表
{\tilde{{\boldsymbol{\lambda}} }}_{i} 和{\tilde{\;{\boldsymbol{\beta}} }}_{i} 的集合;{r_i}({\boldsymbol{x}}) 为第i个时间点下,拟合衰变释热函数与真实值的残差函数定义如下{\boldsymbol{r}}({\tilde {\boldsymbol{\lambda}} _i},{\tilde {\boldsymbol{\beta}} _i}) = \displaystyle\sum\limits_{i \in {\rm{DHP}}} {{{[({{{P}}_{{\mathrm{mix}}}}(t) - {\tilde {{P}}_{{\mathrm{mix}}}}(t,{{\tilde {\boldsymbol{\lambda}} }_i},{{\tilde {\boldsymbol{\beta}} }_i}))/{{{P}}_{{\mathrm{mix}}}}(t)]}^2}} (12) 本文采用Levenberg-Marquardt(L-M)算法[13],搜索衰变热先驱核参数,在该算法中,通过不断更新迭代步,逐步逼近获得最终解,最小逼近步长
\Delta 定义为\Delta = {{\boldsymbol{x}}_{n + 1}} - {{\boldsymbol{x}}_n} = - {({\boldsymbol{J}}_r^{\text{T}}{{\boldsymbol{J}}_r} + a{\boldsymbol{I}})^{ - 1}}{\boldsymbol{J}}_r^{\text{T}}{\boldsymbol{r}} (13) 式中:a为阻尼系数,在迭代过程中通过改变阻尼系数a可调整迭代策略[14],当a较大时,该方法遵循梯度下降法;而当a较小时,则更接近于高斯-牛顿法。
{{\bf{J}}_r} 为Jacobi矩阵{{\bf{J}}_r} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {r_0}}}{{\partial {{{x}}_0}}}}&{}&{ \cdots \cdots }&{}&{\dfrac{{\partial {r_n}}}{{\partial {{{x}}_{{n}}}}}} \end{array}} \right] (14) 搜索流程如图2所示。每次迭代首先计算最小逼近步长
\Delta ,更新搜索向量x后计算残差平方和函数值,若函数值小于精度要求{\varepsilon _{\min }} ,则输出解向量x;否则调整搜索策略继续迭代,直至残差平方和小于精度要求或迭代次数大于{\max _{{\text{iteration}}}} 。采用16组衰变热先驱核得到拟合结果如图3所示。图3(a)为拟合曲线与TFS中的衰变释热曲线的对比,从开始裂变到衰变1010 s(约300年)内,衰变释热率下降了近14个量级,拟合曲线与TFS中的衰变释热曲线吻合较好。图3(b)展示拟合误差随时间的分布,拟合误差的绝对值控制在1%以内,最大相对误差约为0.712%。
由于衰变常数具有全局性,因此,在TFS中搜索到
{\tilde {\boldsymbol{\lambda }}_i} 可作为各衰变热先驱核的衰变常数。此时第i裂变系统的衰变释热拟合函数形式为{\tilde{P}}_{i}(t,\beta )=\displaystyle\sum _{i\in {\mathrm{DHP}}}{\tilde{{\boldsymbol{\lambda}} }}_{i}{\tilde{{\boldsymbol{\beta}} }}_{i}{\text{e}}^{-{\tilde{{\boldsymbol{\lambda }}}}_{i}t} (15) 该问题已接近于线性回归问题,但考虑到数值稳定性,本文仍采用L-M算法。仿照上述思路,在各个裂变系统中搜索残差函数
{\boldsymbol{r}}({\tilde {\boldsymbol{\beta}} _i}) 的最小平方和,最终确定衰变热先驱核在各裂变系统中的裂变产额。2. 压缩流程
本文选取VERA燃料栅元[15]作为特征问题,燃料棒和慢化剂工作温度均为565 K,标记为VERA_1A。为简化建模,省略包壳,燃料棒半径为0.409 6 cm,栅距为1.26 cm。材料组分如表1所示,运行功率密度设定为40 W/gHM,最大燃耗深度为60 GWd/tHM。
表 1 VERA_1A燃料栅元初始材料明细Table 1. Specification of initial load materials in VERA_1Amaterial nuclides atom density/(1024 cm−3) fuel (3.1%) U234 6.11864 E-06U235 7.18132 E-04U236 3.29861 E-06U238 2.21546 E-02O16 4.57642 E-02moderator 565 K O16 2.48112 E-02H1 4.96224 E-02B10 1.07070 E-05B11 4.30971 E-05首先从中子学-燃耗耦合程序出发,记录各燃耗步下通量、截面和原子核密度场并展开衰变热计算,筛选出对衰变热有重要贡献的少数锕系核素。接着引入衰变假设,重新进行燃耗和衰变热计算,并将两次结果对比分析得到NIFP核素。目标核素的选取对基本压缩操作的评价有十分重要的影响,本文参考了一篇评估压水堆中核素重要性的报告[16],筛选出其中对临界安全以及源项计算有重要贡献的核素,并补充对衰变热有重要影响的少数锕系核素和NIFP核素。对精细燃耗库执行基本压缩操作并基于目标核素展开评价,目标核素的选取如表2所示,共计48种核素。
表 2 最终选取的目标核素Table 2. Final selected target nuclidestarget nuclides actinides fission products U-235 U-236 U-238 Kr-85 Sr-90 Y-90 Mo-95 Tc-99 Ru-101 U-239 Np-237 Np-239 Rh-103 Pd-105 Pd-108 Xe-131 Xe-135m Cs-133 Pu-238 Pu-239 Pu-240 Cs-134 Cs-134m Cs-137 Ba-137 Nd-145 Pm-147 Pu-241 Pu-242 Am-241 Pm-148m Pr-141 Nd-143 La-141 Sm-147 Sm-149 Am-243 Cm-242 Cm-244 Sm-150 Sm-151 Sm-152 Eu-153 Eu-154 Eu-157 Cm-246 Gd-155 Pm-150 原始精细燃耗数据库fine_lib数据来源于ENDF/B-VII.0,包含了中子反应截面数据、区分入射中子能量的裂变产额数据以及中子反应分支比和反应热数据。基于对基本压缩操作的重要性二元组分析及评价,经过重要性筛选,即可得到适用于中子学计算的压缩燃耗库。表3展示了设定不同的截断参数产生的压缩燃耗库明细,压缩后燃耗库命名方式为,VERA_XHX,其中第一个X表示燃耗数据库内锕系核素数量,第二个X表示燃耗库内核素总数。经过压缩的数据库规模明显小于精细燃耗库,锕系核素少于40种,核素总数少于180种,随着
{\varepsilon _1},{\varepsilon _2} 的增大,燃耗库内核素总数不断减少。表 3 压缩燃耗数据库明细Table 3. Specification of the compressed burnup libburnup lib {\varepsilon _1},{\varepsilon _2} number of
nuclidesnumber of fission
productsnumber of decay
channelsnumber of neutron
reaction channelsfine_lib 1547 1137 1341 1626 VERA_36H179 1E-6, 1E-3 179 142 138 364 VERA_36H151 5E-6, 1E-3 151 114 122 324 VERA_36H138 1E-5, 1E-3 138 101 119 298 VERA_34H133 1E-5, 1E-2 133 98 113 292 在上述流程的基础上,再将通过拟合参数构造的衰变热先驱核加入压缩燃耗库中,使压缩燃耗库最终适用于衰变热计算。
3. 燃料栅元问题计算分析
本文从有效增殖因子计算、燃耗计算和衰变热计算三个方面对压缩燃耗库展开验证,测试基准题为VERA燃料栅元的四种工况。1A为燃耗库压缩过程中的特征问题工况,1B、1C、1D相较于1A采用不同硼浓度的慢化剂,并且慢化剂工作温度均为600 K,燃料棒温度分别为600 K、900 K、1 200 K。所有问题均采用精细燃耗库与压缩燃耗库进行计算,由于各压缩燃耗库的结果重合度过高,因此不展示VERA_36H151和VERA_36H138的结果。以蒙卡程序OpenMC[17]对各工况进行燃耗计算的结果作为基准解,计算共分14个燃耗步,最大燃耗深度为60 GWd/tHM。表4统计了燃耗深度从0变化至60 GWd/tHM的燃耗过程中,精细燃耗库与压缩燃耗库相较于OpenMC有效增殖因子计算偏差的最大值与均方根水平。可以看到,压缩前后有效增殖因子的计算偏差水平相近。因此,燃耗库的压缩不会引入明显的偏差。
表 4 有效增殖因子计算偏差的最大值与均方根Table 4. Maximum value and root mean square (RMS) value of keff deviationmaximum keff deviation/10−5 RMS of keff deviation/10−5 fine_lib VERA_36H179 VERA_34H133 fine_lib VERA_36H179 VERA_34H133 VERA_1A 202.07 201.79 200.77 93.72 93.64 92.40 VERA_1B 280.45 281.39 280.30 109.65 108.88 107.94 VERA_1C 273.52 274.47 273.39 115.27 114.90 113.87 VERA_1D 290.05 293.17 292.35 123.77 124.16 123.13 不同燃耗深度下,目标核素核子密度相较于基准解的计算相对偏差均方根如表5统计。当采用较高精度的压缩燃耗库计算时,核子密度计算精度未被影响,均方根小于2%;而采用低精度压缩燃耗库时,计算精度会明显降低,误差会随燃耗深度的增加累计。在燃耗步末,燃耗深度达到60 GWd/tHM时,均方根大于5%。偏差来源于目标核素Gd-155,由于Gd-155主要由Gd-154通过中子俘获产生,而Gd-154在压缩过程中被删除,导致Gd-155的计算结果始终偏负,其余核素的相对偏差均方根小于2%,与精细燃耗库相近。
表 5 目标核素核子密度计算误差均方根值Table 5. Root mean square (RMS) value of target nuclides densities deviationRMS/% fine_lib VERA_36H179 VERA_34H133 0.8
GWd/tHM20.0
GWd/tHM60.0
GWd/tHM0.8
GWd/tHM20.0
GWd/tHM60.0
GWd/tHM0.8
GWd/tHM20.0
GWd/tHM60.0
GWd/tHMVERA_1A 1.77 1.69 1.97 1.14 1.12 2.91 1.33 1.35 5.46 VERA_1B 1.74 1.72 2.00 1.37 1.40 3.15 1.52 1.54 5.85 VERA_1C 1.74 1.71 1.98 1.49 1.52 3.24 1.65 1.66 5.86 VERA_1D 1.72 1.74 2.01 1.60 1.63 3.28 1.67 1.69 5.71 若将Gd-154纳入到目标核素中经压缩操作后的燃耗库,可保留Gd-155的计算精度。VERA_34H133的重要性二元组分别为10−5与10−2,包含了133种核素,对Gd-155的计算结果如图4(a)所示。在目标核素中加入Gd-154后,基于相同重要性二元组压缩后的燃耗库计算,Gd-155的计算结果如图4(b)所示。可以看到,改进后燃耗数据库与精细燃耗数据库的计算结果基本吻合。因此,本文提出的方法能有效保留目标核素核子密度的计算精度。
表6对比了VEAR燃料栅元基于精细燃耗数据库与压缩燃耗数据库,进行25 d燃耗计算的计算内存开销与计算时间。当采用精细燃耗数据库计算时,需要约4 MB的内存,而经过压缩后占用的内存约0.5 MB,相比之下缩减了一个数量级。计算时间方面,从原先的160 ms缩短至不足20 ms。由此可见,燃耗数据库的压缩可有效提升计算效率和节约计算储存开销。
表 6 计算内存开销与计算时间Table 6. Calculating memory overhead and computation timecalculating memory overhead/MB computation time/ms fine_lib 4.060 160 VERA_36H179 0.535 16 VERA_34H133 0.407 12 在不同燃耗深度状态点下,展开1010 s(约300年)的衰变热计算,图5(a)和图5(b)分别展示了VERA_1A中在20 GWd/tHM和60 GWd/tHM下的衰变释热曲线。展示了在加入衰变热先驱核后的衰变释热曲线与基准解吻合较好。在106 s时,衰变释热率从大于200 W衰减至不足10 W,在这之间的误差在5%上下浮动,表明了可通过拟合衰变热先驱核,使裂变产物衰变热计算在保证精度的前提下得以简化。
衰变释热率对燃料栅元功率贡献的分析计算结果如表7所示。随着燃耗深度的增加,衰变热对功率的贡献会减小,但始终超过6%。与精细燃耗库结果相比,加入衰变热先驱核的压缩燃耗库计算偏差小于0.5%,最大偏差为0.48%。本文提出的拟合衰变热先驱核方法,在保证精度的前提下简化了裂变产物的衰变热计算。
表 7 衰变热对功率的贡献Table 7. Decay heat contribution to the powerdecay heat contribution/% fine_lib VERA_36H179_16 deviation fine_lib VERA_36H179_16 deviation fine_lib VERA_36H179_16 deviation 0.8 GWd/tHM 20.0 GWd/tHM 60.0 GWd/tHM VERA_1A 6.37 6.00 −0.37 6.14 5.79 −0.35 6.16 5.69 −0.47 VERA_1B 6.39 6.01 −0.38 6.39 6.01 −0.38 6.18 5.70 −0.48 VERA_1C 6.40 6.02 −0.38 6.40 6.02 −0.38 6.16 5.70 −0.46 VERA_1D 6.40 6.02 −0.38 6.40 6.02 −0.38 6.16 5.70 −0.46 4. 结 论
本文提出了一种基于伪衰变核形成的燃耗数据库压缩方法。实现了基本压缩操作的重要性二元组评价和衰变热先驱核参数的拟合搜索,形成两点优势:(1)除了评价中子学计算精度,还考虑了目标核素选取,从而可根据用户需求定制压缩燃耗库;(2)克服了燃耗数据库压缩后衰变热计算精度损失的问题。VERA燃料栅元测试算例表明,经压缩后的燃耗数据库对有效增殖因子的计算不会引入明显偏差;采用精度较高的压缩燃耗库计算,目标核素的核子密度计算相对偏差均方根小于2%;在较低精度的燃耗库中,仅个别核素的计算结果发生偏离;在加入衰变热先驱核后,衰变热计算精度可以有效复原,衰变释热率对燃料栅元功率贡献的计算偏差小于0.5%。本文提出的方法具有大压缩、高精度特点,可有效减小计算开销,保留计算精度。后续研究拟针对以下方面展开:(1)考虑各基本压缩操作之间的干涉效应,进一步提高燃耗数据库压缩精度;(2)赋予衰变热先驱核中子学特性,进一步压缩燃耗库规模。
-
表 1 VERA_1A燃料栅元初始材料明细
Table 1. Specification of initial load materials in VERA_1A
material nuclides atom density/(1024 cm−3) fuel (3.1%) U234 6.11864 E-06U235 7.18132 E-04U236 3.29861 E-06U238 2.21546 E-02O16 4.57642 E-02moderator 565 K O16 2.48112 E-02H1 4.96224 E-02B10 1.07070 E-05B11 4.30971 E-05表 2 最终选取的目标核素
Table 2. Final selected target nuclides
target nuclides actinides fission products U-235 U-236 U-238 Kr-85 Sr-90 Y-90 Mo-95 Tc-99 Ru-101 U-239 Np-237 Np-239 Rh-103 Pd-105 Pd-108 Xe-131 Xe-135m Cs-133 Pu-238 Pu-239 Pu-240 Cs-134 Cs-134m Cs-137 Ba-137 Nd-145 Pm-147 Pu-241 Pu-242 Am-241 Pm-148m Pr-141 Nd-143 La-141 Sm-147 Sm-149 Am-243 Cm-242 Cm-244 Sm-150 Sm-151 Sm-152 Eu-153 Eu-154 Eu-157 Cm-246 Gd-155 Pm-150 表 3 压缩燃耗数据库明细
Table 3. Specification of the compressed burnup lib
burnup lib {\varepsilon _1},{\varepsilon _2} number of
nuclidesnumber of fission
productsnumber of decay
channelsnumber of neutron
reaction channelsfine_lib 1547 1137 1341 1626 VERA_36H179 1E-6, 1E-3 179 142 138 364 VERA_36H151 5E-6, 1E-3 151 114 122 324 VERA_36H138 1E-5, 1E-3 138 101 119 298 VERA_34H133 1E-5, 1E-2 133 98 113 292 表 4 有效增殖因子计算偏差的最大值与均方根
Table 4. Maximum value and root mean square (RMS) value of keff deviation
maximum keff deviation/10−5 RMS of keff deviation/10−5 fine_lib VERA_36H179 VERA_34H133 fine_lib VERA_36H179 VERA_34H133 VERA_1A 202.07 201.79 200.77 93.72 93.64 92.40 VERA_1B 280.45 281.39 280.30 109.65 108.88 107.94 VERA_1C 273.52 274.47 273.39 115.27 114.90 113.87 VERA_1D 290.05 293.17 292.35 123.77 124.16 123.13 表 5 目标核素核子密度计算误差均方根值
Table 5. Root mean square (RMS) value of target nuclides densities deviation
RMS/% fine_lib VERA_36H179 VERA_34H133 0.8
GWd/tHM20.0
GWd/tHM60.0
GWd/tHM0.8
GWd/tHM20.0
GWd/tHM60.0
GWd/tHM0.8
GWd/tHM20.0
GWd/tHM60.0
GWd/tHMVERA_1A 1.77 1.69 1.97 1.14 1.12 2.91 1.33 1.35 5.46 VERA_1B 1.74 1.72 2.00 1.37 1.40 3.15 1.52 1.54 5.85 VERA_1C 1.74 1.71 1.98 1.49 1.52 3.24 1.65 1.66 5.86 VERA_1D 1.72 1.74 2.01 1.60 1.63 3.28 1.67 1.69 5.71 表 6 计算内存开销与计算时间
Table 6. Calculating memory overhead and computation time
calculating memory overhead/MB computation time/ms fine_lib 4.060 160 VERA_36H179 0.535 16 VERA_34H133 0.407 12 表 7 衰变热对功率的贡献
Table 7. Decay heat contribution to the power
decay heat contribution/% fine_lib VERA_36H179_16 deviation fine_lib VERA_36H179_16 deviation fine_lib VERA_36H179_16 deviation 0.8 GWd/tHM 20.0 GWd/tHM 60.0 GWd/tHM VERA_1A 6.37 6.00 −0.37 6.14 5.79 −0.35 6.16 5.69 −0.47 VERA_1B 6.39 6.01 −0.38 6.39 6.01 −0.38 6.18 5.70 −0.48 VERA_1C 6.40 6.02 −0.38 6.40 6.02 −0.38 6.16 5.70 −0.46 VERA_1D 6.40 6.02 −0.38 6.40 6.02 −0.38 6.16 5.70 −0.46 -
[1] Chadwick M B, Obložinský P, Herman M, et al. ENDF/B-VII. 0: next generation evaluated nuclear data library for nuclear science and technology[J]. Nuclear Data Sheets, 2006, 107(12): 2931-3060. doi: 10.1016/j.nds.2006.11.001 [2] Aldama D L, Leszczynski F, Trkov A. WIMS-D library update[R]. STI/PUB-1264, Vienna: IAEA, 2007. [3] 刘燕平, 李世清. 裂变产物中子数据库重要核素的分析[J]. 武汉大学学报(自然科学版), 1990(2):56-58Liu Yanping, Li Shiqing. The analyses of important nuclides for fission products neutron nuclear data library (FPNNDL)[J]. Journal of Wuhan University (Natural Science Edition), 1990(2): 56-58 [4] Katano R, Yamamoto A, Endo T, et al. Generation of simplified burnup chain using contribution matrix of nuclide production[C]//Proceedings of the PHYSOR2014—The Role of Reactor Physics towards a Sustainable Future. 2014. [5] Kajihara T, Tsuji M, Chiba G, et al. Automatic construction of a simplified burn-up chain model by the singular value decomposition[J]. Annals of Nuclear Energy, 2016, 94: 742-749. doi: 10.1016/j.anucene.2016.04.034 [6] Chiba G, Tsuji M, Narabayashi T, et al. Important fission product nuclides identification method for simplified burnup chain construction[J]. Journal of Nuclear Science and Technology, 2015, 52(7/8): 953-960. [7] 胡钰莹, 廖鸿宽, 姚栋, 等. 双约束核素筛选与燃耗链压缩算法研究[J]. 核动力工程, 2022, 43(4):218-222Hu Yuying, Liao Hongkuan, Yao Dong, et al. Research on dual-restricted nuclide selection and burnup chain compression algorithm[J]. Nuclear Power Engineering, 2022, 43(4): 218-222 [8] 黄凯, 吴宏春, 李云召, 等. 基于定量重要性分析的燃耗链压缩方法[J]. 强激光与粒子束, 2017, 29:036002 doi: 10.11884/HPLPB201729.160302Huang Kai, Wu Hongchun, Li Yunzhao, et al. Depletion chain compression method via quantitative significance analysis[J]. High Power Laser and Particle Beams, 2017, 29: 036002 doi: 10.11884/HPLPB201729.160302 [9] 谢仲生. 核反应堆物理分析[M]. 西安: 西安交通大学出版社, 2004: 165-172Xie Zhongsheng. Nuclear reactor physics analysis[M]. Xi'an: Xi'an Jiaotong University Press, 2004: 165-172 [10] Zhang Tengfei, Xiao Wei, Yin Han, et al. VITAS: a multi-purpose simulation code for the solution of neutron transport problems based on variational nodal methods[J]. Annals of Nuclear Energy, 2022, 178: 109335. doi: 10.1016/j.anucene.2022.109335 [11] Pusa M, Leppänen J. Computing the matrix exponential in burnup calculations[J]. Nuclear Science and Engineering, 2010, 164(2): 140-150. doi: 10.13182/NSE09-14 [12] Al-Baali M, Fletcher R. Variational methods for non-linear least-squares[J]. Journal of the Operational Research Society, 1985, 36(5): 405-421. doi: 10.1057/jors.1985.68 [13] Madsen K, Nielsen H B, Tingleff O. Methods for non-linear least squares problems[M]. 2nd ed. Copenhagen: Technical University of Denmark, 2004: 24-28. [14] Gavin H P. The Levenberg-Marquardt algorithm for nonlinear least squares curve-fitting problems[D]. Durham: Duke University, 2019. [15] Godfrey A T. VERA core physics benchmark progression problem specifications[R]. CASL-U-2012-0131-004, 2014. [16] Gauld I C, Ryman J C. Nuclide importance to criticality safety, decay heating, and source terms related to transport and interim storage of high-burnup LWR fuel[R]. NUREG/CR-6700, Oak Ridge: Oak Ridge National Laboratory, 2000. [17] Romano P K, Forget B. The OpenMC Monte Carlo particle transport code[J]. Annals of Nuclear Energy, 2013, 51: 274-281. doi: 10.1016/j.anucene.2012.06.040 -