一种强中子场裂变产物核燃耗的计算方法

摘要

本发明涉及一种强中子场裂变产物核燃耗的计算方法，该方法将裂变产生的产物核的核反应网络拆分为电荷数Z＝22‑72的51个子网，并建立裂变产物核的子网数据文件，根据裂变产物核的子网数据文件，通过隐式龙格‑库塔方法求解各裂变产物核的核密度随时间的变化函数。本发明可以显著的提高强中子场裂变产物核燃耗的计算效率。

权利要求

1.一种强中子场裂变产物核燃耗的计算方法，包括如下步骤：(1)将裂变产生的产物核的核反应网络拆分为电荷数Z＝22至72的51个子网；

(2)建立裂变产物核的子网数据文件；

(3)根据裂变产物核的子网数据文件，求解各裂变产物核的核密度随时间的变化函数，裂变产物核i核密度变化微分方程如下：其中，

表示为第i个产物核素密度随时间变化；

右边第一项为单位时间裂变产生i核的个数，f(t)为裂变率,yi为裂变产额；

右边第二项为i核的消耗项，n(t)为中子注量率,σout为谱平均消失截面，包括该核的(n,γ)、(n,2n)、(n,3n)反应截面；

右边第三项为i核的产生项，分别为其它核j通过(n,γ)、(n,2n)、(n,3n)反应产生i核的贡献，相应的Nj为其它产物核j核密度，σin为相应的平均截面。

2.如权利要求1所述的强中子场裂变产物核燃耗的计算方法，其特征在于：步骤(2)中所述的子网数据文件包括相同电荷的裂变产物核的ID和核数据，其中，对于电荷数Z，质量数A，同核异能态I，其ID＝A×1000+Z×10+I；核数据包括(n,γ)、(n,2n)、(n,3n)的反应截面、独立产额Y。

3.如权利要求1或2所述的强中子场裂变产物核燃耗的计算方法，其特征在于：步骤(3)中通过隐式龙格-库塔方法求解各裂变产物核的核密度随时间的变化函数。

说明书

一种强中子场裂变产物核燃耗的计算方法

技术领域

[0001] 本发明涉及裂变产物核的燃耗计算方法，具体涉及一种强中子场裂变产物核燃耗的计算方法。

背景技术

[0002] 在高密度的强中子场中，裂变产物核可能发生可观的二次反应，使得产物核的产额(或积存量)发生改变，因此需要通过燃耗计算，来确定产额的改变量。

[0003] 目前除了国内钱晶的博士论文编写了FIRENEQ计算程序，尚无发现有公开资料介绍强中子场燃耗的计算方法。与本发明接近的是反应堆的燃耗计算程序。反应堆燃耗计算方法主要有两大类：基于燃耗方程组矩阵求解的数值计算方法和基于单燃耗链的解析方法。如洛斯阿拉莫斯实验室在上世纪80年代开发的Origen2程序，是基于泰勒展开的解析方法。基于数值计算的芬兰VTT中心开发的Serpent系统中的燃耗计算模块、法国开发的燃耗程序MENDEL、欧洲开发的活化计算系统中的核素存量计算程序FISPACT，关于这方面的介绍，可以参看吴明宇的博士论文《基于STEP1.0和MCMG-II的输运-燃耗耦合计算系统的开发与研究》。下面介绍的是与本发明最为接近的基于解析法的CINDER90程序和基于数值方法的FIRENEQ计算程序。

[0004] CINDER90程序通过对核反应网络进行拆分，得到单个核素链，然后通过解析方法得到每个步长内的核素密度的变化。其关键问题之一是反应网络的拆分，当出现循环以及循环嵌套时，拆分比较困难；关键问题二是确保在步长时间内中子通量的近似不变。通过试算，CINDER90的步长不能小于0.1ns，而在这个步长内，强中子场的中子通量变化很大，所以CINDER90不适合于强中子场的燃耗计算，计算结果会出错。FIRENEQ程序采用数值的GEAR算法，对整个核反应网络和局部网络进行计算。裂变产额生的产物核有1300多个，产物核的参与的反应有(n，γ)、(n，2n)、(n，3n)、(n，p)、(n，d)、(n，t)、(n，3He)和(n，α)等，是一个1300x1300的矩阵数值求解，计算一次需要几个小时，甚至需要1天或几天。

发明内容

[0005] 本发明的目的是提供一种强中子场裂变产物核燃耗的计算方法，用于解决强中子场裂变产物核的燃耗修正，并有效地提高计算效率。

[0006] 本发明的技术方案如下：一种强中子场裂变产物核燃耗的计算方法，包括如下步骤：

[0007] (1)将裂变产生的产物核的核反应网络拆分为电荷数Z＝22-72的51个子网；

[0008] (2)建立裂变产物核的子网数据文件；

[0009] (3)根据裂变产物核的子网数据文件，求解各裂变产物核的核密度随时间的变化函数，裂变产物核i核密度变化微分方程如下：

[0010]

[0011] 其中，

[0012] 表示为第i个产物核素密度随时间变化；

[0013] 右边第一项为单位时间裂变产生i核的个数，f(t)为裂变率，yi为裂变产额；

[0014] 右边第二项为i核的消耗项，n(t)为中子注量率，σout为谱平均消失截面，包括该核的(n，γ)、(n，2n)、(n，3n)反应截面；

[0015] 右边第三项为i核的产生项，分别为其它核j通过(n，γ)、(n，2n)、(n，3n)反应产生i核的贡献，相应的Nj为其它产物核j核密度，σin为相应的平均截面。

[0016] 进一步，如上所述的强中子场裂变产物核燃耗的计算方法，步骤(2)中所述的子网数据文件包括相同电荷的裂变产物核的ID和核数据，其中，对于电荷数Z，质量数A，同核异能态I，其ID＝A×1000+Z×10+I；核数据包括(n，γ)、(n，2n)、(n，3n)的反应截面、独立产额Y。

[0017] 进一步，如上所述的强中子场裂变产物核燃耗的计算方法，步骤(3)中通过隐式龙格-库塔方法求解各裂变产物核的核密度随时间的变化函数，所述的隐式龙格-库塔方法采用PERL代码软件包MATH-ODE，用PERL语言编写CABOR程序实现。

[0018] 本发明的有益效果如下：采用本发明的计算方法，建立相应的模型，可以快速地进行计算，计算一次只需要几分钟。速度的提高，得益于整个网络的分割。裂变产物大概有1300多个，求解整个网络，需要求解矩系数矩阵为1300*1300的微分方程组；分割后成51个子网后，每个子网的产物核只有十几个，以Z＝58的子网为例，其核素只有11个，其矩阵为
11*11，约为整个网络的14000分之一。循环所有51个子网约为整网的280分之一。根据这样粗略的估计，时间可以缩短为整个网络事件的1/280。如果求解整个网络需要24小时，那么单个子网的时间约为6秒，所有51个子网的时间累加约为5分钟。可见，本发明可以显著的提高强中子场裂变产物核燃耗的计算效率。

附图说明

[0019] 图1为任意裂变产物核在强中子场下生成与消失的可能反应路径图；

[0020] 图2为由裂变产物核、核反应、核衰变组成的核反应网络示意图，节点表示产物核，箭头表示可能的核反应(参考图1)，图中给出的是质量数范围A＝118-162的网络，虚线带箭头的是Z＝58的子网拆分线；

[0021] 图3为Ce(Z＝58，A＝142-153)同位素反应子网络示意图，中间向右表示(n，γ)反应，该子网包含许多循环，适用于数值求解；

[0022] 图4为本发明的方法用PERL语言编写CABOR程序实现的流程图；

[0023] 图5为本发明的具体实施方案的流程图。

具体实施方式

[0024] 下面结合附图和实施例对本发明进行详细的描述。

[0025] 在强中子场中，i产物核密度变化微分方程为：

[0026]

[0027] 其中， 表示为第i个产物核素密度[单位：/b-cm]随时间变化，右边第一项为单位时间裂变产生i核的个数，f(t)为裂变率[1/s]，yi为裂变产额；右边第二项为i核的消耗项，n(t)为中子注量率[1/cm2-s]，σout为谱平均消失截面，主要包括该核的(n，γ)、(n，2n)、(n，3n)反应截面；右边第三项为i核的产生项，分别为其它核(j)通过(n，γ)、(n，2n)、(n，3n)等反应产生i核的贡献，相应的Nj为其它产物核核密度，σin为相应的平均截面。

[0028] 每个产物核可以通过其它核的衰变或中子反应产生，也可以通过本身的衰变和核反应消失。这些衰变和核反应包括β-衰变、(n，γ)、(n，2n)、(n，3n)、(n，p)、(n，d)、(n，t)、(n，3He)、(n，a)反应等，如图1所示。裂变产生1300多个产物核，其核素和反应组成的核反应网络如图2所示，质量数的范围为A＝66-172，电荷数Z＝22-72。

[0029] 强中子场的持续时间较短，远远小于衰变过程，衰变过程可以忽略，微分方程(1)式去掉了衰变过程，其它中子反应，如(n，p)、(n，d)、(n，t)、(n，3He)、(n，a)等反应截面远远小于(n，γ)、(n，2n)、(n，3n)反应截面之和，因此也可以忽略。按照上述假设，裂变产物核参加的核过程，电荷数不变，相应的核反应网络方程可以拆分成相同电荷数的子网，整网可以拆分为Z＝22-72的51个子网。拆分前的网络范围为A＝72-172，Z＝22-72，图2(示意图)给出了Z＝54-62(Xe-Sm)的部分网络。拆分的方法是在Z＝22.5，23.5，...，71.5的地方横向切开，(图2中带箭头的虚线为Z＝58的子网拆分线)，得到Z＝22，23，...，72的子网，图3为拆分后的Z＝58的子网，并标出了(n，γ)、(n，2n)、(n，3n)的核反应；如图3是拆分后Ce(Z＝58)同位素组成的子网。Ce同位素的反应子网，存在诸多循环的复杂网络，因此很难把该子网拆分成无循环的单线网络，无法用解析的方法给出核密度的随时间的变化函数N(t)。但是，由于该子网只有11个产物核，所以可以用数值解法进行求解，核素少，可以避免速度慢和刚性问题。常用的隐式的龙格-库塔方法可以满足要求。隐式龙格-库塔方法有现成的代码，本发明采用公开免费的PERL代码软件包MATH-ODE(http：//search.cpan.org/dist/MATH-ODE-0.07/)。

[0030] 生成核反应子网，包括数据文件和代码，本发明的技术方案如下(如图4所示)：

[0031] 1)从数据文件library读入相同电荷的裂变产物核的ID和核数据(包括截面和产额)；

[0032] a)对于某一电荷Z，质量数A，同核异能态I，其身份代码为ID＝A*1000+Z*10+I；

[0033] b)核数据包括(n，g)，(n，2n)和(n，3n)截面、独立产额Y；

[0034] 2)根据产物核ID生成子网嵌入式代码Pu239_Z.code，Z表示电荷数，以Z＝58为例，生成的代码如表1行B1-B27所示，信息包括ID，产额，微分代码；

[0035] 3)主代码中利用require(“Pu239_$Z.code”)，引入子网嵌入式代码，并进行求解；

[0036] 4)得到计算结果。

[0037] 基于以上原理，用PERL语言编写了CABOR程序，流程图如图5所示。主要步骤如下：

[0038] A.参考CINDER90的数据格式建立输入数据表，文件名为library；

[0039] B.根据library文件，根据相同电荷的产物核建立子网，生成子网数据文件OD-Z.dat，Z＝51-72.该文件包括反应截面和产额；

[0040] C.生成子网嵌入式代码Pu239_Z.code，Z＝51-72，该代码采用PERL语言，可以嵌入主程序；

[0041] D.执行主程序Cabor.pl，进行计算，代码见表1；

[0042] E.计算生成输出文件sum.out，见表2；

[0043] F.结束。

[0044] 关键语句如表1所示。以求解Z＝58、Ce133-Ce162核素密度为例，表1列出了核心语句和模块。行A5插入子网代码(B1-B27)。行A7-A16初始化微分方程组。行A18-29是求解微分方程组的循环，每个循环完成一个燃耗步长，燃耗时间到达给定的时间即停止。这样就可以得到不同时刻的核密度，然后换算成产额，即可得到产额的变化数据，表2为输出的例子。

[0045] 通过以上方法，实现了网络分割和子网的数值求解，最后得到产额的变化数据。

[0046] 表1 Cabor.pl主要代码(A)和PU239_Z.code嵌入式代码(B)

[0047]

[0048]

[0049] 表2 Ce的计算结结果文件sum.out

[0050]

[0051]

[0052] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若对本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。