核數據處理軟件NECP-Atlas中的光子相關數據計算方法研究

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

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

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

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

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

1 多群緩發光子產生矩陣計算功能

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

多群緩發光子產生矩陣的計算公式可表示為：

(1)

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

(2)

(3)

(4)

(5)

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

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

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

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

(6)

(7)

其中，Qd為評價核數據庫中458反應道給出的平均裂變緩發光子能量。對于未給出MT=458反應道的評價核數據庫，緩發光子產生矩陣未進行以上修正。

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

圖1 多群緩發光子產生矩陣計算流程Fig.1 Workflow to generate delayed photon library

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

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

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

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

(8)

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

(9)

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

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

(10)

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

(11)

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

(12)

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

(13)

(14)

(15)

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

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

光子輻照損傷截面可通過對3種光子-原子反應道產生的離位損傷截面得到：

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

(16)

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

1) 康普頓散射

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

(17)

(18)

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

(19)

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

(20)

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

2) 光電效應

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

E=Eγ-Be

(21)

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

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

(22)

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

(23)

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

(24)

3) 電子對效應

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

(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 數值驗證

3.1 多群緩發光子產生矩陣驗證

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

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

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

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

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

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

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

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

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

a、b——Fe和Au的負電子離位損傷；c、d——Fe和Au的正電子離位損傷圖5 Fe和Au的負電子、正電子離位損傷截面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 結論

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