CRAM在放射性核素存量計(jì)算中的應(yīng)用

張競(jìng)宇1 馬亞棟1 陳義學(xué)1 高強(qiáng)2

?

CRAM在放射性核素存量計(jì)算中的應(yīng)用

張競(jìng)宇馬亞棟陳義學(xué)高強(qiáng)

1（華北電力大學(xué)核科學(xué)與工程學(xué)院 北京 102206） 2（環(huán)境保護(hù)部核與輻射安全中心 北京 100082）

在反應(yīng)堆中，組成材料的穩(wěn)定核素經(jīng)受強(qiáng)中子輻照后，會(huì)被活化成放射性核素。這些核素及其衰變產(chǎn)物對(duì)工作人員的職業(yè)輻照劑量具有重要貢獻(xiàn)。為了更好地進(jìn)行人員的輻射防護(hù)工作，需要對(duì)放射性核素的存量進(jìn)行精確計(jì)算。相對(duì)于核素平衡方程的其它求解方法，切比雪夫有理逼近方法(Chebyshev Rational Approximation Method, CRAM)在計(jì)算精度和效率方面具有綜合性優(yōu)勢(shì)。首先介紹了CRAM的基本理論，隨后選取典型的例題進(jìn)行了測(cè)試驗(yàn)證。與解析解對(duì)比的結(jié)果表明，采用CRAM進(jìn)行中子輻照下的核素活化衰變計(jì)算能夠取得不錯(cuò)的效果，但是用于核素長(zhǎng)期衰變計(jì)算可能導(dǎo)致計(jì)算錯(cuò)誤。針對(duì)此問(wèn)題，將收縮乘方技術(shù)與CRAM相結(jié)合，取得了正確的計(jì)算結(jié)果，拓展了CRAM的適用范圍。

切比雪夫有理逼近方法，放射性核素，存量計(jì)算，收縮乘方

反應(yīng)堆中存在很強(qiáng)的中子輻照，會(huì)導(dǎo)致結(jié)構(gòu)材料的活化而使其具有放射性，當(dāng)工作人員檢查或維修反應(yīng)堆中的管道或設(shè)備時(shí)，不可避免地會(huì)受到這些γ射線的輻照傷害。此外，結(jié)構(gòu)材料在活化過(guò)程中會(huì)產(chǎn)生氫氣和氦氣，這些氣體物質(zhì)會(huì)對(duì)結(jié)構(gòu)材料的強(qiáng)度造成損傷，其中的氚核素還具有很強(qiáng)的穿透性和放射性，容易對(duì)人員和環(huán)境造成危害。因此，需要特別關(guān)注放射性核素的存量計(jì)算問(wèn)題，盡可能提高計(jì)算結(jié)果的準(zhǔn)確性。

活化和衰變計(jì)算中涉及的核素種類多達(dá)上千種，并且不同核素的反應(yīng)截面和半衰期差別很大，因此所對(duì)應(yīng)的核素平衡方程在數(shù)學(xué)上具有大規(guī)模、強(qiáng)剛性的特點(diǎn)。對(duì)于此類問(wèn)題的求解，學(xué)界經(jīng)過(guò)長(zhǎng)期的研究，開(kāi)發(fā)了一系列的方法和程序，主要包括時(shí)間差分方法、線性子鏈方法和矩陣指數(shù)方法。

時(shí)間差分方法需要采用特殊處理才能提高計(jì)算效率。例如FISPACT-2007程序預(yù)先將短壽命核素從系數(shù)矩陣中移出以降低矩陣的病態(tài)性，再采用顯式的指數(shù)歐拉方法基于大步長(zhǎng)求解僅含有長(zhǎng)壽命核素的方程，最后通過(guò)“平衡假設(shè)”近似計(jì)算短壽命核素，該方法的計(jì)算速度較快，但是由于采用了近似處理，在計(jì)算初期會(huì)有一定誤差。

線性子鏈解析(Transmutation Trajectory Analysis, TTA)方法的基本原理是將復(fù)雜反應(yīng)鏈分解成獨(dú)立的線性反應(yīng)鏈，再解析求解，以避免直接處理強(qiáng)剛性方程，目前已用在以CINDER程序?yàn)榇淼亩鄠€(gè)程序中。該方法求得的解析解精度很高，但是反應(yīng)鏈線性化過(guò)程比較繁瑣使得整體效率較低，對(duì)此吳明宇等采用回溯算法自動(dòng)進(jìn)行反應(yīng)鏈搜索，一定程度上提高了計(jì)算效率。

矩陣指數(shù)方法將對(duì)核素平衡方程的求解轉(zhuǎn)化為對(duì)矩陣指數(shù)的計(jì)算，目前計(jì)算矩陣指數(shù)的方法主要基于數(shù)值逼近理論，代表程序有基于泰勒展開(kāi)多項(xiàng)式逼近方法結(jié)合短壽命核素特殊處理的ORIGEN程序、基于切比雪夫有理逼近方法(Chebyshev Rational Approximation Method, CRAM)的SERPENT程序燃耗模塊和DEPTH程序。其中，CRAM憑借無(wú)需對(duì)短壽命核素進(jìn)行單獨(dú)處理、在較大的時(shí)間步長(zhǎng)下仍能得到很高的計(jì)算精度等優(yōu)點(diǎn)而被廣泛推薦。

1 切比雪夫有理逼近方法

對(duì)于含有種核素的材料，在中子輻照下，其核素平衡方程如下：

式(1)可以寫(xiě)成如下的矩陣形式：

式中：為系數(shù)矩陣。

采用矩陣指數(shù)方法求解式(2)，解的形式如下：

根據(jù)文獻(xiàn)[9]的觀察，系數(shù)矩陣的特征值分布在復(fù)數(shù)平面的負(fù)實(shí)軸附近。此時(shí)，CRAM是指數(shù)函數(shù)在(?∞, 0]區(qū)間的最佳有理逼近，如式(4)所示：

考慮到數(shù)值計(jì)算的穩(wěn)定性，通常將式(4)寫(xiě)為以極點(diǎn)和留數(shù)表示的PFD (Partial Fraction Decomposition)形式如下：

式中：表示在?∞處的極限值；α表示留數(shù)；θ表示極點(diǎn)；表示展開(kāi)階數(shù)。

由于極點(diǎn)通常是以共軛對(duì)的形式出現(xiàn)，因此可以利用此性質(zhì)，將式(5)寫(xiě)成如下形式，其計(jì)算量節(jié)省為原來(lái)的一半。

PFD系數(shù)（、α和θ）可以在文獻(xiàn)[11]提供的多項(xiàng)式系數(shù)的基礎(chǔ)上通過(guò)Carathéodory-Fejér方法計(jì)算得到。根據(jù)文獻(xiàn)[13]的研究，當(dāng)階數(shù)≤16時(shí)，Carathéodory-Fejér方法可保持較高的計(jì)算精度。本文選取=16，計(jì)算得到PFD系數(shù)如表1所示。

表1 16階切比雪夫有理逼近方法的系數(shù)

當(dāng)系數(shù)矩陣的元素是按照核素質(zhì)量數(shù)增加的順序排列時(shí)，矩陣具有稀疏結(jié)構(gòu)，其求逆計(jì)算變得相對(duì)容易。

2 測(cè)試?yán)}的計(jì)算結(jié)果

2.1O活化例題的計(jì)算結(jié)果

O活化例題如下：

O (stable)→(n,γ: 6.601 89×10b) →

O (stable)→(n,γ: 1.184 99×10b) →

O (stable)→(n,γ: 5.426 13×10b) →

O (β, 1.577 88×10s)

其中：相關(guān)的反應(yīng)類型和截面已經(jīng)在括號(hào)中予以標(biāo)注；O的初始數(shù)量為10個(gè)；中子通量密度為10cm·s；中子輻照時(shí)間為1 a。

O活化例題的計(jì)算結(jié)果如表2所示。表2中，相對(duì)偏差定義為：（解析解?CRAM計(jì)算結(jié)果）/解析解。其中解析解采用TTA計(jì)算得到，為追求較高的計(jì)算精度，截?cái)鄻?biāo)準(zhǔn)設(shè)置為10。

表216O活化例題計(jì)算結(jié)果

在計(jì)算時(shí)，CRAM的時(shí)間步長(zhǎng)直接取為1 a，所消耗的計(jì)算時(shí)間為0.0125 s，解析方法所消耗的計(jì)算時(shí)間為0.69 s，對(duì)于這種簡(jiǎn)單線性鏈來(lái)說(shuō)，CRAM的計(jì)算效率高于解析方法。結(jié)合表2可以看到，采用CRAM進(jìn)行中子輻照下的核素活化計(jì)算能夠取得不錯(cuò)的效果，其計(jì)算精度和效率都比較高。

2.2 天然Fe活化例題的計(jì)算結(jié)果

天然Fe的初始質(zhì)量為1 kg，富集度分布如下：Fe 5.8%、Fe 91.72%、Fe 2.2%、Fe 0.28%。相關(guān)的反應(yīng)類型、截面和半衰期取自EAF-2007數(shù)據(jù)庫(kù)。中子通量密度為3.641×10cm?s，中子輻照時(shí)間為1 a。天然Fe活化例題中除了活化反應(yīng)外，還有大量的衰變反應(yīng)，因此本例題為活化衰變復(fù)合計(jì)算例題，計(jì)算結(jié)果如表3所示。

表3 天然Fe活化例題計(jì)算結(jié)果

在計(jì)算時(shí)，CRAM的時(shí)間步長(zhǎng)直接取為1 a，所消耗的計(jì)算時(shí)間為2.37 s，解析方法所消耗的計(jì)算時(shí)間為11.84 s，對(duì)于這種復(fù)雜網(wǎng)狀鏈來(lái)說(shuō)，CRAM的計(jì)算效率高于解析方法。結(jié)合表3可以看到，采用CRAM進(jìn)行中子輻照下的核素活化衰變計(jì)算能夠取得不錯(cuò)的效果，其計(jì)算精度和效率都比較高。

2.3Fe衰變例題的計(jì)算結(jié)果

Fe衰變例題如下：

Fe (, 1.548×10s)→Fe (, 5.106×10s)→

Mn (, 1.166 8×10s)→Cr (stable)

其中：相關(guān)的反應(yīng)類型和半衰期()已經(jīng)在括號(hào)中予以標(biāo)注；Fe的初始數(shù)量為10個(gè)。Fe衰變例題的計(jì)算結(jié)果如表4所示。

在計(jì)算時(shí)，CRAM的時(shí)間步長(zhǎng)直接取為衰變時(shí)間。從表4可以看到，當(dāng)衰變時(shí)間在10?10s時(shí)，CRAM的計(jì)算結(jié)果與解析解符合較好。但是當(dāng)衰變時(shí)間大于10s，CRAM的計(jì)算結(jié)果嚴(yán)重偏離解析解。說(shuō)明采用CRAM進(jìn)行核素長(zhǎng)期衰變計(jì)算可能導(dǎo)致計(jì)算錯(cuò)誤。

表453mFe衰變例題計(jì)算結(jié)果

衰變時(shí)間Decay times53Fe原子數(shù)量（解析法）Number of 53Fe atoms (analytic)53Fe原子數(shù)量(CRAM)Number of 53Fe atoms (CRAM)相對(duì)偏差Relative deviation 1.0×1014.34912×10184.34912×10180.00000×100 1.0×1023.35819×10193.35819×10190.00000×100 1.0×1033.52942×10193.52943×10192.83333×10?6 1.0×1041.82500×1014?8.08213×1013?1.44286×100 1.0×1050.00000×100?8.49758×1012∞ 1.0×1060.00000×100?1.42934×1011∞ 1.0×1070.00000×100?8.16992×109∞ 1.0×1080.00000×100?7.56749×108∞

衰變時(shí)間Decay times53Cr原子數(shù)量（解析法）Number of 53Cr atoms (analytic)53Cr原子數(shù)量(CRAM)Number of 53Cr atoms (CRAM)相對(duì)偏差Relative deviation 1.0×1015.93138×1025.93138×1020.00000×100 1.0×1025.21854×1055.21854×1050.00000×100 1.0×1031.84709×1081.84709×1080.00000×100 1.0×1045.37031×1095.37031×1090.00000×100 1.0×1055.88355×10105.88356×10101.69965×10?6 1.0×1065.93488×10115.93488×10110.00000×100 1.0×1075.94001×10125.94001×10120.00000×100 1.0×1085.94052×10135.94052×10130.00000×100

3 收縮乘方技術(shù)的基本原理及應(yīng)用效果

如前所述，式(3)中矩陣的特征值分布在復(fù)數(shù)平面的負(fù)實(shí)軸附近時(shí)，CRAM是矩陣指數(shù)在(?∞, 0]區(qū)間的最佳有理逼近。對(duì)于核素衰變工況，矩陣的特征值會(huì)隨著衰變時(shí)間而變大，即特征值的虛部也會(huì)相應(yīng)變大，使得特征值顯著偏離復(fù)數(shù)平面的負(fù)實(shí)軸。此時(shí)，如圖1所示（詳見(jiàn)文獻(xiàn)[13]），采用CRAM近似矩陣指數(shù)，其計(jì)算精度會(huì)顯著下降，甚至出現(xiàn)計(jì)算錯(cuò)誤。這問(wèn)題可以通過(guò)引入收縮乘方技術(shù)(Scaling and squaring)予以解決。

圖1 16階切比雪夫有理逼近方法(R16,16(At))近似矩陣指數(shù)的計(jì)算精度隨矩陣At特征值的變化

3.1 收縮乘方技術(shù)的基本原理

收縮乘方技術(shù)就是通過(guò)選擇合適的參數(shù)，使得能夠被CRAM準(zhǔn)確而高效地計(jì)算，隨后通過(guò)對(duì)做次乘方運(yùn)算，得到待求的矩陣指數(shù)。

采用收縮乘方技術(shù)后，式(3)可寫(xiě)成如下形式：

收縮乘方技術(shù)的本質(zhì)是將總的計(jì)算時(shí)間劃分為個(gè)時(shí)間步長(zhǎng)。相對(duì)于矩陣，矩陣/特征值減小了倍，相應(yīng)地其虛部也減小了倍，意味著特征值更靠近復(fù)數(shù)平面的負(fù)實(shí)軸，使得CRAM的適用性得到了更好的滿足。

在選擇參數(shù)時(shí)，一個(gè)常用的標(biāo)準(zhǔn)是取最小的，使得||||/≤1，其中||||是矩陣的譜范數(shù)，即||||=[(A)]。在Fe衰變例題中，矩陣的譜范數(shù)||||=0.00641，在譜范數(shù)計(jì)算領(lǐng)域這是一個(gè)很大的值，意味著矩陣是強(qiáng)剛性的，這主要是矩陣中各個(gè)元素的變化速率差別過(guò)大造成的，比如Fe的衰變常數(shù)為4.478×10s，Mn的衰變常數(shù)為5.969×10s，兩者相差了12個(gè)量級(jí)。所以當(dāng)衰變時(shí)間較長(zhǎng)時(shí)，||||會(huì)遠(yuǎn)大于1，此時(shí)需要應(yīng)用收縮乘方技術(shù)，以滿足||||/≤1。

3.2 收縮乘方技術(shù)的應(yīng)用效果

采用收縮乘方技術(shù)后，F(xiàn)e衰變例題的計(jì)算結(jié)果如表5所示。

表5 采用收縮乘方技術(shù)后53mFe衰變例題計(jì)算結(jié)果

衰變時(shí)間Decay times53Fe原子數(shù)量（解析法）Number of 53Fe atoms (analytic)53Fe原子數(shù)量(CRAM)Number of 53Fe atoms (CRAM)相對(duì)偏差Relative deviation 1.0×1014.34912×10184.34912×10180.00000×100 1.0×1023.35819×10193.35819×10190.00000×100 1.0×1033.52942×10193.52935×1019?1.98333×10?5 1.0×1041.82500×10141.82490×1014?5.47945×10?5 1.0×1050.00000×1000.00000×1000.00000×100 1.0×1060.00000×1000.00000×1000.00000×100 1.0×1070.00000×1000.00000×1000.00000×100 1.0×1080.00000×1000.00000×1000.00000×100

（續(xù)表5）

衰變時(shí)間Decay times53Mn原子數(shù)量（解析法）Number of 53Mn atoms (analytic)53Mn原子數(shù)量(CRAM)Number of 53Mn atoms (CRAM)相對(duì)偏差Relative deviation 1.0×1012.98086×10162.98086×10160.00000×100 1.0×1022.51294×10182.51294×10180.00000×100 1.0×1036.35698×10196.35687×1019?1.73038×10?5 1.0×1049.99998×10199.99990×1019?8.00002×10?6 1.0×1051.00000×10209.99998×1019?2.00000×10?6 1.0×1061.00000×10209.99998×1019?2.00000×10?6 1.0×1071.00000×10209.99998×1019?2.00000×10?6 1.0×1089.99999×10199.99998×1019?1.00000×10?6

衰變時(shí)間Decay times53Cr原子數(shù)量（解析法）Number of 53Cr atoms (analytic)53Cr原子數(shù)量(CRAM)Number of 53Cr atoms (CRAM)相對(duì)偏差Relative deviation 1.0×1015.93138×1025.93138×1020.00000×100 1.0×1025.21854×1055.21854×1050.00000×100 1.0×1031.84709×1081.84706×108?1.62418×10?5 1.0×1045.37031×1095.37028×109?5.58627×10?6 1.0×1055.88355×10105.88356×10101.69965×10?6 1.0×1065.93488×10115.93488×10110.00000×100 1.0×1075.94001×10125.94001×10120.00000×100 1.0×1085.94052×10135.94052×10130.00000×100

鑒于矩陣的譜范數(shù)||||=0.00641，為了滿足標(biāo)準(zhǔn)||||/≤1，時(shí)間步長(zhǎng)/被選為100 s。這樣對(duì)于衰變時(shí)間為10s的情況，就需要?jiǎng)澐譃?0個(gè)時(shí)間步長(zhǎng)進(jìn)行計(jì)算，可預(yù)計(jì)的時(shí)間消耗是很大的。但是由表5可以發(fā)現(xiàn)，對(duì)于像Fe、Fe這樣的短壽命核素，在衰變10s時(shí)，原子數(shù)量已經(jīng)減少到0個(gè)，此時(shí)可以從矩陣中將它們移除，相應(yīng)地新的矩陣的譜范數(shù)變?yōu)?.940 58×10，在||||/≤1條件下，時(shí)間步長(zhǎng)可以取到10s，使得計(jì)算效率得以大幅提高。經(jīng)統(tǒng)計(jì)，對(duì)于衰變時(shí)間為10s的情況，解析方法的計(jì)算耗時(shí)為0.53 s，耦合了收縮乘方技術(shù)的CRAM的計(jì)算耗時(shí)為0.62 s，兩者基本相當(dāng)。

此外，從表5可以看到，耦合了收縮乘方技術(shù)后，采用CRAM進(jìn)行核素長(zhǎng)期衰變計(jì)算就可以得到準(zhǔn)確的結(jié)果。

4 結(jié)語(yǔ)

本文介紹了CRAM的基本理論，以及應(yīng)用于求解核素平衡方程的實(shí)現(xiàn)方式。選取O核素和天然Fe元素活化例題以及Fe核素衰變例題對(duì)CRAM進(jìn)行了測(cè)試驗(yàn)證，并與TTA方法得到的解析解進(jìn)行了對(duì)比。

對(duì)于存在中子反應(yīng)的問(wèn)題來(lái)說(shuō)，TTA方法在將閉環(huán)反應(yīng)鏈線性化過(guò)程中，會(huì)截?cái)嗪蜕釛壱徊糠址磻?yīng)鏈，因此從原理上來(lái)說(shuō)會(huì)有一定誤差，但是如果截?cái)鄻?biāo)準(zhǔn)選取得比較嚴(yán)格，閉環(huán)反應(yīng)鏈線性化引入的誤差是很小的，但是會(huì)增加大量的計(jì)算時(shí)間。而CRAM方法不會(huì)對(duì)閉環(huán)反應(yīng)鏈進(jìn)行截?cái)嗪蜕釛?，因此不存在這方面的誤差，其計(jì)算精度比較高，而且可以采用大的計(jì)算步長(zhǎng)，因此計(jì)算效率也比較高。

對(duì)于核素衰變問(wèn)題，不涉及閉環(huán)反應(yīng)鏈，因此TTA方法可以直接給出精確的解析解，計(jì)算效率也比較高。CRAM方法對(duì)于核素短期衰變計(jì)算是適用的，時(shí)間步長(zhǎng)可以直接取為衰變時(shí)間，但是對(duì)于核素長(zhǎng)期衰變計(jì)算，如果時(shí)間步長(zhǎng)直接取為衰變時(shí)間，矩陣的特征值會(huì)偏離復(fù)數(shù)平面的負(fù)實(shí)軸，此時(shí)采用CRAM方法計(jì)算容易出現(xiàn)計(jì)算錯(cuò)誤。因此需要配合使用收縮乘方技術(shù)，使得CRAM方法的適用性得到更好的滿足，收縮乘方技術(shù)的應(yīng)用會(huì)一定程度上降低CRAM方法的計(jì)算效率，但是仍然與TTA方法相當(dāng)。

后續(xù)將對(duì)CRAM方法進(jìn)行進(jìn)一步的測(cè)試驗(yàn)證和工程應(yīng)用。

1 李璐, 張競(jìng)宇, 郭慶洋, 等. 水冷聚變堆主回路活化產(chǎn)物源項(xiàng)計(jì)算分析[J]. 核技術(shù), 2016, 39(11): 110603. DOI: 10.11889/j.0253-3219.2016.hjs.39.110603. LI Lu, ZHANG Jingyu, GUO Qingyang,. Calculation and analysis of activation products source term in water-cooled fusion reactor primary circuit[J]. Nuclear Techniques, 2016, 39(11): 110603. DOI: 10.11889/j.0253-3219.2016.hjs.39.110603.

2 石巍, 曾勤, 李衛(wèi), 等. CFETR第一壁及赤道面外包層中子輻照損傷初步分析[J]. 核技術(shù), 2016, 39(12): 120602. DOI: 10.11889/j.0253-3219.2016.hjs.39.120602. SHI Wei, ZENG Qin, LI Wei,. Primary analysis of radiation damage on first wall and the outboard blanket on equatorial plane for CFETR[J]. Nuclear Techniques, 2016, 39(12): 120602. DOI: 10.11889/j.0253-3219.2016.hjs.39. 120602.

3 Lyu X W, Xia X B, Zhang Z H,. Analysis of tritium production in a 2 MW liquid-fueled molten salt experimental reactor and its environmental impact[J]. Nuclear Science and Techniques, 2016, 27(4): 78. DOI: 10.1007/s41365-016-0100-z.

4 Cetnar J. General solution of bateman equations for nuclear transmutations[J]. Annals of Nuclear Energy, 2006, 33(7): 640?645. DOI: 10.1016/j.anucene.2006.02. 004.

5 吳明宇, 王事喜, 楊勇, 等. 回溯算法在燃耗計(jì)算中的應(yīng)用[J]. 原子能科學(xué)技術(shù), 2013, 47(7): 1127?1132. DOI: 10.7538/yzk.2013.47.07.1127. WU Mingyu, WANG Shixi, YANG Yong,. Application of backtracking algorithm to depletion calculations[J]. Atomic Energy Science and Technology, 2013, 47(7): 1127?1132. DOI: 10.7538/yzk.2013.47.07. 1127.

6 Croff G. Origen2: a versatile computer code for calculating the nuclide compositions and characteristics of nuclear materials[J]. Nuclear Technology, 1983, 62: 335?352. DOI: 10.13182/NT83-1.

7 Lepp?nen J, Pusa M, Viitanen T. The serpent Monte Carlo code: status, development and applications in 2013[J]. Annals of Nuclear Energy, 2015, 82: 142?150. DOI: 10.1016/j.anucene.2014.08.024.

8 She D, Wang K, Yu G L. Development of the point-depletion code DEPTH[J]. Nuclear Engineering and Design, 2013, 258(17): 235?240. DOI: 10.1016/ j.nucengdes.2013.01.007.

9 Pusa M, Leppanen J. Computing the matrix exponential in burnup calculations[J]. Nuclear Science and Engineering, 2010, 164(2): 140?150. DOI: 10.13182/NSE09-14.

10 Cody W J, Meinardus G, Varga R S. Chebyshev rational approximation to exp(?) in [0, +∞] and applications to heat conduction problems[J]. Journal of Approximation Theory, 1969, 2(1): 50?65. DOI: 10.1016/0021-9045 (69)90030-6.

11 Carpenter A J, Ruttan A, Varga R S. Extended numerical computations on the “1/9” conjecture in rational approximation theory[J]. Lecture Notes in Mathematics, 1984, 1105(4): 383?411. DOI: 10.1007/BFb0072427.

12 Trefethen L N, Weideman J A C, Schmelzer T. Talbot quadratures and rational approximations[J]. BIT Numerical Mathematics, 2006, 46(3): 653?670. DOI: 10.1007/s10543-006-0077-9.

13 Pusa M. Higher-order Chebyshev rational approximation method and application to burnup equations[J]. Nuclear Science and Engineering, 2016, 182(3): 297?318. DOI: 10.13182/NSE15-26.

14 Moler C, Loan C V.Nineteen dubious ways to compute the exponential of a matrix[J]. SIAM Review,1978,20(4): 801?836. DOI: 10.1137/1020098.

Application of Chebyshev rational approximation method in inventory calculation of radioactive nuclides

ZHANG JingyuMA YadongCHEN YixueGAO Qiang

1(School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China) 2(Nuclear and Radiation Safety Centre, Ministry of Environmental Protection, Beijing 100082, China)

Background: The material suffering from strong neutron irradiation in the nuclear reactor will be activated to be radioactive nuclides. These nuclides and their decay products contribute a significant part to the occupational radiation exposure (ORE) of personnel. Purpose: For better radiation protection of the workers in nuclear reactor, it is supposed to calculate the inventory of radioactive nuclides accurately. Methods: Compared with other methods for solving the equilibrium equations of nuclides, the Chebyshev rational approximation method (CRAM) has comprehensive advantages on computational accuracy and efficiency. In this paper, the theory of CRAM method is described firstly, and then some typical cases are tested to verify CRAM method. Results & Conclusion: Compared with the analytical solution, CRAM method shows good effect on activation and decay calculation of nuclides under neutron irradiation, but may cause obvious error on long-term decay calculation of nuclides. After coupling with technique of scaling and squaring, CRAM method can derive accurate results for long-term decay calculation of nuclides and its scope of application is extended.

CRAM, Radioactive nuclides, Inventory calculation, Scaling and squaring

TL99

10.11889/j.0253-3219.2017.hjs.40.080502

國(guó)家自然科學(xué)基金(No.11605058)、國(guó)際熱核聚變實(shí)驗(yàn)堆計(jì)劃專項(xiàng)(No.2014GB119000)、中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(No.2017MS041)資助

張競(jìng)宇，男，1984年出生，2013年于清華大學(xué)獲博士學(xué)位，研究領(lǐng)域?yàn)楹朔磻?yīng)堆燃料燃耗/材料活化計(jì)算

高強(qiáng)，E-mail: gaoqiang@chinansc.cn

2017-01-20，

2017-03-24

National Natural Science Foundation of China (No.11605058), National Special Project for Magnetic Confined Nuclear Fusion Energy (No.2014GB119000),the Fundamental Research Funds for the Central Universities(No.2017MS041)

ZHANG Jingyu, male, born in 1984, graduated from Tsinghua University with a doctoral degree in 2013, focusing on the calculation of fuel depletion/material activation

GAO Qiang, E-mail: gaoqiang@chinansc.cn

2017-01-20, accepted date: 2017-03-24