祖鐵軍，湯勇強(qiáng)，曹良志，吳宏春

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

因具有高氫含量和高抗輻照損傷能力，氫化鋯(ZrHx)作為一種高質(zhì)量的慢化劑被廣泛應(yīng)用于TRIGA等反應(yīng)堆中[1]。氫化鋯的熱散射截面，特別是氫化鋯中氫的熱散射截面對反應(yīng)堆中子學(xué)計(jì)算有著重要影響，獲得高精度的氫化鋯熱散射截面是十分必要的。熱散射律數(shù)據(jù)是計(jì)算熱散射截面的基本數(shù)據(jù)。熱中子散射受材料晶體結(jié)構(gòu)的影響較大，不同氫含量的氫化鋯會有不同的晶體結(jié)構(gòu)，因此需根據(jù)氫化鋯實(shí)際晶體結(jié)構(gòu)計(jì)算其熱散射律數(shù)據(jù)[2]。

在基于聲子展開方法的熱散射律數(shù)據(jù)計(jì)算模型[3]中，基本的輸入?yún)?shù)是反映散射材料晶體結(jié)構(gòu)特性的聲子態(tài)密度。對于ZrHx中氫的聲子態(tài)密度，1968年，Slaggie等[4]建立了ZrH2的中心力晶格動力學(xué)模型(CF模型)，采用面心立方晶體結(jié)構(gòu)來近似ZrH2的微四方晶體結(jié)構(gòu)，然后通過擬合熱容數(shù)據(jù)得到聲子態(tài)密度。美國評價(jià)核數(shù)據(jù)庫ENDF/B采用了CF模型。2005年，Mattes等[5]建立了Debye-Gaussian(DG)模型，采用參數(shù)化的公式來描述ZrHx中氫的聲子態(tài)密度，認(rèn)為氫的聲子態(tài)密度由符合Debye分布的聲學(xué)部分和高斯分布的光學(xué)部分組成。歐洲評價(jià)核數(shù)據(jù)庫JEFF采用了DG模型。Zheng等[6]基于DG模型，以TRIGA反應(yīng)堆的實(shí)驗(yàn)結(jié)果為基礎(chǔ)，建立了ZrHx聲子態(tài)密度參數(shù)標(biāo)定方法。Malik等[7]建立了ZrH1.58中氫的三高斯模型。

上述模型的共同特點(diǎn)是采用參數(shù)化公式表示ZrHx中氫的聲子態(tài)密度?；跓嵘⑸浣孛?、反應(yīng)堆的反應(yīng)性等實(shí)驗(yàn)結(jié)果擬合得到聲子態(tài)密度公式中的參數(shù)，不能反映ZrHx的真實(shí)晶體結(jié)構(gòu)，會在熱散射截面計(jì)算中引入誤差。近年來，核材料研究領(lǐng)域通過第一性原理計(jì)算成功得到與實(shí)驗(yàn)吻合的ZrHx的晶體結(jié)構(gòu)[8-9]，計(jì)算中可得到相應(yīng)晶體結(jié)構(gòu)的聲子態(tài)密度。與參數(shù)化聲子態(tài)密度相比，第一性原理計(jì)算得到的聲子態(tài)密度能較好地反映ZrHx的真實(shí)晶體結(jié)構(gòu)。因此，本文利用第一性原理計(jì)算方法計(jì)算ZrHx中氫的聲子態(tài)密度，然后通過核數(shù)據(jù)處理程序NECP-Atlas[10-11]計(jì)算其熱散射律數(shù)據(jù)和熱散射截面，分析不同方法獲得的聲子態(tài)密度對ZrHx熱散射截面的影響以及對含ZrHx的TRIGA反應(yīng)堆反應(yīng)性的影響。

1 熱中子散射截面計(jì)算方法

熱中子散射的雙微分截面σ(E→E′，μ)可表示為：

σ(E→E′,μ)=

(1)

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

(2)

(3)

式中，A為原子與中子的質(zhì)量比。

對于固態(tài)晶體材料，原子在平衡位置附近做微小的振動。在這種情況下，晶體中原子間的相互作用力可認(rèn)為正比于原子的位移，原子間的力可假定為簡諧函數(shù)，稱之為簡諧近似。因此，熱中子散射截面可用聲子展開模型表示為：

(4)

對于ZrHx中的氫，在熱能區(qū)有3種類型的散射：非彈性散射、非相干彈性散射和相干彈性散射。目前的評價(jià)核數(shù)據(jù)庫中忽略ZrHx中氫的相干彈性散射項(xiàng)。因此，為保持與評價(jià)核數(shù)據(jù)庫一致，以下只討論非彈性散射和非相干彈性散射。

1.1 非彈性散射

對于非彈性散射過程，散射系統(tǒng)中至少有1個(gè)聲子產(chǎn)生或者湮滅，因此，保留式(4)中n≥1項(xiàng)，且采用非相干近似[3]，非彈性散射的雙微分截面σine(E→E′，μ)可表示為：

(5)

(6)

(7)

1.2 非相干彈性散射

對于非相干彈性散射，散射系統(tǒng)中沒有聲子產(chǎn)生或湮滅，只保留n=0項(xiàng)，所以非相干彈性散射截面σel，inc(E，μ)計(jì)算表達(dá)式為：

(8)

式中，2wE為德拜-沃勒積分，w由以下公式計(jì)算：

(9)

2 第一性原理計(jì)算及結(jié)果

慢化劑材料ZrHx的氫含量一般在1.4

分別創(chuàng)建δ-ZrH1.510個(gè)原子(Zr4H6)的晶胞和ε-ZrH212個(gè)原子(Zr4H8)的晶胞，如圖1所示；采用第一性原理材料模擬軟件包VASP[12]對以上兩個(gè)晶胞進(jìn)行幾何優(yōu)化和能量評價(jià)，獲得晶胞最低能量時(shí)的晶格常數(shù)；本文計(jì)算的δ-ZrH1.5的晶格常數(shù)為a=b=c=0.476 800 nm；ε-ZrH2的晶格常數(shù)為a=b=0.500 895 nm，c=0.445 270 nm。然后，采用晶格動力學(xué)計(jì)算程序PHONOPY[13]分別將δ-ZrH1.5和ε-ZrH2的晶胞擴(kuò)展為2×2×2的超胞，然后用VASP中的密度泛函微擾理論(DFPT)方法計(jì)算超胞的原子間力常數(shù)。在VASP的計(jì)算中，用投影綴加波(PAW)方法描述電子離子相互作用，采用Perdew-Burke-Ernzerhof(PBE)公式中的廣義梯度近似(GGA)來描述交換和關(guān)聯(lián)勢。布里淵區(qū)k點(diǎn)采用以Γ為中心的Monkhost-Pack方案取樣，k點(diǎn)間距最小值為2π×0.04 ?-1。采用高斯彌散(Gaussian smearing)方法，并且彌散寬度取為0.05 eV。平面波基波的截止能量為500 eV。最后，根據(jù)VASP計(jì)算得到的原子間力常數(shù)，使用PHONOPY程序來計(jì)算聲子態(tài)密度。

圖1 δ-ZrH1.5(a)和ε-ZrH2 (b)的晶胞Fig.1 Unit cell of δ-ZrH1.5 (a) and ε-ZrH2 (b)

圖2給出基于第一性原理計(jì)算獲得的ε-ZrH2聲子態(tài)密度和ENDF/B-Ⅷ.0中ZrH2的CF模型、JEFF-3.3中ZrHx的DG模型的聲子態(tài)密度光學(xué)部分，并與Evans等[14]獲得的ZrH2中氫的聲子態(tài)密度光學(xué)部分的實(shí)驗(yàn)結(jié)果進(jìn)行了比較。可見，采用第一性原理計(jì)算的ε-ZrH2聲子態(tài)密度與實(shí)驗(yàn)結(jié)果吻合較好，很好地復(fù)現(xiàn)了在0.136和0.144 eV處的兩個(gè)光學(xué)峰；在0.15～0.16 eV范圍內(nèi)，ε-ZrH2聲子態(tài)密度與實(shí)驗(yàn)數(shù)據(jù)也有較好的一致性，兩者的聲子態(tài)密度都有一個(gè)明顯的平臺。CF模型在以上兩個(gè)能量也出現(xiàn)了峰，但峰的位置與實(shí)驗(yàn)結(jié)果有偏差。DG模型則只有一個(gè)光滑的峰。

圖2 ZrH2中氫聲子態(tài)密度光學(xué)部分Fig.2 Optical phonon density of state of H in ZrH2

圖3比較了第一性原理計(jì)算獲得的δ-ZrH1.5和ε-ZrH2中氫的聲子態(tài)密度、JEFF-3.3中DG模型、ENDF/B-Ⅷ.0中CF模型的聲子態(tài)密度。ZrHx中氫的聲子態(tài)密度分成兩個(gè)區(qū)域：低能區(qū)的聲學(xué)部分和高能區(qū)的光學(xué)部分。由于聲學(xué)部分?jǐn)?shù)值較小，圖3中將聲學(xué)部分的數(shù)值乘以100。DG模型假設(shè)聲子態(tài)密度由Debye溫度為20 meV的聲學(xué)部分和高斯分布(在0.137 eV處達(dá)到峰值)的光學(xué)部分組成，因此聲子態(tài)密度曲線是平滑變化的。δ-ZrH1.5聲學(xué)部分的結(jié)果在形狀和量級上與CF模型基本吻合，但峰值處的能量低于CF模型，兩種模型出現(xiàn)峰值的能量點(diǎn)分別為0.015 eV和0.017 eV。對于ε-ZrH2的聲學(xué)部分，主峰出現(xiàn)在20 meV處，與DG模型一致。在光學(xué)部分，δ-ZrH1.5的聲子態(tài)密度振蕩大于其他模型。Malik等[7]將ZrHx中氫的聲子態(tài)密度光學(xué)部分被表示為1個(gè)三高斯模型，3個(gè)高斯分布的中心點(diǎn)分別在0.132、0.137和0.151 eV處，而δ-ZrH1.5的計(jì)算結(jié)果顯示在0.135、0.139和0.152 eV處有3個(gè)峰值，這與三高斯模型一致。

圖3 4種模型聲子態(tài)密度的比較Fig.3 Comparison of phonon density of state among four models

3 熱散射截面及TRIGA反應(yīng)堆計(jì)算結(jié)果

3.1 熱散射截面

基于以上獲得聲子態(tài)密度，利用核數(shù)據(jù)處理程序NECP-Atlas中sab_calc模塊[15]計(jì)算了ZrHx中氫的熱散射律數(shù)據(jù)，然后采用therm_calc模塊[11]計(jì)算了熱散射截面。將截面計(jì)算結(jié)果和國際原子能機(jī)構(gòu)核反應(yīng)實(shí)驗(yàn)數(shù)據(jù)庫EXFOR[16]中的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了詳細(xì)對比，由于篇幅原因，下文僅給出部分結(jié)果。

圖4示出能量范圍為1.07×10-3～2.17×10-2eV、5.3×10-2～5.3×10-1eV和熱能區(qū)的ZrH1.5中氫的總散射截面的計(jì)算結(jié)果及實(shí)驗(yàn)結(jié)果[16-17]。由圖4a可看出，在此能量范圍內(nèi)基于不同模型的聲子態(tài)密度計(jì)算的總散射截面無明顯差異。由圖4b可看出，ε-ZrH2、CF模型和DG模型的結(jié)果基本吻合，δ-ZrH1.5與其他模型計(jì)算結(jié)果存在差異。

圖5示出入射能量Ei=0.232 eV、散射角θ=60°、90°和120°時(shí)ZrH2的雙微分散射截面的計(jì)算值及實(shí)驗(yàn)結(jié)果[18]。由于雙微分截面描述的是在某一入射能量下，中子在特定能量、角度出射的概率，對晶體結(jié)構(gòu)較為敏感，所以相對于圖4a、b展示的總散射截面，不同的聲子態(tài)密度計(jì)算模型體現(xiàn)出明顯的差異，且由于δ-ZrH1.5的結(jié)構(gòu)實(shí)驗(yàn)樣品不一致，基于δ-ZrH1.5聲子態(tài)密度計(jì)算的雙微分散射截面有明顯偏差。基于ε-ZrH2的聲子態(tài)密度計(jì)算的雙微分散射截面與實(shí)驗(yàn)值吻合最好，其次是CF模型和DG模型。

能量范圍：a——1.07×10-3～2.17×10-2 eV；b——5.3×10-2～5.3×10-1 eV；c——熱能區(qū)圖4 ZrH1.5中氫在不同能量范圍的總散射截面Fig.4 Total scattering cross section of H in ZrH1.5 in different energy regions

a——Ei=0.232 eV，θ=60°；b——Ei=0.232 eV，θ=90°；c——Ei=0.232 eV，θ=120°圖5 ZrH2的雙微分散射截面Fig.5 Double differential scattering cross section of ZrH2

3.2 TRIGA反應(yīng)堆計(jì)算結(jié)果

利用基于不同聲子態(tài)密度獲得的熱散射截面，對國際臨界安全基準(zhǔn)題庫ICSBEP[19]中含ZrHx的TRIGA反應(yīng)堆基準(zhǔn)題進(jìn)行了計(jì)算，計(jì)算的問題列于表1。計(jì)算時(shí)ZrHx中鋯采用的數(shù)據(jù)為：ENDF/B-Ⅷ.0采用其中給出的ZrHx中鋯的熱散射律數(shù)據(jù)；JEFF-3.3中無鋯的熱散射律數(shù)據(jù)則采用自由氣體模型；本文采用不同晶體結(jié)構(gòu)分別獲得ZrHx中鋯的熱散射律數(shù)據(jù)用于相應(yīng)計(jì)算。研究中發(fā)現(xiàn)ZrHx中Zr的熱散射律數(shù)據(jù)對結(jié)果影響很小。其他所有核素的截面數(shù)據(jù)均來自ENDF/B-Ⅷ.0。采用蒙特卡羅程序MCNP進(jìn)行中子學(xué)計(jì)算，計(jì)算中有效增殖因數(shù)keff的統(tǒng)計(jì)偏差控制在±10 pcm以內(nèi)。計(jì)算的keff與實(shí)驗(yàn)值比較列于表2。表2中：ICT003和ICT013基準(zhǔn)題中ZrHx處于δ相，所以計(jì)算中使用了基于δ-ZrH1.5聲子態(tài)密度的截面；HCM003基準(zhǔn)題中ZrHx處于ε相，采用基于ε-ZrH2聲子態(tài)密度的截面進(jìn)行計(jì)算。

表1 ICSBEP臨界基準(zhǔn)題中氫/鋯原子比Table 1 H/Zr atom ratio in ICSBEP criticality benchmark

由表2可見，CF和DG模型對所有基準(zhǔn)題結(jié)果比較吻合?；诘谝恍栽砺曌討B(tài)密度的熱散射截面對計(jì)算結(jié)果有明顯改善，特別是對于含有δ-ZrH1.5的ICT003和ICT013基準(zhǔn)題。與CF模型相比，采用第一性原理計(jì)算的聲子態(tài)密度可使ICT003_02基準(zhǔn)題偏差絕對值最大減少332 pcm；與DG模型相比，也可使ICT003_02基準(zhǔn)題偏差絕對值最大減少448 pcm。

表2 計(jì)算的keff與實(shí)驗(yàn)值比較Table 2 Comparison between calculated keff and experiment value

對于HCM003系列基準(zhǔn)題，本文計(jì)算結(jié)果與實(shí)驗(yàn)值的偏差控制在100 pcm以內(nèi)。

為進(jìn)一步展示不同模型對中子能譜的影響，對ICT003_02基準(zhǔn)題中A1燃料棒的中子能譜進(jìn)行了對比，如圖6所示。由圖6可見，CF模型和DG模型的中子能譜無明顯差異，而本文結(jié)果相對于前兩者，在0.08～0.16 eV范圍內(nèi)中子通量有所增加，出現(xiàn)了能譜硬化。

圖6 ICT003_02基準(zhǔn)題中A1燃料棒的中子通量Fig.6 Neutron flux of A1 fuel rod in ICT003_02 benchmark

4 結(jié)論

采用第一性原理計(jì)算了氫化鋯的δ-ZrH1.5和ε-ZrH2兩種晶體結(jié)構(gòu)中氫的聲子態(tài)密度，以此計(jì)算了ZrHx中氫的熱散射律數(shù)據(jù)和熱散射截面，然后將熱散射截面應(yīng)用到國際臨界安全基準(zhǔn)題中含ZrHx的TRIGA反應(yīng)堆的中子學(xué)計(jì)算，分析了不同聲子態(tài)密度計(jì)算方法對熱中子散射截面、TRIGA反應(yīng)堆反應(yīng)性的影響。結(jié)果表明：基于實(shí)際晶體結(jié)構(gòu)，由第一性原理計(jì)算得到的聲子態(tài)密度可給出更精確的熱中子散射截面數(shù)據(jù)；與ENDF/B-Ⅷ.0和JEFF-3.3評價(jià)核數(shù)據(jù)庫中模型相比，基于第一性原理獲得的聲子態(tài)密度可明顯提高TRIGA反應(yīng)堆反應(yīng)性的計(jì)算精度。