湯勇強(qiáng)，祖鐵軍，易思宇，徐嘉隆，曹良志，吳宏春

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

熱能區(qū)中子散射數(shù)據(jù)對(duì)反應(yīng)堆設(shè)計(jì)、輻射屏蔽計(jì)算、臨界安全分析以及冷中子源設(shè)計(jì)等的計(jì)算結(jié)果具有重要影響。在熱能區(qū)，慢化劑靶核的化學(xué)鍵作用和熱運(yùn)動(dòng)都會(huì)強(qiáng)烈影響中子的散射行為，由于熱能區(qū)中子的能量與靶核相當(dāng)，中子在與靶核的作用過程中可能從靶核獲得能量或失去能量。與自由靶核原子相比，這些效應(yīng)將影響熱能區(qū)中子的散射截面，以及次級(jí)中子能量和角度分布。評(píng)價(jià)核數(shù)據(jù)庫(ENDF)采用熱散射律(TSL)數(shù)據(jù)描述以上效應(yīng)，由于國(guó)際上缺少熱散射實(shí)驗(yàn)，且熱散射實(shí)驗(yàn)無法對(duì)每種慢化材料都覆蓋整個(gè)熱能區(qū)的入射中子能量范圍和慢化材料的溫度范圍，因此，ENDF中的TSL數(shù)據(jù)需結(jié)合已有實(shí)驗(yàn)，通過物理模型進(jìn)行評(píng)價(jià)。ENDF中含有TSL數(shù)據(jù)的材料極為有限，最新發(fā)布的ENDF/B-Ⅷ.0中，僅24類常用慢化材料以及實(shí)驗(yàn)冷中子源中的冷慢化劑材料含有TSL數(shù)據(jù)[1]，這些材料的TSL數(shù)據(jù)主要是采用美國(guó)洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室核數(shù)據(jù)處理程序NJOY的LEAPR模塊計(jì)算獲得的。

目前ENDF中的TSL數(shù)據(jù)主要存在如下問題：1) 先進(jìn)反應(yīng)堆中使用的慢化劑材料缺少TSL數(shù)據(jù)，如熔鹽堆中的FLiBe；2) 材料的TSL數(shù)據(jù)與其晶體結(jié)構(gòu)具有很強(qiáng)的相關(guān)性，如先進(jìn)反應(yīng)堆內(nèi)的ZrH2材料，其中H核素含量不同使得其晶體結(jié)構(gòu)不同[2]，數(shù)值結(jié)果顯示，如果直接使用評(píng)價(jià)核數(shù)據(jù)庫中給出的ZrH2的TSL數(shù)據(jù)，會(huì)對(duì)計(jì)算結(jié)果造成一定的誤差。針對(duì)以上問題，需根據(jù)材料實(shí)際的晶體結(jié)構(gòu)計(jì)算產(chǎn)生TSL數(shù)據(jù)。目前國(guó)際上的研究多是采用NJOY程序的LEAPR模塊計(jì)算產(chǎn)生所需材料的TSL數(shù)據(jù)，但ENDF/B-Ⅷ.0中TSL的主要貢獻(xiàn)者之一北卡羅來納大學(xué)Hawari教授研究指出，LEAPR中使用的近似會(huì)引入較大的誤差[3]，限制了其使用范圍。

近年來，隨著反應(yīng)堆高保真計(jì)算方法的發(fā)展，核數(shù)據(jù)的精度成為限制反應(yīng)堆設(shè)計(jì)精度提高的主要問題，其中TSL數(shù)據(jù)也是提高核數(shù)據(jù)精度的重要內(nèi)容，因此，需要對(duì)TSL數(shù)據(jù)進(jìn)行深入研究。Hawari教授針對(duì)以上提到的LEAPR問題開發(fā)了FLASSH程序[3]。然而，F(xiàn)LASSH未提及對(duì)非固體靶核TSL的處理能力。

為滿足反應(yīng)堆設(shè)計(jì)對(duì)高精度核數(shù)據(jù)的需求，西安交通大學(xué)NECP實(shí)驗(yàn)室自主研發(fā)了核數(shù)據(jù)處理軟件NECP-Atlas[4]，本文在該程序的基礎(chǔ)上，開發(fā)通用的TSL數(shù)據(jù)計(jì)算模塊，去除傳統(tǒng)方法中的晶體立方近似和原子位置近似，采用基于各向異性位移參數(shù)(ADPs)[5]的TSL數(shù)據(jù)計(jì)算方法，可得到任意固體晶體結(jié)構(gòu)的TSL數(shù)據(jù)；在固體散射靶的模型基礎(chǔ)上，加入擴(kuò)散模型、離散諧振子模型以及舍爾德近似，使其可產(chǎn)生多原子分子固體材料、液體慢化材料的TSL數(shù)據(jù)。

1 理論方法

1.1 熱中子散射公式

熱能區(qū)中子在多原子的散射系統(tǒng)中，由于有多個(gè)原子參與散射，利用量子力學(xué)處理勢(shì)散射的問題時(shí)，需要用多體薛定諤方程來描述該多體系統(tǒng)。因此，一般引入第一玻恩近似和費(fèi)米贗勢(shì)求解多體薛定諤方程，熱能區(qū)中子散射的雙微分截面[6]可表示為：

(σcohS(α,β,T)+σincSs(α,β,T))

(1)

式中：E為入射中子能量；E′為出射中子能量；μ為實(shí)驗(yàn)室坐標(biāo)系下的散射角余弦；T為材料溫度；kB為玻爾茲曼常數(shù)；σcoh為束縛核的相干散射截面；σinc為束縛核的非相干散射截面；α為中子的無量綱動(dòng)量轉(zhuǎn)移量；β為中子的無量綱能量轉(zhuǎn)移量；S(α,β,T)為熱散射函數(shù)；Ss(α,β,T)為自散射函數(shù)。α和β分別表征中子散射前后角度和能量的變化量，可分別表示為：

(2)

(3)

熱散射函數(shù)S(α,β,T)可分成兩部分：

S(α,β,T)=Ss(α,β,T)+Sd(α,β,T)

(4)

其中：自散射函數(shù)Ss(α,β,T)計(jì)入非相干散射效應(yīng)；相互散射函數(shù)Sd(α,β,T)加上自散射函數(shù)，即S(α,β,T)計(jì)入相干散射效應(yīng)。

引入簡(jiǎn)諧近似，即原子偏離平衡位置的位移較小的情況下，晶體中原子間的相互作用可近似正比于原子的位移。可通過對(duì)散射函數(shù)的展開將熱中子散射律分解為彈性和非彈兩部分。對(duì)于晶體，這種展開方式稱為聲子展開。散射函數(shù)用聲子展開為：

(5)

(6)

其中，上標(biāo)0表示0聲子產(chǎn)生(或湮滅)，代表彈性散射，中子與靶核散射前后能量沒有改變；帶有其他更大數(shù)值上標(biāo)的項(xiàng)表示1，2，3等聲子項(xiàng)，表示有聲子產(chǎn)生(或湮滅)，中子散射后獲得能量或失去能量，代表非彈性散射項(xiàng)。熱能區(qū)中子散射可分為4部分：相干彈性散射、相干非彈性散射、非相干彈性散射和非相干非彈性散射。

通用形式的雙微分截面可進(jìn)一步表示為：

(7)

1.2 固體熱中子散射律數(shù)據(jù)計(jì)算模型

固體材料是由大量粒子(原子、離子等)所組成的具有強(qiáng)相互作用的集體。在一定溫度下，晶體格點(diǎn)粒子在平衡位置附近小振幅振動(dòng)，稱為晶格振動(dòng)。晶格振動(dòng)是由不同頻率的簡(jiǎn)正振動(dòng)疊加而成的一種集體振動(dòng)。對(duì)于給定的晶體，總的振動(dòng)模式數(shù)目是一定的，按振動(dòng)頻率有一分布，用晶格振動(dòng)模式密度(或頻率分布函數(shù))描述。從量子力學(xué)的觀點(diǎn)看，表征原子集體振動(dòng)的能量是量子化的，每個(gè)振動(dòng)模式能量的最小單位稱為聲子，因此晶格振動(dòng)模式密度又稱聲子態(tài)密度。根據(jù)聲子態(tài)密度，便可對(duì)固體慢化材料的TSL數(shù)據(jù)進(jìn)行研究[7]。

1) 固體非彈性散射

中子與所有慢化材料的散射過程中都存在非彈性散射，包括相干部分和非相干部分。一般來說，相互散射律部分較弱，引入非相干近似，即忽略相干散射項(xiàng)中的相互散射律部分，于是，聲子展開的非彈性散射可表示為：

σine(E→E′,μ,T)=

(8)

(9)

其中，ρ(β)為聲子態(tài)密度，是計(jì)算TSL數(shù)據(jù)的重要輸入?yún)?shù)，可通過實(shí)驗(yàn)方式或第一性原理計(jì)算獲得。

2) 固體非相干彈性散射

(10)

其中，w(T)為德拜-沃勒積分，表示評(píng)價(jià)數(shù)據(jù)庫中描述非相干彈性散射的TSL數(shù)據(jù)。

(11)

3) 固體相干彈性散射

(12)

計(jì)算相干彈性散射截面的關(guān)鍵問題是得到S0的計(jì)算表達(dá)式，其介紹詳見文獻(xiàn)[6]。從第一性原理出發(fā)，推導(dǎo)得到了便于計(jì)算的相干彈性散射截面表達(dá)式：

(13)

式中：N為晶胞中的原子數(shù)；v0為晶胞的體積；κ為倒易晶格矢量；τ為中子散射過程中的動(dòng)量轉(zhuǎn)移矢量；dj和σcoh,j分別為晶胞中第j個(gè)原子的位置和束縛態(tài)相干散射截面；e-2Wj(τ,T)為第j個(gè)原子的德拜-沃勒因子；Wj(τ,T)為德拜-沃勒系數(shù)。

傳統(tǒng)方法中，為簡(jiǎn)化德拜-沃勒因子的計(jì)算，采用了2個(gè)近似。首先是晶體立方近似，即假設(shè)在散射系統(tǒng)中每個(gè)原子的位移都是各向同性的，對(duì)于各原子可通過使用統(tǒng)一聲子態(tài)密度來計(jì)算德拜-沃勒系數(shù)。于是，德拜-沃勒系數(shù)可用統(tǒng)一的表達(dá)式表示：

(14)

其中，Mj為晶胞中第j個(gè)原子的質(zhì)量。接下來引入第2個(gè)近似，即原子位置近似。該近似假設(shè)晶胞中原子的位移與原子所處位置和原子類型無關(guān)。因此，德拜-沃勒系數(shù)可進(jìn)一步簡(jiǎn)化為W(τ,T)。采用傳統(tǒng)方法的相干彈性散射截面的最終形式為：

(15)

然而，對(duì)于非立方晶體，原子在不同位置不同方向所受的作用力是不相同的，這些力是相關(guān)的，原子在一個(gè)方向上的運(yùn)動(dòng)會(huì)導(dǎo)致在另一個(gè)方向上的力發(fā)生變化。原子的各向異性熱運(yùn)動(dòng)可用ADPs均方位移矩陣U(T)表示，U(T)是一3×3的實(shí)對(duì)稱矩陣，表示原子在平衡位置附近位移的均方振幅，具有長(zhǎng)度平方量綱。嚴(yán)格的通用固體相干彈性TSL數(shù)據(jù)可由式(16)獲得，稱為各向異性位移參數(shù)法[8-9]。

(16)

2π2τTUj(T)τ

(17)

傳統(tǒng)方法中，德拜-沃勒系數(shù)Wj(τ,T)僅與倒格子向量τ的長(zhǎng)度和材料溫度有關(guān)。采用ADPs方法的Wj(τ,T)不僅與倒格子向量τ的長(zhǎng)度有關(guān)，還與其方向有關(guān)。至此，嚴(yán)格的相干彈性散射截面計(jì)算方法推導(dǎo)完成。由于ADPs方法是基于不同原子均方位移計(jì)算的，因此，該方法適用于任何固體材料的任何晶體結(jié)構(gòu)。

1.3 液體或氣體的散射律計(jì)算模型

(18)

1) 有效寬度模型

有效寬度模型一般用于描述主散射體的擴(kuò)散項(xiàng)，如H2O中的H，其散射律函數(shù)解析表達(dá)式為：

(19)

其中：K1為第2類修正貝塞爾函數(shù)；ωt為平動(dòng)權(quán)重；c為無量綱擴(kuò)散常數(shù)。

2) 自由氣體模型

自由氣體模型一般用于描述次散射體的擴(kuò)散項(xiàng)，如H2O中的O，其表達(dá)式為：

(20)

1.4 多原子分子內(nèi)振動(dòng)模式

無論是固體還是液體散射靶，若是由多原子分子構(gòu)成，此時(shí)分子內(nèi)部的振動(dòng)不可忽略，需引入離散諧振子來描述多原子分子的內(nèi)部振動(dòng)模式。如甲烷中每個(gè)甲烷分子包含4個(gè)H原子和1個(gè)C原子，其分子內(nèi)部具有4個(gè)特征能量簡(jiǎn)正振動(dòng)，分別為162、190、361、374 meV。

(21)

其中：ωi和βi為第i個(gè)諧振子的權(quán)重和特征能量參數(shù)；In(x)為第1類修正貝塞爾函數(shù)。通過與式(18)相似的卷積公式，即可將離散振動(dòng)考慮到熱中子散射律中。

1.5 液體分子間的相干效應(yīng)

前文提到，對(duì)于非彈性散射一般可忽略其中相干散射中的相互散射律部分，然而對(duì)于束縛核相干散射截面σcoh占總束縛核散射截面σb(σb=σcoh+σinc)比例較大的核素，相干散射貢獻(xiàn)部分無法滿足非相干近似條件，需引入舍爾德修正因子s(κ)來考慮非彈性散射中的相互散射律部分。如Márquez-Damián提出的重水CAB模型[11]中，由于σD,coh/σD,b=5.592/7.64=0.731 94，在計(jì)算D2O中D的非彈性散射時(shí)，需考慮相干散射中相互散射律部分的影響。采用舍爾德近似[12]，相干TSLS(α,β,T)可由非相干TSLSs(α,β,T)來計(jì)算。

S(α,β,T)=Ss(α/s(κ),β,T)*s(κ)

(22)

(23)

非彈性散射截面為：

(24)

2 程序開發(fā)與數(shù)值結(jié)果

2.1 程序概述

基于以上熱中子散射理論，在核數(shù)據(jù)處理程序NECP-Atlas中開發(fā)了sab_calc模塊，用于產(chǎn)生慢化材料TSL數(shù)據(jù)，為反應(yīng)堆設(shè)計(jì)提供高精度的TSL數(shù)據(jù)。sab_calc模塊包含NJOY中LEAPR模塊的所有功能，且在相干彈性散射中去除傳統(tǒng)方法中的晶體立方近似和原子位置近似，采用ADPs方法，得到高精度的相干彈性TSL數(shù)據(jù)。利用sab_calc模塊產(chǎn)生的TSL數(shù)據(jù)，借助NECP-Atlas程序中其他已有的模塊，即可產(chǎn)生中子學(xué)計(jì)算所需的熱散射截面數(shù)據(jù)。

產(chǎn)生的TSL數(shù)據(jù)包括描述非彈性散射的Stot(α,β,T)、非相干彈性散射的德拜-沃勒積分w(T)，以及相干彈性散射的Si(E,T)。目前，國(guó)際上的ENDF均采用ENDF-6格式[14]存儲(chǔ)數(shù)據(jù)，為使所開發(fā)程序輸出的數(shù)據(jù)與ENDF一致，程序根據(jù)慢化材料的類型，將計(jì)算得到的TSL數(shù)據(jù)按ENDF-6格式進(jìn)行儲(chǔ)存。該格式按照主散射體的非彈性散射類別、主散射體的彈性散射是否為相干散射、是否有次散射體分成表1所列類型。

2.2 非彈性散射截面的驗(yàn)證

液態(tài)重水的非彈性散射極為復(fù)雜，包括連續(xù)譜項(xiàng)、擴(kuò)散項(xiàng)、離散振蕩頻譜和舍爾德近似等的處理，因此，本文以液態(tài)重水為例，將NECP-Atlas產(chǎn)生的非彈性散射TSL數(shù)據(jù)和熱散射截面與NJOY2016進(jìn)行對(duì)比，驗(yàn)證程序的正確性。

圖1為NECP-Atlas中sab_calc計(jì)算的D2O中D在293.6 K下的非彈性TSL數(shù)據(jù)，每條曲線表示1個(gè)α值下Stot(α,β,T)隨β的變化曲線?？煽吹剑讦梁苄r(shí)，曲線上出現(xiàn)了幾個(gè)峰值，證實(shí)了水分子的擴(kuò)散運(yùn)動(dòng)和水分子內(nèi)的振動(dòng)運(yùn)動(dòng)。圖2為NECP-Atlas采用D2O中D在293.6 K的TSL數(shù)據(jù)計(jì)算的非彈性散射概率密度及其與NJOY2016的相對(duì)偏差，圖中每條曲線代表1個(gè)入射能量。從圖2b可看出，NECP-Atlas與NJOY2016產(chǎn)生的散射概率密度符合良好，在小部分能區(qū)相對(duì)偏差略偏大，主要在準(zhǔn)彈性散射峰附近。準(zhǔn)彈性散射指中子散射前后能量的變化遠(yuǎn)小于中子本身的能量。相對(duì)偏差較大是因?yàn)樵跍?zhǔn)彈性散射峰附近，散射概率變化劇烈，需要更密集的出射能量網(wǎng)格才能滿足線性化收斂準(zhǔn)則，在NJOY2016程序內(nèi)存中對(duì)出射能量網(wǎng)格采用8位有效數(shù)字，而NJOY2016輸出的文件中，有效位數(shù)最多僅為7位，輸出過程中會(huì)產(chǎn)生偏差，由于輸出數(shù)據(jù)有效位數(shù)的限制，在程序內(nèi)部采取更高有效位數(shù)實(shí)際上是無意義的，NECP-Atlas在判斷收斂時(shí)，考慮了輸出數(shù)據(jù)有效位數(shù)的限制，避免了以上問題。以上不同的收斂策略導(dǎo)致NECP-Atlas與NJOY2016得到的散射概率相對(duì)偏差較大。D2O中D在293.6 K、入射能量0.025 3 eV的幾個(gè)出射能量點(diǎn)及散射概率列于表2?？煽吹剑贜JOY2016內(nèi)存中有5組有效數(shù)據(jù)，而在輸出文件中有效數(shù)據(jù)僅為(0.025 299 94，370.374 3)和(0.025 299 95，371.729 0)，證實(shí)在輸出過程中產(chǎn)生了偏差。為進(jìn)一步證明，圖3示出了兩者基于相同能量網(wǎng)格D2O中D在293.6 K、入射能量0.025 3 eV的非彈性散射概率的相對(duì)偏差，可見相對(duì)偏差均小于10-4，兩者結(jié)果吻合很好，證明由不同收斂策略產(chǎn)生的不同能量網(wǎng)格會(huì)帶來偏差。

表1 ENDF-6格式中熱散射律數(shù)據(jù)類型Table 1 Type of thermal scattering law data in ENDF-6 format

圖1 D2O中的D在293.6 K下的散射律Fig.1 TSL for D in D2O at 293.6 K

圖2 NECP-Atlas計(jì)算的D2O中D在293.6 K的非彈性散射概率密度(a)以及與NJOY2016的相對(duì)偏差(b)Fig.2 Thermal inelastic scattering probability for D in D2O at 293.6 K calculated by NECP-Atlas (a) and comparison with NJOY2016 (b)

表2 NJOY2016程序內(nèi)存和輸出文件數(shù)據(jù)對(duì)比Table 2 Comparison of values calculated by NJOY2016 in memory and output file

2.3 非相干彈性散射截面的驗(yàn)證

表3為本文計(jì)算的LiH中H的德拜-沃勒積分和有效溫度與文獻(xiàn)值的對(duì)比，可看出，本文計(jì)算的數(shù)據(jù)與文獻(xiàn)值符合良好。

2.4 相干彈性散射截面的驗(yàn)證

采用傳統(tǒng)的晶體立方近似方法和ADPs方法分別計(jì)算了金屬Be的相干彈性TSL數(shù)據(jù)(圖4a)，可看到，sab_calc采用傳統(tǒng)方法與LEAPR結(jié)果符合很好，采用ADPs方法會(huì)導(dǎo)致TSL數(shù)據(jù)產(chǎn)生最大10%的偏差。這種偏差會(huì)導(dǎo)致由TSL數(shù)據(jù)計(jì)算的相干彈性散射截面產(chǎn)生最大10%的偏差(圖4b)。

圖3 基于相同能量網(wǎng)格非彈性散射概率的相對(duì)偏差Fig.3 Comparison of calculated thermal inelastic scattering probability based on the same secondary energy grids

表3 LiH中H的德拜-沃勒積分和有效溫度Table 3 Debye-Waller integral and effective temperature for H in LiH

圖4 296 K下金屬Be的相干彈性TSL數(shù)據(jù)(a)和相干彈性散射截面(b)Fig.4 Coherent elastic TSL data (a) and coherent elastic scattering data (b) for metal-Be

為檢驗(yàn)以上相干彈性散射數(shù)據(jù)的偏差對(duì)反應(yīng)堆中子學(xué)計(jì)算結(jié)果造成的影響，使用NECP-Atlas程序基于ENDF/B-Ⅷ.0中子庫和計(jì)算的TSL數(shù)據(jù)制作蒙特卡羅程序使用的ACE格式數(shù)據(jù)庫，并采用MCNP程序計(jì)算了27個(gè)以金屬Be為慢化劑的國(guó)際臨界基準(zhǔn)題(ICSBEP)[17]HEU-MET-THERM-026(HMT026)。圖5示出了對(duì)金屬Be分別基于LEAPR、晶體立方近似方法的sab_calc和ADPs方法的sab_calc制作的熱散射庫，計(jì)算得到的基準(zhǔn)題keff，所有的統(tǒng)計(jì)偏差為24 pcm。結(jié)果表明，晶體立方近似方法的sab_calc結(jié)果和基于相同假設(shè)的LEAPR結(jié)果吻合，證明了sab_calc程序開發(fā)正確；其次，相對(duì)于LEAPR，sab_calc的ADPs方法對(duì)27個(gè)基準(zhǔn)題均有改進(jìn)，與實(shí)驗(yàn)基準(zhǔn)結(jié)果的偏差平均降低60 pcm，最大偏差為HMT026_19基準(zhǔn)題的135 pcm。

圖5 金屬Be采用ADPs方法的臨界基準(zhǔn)題keffFig.5 Calculated keff of critical benchmarks with metal beryllium using ADPs method

3 結(jié)論

基于熱中子散射理論，本文在核數(shù)據(jù)處理程序NECP-Atlas中開發(fā)了TSL數(shù)據(jù)產(chǎn)生模塊sab_calc，該模塊可產(chǎn)生任意晶體結(jié)構(gòu)的固體材料、部分液體材料的TSL數(shù)據(jù)，通過數(shù)值結(jié)果證明了程序的正確性。該程序采用ADPs方法，提高了相干彈性TSL數(shù)據(jù)的精度，基于ICSBEP基準(zhǔn)題的計(jì)算結(jié)果表明，采用ADPs方法獲得的金屬Be的TSL數(shù)據(jù)，與傳統(tǒng)晶體立方近似相比，會(huì)使keff的計(jì)算結(jié)果與實(shí)驗(yàn)基準(zhǔn)值的偏差平均降低約60 pcm。