楊洪新 張競宇 李平烈 趙曉澤 譚 琳 肖懿鑫

（華北電力大學(xué)核科學(xué)與工程學(xué)院 北京102206）

在反應(yīng)堆物理的領(lǐng)域中，燃耗計算主要用來分析重核和裂變產(chǎn)物同位素的原子核密度隨時間的變化趨勢以及它們對輸運計算（如中子通量密度分布和反應(yīng)性）的影響［1］。通過對燃耗的計算分析可以得到核素成分與原子核密度隨燃耗的變化以及中子通量密度的密度分布等以供進(jìn)一步研究核燃料的轉(zhuǎn)化及循環(huán)、確定堆芯壽期等，為反應(yīng)堆堆芯燃耗管理、燃料增殖以及燃耗信任制技術(shù)應(yīng)用和乏燃料后處理的研究提供基礎(chǔ)［2］。在反應(yīng)堆源項分析領(lǐng)域，活化計算主要用來分析材料受中子輻照后產(chǎn)生的放射性核素成分和活度，其計算模型和方法與燃耗計算是一致的，僅僅是應(yīng)用領(lǐng)域不同，因此不需特意區(qū)別兩者的差別。對產(chǎn)生的核素進(jìn)行徑跡追蹤以研究核素的產(chǎn)生路徑及每個路徑的貢獻(xiàn)統(tǒng)計，有助于了解核素的轉(zhuǎn)化途徑以辨識重要核素或核反應(yīng)，是反應(yīng)堆物理和源項分析中十分重要的內(nèi)容。

1 國內(nèi)外研究現(xiàn)狀

現(xiàn)有的反應(yīng)堆燃耗/活化計算方法一般可以分為兩類，分別是基于單條燃耗鏈求解的線性子鏈方法（Transmutation Trajectory Analysis，TTA）和基于燃耗方程矩陣求解的矩陣指數(shù)方法［3］。其中CINDER［4］點燃耗程序和MCB［5］等蒙特卡羅程序的燃耗計算模塊使用TTA方法，目前TTA方法已經(jīng)可以實現(xiàn)以常數(shù)項和基于拉普拉斯變換的多項式項為非齊次項的燃耗方程計算［6］。ORIGEN［7］點燃耗程序則使用多項式逼近方法。SERPENT［8］蒙特卡羅程序燃耗計算模塊首先使用了切比雪夫有理近似方法（Chebyshev Rational Approximation Method，CRAM）。歐洲活化程序FISPACT在2007版本中使用了歐拉指數(shù)方法，在后續(xù)版本中使用了GEAR 方法［9］。國內(nèi)計算程序主要有清華大學(xué)開發(fā)的NUIT程序和華北電力大學(xué)核反應(yīng)堆源項研究組開發(fā)的ABURN程序。ABURN程序基于Fortran語言開發(fā)，包含了TTA 以及CRAM 等算法，兼容EAF-2007 和ORIGEN 數(shù)據(jù)庫，程序可實現(xiàn)活化計算、衰變計算、脈沖工況計算、敏感性及不確定性分析等功能。

表1列出了目前典型的可用于復(fù)雜核素系統(tǒng)計算的通用點燃耗程序的基本情況。

表1 通用的點燃耗程序Table 1 General ignition consumption program

目前國內(nèi)所開發(fā)的程序大多基于數(shù)值算法，只提供了核素成分和核素存量的計算結(jié)果，并沒有進(jìn)行核素轉(zhuǎn)變過程的反應(yīng)鏈精細(xì)追蹤，而對于輻射防護(hù)工作者來說，確定放射性核素的來源并制定相應(yīng)的防護(hù)措施是十分重要的。因此本文基于ABURN程序開發(fā)了徑跡追蹤和貢獻(xiàn)統(tǒng)計的計算模塊，可在實際計算中考察某核素的產(chǎn)生路徑及貢獻(xiàn)百分比，以進(jìn)一步滿足輻射防護(hù)領(lǐng)域的工程應(yīng)用要求。

2 活化計算

2.1 計算方法

TTA的基本原理是將復(fù)雜的燃耗鏈拆解成一條條獨立的線性鏈，然后再對每一條燃耗鏈進(jìn)行解析求解，最后綜合各線性鏈中同種核素的計算結(jié)果，即可得到在復(fù)雜網(wǎng)狀燃耗鏈中的對應(yīng)某一時刻的核素存量［10］。TTA方法在理論上可以得到每條鏈的精確解，能夠?qū)μ囟ê怂剡M(jìn)行徑跡追蹤和貢獻(xiàn)統(tǒng)計研究，所以選取TTA方法對核素徑跡追蹤進(jìn)行研究，給出指定核素的產(chǎn)生路徑并給出每一條產(chǎn)生該核素的路徑的貢獻(xiàn)統(tǒng)計。采用TTA 方法求解也可以有效避免矩陣剛性問題引起的數(shù)值不穩(wěn)定性，并且在本質(zhì)上是求出了各條反應(yīng)鏈的解析解，所以結(jié)果也更加準(zhǔn)確。

2.2 線性子鏈方法

燃耗計算即在給定的燃耗步長內(nèi)求解核素原子核密度隨時間的變化情況，有燃耗方程：

在給定燃耗區(qū)，給定燃耗步長內(nèi)，式（1）可以簡化為：

活化計算主要研究的問題是核反應(yīng)堆的材料受到中子輻照時的核素成分隨時間的變化，活化方程的表達(dá)式和燃耗方程比較類似，區(qū)別在于通常活化方程不考慮裂變反應(yīng)，其主要關(guān)注核素的衰變和活化反應(yīng)，因此有活化方程：

式中：Ni(t)表示第i種核素在t時刻的存量，cm-3；σi表示第i種核素的中子誘發(fā)反應(yīng)微觀截面，cm2；σij表示第j種核素活化生成第i種核素的中子誘發(fā)微觀截面，cm2；λi表示第i種核素的衰變常量，s-1；λij表示第j種核素活化生成第i種核素的衰變常量，s-1；φ表示中子通量密度，cm-2·s-1。

對于一條含有n個核素的核素鏈，假設(shè)只有第一個核素具有初始核素密度N1(0)，利用TTA 方法可以計算得到第k個核素經(jīng)過時間t后的核素存量為：

3 徑跡追蹤貢獻(xiàn)統(tǒng)計研究

線性子鏈方法將復(fù)雜的網(wǎng)狀燃耗鏈拆解成一條條簡單的鏈進(jìn)行計算，因此實現(xiàn)徑跡貢獻(xiàn)統(tǒng)計功能，需要利用回溯法對計算中的核素和反應(yīng)信息進(jìn)行搜索，要遍歷整個解空間搜索所有核素可能的反應(yīng)路徑，并使用截斷誤差進(jìn)行截斷［11］。

3.1 回溯法

回溯算法是一種類似枚舉的搜索嘗試過程，按照深度優(yōu)先的準(zhǔn)則進(jìn)行搜索，從一個頭節(jié)點開始，向核反應(yīng)的縱深方向搜索，直到搜索到某一個核素不再發(fā)生其他反應(yīng)或滿足預(yù)先設(shè)置的截斷標(biāo)準(zhǔn)時，回溯到反應(yīng)鏈的上一個核素檢驗有無其他未反應(yīng)的反應(yīng)道，若有則以這個核素為新節(jié)點繼續(xù)往縱深方向搜索，若無則再回溯到該核素的上一個核素繼續(xù)進(jìn)行檢驗，直到回溯到反應(yīng)鏈的頭節(jié)點，對所有初始核素完成相關(guān)反應(yīng)的搜索［12］。

回溯法可以在搜索過程中遍歷整個解空間，但是真實的反應(yīng)堆級別的反應(yīng)網(wǎng)是非常復(fù)雜的，并且如果出現(xiàn)閉環(huán)就會導(dǎo)致反應(yīng)鏈出現(xiàn)無盡循環(huán)，因此需要制定一個截斷標(biāo)準(zhǔn)，在滿足條件的情況下對反應(yīng)鏈進(jìn)行截止，剔除不重要的反應(yīng)以提高計算效率［13］。

定義子鏈中每個核素的截斷比例Pk(t)，即t時刻核素鏈中當(dāng)前核素的原子核密度占核素鏈母核初始密度的比例份額［14］：

再根據(jù)線性子鏈的截斷標(biāo)準(zhǔn)判斷是否截斷：

3.2 徑跡貢獻(xiàn)統(tǒng)計研究方案

假設(shè)存在如下核素A生成核素F的反應(yīng)網(wǎng)絡(luò)：

根據(jù)線性鏈方法，可以將復(fù)雜的反應(yīng)網(wǎng)絡(luò)拆解成一條條單獨的鏈進(jìn)行解析計算，則核素A生成核素F的反應(yīng)網(wǎng)絡(luò)可以拆解為：

在圖1顯示的反應(yīng)網(wǎng)絡(luò)中，初始核素A共有4條路徑生成核素F，在此反應(yīng)網(wǎng)絡(luò)中，各核素是相互獨立的，因此共有4 組F核素輸出。在計算中，核素F的總量可以表示為：

圖1 反應(yīng)網(wǎng)Fig.1 Reaction network

如圖2 所示，即在進(jìn)行活化計算時展示的初始核素A生成核素F的路徑，其每條路徑的貢獻(xiàn)百分比ηi即此條鏈生成的F核素量所占F核素總量的份額，即：

圖2 反應(yīng)鏈Fig.2 Reaction chain

ABURN 程序的徑跡追蹤與貢獻(xiàn)統(tǒng)計功能是基于線性子鏈方法結(jié)合回溯算法得到的，理論上可以對反應(yīng)中涉及到的所有核素進(jìn)行回溯追蹤并將其路徑及貢獻(xiàn)百分比表示出來，因此具備較高的計算精度。

4 算例測試

本文基于TTA 方法自主開發(fā)了徑跡貢獻(xiàn)統(tǒng)計計算模塊，并將其集成到華北電力大學(xué)源項研究組自主開發(fā)的活化程序ABURN 中。為了對ABURN程序的徑跡追蹤計算結(jié)果進(jìn)行驗證，本文選取IAEA發(fā)布的天然鐵活化基準(zhǔn)題、50Cr 活化基準(zhǔn)題［16］以及UKAEA 發(fā)布的95Zr 衰變基準(zhǔn)題［17］，ABURN 程序和歐洲活化計算程序FISPACT 均使用EAF-2007 數(shù)據(jù)庫進(jìn)行核素徑跡貢獻(xiàn)統(tǒng)計計算結(jié)果對比。

4.1 天然鐵活化基準(zhǔn)題

在天然鐵的活化計算中，初始核素包括鐵的4種同位素，包括54Fe、56Fe、57Fe以及58Fe，能群數(shù)為100群，中子通量密度為3.640 99×1015cm-2·s-1，計算時間為1 a，截斷參數(shù)選取1.0×10-29，具體參數(shù)如表2 所示。計算中能群數(shù)為100，能群分布如表3所示。

表2 天然鐵活化計算輸入?yún)?shù)Table 2 Input parameters for natural Fe activation calculation

表3 能群分布Table 3 Energy group distribution

續(xù)表

表4 比較了活化反應(yīng)產(chǎn)生的各核素原子核數(shù)目，其中A表示ABURN 程序的計算結(jié)果，F(xiàn)表示FISPACT 程序的計算結(jié)果，兩者的差異用相對誤差δ來表示，相對誤差δ計算公式為：

天然鐵活化計算結(jié)果為：

在計算結(jié)果中選取52V、51Cr、54Mn、55Fe 以及60Co進(jìn)行徑跡追蹤研究（表5）。

表5 天然鐵徑跡追蹤結(jié)果對比Table 5 Pathway analysis of Fe activation calculation

對上述結(jié)果進(jìn)行對比分析可以發(fā)現(xiàn)：對于核素原子核數(shù)目的計算結(jié)果，ABURN 程序與FISPACT程序的相對誤差絕對值均小于0.5%且大部分小于0.1%，說明ABURN程序的計算結(jié)果精度較高；對于目標(biāo)核素的路徑及各路徑貢獻(xiàn)百分比的計算結(jié)果，ABURN 程序與FISPACT 程序基本一致，說明ABURN徑跡追蹤的計算結(jié)果較為準(zhǔn)確。

4.2 95Zr衰變基準(zhǔn)題

在95Zr 衰變計算中，初始核素為95Zr，初始核素數(shù)目為1×1020，計算時間為5.0×107s，截斷參數(shù)選取1.0×10-29，具體參數(shù)如表6 所示。95Zr 計算結(jié)果見表7。對各核素進(jìn)行徑跡追蹤研究，結(jié)果如表8所示。

表6 95Zr衰變計算輸入?yún)?shù)Table 6 Input parameters of95Zr decay calculation

表7 95Zr衰變計算結(jié)果對比Table 7 Comparison of95Zr decay calculation results

表8 95Zr計算徑跡追蹤結(jié)果對比Table 8 Pathway analysis of95Zr decay calculation

對上述結(jié)果進(jìn)行對比分析，可以發(fā)現(xiàn)：對于核素原子核數(shù)目的計算結(jié)果，ABURN 程序與FISPACT程序的相對誤差絕對值均小于0.01%，說明ABURN程序的計算結(jié)果精度較高；對于目標(biāo)核素的路徑及各路徑貢獻(xiàn)百分比的計算結(jié)果，ABURN 程序與衰變鏈理論計算結(jié)果基本一致，說明ABURN 徑跡追蹤的計算結(jié)果較為準(zhǔn)確。

4.3 50Cr先活化后衰變基準(zhǔn)題

在50Cr 的先活化后衰變基準(zhǔn)題計算中，初始核素為50Cr，能群數(shù)100，能群分布如表3 所示，中子通量密度為3.640 99×1015cm-2·s-1，計算時間為先活化計算1 a，再衰變計算0.5 a，截斷參數(shù)選取1.0×10-29，具體參數(shù)如表9 所示。50Cr 先活化后衰變計算結(jié)果見表10。在計算結(jié)果中選取39Ar、50V、48Ca以及54Mn進(jìn)行徑跡追蹤研究（表11）。

表9 50Cr先活化后衰變計算輸入?yún)?shù)Table 9 Input parameters of50Cr activation and decay calculation

表10 50Cr先活化后衰變計算結(jié)果對比Table 10 Comparison of50Cr activation and decay calculation results

表11 50Cr先活化后衰變計算徑跡追蹤結(jié)果對比Table 11 Pathway analysis of50Cr activation and decay calculation

對上述結(jié)果進(jìn)行對比分析，可以發(fā)現(xiàn)：對于核素原子核數(shù)目的計算結(jié)果，ABURN 程序與FISPACT程序的相對誤差絕對值均小于0.1%，說明ABURN程序的計算結(jié)果精度較高；對于目標(biāo)核素的路徑及各路徑貢獻(xiàn)百分比的計算結(jié)果，ABURN 程序與衰變鏈理論計算結(jié)果基本一致，說明ABURN 徑跡追蹤的計算結(jié)果較為準(zhǔn)確。

5 結(jié)語

對活化產(chǎn)生的核素進(jìn)行徑跡追蹤貢獻(xiàn)統(tǒng)計以研究在實際反應(yīng)中核素的產(chǎn)生路徑及各路徑的貢獻(xiàn)百分比，是反應(yīng)堆物理和源項分析中十分重要的內(nèi)容。本文采用線性子鏈方法，結(jié)合回溯法，實現(xiàn)在活化計算中對指定核素的產(chǎn)生路徑進(jìn)行追蹤分析，并給出每條路徑的貢獻(xiàn)百分比，并與國外廣泛使用的同類軟件FISPACT進(jìn)行了對比驗證，分別驗證了天然鐵活化基準(zhǔn)題，95Zr 衰變基準(zhǔn)題，50Cr 先活化后衰變基準(zhǔn)題等三組基準(zhǔn)題，對其中部分核素進(jìn)行了徑跡追蹤貢獻(xiàn)統(tǒng)計研究，核素的路徑一致，各路徑的貢獻(xiàn)百分比也基本一致。

經(jīng)過測試驗證，本文自主編寫的徑跡追蹤程序模塊的計算結(jié)果與FISPACT程序的基本一致，可以認(rèn)為本文徑跡追蹤方法和程序是可行的，從而為活化計算中重要核素或核反應(yīng)的篩選提供了有力的仿真工具。

作者貢獻(xiàn)聲明楊洪新：實施研究，采集數(shù)據(jù)，分析解釋數(shù)據(jù)，起草文章；張競宇：指導(dǎo)，技術(shù)材料支持，對文章的知識性內(nèi)容進(jìn)行審閱；李平烈：技術(shù)材料支持，統(tǒng)計分析，內(nèi)容校核；趙曉澤：采集數(shù)據(jù)，分析數(shù)據(jù)；譚琳：采集數(shù)據(jù)，分析數(shù)據(jù)；肖懿鑫：采集數(shù)據(jù)，分析數(shù)據(jù)。