祖鐵軍，湯勇強，曹良志，吳宏春

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

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

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

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

1 熱中子散射截面計算方法

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

σ(E→E′,μ)=

(1)

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

(2)

(3)

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

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

(4)

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

1.1 非彈性散射

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

(5)

(6)

(7)

1.2 非相干彈性散射

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

(8)

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

(9)

2 第一性原理計算及結果

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

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

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

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

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

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

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

3 熱散射截面及TRIGA反應堆計算結果

3.1 熱散射截面

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

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

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

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

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

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

表2 計算的keff與實驗值比較Table 2 Comparison between calculated keff and experiment value

對于HCM003系列基準題，本文計算結果與實驗值的偏差控制在100 pcm以內(nèi)。

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

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

4 結論

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