U1-xThxO2混合燃料力學性能的分子動力學模擬*

辛勇 包宏偉 孫志鵬 張吉斌 劉仕超 郭子萱王浩煜 馬飛? 李垣明?

1) (中國核動力研究設計院, 核反應堆系統設計技術重點實驗室, 成都 610213)

2) (西安交通大學材料科學與工程學院, 金屬材料強度國家重點實驗室, 西安 710049)

1 引 言

二氧化鈾(UO2)陶瓷燃料是目前商用反應堆的主要核燃料, 包括壓水堆、輕水堆及重水堆等,具有熔點高、抗輻照性能好、與包殼和冷卻劑相容性好等優點[1,2].然而, 隨著反應堆向長周期、高燃耗的方向發展, 工況條件愈加苛刻, UO2核燃料在高溫、氣態和固態核裂變產物及高能粒子沖擊等多種因素的作用下, 將導致燃料內部熱應力升高、晶格畸變甚至相變、腐蝕等作用, 嚴重限制燃料燃耗及安全性[3,4].為提高核動力效率以及反應堆安全性能, 尋求高熱力學穩定性、抗輻照腫脹性能優異的核燃料至關重要.

目前各國常采用的改進方法為通過摻雜形成混合氧化物燃料(MOX)[5].釷(Th)元素, 在地殼中的含量是U的3至4倍, 不易裂變, 與UO2無限互溶形成U1—xThxO2固溶體混合燃料, 具有更高的熔點、熱導率和抗氧化性能[6-8].目前至少5%核燃料為U1—xThxO2MOX, 廣泛存在于核燃料循環的各個階段.核燃料循環過程中, 高的熱應力導致燃料芯塊變形甚至斷裂, 造成裂變氣體釋放, 使熱導率降低, 大大提高了芯塊包殼界面的失效概率,嚴重影響反應堆安全運行[9,10].然而, 受限于實驗條件, 目前關于U1—xThxO2混合燃料熱力學行為的研究仍具有較大挑戰.利用分子動力學(molecular dynamics, MD)模擬, 可從原子尺度準確描述基于UO2燃料的結構演變與熱力學特性.Xiao等[11,12]研究了U1—xThxO2混合燃料的晶格常數、熱膨脹系數、熱導率、熱容及焓在300—2100 K溫度范圍內的變化規律.此外, 借助MD模擬可研究UO2基混合燃料中晶界遷移、孔洞演化、氣體擴散及級聯碰撞過程[13-16].近年來, Cooper等[17,18]成功開發了一種混合多體勢, 精確描述了UO2基混合燃料的熱力學特性, 例如, 基于該勢函數, Balboa等[19]計算了純UO2的晶格常數隨溫度變化規律,與實驗結果高度吻合; Rahman等[20,21]成功計算了U0.5Th0.5O2體系的熱導率、剪切模量、楊氏模量及壓縮系數等物理參量; Ghosh等[22]則利用固液兩相法成功計算了(U,Th)O2體系的熔點, 發現UO2和ThO2的熔點分別在3650—3675 K 和3050—3075 K 之間; Palomares等[23]則模擬了600 ℃下UO2在氧富集條件下的結構穩定性, 計算得到的徑向分布函數與實驗測定的結果高度吻合.

實際反應堆條件下, UO2基燃料芯塊缺乏宏觀塑性, 承受較高溫度和較大溫度梯度下, 裂紋快速形核擴展, 發生脆性斷裂來釋放應力, 此外, 斷裂的芯塊碎片與包殼材料接觸, 使裂變氣體釋放,碘、鎘等裂變產物與包殼材料相互作用, 使其脆化,進而造成嚴重事故[24].因此, UO2基燃料芯塊的斷裂行為在過去半個多世紀以來受到廣泛關注.UO2基燃料芯塊具有復雜斷裂行為, 對溫度、晶體組織、加載方式等均很敏感[25], 但大多數研究多著眼對其斷裂過程的觀察表征.相關實驗與理論研究發現純UO2單晶的韌脆轉變溫度在1300—1900 K之間, 而輻照條件下裂變氣體的擴散通常在溫度低于1000 K時即可發生[26].Mo等[27]利用放電等離子體燒結(spark plasma sintering, SPS)法制備了UO2塊體, 利用先進同步微X衍射(micro-region X-ray diffraction, mXRD)技術測量了裂紋尖端區域的高晶格應力.Desai等[28]、Tian等[29]和Zhang等[30]利用MD模擬觀察到了多晶UO2中裂紋的形成與快速擴展, 最終導致沿晶脆斷.盡管圍繞U1—xThxO2混合燃料的基本熱力學特性有了一些MD 模擬結果, 但對UO2基燃料在一定溫度范圍內的變形斷裂行為的研究還比較欠缺, 裂紋的形核與擴展機理仍不清楚, 尤其是對其塑性仍有一定爭議[31].特別的是, 高壓試驗、密度泛函理論計算和MD模擬均證明, UO2在高壓或應力誘導下發生由立方結構的螢石相(fluorite phase)至斜方結構的氯鉛礦相(cotunnite phase)的相變[32,33].Zhang等[30]、Fossati等[34]及Balboa等[19]甚至在單晶UO2晶體中的裂紋附近觀察到了更為特殊的金紅石相(rutile phase)或scrutinyite相.上述相變勢必影響UO2基燃料的斷裂特性.本文采用MD模擬, 系統地研究了U1—xThxO2混合燃料的力學行為, 預測可能發生的相變, 考察溫度、摻雜濃度等因素對彈性模量、斷裂強度、斷裂應變等力學性能的影響.

2 計算方法與模型

MD模擬中勢函數的選取至關重要, 決定了計算結果的精確與否.選取Cooper等[17,18]開發的混合多體勢來描述U, Th與O之間的關系, 該混合勢耦合了多體的原子嵌入勢(EAM)及對勢作用,表達式為

其中第一項為原子i與j之間的對勢作用, 第二項為EAM勢.對勢表達式為

式中φC(rij) 為長程庫侖力, 通過基于多粒子網格(particle-particle particle-mesh, PPPM)的埃瓦爾得法計算,

φB(rij) 與φM(rij) 分別為短程的Buckingham勢與Morse勢, 其表達式為

其中Aαβ,Cαβ和Dαβ為經驗參數.該混合勢已廣泛應用于U1—xThxO2混合燃料體系熱力學特性的計算, 精確描述了晶格常數、熱導率、熱容、焓等物理量的變化規律.Cooper等[17]給出了通用于MD模擬的勢函數文件, 為保證計算精度, 本文所設定的截斷半徑為11.0, PPPM參數為10—5.

室溫下UO2與ThO2均為立方螢石結構, 晶格常數非常接近, 分別為5.4704 ?和5.5975 ?, 其中U和Th離子占據面心立方格點位置[35].采用三維周期邊界條件模擬塊體UO2混合燃料的力學特性.鑒于長程庫侖需耗費大量計算資源, 首先對比了兩種尺寸模型(7.75 nm × 8.05 nm × 7.65 nm,15.5 nm × 16.0 nm × 15.1 nm)的差異性, 發現周期邊界條件下尺寸并未影響UO2基混合燃料的基本力學特性, 因此, 選取較小尺寸模型進行研究.如圖1所示為7.75 nm × 8.05 nm × 7.65 nm 的U1—xThxO2混合燃料模型, 考慮兩種晶體取向的模型,x,y和z坐標軸方向分別為和[111]或[100], [010]和[001].總原子數為34560, 其中陽離子(U, Th)總數為11520, O離子總數為23040,Th離子以隨機替位摻雜形式取代U離子, 其濃度在0.05至0.55范圍內變化.在NPT系綜下利用Nosé-Hoover熱浴法對構建好的初始模型進行結構弛豫, 溫度設定為300 K, 壓強為0 Pa, 弛豫時間為100 ps, 時間步長為1 fs.然后每隔100 K將溫度逐漸升至1200 K, 每次升溫后弛豫100 ps.然后在給定溫度下進行沿z軸方向的單軸拉伸加載, 如圖1所示, 固定上下兩端5原子層為加載層, 施加單軸拉伸加載, 應變速率保持在10—9s—1.利用維里應力準則計算體系中每個原子的應力, 表達式為[36]

圖1 U1—xThxO2 混合燃料模型及單軸拉伸加載示意圖Fig.1.Atomic model of U1—xThxO2 upon uniaxial tensile loading.

其中V為體積,N為原子總數,為α原子的第i個分速度,為α原子的質量,γαβ為α和β原子之間的距離,U′為勢能函數.加載過程中整個系統的應力由以下公式給出:

式中,S為原子應力.應變由以下公式給出:

式中,L(Nsteps) 為體系瞬時長度,L(0) 為初始長度.彈性模量代表著使原子離開平衡位置的難易程度,是表征原子間結合力強弱的物理量.根據(7)式和(8)式, 繪制應力-應變曲線, 在線彈性應變范圍內進行線性擬合得到彈性模量與泊松比隨溫度的變化規律, 計算公式為

式中 ΔU為勢能變化量,ε為沿某一方向的應變.繼續加載直至體系失效, 得到斷裂強度與溫度的關系.所有模擬工作均使用開源軟件LAMMPS[37]完成, 原子結構通過開源軟件OVITO進行可視化[38].

3 結果與討論

3.1 力致相變

以UO2與U0.75Th0.25O2為例, 圖2(a)給出了室溫300 K時沿[001]和[111]晶向拉伸時的應力-應變曲線.若z軸沿[001]晶向, 當應變增加至0.08附近, 應力發生突降, 當應變增至0.13左右時, 應力重新開始線性增加, 直至斷裂, 而若z軸沿[111]晶向, 應力隨應變線性增加至最高點后急劇下降,最終斷裂.圖2(b)和圖2(c)分別為沿[001]和[111]拉伸過程中不同應變下的徑向分布函數.z軸沿[001]晶向時, 無論是對UO2還是U0.75Th0.25O2,應變超過0.08后, O—O鍵特征峰逐漸消失, 意味著兩個O原子之間距離變遠, U—O和Th—O鍵特征峰藍移, 即U與O和Th與O之間距離縮短, 暗示相變的發生.而沿[111]晶向時, U—O或Th—O鍵及O—O鍵特征峰, 均只隨應變的增加而紅移, 證明在拉伸加載時鍵長伸長.圖2(d)為U0.75Th0.25O2沿[001]晶向拉伸時的原子結構演化, 插圖為原子結構放大圖, 與應力-應變曲線首次下降區域對應, 觀察到由初始的螢石結構至scrutinyite結構的相變, 該結構對稱性較低, 與Fossati等[34]的MD 模擬結果相吻合, 相變極大地增強了UO2混合燃料的塑性.如圖2(e)所示, 沿[111]晶向拉伸時原子結構僅發生變形, 未觀察到相變.雖然實驗上已經證明UO2具有多種相, 但其僅能在高壓等極端條件下存在, 目前還未確切觀測到單軸應力即可誘導的相變.這種局部高應力誘發的固體相變, 類似于馬氏體相變, 僅僅原子相對位置發生改變, 由于Th離子以替位摻雜的形式存在,并未改變U1—xThxO2混合燃料的物理化學環境, 因此相變特征對摻雜濃度不敏感.溫度是影響相變的另一重要因素, 隨溫度的升高, 熱振動勢必會改變相變勢壘.鑒于實驗中難以觀察到這些特殊相變,且實際中U1—xThxO2混合燃料大多發生脆性斷裂.因此, 本文僅研究了室溫下的相變行為, 未來還需進一步研究相變的溫度依賴性, 為實驗設計提供理論預測.

圖2 UO2與U0.75Th0.25O2沿[001]和[111]方向單軸拉伸加載過程 (a) 應力-應變曲線; (b) 沿[001]拉伸過程中的徑向分布函數;(c) 沿[111]拉伸過程中的徑向分布函數; (d) U0.75Th0.25O2沿[001]拉伸的原子結構演化; (e) U0.75Th0.25O2沿[111]拉伸時的原子結構演化Fig.2.Tensile behaviors of UO2 and U0.75Th0.25O2 along [001] and [111] direction: (a) Stress-strain curves; (b) radial distribution function (RDF, g(r)) along [001] direction; (c) RDF along [111] direction; (d) the atomic structure evolution of U0.75Th0.25O2 along[001] direction; (e) the atomic structure evolution of U0.75Th0.25O2 along [111] direction.

3.2 濃度與溫度對力學性能的影響

鑒于沿[111]晶向拉伸時未發生相變, 以其為典型代表考察摻雜濃度和溫度對U1—xThxO2混合燃料力學性能及斷裂行為的影響.圖3(a)—(c)分別給出了UO2, U0.95Th0.05O2, 以及U0.45Th0.55O2混合燃料在300—1200 K溫度范圍內單軸拉伸加載時的應力-應變曲線.針對不同摻雜濃度的混合燃料, 在室溫下進行加載時, 應力隨應變的增加而增加, 當達到臨界值后迅速下降, 對應脆性斷裂,該點對應體系的斷裂應力和應變.隨著溫度的升高, 當應力增加到臨界值后發生一定程度的流動,然后突降, 即UO2晶體發生斷裂.對線彈性范圍(應變在0—0.04之間)的應力-應變曲線進行線性擬合, 可得到彈性模量.圖3(d)—(f)分別給出了對應的三種體系彈性模量隨溫度的變化.當溫度由300 K升高至1200 K時, UO2的彈性模量由140.3 GPa迅速降至118.2 GPa; U0.95Th0.05O2的彈性模量由139.93 GPa線性減小至121.04 GPa;U0.45Th0.55O2的彈性模量由140.14 GPa減小至121.64 GPa.圖3(g)—(i)顯示出對應的斷裂應力隨溫度升高而減小的變化趨勢, 圖3(j)—(l)顯示出斷裂應變隨溫度的升高而增加的變化趨勢.當溫度由300 K升高至1200 K時, UO2的斷裂應力由9.53 GPa線性降低至7.23 GPa, 斷裂應變則由0.108增加至0.15; U0.95Th0.05O2的斷裂應力由9.96 GPa減至7.46 GPa, 斷裂應變則由0.105增至0.141; U0.45Th0.55O2的斷裂應力由10.146 GPa減至8.736 GPa, 斷裂應變則由0.102增至0.145.

圖3 溫度對U1—xThxO2 混合燃料沿[111]方向單軸拉伸時的力學性能的影響 (a), (b), (c) UO2, U0.95Th0.05O2及U0.45Th0.55O2不同溫度下的應力-應變曲線; (d), (e), (f) 對應的彈性模量隨溫度的變化; (g), (h), (i) 對應的斷裂應力隨溫度的變化; (j), (k), (l) 對應的斷裂應變隨溫度的變化Fig.3.Effect of temperature on the mechanical properties of U1—xThxO2 loaded along the [111] direction: (a), (b), (c) Stress-strain curves of UO2, U0.95Th0.05O2 and U0.45Th0.55O2; (d), (e), (f) the corresponding elastic modulus as a function of temperature; (g), (h),(i) the corresponding fracture stress as a function of temperature; (j), (k), (l) the corresponding fracture strain as a function of temperature.

總體來講, UxTh1—xO2混合燃料體系的力學性能強烈依賴于溫度, 不論摻雜濃度有多大, 彈性模量均隨溫度的升高而減小, 根據經典的固體物理理論, 晶體的彈性模量可表述為[39]

式中E0為0 K時的彈性模量;T0為愛因斯坦溫度;a0為材料相關的常數, 由Mie-Gruneisen描述:

其中r和δT為 Gruneisen 和 Anderson-Gruneisen參數, 與溫度相關;V0為原子的本征體積.根據(10)式和(11)式, 可得彈性模量與溫度呈負相關關系.針對UxTh1—xO2混合燃料體系, 其熔點均超過3000 K, 一般認為溫度低于1200 K時表現脆性, 僅當溫度高于1200 K時, 才會發生塑性變形, 我們的結果也佐證了這一觀點.在核裂變反應堆中燃料的服役溫度通常在1200 K左右, 此時UxTh1—xO2混合燃料在受到機械應力的情況下易發生脆性斷裂, 進而導致部件失效, 造成安全事故.

圖4詳細分析了不同溫度下U1—xThxO2混合燃料力學性能隨摻雜濃度的變化趨勢.圖4(a)給出了300 K時不同摻雜濃度對應的應力-應變曲線, 摻雜濃度的增加未引起應力-應變曲線顯著變化.圖4(b)為不同溫度條件下彈性模量隨摻雜濃度的變化趨勢, 當摻雜濃度小于0.1時, 彈性模量隨摻雜濃度的增加減小, 而摻雜濃度高于0.1時,彈性模量呈增加趨勢.如圖4(c)和圖4(d)所示,斷裂應力隨摻雜濃度的增加而提升, 但斷裂應變隨濃度的增加而減小.

圖4 Th摻雜濃度對U1—xThxO2 混合燃料沿[111]方向單軸拉伸時的力學性能的影響 (a) 300 K下不同摻雜濃度時的應力-應變曲線; (b) 不同溫度下彈性模量隨摻雜濃度的變化; (c) 不同溫度下斷裂應力隨摻雜濃度的變化; (d) 不同溫度下斷裂應變隨摻雜濃度的變化Fig.4.Effect of Th concentration on the mechanical properties of U1—xThxO2 loaded along the [111] direction: (a) Stress-strain curves; (b) the elastic modulus as a function of Th concentration; (c) the fracture stress as a function of Th concentration; (d) the fracture strain as a function of Th concentration.

不同溫度和摻雜濃度下, U1—xThxO2混合燃料的應力在達到最高值后均發生突降, 預示其發生脆性斷裂.圖5以UO2與U0.45Th0.55O2為例, 給出了其在室溫和高溫下拉伸加載過程中的原子結構演化, 其中原子顏色由其切應變標定, 紅色和綠色代表較大應變, 藍色代表較小應變, 計算公式為

式中,h為原子i的局部Lagrangian應變矩陣.針對不同的UxTh1—xO2混合燃料, 在承受拉伸加載時, 彈性應變范圍內, 離子鍵長變形伸長, 在達到斷裂應變后, 在切應變較大區域內U—O或Th—O離子鍵破壞, 形成初始裂紋, 然后迅速擴展直至最終斷裂.斷裂面與拉伸方向呈45°左右, 滿足經典的脆性材料斷裂準則.利用位錯拾取算法(dislocation extraction algorithm, DXA)對斷裂時的原子結構進行表征, 如圖5中最后一列所示,灰色曲面為裂紋表面, 裂紋尖端附近只有一根1/2〈110〉全位錯, 意味著基本沒有塑性變形, 證明了U1—xThxO2混合燃料的脆性斷裂特征, 與Xiao等[12]、Balboa等[19]和Zhang等[30]的MD模擬結果及Mo等[27]的納米壓入結果相吻合.

圖5 U1—xThxO2 混合燃料沿[111]方向單軸拉伸時的原子結構演化 (a) UO2, 300 K; (b) UO2, 1000 K; (c) U0.45Th0.55O2, 300 K;(d) U0.45Th0.55O2, 1200 KFig.5.Typical atomic structure evolution of U1—xThxO2 upon tensile loading along [111] direction: (a) UO2, 300 K; (b) UO2, 1000 K;(c) U0.45Th0.55O2, 300 K; (d) U0.45Th0.55O2, 1200 K.

實際當中制備的UO2燃料芯塊通常為多晶,晶界在U1—xThxO2混合燃料力學性能與斷裂特性中扮演重要角色, 本文以UO2, U0.75Th0.25O2和U0.50Th0.50O2為例, 考察了溫度與摻雜濃度對多晶混合燃料體系力學特性的影響.基于Voronoi算法構建了三維方向尺寸均為20 nm的多晶模型, 平均晶粒尺寸為8.5 nm.為使晶界充分弛豫, 將多晶模型在較高溫度(1200 K)下弛豫100 ps, 再退火至300 K繼續弛豫100 ps.圖6(a)為在不同溫度下三種多晶模型沿z軸拉伸時的應力-應變曲線,應力隨應變的增加達到峰值后發生緩慢下降, 保持了一定程度的塑性流動, 與文獻[12,29]的MD模擬結果吻合.圖6(b)和圖6(c)分別顯示了彈性模量與斷裂應力隨溫度的變化趨勢.溫度升高,彈性模量與斷裂應力均減小, 與單晶中的變化趨勢一致, 特別的是, 多晶中彈性模量對摻雜濃度不敏感, 而斷裂應力卻隨摻雜濃度的增加而提高.以U0.75Th0.25O2為例, 圖6(d)為300 K時不同應變下的原子結構, 晶界處的離子所受切應變較大, 在拉伸應變達到臨界值后率先斷裂, 即發生了沿晶斷裂.圖6(e)為拉伸過程中某三角晶界附近的原子結構演化過程, 可清楚地觀察到晶界處的離子鍵的破壞, 然后沿晶界方向持續擴展, 直至最終斷裂.Feng等[40]873和1973 K高溫下將UO2燒結成高密度的板狀塊材, 通過聚焦離子束切割加工成微米單晶或多晶材料, 然后在樣品表面預先切割出“V”形缺口, 進行原位單軸拉伸, 發現所有樣品均呈脆性斷裂特征, 與本文的MD模擬結果高度吻合.

圖6 多晶U1—xThxO2 混合燃料單軸拉伸力學特性 (a) 三種摻雜濃度下不同溫度時的應力-應變曲線; (b) 彈性模量隨溫度的變化; (c) 斷裂應力隨溫度的變化; (d) 斷裂過程中的三維結構; (e) 不同應變下三角晶界區域的放大圖, 原子顏色由其應變標定Fig.6.Mechanical behaviors of polycrystalline U1—xThxO2: (a) Stress-strain curves for different temperature and Th concentration;(b) the elastic modulus as a function of temperature; (c) the fracture stress as a function of temperatures; (d) the three-dimensional atomic structure; (e) the atomic structure evolution around a triple grain boundary.

4 結 論

本文利用分子動力學模擬, 系統研究了溫度與摻雜濃度對U1—xThxO2混合燃料結構與力學特性的影響.結果表明:

1) U1—xThxO2沿[001]晶向單軸拉伸, 可發生由初始的面心立方結構的螢石相至scrutinyite相的相變, 新相具有低對稱結構, 該結構已被第一性原理計算證實.

2) UxTh1—xO2混合燃料體系的力學性能強烈依賴于溫度, 不論摻雜濃度有多大, 彈性模量均隨溫度的升高而減小, 斷裂應力均隨溫度的升高而減小, 但斷裂應變卻呈增加趨勢, 符合經典固體物理理論; 當摻雜濃度小于0.1時, 彈性模量呈下降趨勢,而摻雜濃度高于0.1時, 彈性模量呈增加趨勢; 斷裂應力隨摻雜濃度的增加而增加, 斷裂應變則減小;

3) 溫度在300—1200 K范圍內, 不同摻雜濃度條件下混合燃料體系均表現脆性斷裂特性, 多晶樣品中發生脆性沿晶斷裂.

本文發現了UxTh1—xO2混合燃料體系具有較優異的熱力學特性, 可作為新型核燃料的備選, 但未來還需針對其中缺陷, 如氧空位、位錯與晶界等的復雜作用展開研究.