核數(shù)據(jù)處理軟件NECP-Atlas中的光子相關(guān)數(shù)據(jù)計算方法研究

祖鐵軍，徐 寧，尹 文，曹良志，吳宏春

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

核數(shù)據(jù)處理軟件是聯(lián)系評價核數(shù)據(jù)與反應(yīng)堆核設(shè)計的紐帶，可為反應(yīng)堆核設(shè)計提供基礎(chǔ)核數(shù)據(jù)庫(應(yīng)用核數(shù)據(jù)庫)，是關(guān)鍵的反應(yīng)堆核設(shè)計軟件之一，其精度將直接影響最終核設(shè)計結(jié)果的精度。長期以來，國內(nèi)外主要使用美國洛斯阿拉莫斯國家實驗室的NJOY[1]、國際原子能機(jī)構(gòu)的PREPRO[2]等程序制作應(yīng)用核數(shù)據(jù)庫，由此會產(chǎn)生兩方面問題：一是這些軟件主要針對中子輸運(yùn)計算的需求開發(fā)，能提供的數(shù)據(jù)類型有限，不能滿足反應(yīng)堆核設(shè)計各方面的需求；二是核數(shù)據(jù)處理理論模型、數(shù)據(jù)庫制作方法研究基本停滯，并未隨核設(shè)計方法的發(fā)展得到相應(yīng)的發(fā)展，成為核反應(yīng)堆設(shè)計精度提高的障礙。近年來，隨著核反應(yīng)堆高保真數(shù)值模擬技術(shù)的發(fā)展，國際上對應(yīng)用核數(shù)據(jù)庫制作方法開展了一定研究[3]，這些研究僅是在NJOY程序基礎(chǔ)上對應(yīng)用數(shù)據(jù)進(jìn)行的針對性的改進(jìn)。2015年國際原子能機(jī)構(gòu)成立了核數(shù)據(jù)處理項目[4]，旨在推動核數(shù)據(jù)處理方法及軟件的發(fā)展，國際上一批新的核數(shù)據(jù)處理軟件參與了該項目，如中國核數(shù)據(jù)中心的Ruler、西安交通大學(xué)的NECP-Atlas、日本原子能機(jī)構(gòu)的FRENDY、法國原子能委員會的GALIEE等。

NECP-Atlas是西安交通大學(xué)研發(fā)的核數(shù)據(jù)處理軟件[5]，旨在建立核數(shù)據(jù)處理研究平臺，開展高精度的核數(shù)據(jù)處理方法研究，實現(xiàn)核數(shù)據(jù)處理軟件的自主化。該軟件自2015年開始研發(fā)，目前已具備豐富的核數(shù)據(jù)處理能力，可為反應(yīng)堆核設(shè)計提供不同類型的核數(shù)據(jù)，包括中子輸運(yùn)計算的反應(yīng)截面數(shù)據(jù)庫[6]、燃耗數(shù)據(jù)庫[7-8]、屏蔽數(shù)據(jù)庫[9]、反應(yīng)截面的多群協(xié)方差數(shù)據(jù)[10]、裂變產(chǎn)額的協(xié)方差數(shù)據(jù)[7-8]、釋熱與輻照損傷數(shù)據(jù)[11]、熱散射律數(shù)據(jù)[12]、活化及源項數(shù)據(jù)庫[13]等；開展了先進(jìn)的核數(shù)據(jù)處理方法研究，如共振彈性散射效應(yīng)處理方法[14-15]、熱能散射數(shù)據(jù)計算方法[16]等，提高了應(yīng)用數(shù)據(jù)庫的計算精度；基于NECP-Atlas軟件研制了國產(chǎn)化的應(yīng)用核數(shù)據(jù)庫，將我國最新的評價核數(shù)據(jù)庫CENDL-3.2應(yīng)用于核反應(yīng)堆物理設(shè)計及屏蔽設(shè)計軟件[17-18]。

光子是核反應(yīng)堆內(nèi)重要的輻射粒子，是堆內(nèi)釋熱率、材料輻照損傷等的重要來源。目前國際上現(xiàn)有的核數(shù)據(jù)處理軟件僅可給出中子的非彈性散射、輻射俘獲、裂變等核反應(yīng)釋放的瞬發(fā)光子，不能給出裂變產(chǎn)物衰變過程中釋放的光子，這將影響反應(yīng)堆內(nèi)功率或釋熱率的計算精度[19]。對于輻照損傷截面數(shù)據(jù)，目前核數(shù)據(jù)處理軟件僅考慮中子入射產(chǎn)生的輻照損傷，不能提供光子引起的輻照損傷截面數(shù)據(jù)，而Remec等[20]和Alexander等[21]的研究發(fā)現(xiàn)，光子的離位損傷是高通量核素生產(chǎn)堆、沸水堆壓力容器材料輻照損傷的重要組成部分。針對以上問題，本文擬進(jìn)行光子產(chǎn)生數(shù)據(jù)、光子輻照損傷數(shù)據(jù)計算方法研究，以完善NECP-Atlas的功能。

1 多群緩發(fā)光子產(chǎn)生矩陣計算功能

光子為中性粒子，與中子類似可在堆芯內(nèi)輸運(yùn)一定距離，在輸運(yùn)過程中將其能量沉積在發(fā)生核反應(yīng)的位置。核數(shù)據(jù)處理軟件需要為光子輸運(yùn)計算提供光子產(chǎn)生、光子與原子反應(yīng)相關(guān)的截面數(shù)據(jù)。堆芯內(nèi)光子按產(chǎn)生方式可分為2類：第1類是在中子引發(fā)的非彈性散射、輻射俘獲、裂變等反應(yīng)發(fā)生的瞬間釋放的光子，稱為瞬發(fā)光子；第2類是裂變產(chǎn)物衰變過程中釋放的光子，稱為緩發(fā)光子。傳統(tǒng)的核數(shù)據(jù)處理程序可根據(jù)中子反應(yīng)核數(shù)據(jù)計算中子引發(fā)核反應(yīng)產(chǎn)生瞬發(fā)光子的截面，但未考慮緩發(fā)光子的產(chǎn)生數(shù)據(jù)。對于燃料組件，光子對總釋熱率的貢獻(xiàn)約占10%[19]，其中緩發(fā)光子貢獻(xiàn)約為30%[22]，所以一般忽略緩發(fā)光子對釋熱的貢獻(xiàn)。但鈉冷快堆等反應(yīng)堆中包含了鋼反射層、控制組件等不含裂變材料的組件，這些組件內(nèi)的釋熱率主要來自于光子釋熱，文獻(xiàn)[19]中指出，快堆反射層內(nèi)光子對釋熱率的貢獻(xiàn)達(dá)90%，因此，需要精確的光子產(chǎn)生數(shù)據(jù)。NECP-Atlas根據(jù)評價數(shù)據(jù)庫提供的裂變產(chǎn)額、衰變子庫提供的數(shù)據(jù)，計算產(chǎn)生多群緩發(fā)光子產(chǎn)生矩陣，其計算理論如下。

多群緩發(fā)光子產(chǎn)生矩陣的計算公式可表示為：

(1)

其中：yi(h→g)為裂變核素i的第h群中子對第g群緩發(fā)光子的產(chǎn)額；Eh和Eh-1為中子能群邊界；p為裂變產(chǎn)物編號；yi→p為裂變核素i產(chǎn)生裂變產(chǎn)物p的產(chǎn)額，由裂變產(chǎn)額子庫給出；Ein為入射中子能量；Dg為多群衰變光子產(chǎn)額，根據(jù)衰變子庫中核素p衰變光子能譜的形式分別采用式(2)、(3)進(jìn)行計算，前者用于光子能譜為離散能級的情況，后者用于連續(xù)能譜情況。

(2)

(3)

(4)

(5)

其中，〈Ek〉為第k種衰變能譜的平均衰變光子能量。以JEFF-2.2評價核數(shù)據(jù)庫為例，部分核素的Ak值列于表1。由表1可看出，Ak值基本在1.0附近，表明了上述對衰變光子數(shù)據(jù)處理的準(zhǔn)確性。

表1 部分核素的Ak值Table1 Ak values of some nuclides

在評價核數(shù)據(jù)庫的458反應(yīng)道(MT=458)給出了核素裂變后釋放的衰變光子能量，為保證以上通過采用衰變數(shù)據(jù)計算的光子能量與458反應(yīng)道給出的總能量一致，將式(1)計算的多群緩發(fā)光子產(chǎn)生矩陣做進(jìn)一步處理：

Yi(h→g)=f(g)·yi(h→g)

(6)

(7)

其中，Qd為評價核數(shù)據(jù)庫中458反應(yīng)道給出的平均裂變緩發(fā)光子能量。對于未給出MT=458反應(yīng)道的評價核數(shù)據(jù)庫，緩發(fā)光子產(chǎn)生矩陣未進(jìn)行以上修正。

多群緩發(fā)光子產(chǎn)生矩陣計算流程如圖1所示。首先利用評價核數(shù)據(jù)庫衰變子庫給出的衰變光子能譜，基于式(2)、(3)計算多群衰變光子產(chǎn)額Dg；然后利用裂變產(chǎn)額子庫給出的獨(dú)立裂變產(chǎn)額，基于式(1)計算多群緩發(fā)光子產(chǎn)生矩陣yi(h→g)；最后為保證緩發(fā)光子釋熱計算的準(zhǔn)確性，利用中子子庫中給出的平均裂變緩發(fā)光子能量Qd對計算所得多群緩發(fā)光子產(chǎn)生矩陣進(jìn)行修正，得到最終的多群緩發(fā)光子產(chǎn)生矩陣Yi(h→g)。若評價核數(shù)據(jù)庫中未給出平均裂變緩發(fā)光子能量，則不對yi(h→g)進(jìn)行修正。

圖1 多群緩發(fā)光子產(chǎn)生矩陣計算流程Fig.1 Workflow to generate delayed photon library

2 光子離位輻照損傷截面計算功能

光子在輸運(yùn)過程中與堆內(nèi)材料原子發(fā)生反應(yīng)，包括康普頓散射、電子對效應(yīng)、光電效應(yīng)等。這些反應(yīng)將產(chǎn)生次級電子，電子與原子核發(fā)生碰撞后會造成原子離位損傷。因此，計算光子的離位損傷截面首先需建立電子離位輻照損傷截面計算模型，然后在此基礎(chǔ)上建立光子離位輻照損傷截面計算模型。需要指出的是，此處建立的電子離位輻照損傷截面計算方法，可直接單獨(dú)用于電子入射時產(chǎn)生的離位輻照損傷截面的計算。

2.1 電子離位損傷截面計算方法

電子與原子核發(fā)生碰撞的離位輻照損傷截面的計算公式為：

(8)

其中：υ(T)為材料的離位損傷函數(shù)，可采用K-P模型、NRT模型或ARC模型[11]計算獲得；T為反沖核的能量；Tmax為電子能產(chǎn)生的反沖核能量的最大值。

(9)

其中：E為電子能量；M為材料原子核質(zhì)量；m0為電子靜止質(zhì)量；c為光速。

式(8)中，dσ/dT為電子的微分散射截面，可表示為Mott截面形式：

(10)

σMott=σRR(θ,Z,E)

(11)

其中，σR為盧瑟福微分散射截面，可表示為：

(12)

對于R(θ,Z,E)，文獻(xiàn)[23]針對1 keV～900 MeV的電子與1≤Z≤118的元素反應(yīng)，提出一種插值多項式形式，如式(13)所示：

(13)

(14)

(15)

參照文獻(xiàn)[23]針對正負(fù)電子分別擬合獲得系數(shù)bk,j(Z)。最終，通過以上方法可計算不同入射電子能量E、不同散射角度θ下的電子微分散射截面。

2.2 光子離位損傷截面計算方法

光子輻照損傷截面可通過對3種光子-原子反應(yīng)道產(chǎn)生的離位損傷截面得到：

σt(Eγ)=σCS(Eγ)+σPE(Eγ)+σPP(Eγ)

(16)

其中：σCS(Eγ)、σPE(Eγ)、σPP(Eγ)分別為康普頓散射、光電效應(yīng)和電子對效應(yīng)對總光子離位損傷截面的貢獻(xiàn)；Eγ為入射光子能量。為獲得光子的離位損傷截面，需分別計算這3個分反應(yīng)道的離位損傷截面。

1) 康普頓散射

對于康普頓散射，其離位損傷截面可表示為：

(17)

(18)

式(17)中的積分上限可表示為：

(19)

式(17)中電子離位損傷函數(shù)可通過下式積分得到：

(20)

其中：NV為材料的原子核密度；σD(T)為電子離位損傷截面，通過式(8)計算；S(T)為電子阻止本領(lǐng)，本文采用ICRU1984電子組織本領(lǐng)數(shù)據(jù)庫[25]計算，對于該數(shù)據(jù)庫中未給出正電子組織本領(lǐng)的材料采用文獻(xiàn)[26]中的方法近似計算。

2) 光電效應(yīng)

對于光電效應(yīng)，電子的動能可通過以下公式計算：

E=Eγ-Be

(21)

其中，Be為電子結(jié)合能，一般是幾百eV的量級，相比于光子能量這個數(shù)值很小，可忽略不計。因此，對于一給定的入射光子能量，出射電子能量是一固定值，光電效應(yīng)的光子輻照損傷截面可表示為：

σPE(Eγ)=σPE(E)n(E)

(22)

其中：n(E)為電子離位損傷函數(shù)，與式(8)中υ(T)的計算方法相同；σPE(E)為光電效應(yīng)反應(yīng)截面，可用Hall公式[27]表示：

(23)

其中：α=Z/137；γ為Lorentzian系數(shù)，可表示為：

(24)

3) 電子對效應(yīng)

對于電子對效應(yīng)，離位損傷截面可通過下式計算：

(25)

(26)

其中：σ∞=5.18×10-28cm2；變量F(s)通過下式計算：

(27)

其中：l=2；n=8。其他變量的計算公式如下：

s=E/(Eγ-2m0c2)

(28)

u=ln(Eγ/m0c2)

(29)

g(u)=-0.183 5u3+1.653u2-

2.154 3u+0.761 4

(30)

h(u)=0.219 3u+0.182 5

(31)

3 數(shù)值驗證

3.1 多群緩發(fā)光子產(chǎn)生矩陣驗證

本文針對美國阿貢國家實驗室的EBR-Ⅱ鈉冷快堆，采用先進(jìn)反應(yīng)堆計算程序SARAX基于歐洲快堆計算程序ERANOS的數(shù)據(jù)庫進(jìn)行光子釋熱率計算。EBR-Ⅱ堆芯采用六邊形柵元布置，由16圈共637個組件構(gòu)成，其中組件的對邊距為5.892 9 cm。根據(jù)組件的分布，可將整個堆芯分成3部分，由內(nèi)到外分別為燃料區(qū)、反射層區(qū)和外增殖區(qū)。在燃料區(qū)內(nèi)除燃料組件外，還布置有8個控制棒組件和2個安全棒組件。ERANOS數(shù)據(jù)庫主要基于JEFF-2.2[28]研制，包含多群瞬發(fā)與緩發(fā)光子產(chǎn)生矩陣。本文采用相同的評價核數(shù)據(jù)庫制作多群瞬發(fā)和緩發(fā)光子數(shù)據(jù)庫，并與基于ERANOS數(shù)據(jù)庫的計算結(jié)果進(jìn)行對比。采用本文方法獲得的考慮緩發(fā)光子產(chǎn)生矩陣的數(shù)據(jù)庫計算的堆內(nèi)不同組件的光子功率及ERANOS的數(shù)據(jù)庫的結(jié)果及比較示于圖2。由圖2可見，兩者的最大相對偏差為0.93%，表明了本文方法的正確性。為體現(xiàn)緩發(fā)光子對光子功率計算結(jié)果的影響，計算了不考慮緩發(fā)光子和考慮緩發(fā)光子計算得到的光子功率的相對偏差，結(jié)果示于圖3，其中，燃料組件的最大偏差為32.91%，控制棒組件和反射層組件最大相對偏差為32.58%、20.41%。

圖2 組件光子功率計算結(jié)果Fig.2 Calculated result of photon heat for assembly

圖3 考慮與不考慮緩發(fā)光子對組件光子功率計算結(jié)果的影響Fig.3 Photon heat for assemby with and ignoring delayed photon contribution

為考慮緩發(fā)光子對商用壓水堆光子通量的影響，利用壓水堆燃料管理程序NECP-Bamboo[29]對VERA-2D基準(zhǔn)題的光子通量進(jìn)行計算，中子能群結(jié)構(gòu)選用WIMS-69群能群結(jié)構(gòu)，光子能群結(jié)構(gòu)選用美國洛斯阿拉莫斯國家實驗室給定的48群能群結(jié)構(gòu)[1]。不考慮緩發(fā)光子和考慮緩發(fā)光子計算得到的組件光子總通量示于圖4。從圖4可見，兩者最大相對偏差為11.68%。

圖4 VERA-2D組件光子總通量計算結(jié)果Fig.4 Photon flux calculation result for VERA-2D

3.2 電子離位輻照損傷截面的驗證

電子離位輻照損傷截面的計算是光子離位輻照損傷截面計算的基礎(chǔ)，因此，首先對其計算精度進(jìn)行驗證。選取美國橡樹嶺國家實驗室Oen等[30]的正、負(fù)電子離位損傷截面作為參考解，該實驗室采用K-P模型[5]用于材料原子離位損傷函數(shù)計算，NECP-Atlas計算中采用相同的模型。負(fù)電子、正電子引起的Fe、Au的離位輻照損傷截面示于圖5。由圖5可見，NECP-Atlas獲得的結(jié)果與Oen的結(jié)果符合良好。

3.3 光子離位輻照損傷截面驗證

選取Fukuya等[26]的結(jié)果作為參考解驗證NECP-Atlas計算的光子離位損傷截面的準(zhǔn)確性，結(jié)果示于圖6。由圖6可看出，NECP-Atlas計算得到的光子離位損傷截面與參考解吻合良好，證明了程序的準(zhǔn)確性。

a、b——Fe和Au的負(fù)電子離位損傷；c、d——Fe和Au的正電子離位損傷圖5 Fe和Au的負(fù)電子、正電子離位損傷截面Fig.5 Displacement cross section of Fe and Au caused by electron and positron

圖6 Fe的光子離位損傷截面Fig.6 Gamma displacement damage cross section of Fe

4 結(jié)論

本文基于核數(shù)據(jù)處理程序NECP-Atlas開展了光子相關(guān)數(shù)據(jù)計算方法研究，在程序中建立了完善的光子數(shù)據(jù)計算方法，除傳統(tǒng)核數(shù)據(jù)處理程序可產(chǎn)生的中子核反應(yīng)釋放的瞬發(fā)光子產(chǎn)生截面、光子與原子的反應(yīng)截面，針對NECP-Atlas程序新開發(fā)了裂變產(chǎn)物衰變釋放的緩發(fā)光子多群產(chǎn)生矩陣、光子離位輻照損傷截面計算功能。數(shù)值驗證結(jié)果表明，鈉冷快堆中緩發(fā)光子對非燃料組件的釋熱將產(chǎn)生顯著影響；針對光子離位輻照損傷截面的數(shù)值驗證證明了程序的正確性。此外，由于光子離位輻照損傷截面的計算需以電子離位輻照損傷截面的計算為基礎(chǔ)，因此，形成了電子離位輻照損傷截面的計算功能，并對該功能進(jìn)行了驗證。