氮化鈾熱中子截面的第一性原理計(jì)算

王立鵬 江新標(biāo) 吳宏春 樊慧慶

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

2)(西北核技術(shù)研究所,西安 710024)

3)(西北工業(yè)大學(xué)材料學(xué)院,西安 710072)

氮化鈾(UN)因其較好的熱物性和耐事故容錯性成為先進(jìn)動力堆的候選燃料,但目前熱能區(qū)缺少可靠的UN熱中子截面數(shù)據(jù),這對于熱中子反應(yīng)堆物理計(jì)算是很不利的.本文基于量子力學(xué)的第一性原理,利用VASP/PHONON軟件模擬計(jì)算了UN的聲子態(tài)密度,以此為積分得到UN的定容比熱容,并基于新制作的聲子態(tài)密度,采用核截面處理程序NJOY/LEAPR,利用熱中子散射理論,得到UN的S(α,β)數(shù)據(jù),進(jìn)而研究UN的熱中子散射截面,并與傳統(tǒng)壓水堆的二氧化鈾(UO2)進(jìn)行對比.結(jié)果表明:優(yōu)化的晶格參數(shù)與數(shù)據(jù)庫符合較好,UN聲子態(tài)密度的聲子項(xiàng)和光子項(xiàng)較UO2的分隔更加明顯,定容比熱容計(jì)算結(jié)果與實(shí)驗(yàn)值一致,基于該聲子態(tài)密度計(jì)算得到的UN中238U的非彈性散射和彈性散射截面比相同溫度下UO2中238U小,UN中N僅考慮了非相干散射部分,隨著溫度升高,UN彈性散射截面變小,非彈性散射變大,并在高能段趨于自由核散射截面.本文的研究結(jié)果填補(bǔ)了UN熱中子截面數(shù)據(jù)的缺失,為下一步系統(tǒng)研究UN燃料在輕水堆中的中子學(xué)性能奠定了基礎(chǔ).

1 引 言

相較于二氧化鈾(UO2)燃料,氮化鈾(UN)具有鈾密度高、熔點(diǎn)高、熱導(dǎo)率高、熱膨脹系數(shù)低、輻照穩(wěn)定性好、裂變氣體釋放率低等優(yōu)點(diǎn),成為先進(jìn)動力堆的候選燃料,也是新型的耐事故容錯燃料,具有較好的發(fā)展前景[1,2].熱中子反應(yīng)堆設(shè)計(jì)計(jì)算使用的核材料的熱中子截面對反應(yīng)堆臨界安全特性、中子能譜等都會產(chǎn)生較大影響,需要較為可靠的熱中子截面數(shù)據(jù)才能準(zhǔn)確地計(jì)算出堆芯的物理參數(shù).目前國際上已建立的多種常用慢化劑材料的熱中子截面產(chǎn)生技術(shù)和處理方法大多數(shù)基于半經(jīng)驗(yàn)提出的簡化模型,再引入較多近似從而得到熱中子截面數(shù)據(jù).近年來,國外發(fā)展的先進(jìn)模擬與仿真技術(shù),依托高性能計(jì)算能力,以“第一性原理”的物理學(xué)為基礎(chǔ)建立模型,取代經(jīng)驗(yàn)公式,已經(jīng)實(shí)現(xiàn)了許多復(fù)雜系統(tǒng)的預(yù)測性模擬,包括熱中子截面的模擬計(jì)算[3?5].以往的熱中子截面庫僅給出了傳統(tǒng)的UO2燃料的熱中子截面數(shù)據(jù),比如MCNP5自帶的熱中子截面ENDF70SAB數(shù)據(jù)文件,缺少UN燃料的熱中子數(shù)據(jù)[6],這樣就會給以UN為燃料的先進(jìn)輕水堆或者其他類型熱中子反應(yīng)堆的物理計(jì)算帶來較大誤差.最新版的ENDF/B VIII.0中,美國北卡羅來納州立大學(xué)的LEIP實(shí)驗(yàn)室利用“in-house code”制作了UN的TSL(thermal scattering library)文件庫[7?9],相干彈性散射部分歸并到U核內(nèi),不相干彈性散射歸并到14N核內(nèi),并對NJOY程序[10,11]中LEAPR模塊進(jìn)行了大幅度修改,其中相干彈性散射部分采用了Cubic approximation和Exact Debye-Waller Matrix方法,僅依靠現(xiàn)有TSL庫,沒有相關(guān)程序,無法對其進(jìn)行更加深入的研究.因此,本文基于UN熱中子數(shù)據(jù)的迫切需求,研究了UN精確聲子態(tài)密度的產(chǎn)生方法和驗(yàn)證,對NJOY現(xiàn)有模型稍加改動,使其適用于對UN的研究,以補(bǔ)充UN燃料的熱中子截面數(shù)據(jù).

熱能區(qū)的中子散射截面不單純與中子能量有關(guān),還與散射介質(zhì)的溫度及物理、化學(xué)性質(zhì)有關(guān),它反映的是材料自身的聲子態(tài)密度特征,熱中子截面產(chǎn)生技術(shù)與處理方法是一項(xiàng)多學(xué)科相互交叉的難題.在熱中子反應(yīng)堆中,對于能量低于1 eV的中子,由于中子能量與散射核的熱能是可比較的,不能簡單認(rèn)為靶核是靜止的.由于熱中子與運(yùn)動的靶核發(fā)生碰撞時(shí),能從靶核獲得能量,所以,出射中子的能量可能大于入射中子的能量,這即是熱能區(qū)中子的向上散射現(xiàn)象.在分子或固體中,散射核與鄰近核之間存在著相互作用,原子核處于束縛狀態(tài),與中子發(fā)生碰撞時(shí)不能自由地反沖.由于較低能量中子的德布羅意波長可與分子或晶體內(nèi)核的間距相比較,與不同核發(fā)生散射的中子間可能發(fā)生干涉效應(yīng)[12,13].一般用熱散射S(α,β)來表示熱中子截面,ENDF格式評價(jià)庫中有對熱中子截面專門的TSL文件描述[14],特定的評價(jià)數(shù)據(jù)庫僅給出了少數(shù)常用慢化劑材料的TSL文件,一般用戶無法根據(jù)自己的需要處理評價(jià)庫以外的材料,要制作UN的熱中子散射截面數(shù)據(jù)就需要UN準(zhǔn)確的聲子態(tài)密度和可靠的熱中子截面處理方法.

2 熱中子散射理論

熱中子相干散射和非相干散射包括彈性散射和非彈性散射部分,彈性散射不會帶來系統(tǒng)能量的變化,但是由于熱中子能量與晶格振動能態(tài)相當(dāng),認(rèn)為熱中子彈性散射是整個晶格的散射,這樣靶核的有效質(zhì)量將會很大,系統(tǒng)在中子散射過程中不會失去能量.相反,熱中子的非彈性散射會帶來能量的損失或者增加,這主要是由碰撞核處于激發(fā)態(tài)引起的,即非彈性散射會伴隨原子的一個或多個聲子的發(fā)射或吸收.雖然熱中子的非彈性散射不會引起整個靶核本身的激發(fā)態(tài),但可以引起分子(或晶體)中原子的一個或幾個振動量子態(tài)的改變.

熱中子截面通?？梢苑譃?部分[15].

1)非彈性散射:包括相干和非相干非彈性散射,對所有物質(zhì)都重要,用熱散射律表示.非彈性散射的雙微分散射截面是散射律S(α,β)的函數(shù):

其中E和E′是實(shí)驗(yàn)系中入射和散射中子的能量,?是實(shí)驗(yàn)系下的散射角度,σcoh是材料的束縛核的相干散射截面,σinc是材料的束縛核的非相干散射截面,α為動量轉(zhuǎn)移量,β為能量轉(zhuǎn)移量,S(α,β)為熱中子散射律.表達(dá)式分別為

其中A表示散射核的原子質(zhì)量M和中子質(zhì)量的比值;κ為散射矢量;kT是溫度,單位是eV;~為普朗克常數(shù);ε為能量變化量;SS(α,β)為不考慮相互作用的自散射律,Sd(α,β)為考慮相干效應(yīng)的分立散射律.由(1)式可以看出,要求非彈性散射截面需要知道束縛核的截面以及相應(yīng)的散射律.以往在熱中子散射計(jì)算中通常引入“非相干近似”,比如NJOY程序中的LEAPR模塊,即認(rèn)為(4)式中的Sd(α,β)為零,這樣(1)式變?yōu)?/p>

對于引入非相干近似后的 (5)式,認(rèn)為SS(α,β)是高斯分布的函數(shù),這時(shí),熱中子散射律SS(α,β)可以寫成

其中的ρ(β)是以β為單位的聲子態(tài)密度,并且是歸一化的.

當(dāng)α和|β|值較大時(shí)(α > αsw,|β|> βsw), 即在入射能量較高時(shí),α和β已經(jīng)超出了S(α,β)范圍,這時(shí)就需要引入短時(shí)間碰撞近似(SCT)[10,11,15],從而熱中子的非彈性散射表達(dá)變?yōu)?/p>

進(jìn)而得到“短時(shí)間碰撞下”的雙微分截面為

2)非相干彈性散射:對含氫的固態(tài)物質(zhì)重要,如ZrHx、聚乙烯、固態(tài)輕水等.在含氫固體中,有一種由于零階聲子項(xiàng)產(chǎn)生的彈性散射(不造成能量損失),稱為“非相干彈性散射”,對于非相干彈性散射,其雙微分散射截面為

μi是第i個等概率區(qū)間的散射角余弦的上限,ˉμi是該區(qū)間內(nèi)平均散射角余弦值,N是等概率區(qū)間的個數(shù),并且μ0=?1.

3)相干彈性散射:對晶體重要,如石墨、鈹和UO2等.在包含相干散射的固體中,組成固體的晶體不同平面的原子會發(fā)生干涉散射,在ENDF庫中,這一過程被稱為“相干彈性散射”,因?yàn)闆]有能量的損失.

多晶材料的雙微分相干彈性散射表達(dá)式如下:

其中σb是材料有效的束縛核相干散射截面,Ei是所謂的布拉格閾值,fi是對應(yīng)的晶格學(xué)結(jié)構(gòu)因子,不同晶格結(jié)構(gòu)的fi不同.

3 計(jì)算方法與建模

氮化鈾燃料的原子結(jié)構(gòu)如圖1所示,它與堿金屬鹵化物NaCl,KCl和MgO的結(jié)構(gòu)類似,這種結(jié)構(gòu)被稱之為面心立方結(jié)構(gòu)(FCC),它的原胞結(jié)構(gòu)中只包含兩個原子,每個原子被6個其他原子所包圍.UN燃料結(jié)合了金屬燃料和氧化物燃料的雙重優(yōu)點(diǎn),既有像金屬燃料一樣的高熱導(dǎo)率和高密度,又有像氧化物燃料一樣的高熔點(diǎn)和較高的結(jié)構(gòu)完整性.

圖1 UN原子結(jié)構(gòu)圖Fig.1.Atom structure of UN.

利用美國MaterialDesign公司研制的MedeA[16]材料計(jì)算平臺,通過調(diào)用平臺下的VASP和PHONON軟件完成UN材料的聲子態(tài)密度計(jì)算.本文首先計(jì)算了UN在基態(tài)時(shí)的能量與結(jié)構(gòu),交換關(guān)聯(lián)函數(shù)使用Perdew在1991年提出的梯度密度修正近似GGA(Gradient Corrected Approximation)[17],U和N的截?cái)嗄芰窟x取為400 eV,計(jì)算采用周期性超晶格方法,UN立方相的布里淵區(qū)積分在6×6×6的Monkhost-Pack格子中進(jìn)行.根據(jù)體系的周期性,移動每個原子的位置,計(jì)算出原胞中所有原子受力,利用第一性原理線性響應(yīng)理論計(jì)算Hellmann-Feynman(HF)力常數(shù),進(jìn)而得到UN的聲子譜,流程如圖2所示.

圖2 VASP/PHONON產(chǎn)生聲子態(tài)密度的流程圖Fig.2.Flowchart of phonon density of states generation in VASP/PHONON code.

熱力學(xué)函數(shù)內(nèi)能(E)、熵(S)、自由能(F)、恒容熱容(Cv)與聲子譜密切相關(guān),一旦得到了UN的聲子譜,這些熱力學(xué)參數(shù)就可以在簡諧近似的模型下確定,因?yàn)檫@些參數(shù)主要采用了聲子譜作為它們積分的權(quán)重譜,其中Cv是比較重要的參數(shù),因?yàn)樗梢栽趯?shí)驗(yàn)中準(zhǔn)確測量,通常用它來檢驗(yàn)聲子譜計(jì)算的準(zhǔn)確性,在PHONON/VASP程序的計(jì)算中,比熱容有如下表示:

其中r是晶胞內(nèi)的自由度個數(shù),kB是玻爾茲曼常數(shù),T是溫度,~是簡化普朗克常數(shù).

因?yàn)?4N在熱能區(qū)具有較高的熱中子吸收截面,而且還會產(chǎn)生半衰期較長的14C,所以在UN的實(shí)際應(yīng)用中,多采用高富集度的U15N.UN中U的熱中子分析采用和UO2一致的方法,僅處理UN中U的238U核素,表1所列為UN中238U,14N和15N各核素原子量以及各自束縛核的相干和非相干散射截面[18],其中自由散射截面由下式得到:

從表1中數(shù)據(jù)可以看出,238U采用相干彈性散射截面和非相干非彈性散射,14N和15N采用非相干彈性散射截面是比較合理的.因此,在NJOY/LEAPR的計(jì)算中,僅需要對UN中238U的相干彈性散射截面的晶格學(xué)結(jié)構(gòu)因子進(jìn)行重新構(gòu)置,其他參數(shù)按照表中數(shù)據(jù)和聲子態(tài)密度進(jìn)行輸入,而對于非彈性散射截面部分,則采用程序自帶的“非相干近似”忽略其中的相干部分.

表1 UN各核素束縛核的不同反應(yīng)截面Table 1.Dif f erent cross sections of isotopes in UN.

4 計(jì)算結(jié)果與討論

4.1 晶格參數(shù)的比較

UN的幾何結(jié)構(gòu)優(yōu)化采用平面波贗勢程序VASP進(jìn)行,初始參數(shù)選取Pearson數(shù)據(jù)庫中的實(shí)驗(yàn)值,對晶格參數(shù)和原子位置進(jìn)行模擬,利用VASP程序包中的structure optimization選項(xiàng)將UN結(jié)構(gòu)松弛到它的最低能量狀態(tài),電子能量收斂的標(biāo)準(zhǔn)選取為1×10?5eV,平面波截?cái)嗄転?00 eV,分別用廣義梯度近似GGA和局域密度泛函LDA進(jìn)行計(jì)算,結(jié)果如表2所列.可以看出,采用GGA的贗勢優(yōu)化得到的晶格參數(shù)更符合真實(shí)值.因此,接下來基于GGA贗勢優(yōu)化得到的結(jié)構(gòu)進(jìn)行聲子譜的計(jì)算.

表2 晶格參數(shù)的對比Table 2.Comparison in lattice parameters.

4.2 聲子譜的比較

UN聲子色散關(guān)系在布里淵區(qū)不同方向的計(jì)算結(jié)果如圖3所示,較低的分支為聲子項(xiàng),較高的分支為光子項(xiàng).聲子態(tài)密度的計(jì)算結(jié)果如圖4所示,與聲子色散關(guān)系結(jié)論一致,低頻的聲學(xué)支主要是體系的整體運(yùn)動,光學(xué)支是原子間的相互運(yùn)動,而且UN的聲子項(xiàng)和光子項(xiàng)分隔較為明顯.將UO2的分聲子態(tài)密度與UN進(jìn)行比較[19],如圖5和圖6所示.由圖5可以看出兩者的U分聲子態(tài)密度相差不多,而由圖6可以看出UO2中O的分聲子態(tài)密度比UN中N的分聲子態(tài)密度作用范圍更廣,但絕對值有所降低,同時(shí)UO2的聲子項(xiàng)和光子項(xiàng)分隔不明顯,說明聲子態(tài)密度不僅與核素相關(guān)還與其晶體結(jié)構(gòu)密切相關(guān).

圖3 UN聲子色散關(guān)系Fig.3.Phonon dispersion of UN.

圖4 UN聲子態(tài)密度圖Fig.4.Phonon density of states in UN.

圖5 UN和UO2中U的分聲子態(tài)密度Fig.5.Phonon spectrum of U in UN and UO2.

圖6 UN中N和UO2中O的分聲子態(tài)密度Fig.6.Phonon spectrum of N in UN and O in UO2.

4.3 比熱容的比較

圖7 所示為本文利用VASP+PHONON軟件計(jì)算的UN比熱容隨溫度的變化,同時(shí)與文獻(xiàn)[20]中分子動力學(xué)(molecular dynamics,MD)的模擬結(jié)果以及文獻(xiàn)[21]中的實(shí)驗(yàn)結(jié)果進(jìn)行了比較.從圖7可以看出,本文模擬的Cv較文獻(xiàn)[20]中MD模擬的值更接近實(shí)驗(yàn)結(jié)果,同時(shí),高溫下三者都趨近于佩蒂特杜隆極限.通過比較可以看出,本文計(jì)算得到的UN聲子態(tài)密度是較為準(zhǔn)確的.

圖7 UN比熱容隨溫度的變化Fig.7.Heat capacity changes of UN with temperature.

4.4 熱中子散射律、Debye Waller積分和Teff

UN中N和U的S(α,β)隨β的變化情況如圖8所示,圖中給出了α=0.5,α=1.33和α=10的計(jì)算結(jié)果.從圖8(a)中可以看出,在α較小時(shí),S(α,β)隨β變化明顯,表現(xiàn)出顯著的振動,這種振動是聲子譜導(dǎo)致的,隨著α值的增大,β的變化范圍慢慢增大,同時(shí)振動減弱,這就要求在較大動量變化情況下進(jìn)行計(jì)算時(shí),β的取值范圍要盡可能大一些,這樣能更好地反映振子情況.從圖8(b)可以看出,UN中U的S(α,β)隨β變化的振蕩特性不明顯,這主要是因?yàn)閁原子質(zhì)量比較大,散射中的最大能量損失很小,這種情況下就沒有必要擴(kuò)展β到比較高的能量.

圖8 (a)UN中N的S(α,β)隨β的變化情況;(b)UN中U的S(α,β)隨β的變化情況Fig.8. (a)S(α,β)of N in UN changes with β;(b)S(α,β)of U in UN changes with β.

圖9 (a)293.6 K溫度下UO2中16O和UN中15N(14N)的非彈性散射和彈性散射截面;(b)293.6 K溫度下UO2中U和UN中U的非彈性散射和彈性散射截面Fig.9.(a)Inelastic and elastic cross sections of16O in UO2and15N(14N)in UN at 293.6 K;(b)inelastic and elastic cross sections of U in UO2and UN at 293.6 K.

由(8)式—(10)式可以得到,較高的入射能量超出了S(α,β)范圍,引入SCT后,熱中子的非彈性散射由Teff決定,由 (11)式—(16)式可以得到,非相干彈性散射截面由Debye Waller積分γ(0)決定.表3給出為UN中238U和14N(15N)的Debye Waller積分和Teff數(shù)值,從表3可以看出有效溫度要略高于實(shí)際溫度.

表3 Debye Waller積分和TeffTable 3.Debye Waller integral and Teffparameters.

4.5 熱中子截面

為了比較UO2和UN熱中子截面的差異,分別對比了293.6 K溫度下UO2中16O和UN中14N,15N以及UO2中U和UN中U的非彈性散射截面和彈性散射截面,如圖9所示.同時(shí)給出了不同溫度下UN中14N以及UN中U非彈性散射截面和彈性散射截面,如圖10所示.由圖9和圖10可以看出,UO2中16O考慮了相干彈性散射部分,UN中14N和15N忽略了彈性散射的相干部分,并且14N的非彈性散射和彈性散射截面略高于15N,14N的非彈性散射截面與16O較為接近,相同溫度下UN中238U非彈性散射和彈性散射截面比UO2中238U偏小.隨著溫度的升高,UN中N和U的非彈性散射截面升高,中子與UN作用更激發(fā)晶格態(tài),從而獲得聲子;相反,彈性散射截面隨溫度的升高是降低的,當(dāng)能量為1 eV左右時(shí)彈性散射截面基本為零,總截面主要是非彈性散射截面的貢獻(xiàn),等于自由核散射截面,與自由氣體模型一致.

圖10 非彈性散射和彈性散射截面隨溫度的變化 (a)UN中14N;(b)UN中UFig.10.Inelastic and elastic cross sections changes with temperatures:(a)N in UN;(b)U in UN.

5 結(jié)論與展望

本文采用第一性原理,基于熱中子散射理論分析制作了氮化鈾的熱中子截面,提出合適的UN熱中子截面的預(yù)測方法.從UN晶格參數(shù)出發(fā),利用第一性原理晶格動力學(xué)的直接方法得到了UN聲子態(tài)密度,補(bǔ)充了熱能區(qū)多溫度點(diǎn)的UN熱中子散射截面,填補(bǔ)UN熱能區(qū)數(shù)據(jù)的缺失;分別將影響較大的中子熱散射律、熱中子彈性和非彈性散射截面與傳統(tǒng)UO2進(jìn)行對比,相同溫度下,UN中238U非彈性散射和彈性散射截面比UO2小,UN中N忽略了相干散射部分;隨著溫度升高,UN彈性散射截面變小,非彈性散射變大;彈性散射截面在高能段等于零,非彈性散射截面在低能區(qū)主要通過聲子的吸收獲得能量,截面的變化符合1/v規(guī)律,在中能區(qū),通過與UN核碰撞產(chǎn)生或發(fā)射聲子,隨著能量的升高截面增加,在高能區(qū)符合自由核模型.本文的研究結(jié)論揭示了UN在熱中子反應(yīng)堆中的熱中子散射特性,為以UN為燃料的熱中子反應(yīng)堆的研究奠定了基礎(chǔ).