基于衰变链计算的裂变产物缓发γ能谱分析

时间：2024-07-28

田晨扬，过惠平，欧阳晓平，许鹏，张立波，白德平，田庆云，吕汶辉

(1.火箭军工程大学，陕西西安 710025；2.96607部队，陕西宝鸡 721000；3.西北核技术研究所，陕西西安 710000)

在核军控领域，探测未知材料中235U是否存在以及存在多少数量，一直是一技术难题[1-3]。为获取铀材料中235U的含量，需对铀裂变产物的缓发γ能谱进行分析[4-6]。对缓发γ能谱相关的测量与分析成果可应用于许多领域，如核武器核查、核材料安全管控、海关货物、核反应堆工程及退役核设施的辐射监测等领域[7-10]。

在军控核查技术中，当采用被动法探测铀材料的235U质量丰度时，由于235U的185.7 keV特征γ射线能量较低，很容易被高原子序数物质和铀的自吸收所屏蔽，导致测量十分困难[1-5]。本文采用238U裂变阈以下的热中子主动诱发铀材料裂变，结合裂变产物核素的放射性衰变数据和裂变产额数据，计算分析裂变产物发出的缓发γ射线，得到缓发γ能谱。最终，从对应缓发γ能谱选择特征峰来表征235U的存在性及丰度，确定被测材料是否含浓缩铀。在实际核查中，只要探测到特征峰，即说明被测物质中有235U存在；只要特征峰的计数达到特定数值范围，即可判断所测的铀材料是浓缩铀。

1 铀裂变产物的衰变链计算推导

当铀裂变产物放射出瞬发中子和γ射线后，由于具有较高的中子-质子比，继续经历一系列β衰变，最后变为稳定的核素[11-13]。通过β衰变，一个核素转变为另一个核素的激发态，当激发态跃迁到基态时，放射出缓发γ射线。通过β衰变和电子俘获，核素原子序数或增或减，但质量数保持不变[14-16]。这些质量数彼此相等的核素构成一个衰变链，铀裂变产物构成若干条衰变链，质量数从66到172不等。

图1示出1条简单衰变链和1条复杂衰变链，其中，115Ag具有特定的母体和复杂的分支衰变。

1.1 物理模型

在用热中子照射材料时，照射结束时刻称为零时[18-19]。为得到缓发γ能谱，需确定零时前后各级核素数目的变化规律。

根据文献[9]，可将1个复杂衰变链划分为几个线性衰变链，在线性衰变链中，每个核素只有1个母核，没有分支衰变。对于某个核素，通过叠加它在每个线性衰变链上的分支数目，可得到总数目。为简化物理模型，将衰变链分为基态衰变链和激发态衰变链，前者包含1个处于基态的母核，而后者包含1个处于激发态的母核。

图1 典型衰变链[17]Fig.1 Typical decay chain[17]

图1a中基态衰变链和激发态衰变链分别如下：

其中，IT表示同质异能跃迁。

假设样品材料中234U、235U、236U和238U的质量分数分别为p1、p2、p3和p4，不考虑样品形状，某种铀同位素的初始数目为Ni，0，则有：

(1)

其中：m为样品质量；M和NA分别为靶核的原子质量和阿伏伽德罗常数。

以235U为例，当00,0≤k0，0≤k

假设入射中子的能量单一，靶核的核反应率如下：

N=N0e-ut

(2)

1.2 零时前核素数目的变化规律

1) 基态衰变链的核素数目

推导零时前基态衰变链的核素数目。根据以上假设，母核的衰变率如下：

(3)

N1=φσfN0b10(e-ut-e-λ1t)

(4)

第1个子核的数目由235U裂变、前级衰变和自身衰变3种方式贡献：

N2=φσfN0(b20(e-ut-e-λ2t)+

b21(e-λ2t-e-λ1t))

(5)

(6)

(7)

根据数学归纳法，可得各级核素的数目：

Nj=φσfN0(bj0(e-ut-e-λjt)+

bj1(e-λjt-e-λ1t)+…+bjk(e-λjt-e-λkt)+

…+bj,j-1(e-λjt-e-λj-1t))j>1

(8)

其中：

(9)

(10)

(11)

式(8)中，令t=t0，可得基态衰变链各级核素的零时刻数目。

2) 激发态衰变链的核素数目

(12)

(13)

(14)

(15)

(16)

(17)

利用同样方法，可得：

j>1

(18)

其中：

(19)

(20)

(21)

式(18)中，令t=t0，可得激发态衰变链各级核素的零时刻数目。

1.3 零时后核素数目的变化规律

在衰变链中，每个核素数目的增加和减少过程可用Bateman方程[17]表示：

(23)

根据Bateman方程，假设前3个核素A、B、C的衰变常量分别为λ1、λ2和λ3，在t0时刻的数目分别为N1(t0)、N2(t0)和N3(t0)，在t(t>t0)时刻的数目分别为N1(t)、N2(t)和N3(t)。

1) 核素A、B、C数目的变化规律

随着时间的增加，核素A呈指数衰减：

N1(t)=c10e-λ1(t-t0)

c10=N1(t0)

(24)

核素B以速率r1λ1N1由核素A产生，同时以速率λ2N2生成核素C。

N2(t)=c20e-λ1(t-t0)+c21e-λ2(t-t0)

(25)

因此，零时后衰变链上核素j的数目Nj(t)如下：

Nj(t)=cj0e-λ1(t-t0)+…+

cjke-λk+1(t-t0)+…+cj,j-1e-λj(t-t0)

(26)

(27)

(28)

(29)

2) 衰变链上任意核素的缓发γ射线数目

(1) 衰变校正

在测量时间t2-t1过程中，令任意核素j发生衰变的原子数目为nj(j=1，2，3，…)。对于母核，有：

N1(t0)(e-λ1t1-e-λ1t2)

(30)

对于子核，有：

(31)

(32)

同一个核素可能属于多个线性衰变链，所以，应对每个线性衰变链的衰变数目进行叠加。

(2) γ 发射概率校正

2 缓发γ能谱的Monte Carlo计算

2.1 主程序结构

使用Monte Carlo程序MCNP建立缓发γ能谱的物理模型，利用衰变数据和裂变产额数据[20]进行具体计算。

主程序结构如图2所示。

1) 输入计算参数，读取裂变产物质量数，从文件CHAIN.txt读取每个核素的信息，包括质子数、衰变模式、分支比和半衰期等。

2) 在1个衰变链中，如果1个核素的质子数最小，则以它为母核建立1个β-基态链；如果1个核素的质子数最大，则以它为母核建立1个EC基态链。如果核素以激发态存在，则以它为母核建立激发态链。

3) 根据入射中子类型(热中子、裂变中子或高能中子)，读取每个铀同位素(234U、235U、236U、238U)的裂变产额数据。

4) 从文件DECAY.txt读取γ能量和发射概率数据，利用前文推导的公式，计算对应衰变链上每个核素的衰变原子数目，获得不同能量的缓发γ射线计数。

图2 主程序结构Fig.2 Structure of main program

5) 按照质量数A从66到172循环操作以上步骤，根据特定算法，对所有裂变产物的缓发γ射线数据进行排序，将γ能谱数据输出到文件SPECT.txt中。

2.2 实验与计算结果对比

利用脉冲反应堆照射铀材料(235U含量90%)，照射时间为10 min，235U靶质量为0.3 mg，热中子通量为1×1012cm-2·s-1。利用跑兔装置将照射后的样品直接打入铅屏蔽体后进行测量，铅屏蔽体的厚度为10 cm，内衬的铜和镉板厚度均为1 mm，有机玻璃厚度为5 mm。采用GMX10180同轴N型HPGe探测器，样品在距探测器表面35 cm处测量1 h，死时间比例为35.95%。测量得到的缓发γ能谱如图3所示。

经过测量，选择几个典型的裂变产物获取相应的γ计数，根据相同的测量条件，使用本文所建程序进行模拟计算，缓发γ计数的实验结果与模拟计算结果的对比列于表1。从表1可知，实验结果与模拟计算结果的相对偏差在-2.25%～-5.67%之间。模拟计算结果较实验结果偏小，原因如下：1) 对于某一裂变产物特定能量的γ射线峰，实验测得的计数为此能峰的净面积，这就可能将不同核素的能量很接近的缓发γ射线一起包含在内，而理论计算值只是对应单一核素特定能量γ射线计数；2) 在建立衰变链时，理论计算没有考虑分支比很小的母核基态向子核激发态衰变的线性链，导致处在这种链上的核素数目少于实际数目，使计算值小于实验值；3) 在推导公式时，未考虑瞬发中子的诱发裂变、裂变产物核素发生(n,γ)反应等情况，这些情况也产生缓发γ射线。

图3 缓发γ能谱测量结果Fig.3 Measurement result of delayed gamma spectrum

核素名称缓发γ射线能量/keV缓发γ计数实验结果计算结果相对偏差/%93Sr590.006.76×1056.58×105-2.6894Y918.001.26×1061.23×106-2.25134I1 072.552.21×1052.14×105-3.1489Rb1 248.102.21×1051.59×105-3.5689Rb2 196.001.48×1041.41×104-4.25138Cs1 435.865.27×1055.12×105-2.81138Xe1 768.003.67×1043.54×104-3.75142La2 397.802.58×1042.48×104-3.97142La2 542.701.76×1041.69×104-4.13142La2 971.004.36×1034.17×103-4.6995Y2 632.003.72×1033.53×103-4.94133Sb2 755.001.70×1031.60×103-5.5190Rb3 383.001.38×1031.30×103-5.67