趙中華 渠廣昊 姚佳池 閔道敏翟鵬飛 劉杰 李盛濤?

1) (西安交通大學電氣工程學院, 電力設(shè)備電氣絕緣國家重點實驗室, 西安 710049)

2) (中國科學院近代物理研究所, 蘭州 730000)

ZrO2陶瓷耐高溫、耐腐蝕、抗輻照性能強, 是極具前景的反應堆惰性基質(zhì)燃料和錒系元素固化材料.本文聯(lián)合使用熱峰模型和分子動力學方法, 模擬了核輻射環(huán)境下ZrO2的相變過程: 基于熱峰模型, 從快速重離子注入后能量沉積和傳導的多物理過程出發(fā), 建立熱擴散方程, 求得ZrO2晶格溫度時空演變特性; 然后運用分子動力學方法模擬了該熱峰作用下, 單斜ZrO2相變的微觀物理過程.研究發(fā)現(xiàn), 電子能損為30 keV·nm–1的單一快速重離子注入后, ZrO2中心產(chǎn)生一個半徑為7 nm的柱形徑跡, 徑跡中心晶格迅速熔融, Zr原子配位數(shù)由7降至4—6, 2 ps時開始結(jié)晶并形成空洞, 空洞周圍為非晶區(qū), 非晶區(qū)外Zr原子配位數(shù)變?yōu)?, 同時X射線衍射(X-ray diffraction, XRD)計算和分析結(jié)果確認發(fā)生了單斜相向四方相的轉(zhuǎn)變.隨著熱峰能量向周圍傳遞, 相變區(qū)逐漸擴大.經(jīng)熱峰計算和分子動力學模擬, 輻照誘導ZrO2由單斜相轉(zhuǎn)為四方相的快速重離子的電子能損閾值為21 keV·nm–1.

1 引 言

隨著核能的不斷開發(fā)利用, 高放射性核廢料持續(xù)增多, 如何有效地提高核燃料利用率并減少放射性核廢料, 已成為核能可持續(xù)發(fā)展所必須解決的關(guān)鍵問題.近年來研究發(fā)現(xiàn), 在反應堆以及加速器驅(qū)動的次臨界系統(tǒng)中, 采用惰性基質(zhì)燃料能夠有效地促進钚的回收利用和次錒系元素的嬗變, 從源頭上降低核廢物中長壽命放射性元素的含量[1,2].氧化鋯(ZrO2)陶瓷耐高溫、耐腐蝕、抗輻照性能強, 是惰性基質(zhì)燃料的候選材料之一[3,4].然而, ZrO2在不同溫度下呈現(xiàn)三種不同的晶相, 1170 ℃以下為單斜相, 1170—2370 ℃之間為四方相, 2370—2715 ℃之間為立方相.ZrO2在加熱冷卻過程中,可發(fā)生單斜至四方相變(1100—1250 ℃), 以及四方至單斜相變(1050—700 ℃), 后者的剪切應變?yōu)?6%—18%, 導致體積增加約5%[5].核反應堆的強輻射環(huán)境可誘導ZrO2發(fā)生相變, 這將對ZrO2在核工業(yè)方面的應用產(chǎn)生極大影響, 例如惰性基質(zhì)燃料相變引起的體積和應力變化將可能損壞燃料棒結(jié)構(gòu), 損傷核反應裝置, 帶來核泄漏風險[1,2,6].因此, 研究ZrO2在高能粒子輻照下的相變行為具有重要工程意義.

近20年來, 利用加速器技術(shù)對ZrO2在不同種類和能量的快速重離子作用下的抗輻照性能進行了廣泛研究[7?12].研究發(fā)現(xiàn), 單斜ZrO2晶體在快速重離子輻照下可以轉(zhuǎn)化為四方相, 并且能夠在室溫下穩(wěn)定存在.當前主要采用熱峰模型對快速重離子誘導材料相變現(xiàn)象進行機理解釋[13?15], 即快速重離子輻照促進材料局部溫度升高, 超過材料的相變溫度, 從而誘導相變.但熱峰模型只能定量計算輻照后材料內(nèi)部的溫度分布, 不能直觀地反映材料微觀結(jié)構(gòu)變化[16].分子動力學能夠模擬幾百皮秒和幾十納米時空范圍內(nèi)材料微觀結(jié)構(gòu)演化, 已廣泛應用于材料的低能輻照損傷研究中[17,18].然而快速重離子輻照材料時, 能量損失過程由電子非彈性碰撞主導.分子動力學無法描述電子與電子以及電子與聲子之間的相互作用, 因此不能模擬高能粒子對材料的輻照損傷作用.

本文結(jié)合熱峰模型和分子動力學兩種方法的優(yōu)勢, 對快速重離子輻照下, ZrO2由單斜相轉(zhuǎn)化為四方相的微觀過程進行模擬, 并計算誘導ZrO2相變的快速重離子的電子能損閾值.首先, 利用熱峰模型描述快速重離子入射后, ZrO2內(nèi)部電子激發(fā)、電子-聲子能量耦合、電子-電子以及聲子-聲子能量傳遞過程, 計算得出晶格溫度; 然后, 根據(jù)晶格溫度標定體系原子動能, 將其作為分子動力學模擬初始條件, 研究ZrO2微觀結(jié)構(gòu)演化, 并通過X射線衍射(X-ray diffraction, XRD)計算確定ZrO2相變類型; 隨后, 根據(jù)Zr原子配位數(shù)變化具體分析ZrO2相變過程; 最后, 改變?nèi)肷潆x子的電子能損,進行多次熱峰計算和分子動力學模擬, 確定誘導ZrO2發(fā)生相變的快速重離子的電子能損閾值.

2 計算模型

2.1 熱峰模型

熱峰模型認為快速重離子穿過靶材料時, 首先以離子運動路徑為軸線形成圓柱形的電離核心區(qū)域, 稱為徑跡.在大約10–16至10–15s時間內(nèi)將能量傳遞給靶材電子系統(tǒng); 隨后, 徑跡中心電子通過熱傳導將部分能量傳遞給徑跡周圍電子; 同時, 通過電子-聲子耦合作用將部分能量傳遞給原子; 徑跡中心原子獲得能量后, 通過晶格振動將能量傳遞給徑跡周圍原子.熱峰模型的熱過程可以通過以下兩個耦合微分方程描述[19?21], 分別對應電子和原子系統(tǒng).

式中,Te和Ta,Ce和Ca,Ke和Ka分別為電子和原子的溫度、比熱容和熱導率,Ce取1 J·cm–3·K–1,Ke取2 J·cm–1·s–1·K–1[22],Ca[23,24]和Ka[24?26]為 隨溫度變化的量, 其表達式分別為

ρ為晶體質(zhì)量密度, 取5.68 g·cm–3;g為電子-聲子耦合系數(shù), 由電子平均自由程λ決定,g=Ke/λ,對于ZrO2,λ= 4.5 nm[14];A(r,t)為入射離子在半徑r, 時間t時沉積在電子系統(tǒng)的能量密度, 其表達式為

式中,D(r)為初始沉積能量空間分布, 指數(shù)因子α= 1/τ,τ為電子能量沉積時間, 取10–15s, 可調(diào)參數(shù)A表達式為

式中,rm為電離徑跡半徑,Se為入射離子的電子能損值, 對于給定能量和種類的入射離子, 其電子能損值可由SRIM軟件包[27]求得.

數(shù)值求解熱峰方程, 獲取ZrO2原子溫度的時空分布.

2.2 分子動力學模擬

建立如圖1所示的單斜ZrO2原胞模型, 將其擴展為40a0× 10b0× 40c0的超晶胞模型, 模型尺寸為20.5816 nm × 5.2075 nm × 21.2428 nm, 共包 含1.92 × 106個 原 子.其 中a0= 5.1454 ?,b0= 5.2075 ?,c0= 5.3107 ?, 為單斜ZrO2的晶格常數(shù).

圖1 單斜ZrO2原胞結(jié)構(gòu)模型Fig.1.The crystal structure of monoclinic ZrO2.

計算采用LAMMPS軟件[28], ZrO2原子間相互作用力采用Buckingham勢函數(shù)[29,30]描述, 勢能表達式為

式中, 第一項為長程庫倫勢, 后兩項為短程Buckingham勢;qi和qj分別為原子i和j的電荷,rij為原子間距,A,ρ,C為可調(diào)參數(shù).在快速重離子輻照下, 熱峰中心區(qū)域原子處于非平衡態(tài), 原子間距小于平衡距離, 相互作用力為較強的斥力.Buckingham勢能夠準確模擬平衡態(tài)下ZrO2原子間相互作用, 然而, 由其描述的原子間近核作用力為相互吸引力, 與實際不符.因此, 采用普適的Ziegler-Biersack-Littmark (ZBL)庫倫屏蔽勢[27]對Buckingham勢的近核力進行修正, 并用Fermi函數(shù)將兩種勢函數(shù)平滑連接[31]:

式中,rij為原子間距;rc,AF為可調(diào)參數(shù),rc決定兩種勢函數(shù)間連接點的選取,AF控制連接處的平滑程度.

模擬過程中采用三維周期性邊界條件, 截斷半徑取1.6 nm.采用自適應時間步長[32], 即根據(jù)原子速度和受力情況, 實時更新積分步長, 限制每個原子在單個步長內(nèi)最大移動距離為0.5 ?, 最大動能變化為100 eV, 并設(shè)置最大允許步長為1 fs.

熱峰注入前, 采用NVT系綜進行10 ps弛豫,使用Nose-Hoover恒溫器進行300 K熱浴.隨后,如圖2所示, 將熱峰沿[010]方向注入, 根據(jù)熱峰模型計算得出的晶格溫度標定原子動能.采用NVE系綜對弛豫階段進行模擬, 并將X和Z方向3 ?厚邊界層與300 K Berendsen恒溫器連接, 模擬能量通過熱傳導耗散到環(huán)境的過程.

圖2 熱峰注入靶材料示意圖(Se = 30 keV·nm–1)Fig.2.Diagram of thermal spike injection into target material (Se = 30 keV·nm–1).

3 結(jié)果與討論

3.1 熱峰作用下的晶格溫度時空演化

圖3展示了不同位置和不同時刻的晶格溫度變化趨勢, 入射離子的電子能損為30 keV·nm–1,為實驗值(10—50 keV·nm–1)的中位數(shù).

圖3 晶格溫度分布曲線 (a) 不同位置晶格溫度隨時間演化; (b) 不同時刻晶格溫度的徑向分布Fig.3.Distributions of lattice temperature: (a) Evolution of lattice temperature at different radial distances; (b) radial distribution of lattice temperature at different times.

由圖3可知, 在電子能損為30 keV·nm–1的快速重離子輻照下, ZrO2徑跡周圍晶格溫度迅速升高, 60 fs時各處溫度基本已達最大值, 徑跡中心溫度最高接近12000 K, 遠超ZrO2熔點(2715 ℃).徑跡中心溫度達到峰值后迅速降低, 距徑跡中心較遠的晶格溫度達到峰值后, 先維持一段時間再緩慢下降, 徑跡內(nèi)各處溫差隨時間減小, 60 ps時各處溫度已接近環(huán)境溫度300 K.熱峰作用下, 距徑跡中心4 nm內(nèi)晶格溫度會超過熔點(2715 ℃),6 nm內(nèi)晶格溫度會超過單斜相至四方相的相變溫度(1170 ℃).

3.2 熱峰作用下的晶格溫度時空演化

由熱峰模型得出的晶格溫度徑向分布, 反映了快速重離子輻照后, 沉積在不同位置處的原子動能.根據(jù)計算結(jié)果, 60 fs時徑跡內(nèi)各處溫度基本已達最大值, 故取此時刻晶格溫度, 標定不同位置處原子動能.經(jīng)分子動力學模擬, 得出不同時刻ZrO2微觀結(jié)構(gòu)變化如圖4所示.

圖4 不同時刻ZrO2原子結(jié)構(gòu)變化.所有原子沿熱峰注入方向投影在(010)平面上Fig.4.Atomic structure of ZrO2 at different times.All the atoms are projected on (010) plane along the injection direction of thermal spike.

由圖4可知, 0 ps時體系原子處于熱平衡態(tài).熱峰注入后, 中心區(qū)域原子能量較高, 受熱激發(fā)離開初始晶格位置, 并與周圍原子發(fā)生劇烈碰撞.0.15 ps時, 產(chǎn)生一個熔融的柱形徑跡, 距離熱峰中心3 nm內(nèi)原子無序程度較高, 3至7 nm內(nèi)原子無序程度較低.2 ps時徑跡中心開始形成空洞, 空洞外圍為高度無序的非晶區(qū), 非晶區(qū)周邊原子熔融重排, 形成了兩種不同的晶區(qū)(晶區(qū)一和晶區(qū)二).在(010)投影平面上, 晶區(qū)二內(nèi)Zr, O原子空間分布均勻, 且Zr原子位于臨近4個O原子的對稱中心, O原子位于臨近4個Zr原子的對稱中心; 晶區(qū)一內(nèi)Zr與Zr原子間距較遠, O與O原子間距較近, 且O原子不處于臨近Zr原子的對稱中心位置.由4, 6和100 ps的微觀結(jié)構(gòu)可知, 隨著能量由徑跡中心向四周傳遞, 重結(jié)晶過程逐漸發(fā)展至整個體系, 空洞、晶區(qū)一和晶區(qū)二的尺寸逐漸增大, 100 ps時體系已達平衡態(tài).

XRD圖譜能夠反映材料成分、形態(tài)以及內(nèi)部原子排列結(jié)構(gòu)等信息, 廣泛應用于物質(zhì)的晶體結(jié)構(gòu)表征.取100 ps時晶區(qū)二內(nèi)部分原子坐標, 進行分子動力學計算, 獲得輻照后ZrO2晶體XRD圖譜[33],并與輻照前XRD圖譜以及單斜(m-ZrO2)、四方(t-ZrO2)標準譜進行對比, 結(jié)果如圖5所示.

圖5 輻照前后以及標準ZrO2 XRD圖譜Fig.5.XRD patterns of pristine, ion-irradiated and standard ZrO2.

由圖5可知, 輻照前XRD圖譜在28.2°, 31.5°,34.2°, 35.3°, 49.3°和50.1°出現(xiàn)衍射峰, 分別對應m-ZrO2的( 111 ), (111), (002), (200), (022), (220)晶面衍射(JCPDS.No 83-0936), 表明輻照前ZrO2為單斜相.經(jīng)快速重離子輻照后, 單斜相特有的特征峰消失, 并且在30.3°, 34.6°, 50.3°, 50.8°, 59.4°,60.3°出現(xiàn)新的特征峰, 與t-ZrO2的(101), (002),(112), (200), (103), (211)衍射晶面良好對應(JCPDS.No 88-1007), 證明在快速重離子輻照作用下, ZrO2由單斜相轉(zhuǎn)為了四方相.

ZrO2晶體中, 單斜相的Zr原子配位數(shù)為7,四方相的Zr原子配位數(shù)為8, 因此可根據(jù)Zr原子配位數(shù)的變化情況研究輻照條件下ZrO2相變過程.圖6展示了不同配位數(shù)的Zr原子空間分布隨時間演化過程.

圖6 不同時刻不同配位數(shù)Zr原子空間分布.所有原子沿離子入射方向投影在(010)平面上, 顏色不同代表配位數(shù)不同F(xiàn)ig.6.Spatial distribution of Zr atoms with different coordination number at various times.All the atoms are projected on (010) plane along the injection direction of thermal spike.Different coordination numbers are distinguished by the atomic color.

如圖6所示, 初始時刻所有Zr原子配位數(shù)均為7, 表明ZrO2為單斜相.0.15 ps時, 距熱峰中心7 nm范圍內(nèi)形成一個熔融的柱形徑跡, 其中28%的Zr原子配位數(shù)降為6, 4.4%的Zr原子配位數(shù)降為5, 距徑跡中心3 nm范圍內(nèi), 47.8%的Zr原子配位數(shù)降為6, 14.8%的Zr原子配位數(shù)降為5, 0.7%的Zr原子配位數(shù)降為4, 表明越靠近徑跡中心, 原子無序程度越高.隨著能量向四周傳遞,徑跡中心溫度逐漸降低, 原子發(fā)生結(jié)晶重排, 2 ps時, 徑跡中心出現(xiàn)空洞, 其周邊由配位數(shù)為3或4的Zr原子組成非晶區(qū), 非晶區(qū)兩側(cè)生長出兩種不同的晶相, 分別為配位數(shù)為8的四方相以及配位數(shù)為6的低配位相.由4, 6和100 ps微觀結(jié)構(gòu)可知, 隨著能量繼續(xù)由中心向四周傳遞, 四方相以及低配位相逐漸發(fā)展至整個體系, 同時空洞范圍不斷擴大.4和6 ps時, 四方相區(qū)域內(nèi)仍有部分Zr原子配位數(shù)為7, 100 ps時, 四方相區(qū)域內(nèi)所有Zr原子配位數(shù)均為8, 表明隨著時間推移, 相變區(qū)域變得更加致密有序.100 ps時, 空洞周邊部分低配位相逐漸轉(zhuǎn)變?yōu)樗姆较? 將兩側(cè)相變區(qū)連接起來.本文模擬所得到的空洞, 在實驗中尚未得到驗證.然而, 單斜ZrO2密度為5.65 g/cm3, 四方ZrO2密度為6.10 g/cm3, 一定質(zhì)量的ZrO2由單斜相轉(zhuǎn)為四方相后, 體積減小約7.4%.本文分子動力學模擬中采用的系綜為NVE系綜, 模擬過程中體系體積固定, 因此, 在NVE系綜下進行ZrO2由單斜相至四方相的相變模擬, 必然導致空洞的產(chǎn)生.此外, 在其他材料(如Ge[31], Fe-Cr[34]合金等)的輻照損傷研究中, 曾觀察到空洞的產(chǎn)生.

3.3 輻照誘導單斜ZrO2相變的電子能損閾值

為獲得輻照誘導ZrO2由單斜相轉(zhuǎn)為四方相的快速重離子的電子能損閾值, 在電子能損為15—30 keV·nm–1范圍內(nèi)進行了多次熱峰計算和分子動力學模擬.圖7展示了Zr原子配位數(shù)隨電子能損值變化趨勢.

如圖7所示, 電子能損值為15—20 keV·nm–1時, 配位數(shù)為7的Zr原子個數(shù)隨電子能損值增加而減少, 未出現(xiàn)配位數(shù)為8的Zr原子, 表明ZrO2仍保持單斜相.當電子能損值為21 keV·nm–1時,72.4% Zr原子配位數(shù)由7變?yōu)?, 表明ZrO2已經(jīng)由單斜相轉(zhuǎn)為四方相.由此得出輻照誘導單斜ZrO2相變的快速重離子的電子能損閾值為21 keV·nm–1.電子能損值為21—30 keV·nm-1時, 隨著電子能損值增加, 配位數(shù)為8的Zr原子個數(shù)先增加后減小,占Zr原子總數(shù)比例在71.5%—76%之間.

圖7 配位數(shù)隨電子能損值變化Fig.7.The change of coordination number with electronic energy loss.

當入射離子電子能損分別為20和21 keV·nm–1時, 輻照后不同配位數(shù)的Zr原子空間分布如圖8所示.

圖8 不同配位數(shù)Zr原子空間分布 (a) 電子能損值為20 keV·nm–1; (b)電子能損值為21 keV·nm–1.所有原子沿離子入射方向投影在(010)平面上, 顏色不同代表配位數(shù)不同F(xiàn)ig.8.Spatial distribution of Zr atoms with different coordination number: (a) Se = 20 keV nm–1; (b) Se = 21 keV nm–1.All the atoms are projected on (010) plane along the injection direction of thermal spike.Different coordination numbers are distinguished by the atomic color.

如圖8所示, 電子能損值為20 keV·nm–1時,輻照未能造成ZrO2由單斜相至四方相的轉(zhuǎn)變.在熱峰作用下, 部分原子離開其正常晶格位置, 成為點缺陷, 點缺陷的引入使得部分Zr原子配位數(shù)降為6或5, 點缺陷在空間上均勻分布.電子能損值低于20 keV·nm–1時, 不同配位數(shù)Zr原子的空間分布圖與圖8(a)相似.當快速重離子電子能損值為21 keV·nm–1時, 72.4% Zr原子配位數(shù)由7變?yōu)?, 并且在空間上形成連續(xù)致密的排列結(jié)構(gòu), 表明ZrO2已經(jīng)由單斜相轉(zhuǎn)為四方相.當電子能損值低于27 keV·nm–1時, 低配位相和四方相中心未形成連續(xù)空洞, 與圖8(b)相似, 這是由于熱峰溫度較低, 中心原子動能較小, 不足以引起足夠劇烈的碰撞使得熱峰中心所有原子運動到足夠遠的區(qū)域, 隨著系統(tǒng)快速淬火, 原子迅速冷卻, 形成非連續(xù)的空洞.當電子能損值高于27 keV·nm–1時, 低配位相和四方相中心形成連續(xù)空洞, 不同配位數(shù)Zr原子空間分布與圖6(f)相似.

根據(jù)熱峰模型[13,35], 快速重離子輻照后, 材料相變區(qū)域面積σ與入射離子的電子能損Se存在如下關(guān)系:

式中,a0為熱峰溫度的初始分布半徑;Set為入射離子的電子能損閾值.

經(jīng)熱峰模型計算和分子動力學模擬得出, 輻照致單斜ZrO2相變的快速重離子的電子能損閾值為21 keV·nm–1, 與實驗值12 keV·nm–1存在一定誤差[7?10], 這主要與輻照劑量以及勢函數(shù)的準確性有關(guān).模擬條件為單離子入射, 對應的輻照劑量為2.37 × 1011ions·cm–2, 而實驗中通常輻照劑量較大, 為1012—1015ions·cm–2, 可能存在同一區(qū)域被多個離子輻照的情況.此外, 分子動力學模擬中采用普適的Ziegler-Biersack-Littmark (ZBL)庫倫屏蔽勢修正的Buckingham勢描述ZrO2原子間相互作用力, 其與真實的原子間作用力相比存在一定差異, 這將導致模擬結(jié)果與實驗結(jié)果存在一定偏差.未來可開發(fā)更為準確的勢函數(shù), 開展多個離子同時入射以及多個離子先后入射的研究, 模擬高劑量快速重離子輻照對材料的損傷作用.

4 結(jié) 論

本文結(jié)合了熱峰模型和分子動力學的優(yōu)勢, 模擬了快速重離子輻照下, ZrO2由單斜相轉(zhuǎn)為四方相的微觀物理過程, 并計算了輻照誘導ZrO2相變的快速重離子的電子能損閾值.首先, 基于熱峰模型, 考慮快速重離子輻照單斜ZrO2時, 電子激發(fā)、電子-聲子能量耦合以及電子-電子、聲子-聲子能量傳導多物理過程, 建立熱擴散方程, 求得ZrO2晶格溫度時空演變特性.然后, 運用分子動力學方法模擬了熱峰作用下ZrO2相變的微觀物理過程.熱峰注入后, ZrO2中心形成一個半徑為7 nm的徑跡, 徑跡中心晶格迅速熔融; 2 ps時熔融區(qū)逐漸冷卻重結(jié)晶, 其中心開始形成空洞, 空洞外圍為高度無序的非晶區(qū), 非晶區(qū)外圍生成四方相和低配位相; 隨著能量由徑跡中心向四周傳遞, 四方相和低配位相逐漸發(fā)展至整個體系, 同時空洞尺寸隨之增大.經(jīng)熱峰計算和分子動力學模擬, 輻照誘導單斜ZrO2相變的快速重離子的電子能損閾值為21 keV·nm–1, 高于實驗值12 keV·nm-1, 主要原因為: 1) 入射離子劑量的模擬值小于實驗值; 2) 分子動力學中采用勢函數(shù)描述的原子間相互作用力與真實情況存在一定誤差, 未來可開發(fā)更精確的勢函數(shù)描述ZrO2原子間相互作用力, 開展高劑量快速重離子輻照下, 材料輻照損傷的模擬研究.本文建立的研究方法可有效地模擬高能離子對材料輻照損傷的微觀物理過程, 將來可用于其他惰性基質(zhì)燃料(如YSZ/ZrO2-PuO2等)的高能輻照損傷研究.