TRISO顆粒SiC層輻照行為與力學性能的分子動力學模擬

嚴澤凡,劉澤兵,田 宇,劉榮正,劉 兵,邵友林,唐亞平,劉馬林

(清華大學 核能與新能源技術研究院,北京 100084)

碳化硅(SiC)材料具有優異的力學性能、熱學性能和抗輻照性能,是一類非常重要的核材料,其在先進核能系統中有廣泛應用。目前SiC已被作為高溫氣冷堆所使用的三元結構各向同性(tri-structural isotropic, TRISO)包覆核燃料顆粒的包覆層材料[1-2],通常使用流化床-化學氣相沉積(fluidized bed-chemical vapor deposition, FB-CVD)方法包覆制備[3]。在TRISO顆粒中,SiC作為主要的裂變產物阻擋層和承壓層,可保證燃料顆粒的結構穩定性[4-5]。因此,SiC層對TRISO顆粒的安全性能有重要影響,有必要對SiC層的輻照行為和輻照下力學性能變化進行研究。

在輻照實驗中,通常難從微觀角度細致觀察樣品結構在輻照過程中的演化情況。采用分子動力學(molecular dynamics, MD)模擬可精確描述輻照過程中的缺陷和微觀結構演變,有助于分析材料的輻照行為。目前已有大量學者對SiC材料的輻照行為進行了MD模擬研究[6-9],但目前相關研究集中在SiC單晶的單次級聯或級聯交疊,針對TRISO顆粒復雜的SiC層微觀結構的MD級聯交疊模擬有待深入研究。輻照會引起SiC層力學性能的改變,可通過納米壓痕方法獲取力學性能和力學行為的信息[10],從而幫助探究SiC層因輻照發生力學性能變化的原因。MD模擬對開展納米壓痕研究有著類似的優勢,有助于分析材料的力學性能。針對SiC材料的納米壓痕的MD模擬研究有很多文獻報道[11-13],但研究對象大多為單一類型的SiC單晶或多晶,缺乏對輻照后的TRISO顆粒SiC層進行MD納米壓痕研究。因此,本文擬采用MD模擬詳細研究TRISO顆粒SiC層在級聯交疊過程中的輻照行為以及力學性能。

本文將首先通過MD模擬計算SiC在輻照中的體積腫脹理論值,以及納米壓痕下的力學性能理論值,并與實驗值比較,以證明使用Tersoff/ZBL勢進行SiC層的級聯交疊模擬,以及使用Vashishta勢和Tersoff勢進行SiC層的納米壓痕MD模擬的適用性。為研究SiC層的輻照行為,考察SiC層在級聯交疊過程中的輻照行為,通過腫脹程度、密度、原子結構類型、點缺陷演化等參量對輻照行為進行定量化分析,并計算不同輻照劑量下SiC層的力學性能。為進一步研究SiC層在輻照前后的力學性能變化原因,通過載荷-深度曲線、應力應變等參量對SiC層的力學行為進行定量化分析,以在微觀層次解釋輻照對SiC層力學性能的影響。

1 計算方法

1.1 模擬體系的提出和構建

本文研究團隊對TRISO顆粒制備工藝有著長期的實驗探索和優化。通過對TRISO顆粒的長期實驗研究發現,不同包覆條件下制備的SiC層的微觀結構不同,當包覆濃度適中、顆粒較多時制備的SiC層多為等軸狀多晶,當包覆濃度極低、顆粒較少時制備的SiC層多為長軸狀多晶。實驗中獲得的不同SiC層電子背散射衍射(electron back scatter diffraction, EBSD)圖像如圖1所示。

圖1 TRISO顆粒SiC層的分子動力學模擬基礎

采用LAMMPS軟件[14]進行MD模擬,后處理主要采用OVITO軟件[15]。本文使用的各類SiC層模型如圖1所示。圖1中的模型采用識別金剛石結構(identify diamond structure, IDS)方法[16]顯示。

SiC單晶用于驗證MD方法,在實際包覆中難以出現,而等軸狀多晶和長軸狀多晶在實際包覆中常見,故選取模型Ⅱ、Ⅲ進行輻照模擬研究。模型先在300 K、0 Pa的條件下用NPT系綜進行10 ps的熱平衡。然后從中隨機選擇1個Si原子作為初級撞出原子(primary knock-on atom, PKA),賦予其5 keV的動能,使其引發碰撞級聯并造成輻照損傷。為緩解輻照引起的劇烈升溫,在模型邊緣設置厚度為0.3 nm、溫度為300 K的恒溫區,使用NVT系綜控溫。剩余區域為使用NVE系綜的牛頓區,使其運動受牛頓第二定律控制,模擬原子間因碰撞產生的動能傳遞。模擬體系在x、y、z方向上均為周期性邊界條件。一次輻照從賦予PKA動能開始,模擬體系在經過20 ps的碰撞級聯后,在300 K、0 Pa的條件下用NPT系綜對模型熱平衡20 ps。級聯交疊通過不斷重復以上過程實現,本文共連續模擬了1 000次級聯,輻照劑量達0.444 dpa。

進行力學性能研究的模擬體系包括SiC層和剛體球形金剛石壓頭。在壓痕開始前,先將模擬體系能量最小化,然后在300 K、0 Pa的條件下用NPT系綜熱平衡。在壓痕過程中,SiC層分固定區、恒溫區和牛頓區。固定區用于固定邊界原子,消除壓痕過程中SiC層模型的剛體運動。恒溫區在NVE系綜下采用速度標定法將SiC層模型的平衡溫度維持在300 K。牛頓區采用NVE系綜,使其運動受牛頓第二定律控制。壓頭先沿z軸向下加載,并在壓入SiC層表面與壓頭半徑相同的深度后,以大小相等、方向相反的速度卸載。模擬體系在x、y軸方向上為周期性邊界條件,在z軸方向上為固定邊界條件。模擬過程的時間步長為0.001 ps。納米壓痕過程的模擬體系設置如圖2所示。

圖2 納米壓痕過程的模擬體系設置

1.2 MD方法和勢函數

在MD模擬研究中,勢函數被用于描述原子間的相互作用,是模擬的基礎和關鍵。本文在SiC層輻照行為研究中采用Tersoff/ZBL勢描述SiC層中的原子相互作用;在SiC層力學性能研究中采用Vashishta勢描述SiC層中的原子相互作用,采用Erhart和Albe[17]改良的Tersoff勢描述SiC層與金剛石之間的原子相互作用以及金剛石中的原子相互作用。

Tersoff/ZBL勢被廣泛用于SiC輻照的MD模擬研究[18-19],其基本形式如下所示[20]:

VTersoff/ZBL(rij)=(1-fF(rij))VZBL(rij)+

fF(rij)VTersoff(rij)

(1)

其中:VTersoff/ZBL為Tersoff/ZBL勢;VZBL為ZBL勢,用于修正Tersoff勢在描述近程相互作用時的不足;fF為一類費米函數,用于平滑地連接ZBL勢和Tersoff勢;rij為原子i和j之間的距離。

Tersoff勢的經典形式如下所示[21]:

VTersoff(rij)=fC(rij)(fR(rij)+bijfA(rij))

(2)

其中:VTersoff為Tersoff勢;fR為二體勢,代表排斥作用;fA為與鍵合相關的吸引作用;bij為連接原子i和j的鍵級,代表了局部鍵合并確定勢對鍵角的依賴性;fA與bij的乘積為三體勢;fC為截止函數。

Vashishta勢被廣泛用于SiC材料納米壓痕的MD模擬研究[22-23],其基本形式如下所示[24]:

(3)

(4)

(5)

1.3 模擬方法的驗證

為驗證Vashishta勢和Tersoff勢描述SiC層力學行為的可靠性,根據Oliver&Pharr方法,計算出模型Ⅰ的楊氏模量,并與實驗值進行對比,匯總列于表1。同時計算出模型Ⅰ～Ⅲ的硬度理論值,并與不同晶粒尺寸SiC多晶硬度的實驗值匯總示于圖3所示,此處的晶粒尺寸為多晶中晶粒的平均直徑,可看出,通過MD模擬獲得的力學性能理論值與實驗測量值較接近。

表1 SiC單晶楊氏模量的理論值和實驗值對比

圖3 模型Ⅰ～Ⅲ的硬度理論值與SiC多晶硬度實驗值的對比[27]

為進一步驗證Tersoff/ZBL勢描述SiC層輻照過程的可靠性,對模型Ⅱ和模型Ⅲ在輻照過程中的體積腫脹程度進行計算,并與輻照得到的實驗值進行對比,如圖4所示,可看出,本文體積腫脹程度的MD計算結果與輻照實驗結果基本一致。綜上,本文采用的勢能函數和相關參數可用于SiC輻照行為和力學性能模擬。

圖4 模型Ⅱ、Ⅲ在輻照條件下的體積腫脹程度理論值與實驗值對比[28]

2 結果與討論

2.1 SiC輻照腫脹分析

為研究SiC層的輻照行為,給出SiC層的輻照腫脹與非晶化的演化過程,并進行定量化分析。模型Ⅱ、Ⅲ在輻照過程中的體積腫脹程度與密度如圖5所示。從圖5可看出,兩種模型在輻照過程中發生了體積膨脹、密度下降的現象。兩種晶體模型的變化過程基本一致,說明長軸晶和等軸晶受輻照影響規律類似。SiC層的體積腫脹程度先迅速升高,在輻照劑量0.2 dpa左右開始放緩,并逐漸達到穩定;密度的變化趨勢與體積腫脹程度類似。

圖5 模型Ⅱ、Ⅲ在輻照過程中的體積腫脹程度與密度

2.2 SiC輻照非晶化分析

通過對模型Ⅱ、Ⅲ在輻照過程中各類結構的原子比例進行統計有助于了解輻照過程中的缺陷和非晶化的演化情況,如圖6所示。其中金剛石第一近鄰結構和金剛石第二近鄰結構可看作金剛石結構與非晶結構之間的中間態,前者更接近晶體結構,后者更接近非晶結構。兩種模型在輻照過程中各類結構的原子比例基本一致。在輻照初期,晶體結構(立方金剛石)迅速減少,中間態結構(立方金剛石第一近鄰、立方金剛石第二近鄰)和非晶結構原子迅速增加。在輻照劑量達到0.05 dpa左右時,中間態結構原子的比例達最大,之后與晶體結構原子一起緩慢減少,共同轉變為非晶結構原子。在輻照劑量達0.25 dpa左右時,晶體結構與中間態結構的原子比例基本降為0,而非晶結構的原子比例也達最大,此后一直穩定不變。值得注意的是,中間態結構中立方金剛石第二近鄰結構的原子比例相較于立方金剛石第一近鄰結構有一定滯后性,說明前者在一定程度上是由后者轉化而來。

圖6 模型Ⅱ、Ⅲ在輻照過程中的各類原子比例

圖7展示了模型Ⅱ、Ⅲ在輻照過程中的中間態結構和非晶結構原子的演化過程,其中用紅色標記了初始結構中的晶界原子,并隱去了晶體結構原子,中間態結構原子和非晶結構原子的顏色標注與圖1中一致。可看出,在輻照過程中兩種模型的中間態結構原子均傾向于先以團簇形式在晶界附近生成。在中間態結構原子團簇不斷長大并充滿晶粒內部后,非晶結構原子同樣傾向于以團簇形式出現在晶界附近,并不斷擴張,直至完全充滿晶粒,使SiC層完全非晶化。這些現象與文獻[29-30]中的輻照過程點缺陷結果基本一致。如Jin等[29]研究了晶界對SiC輻照缺陷的影響。研究表明,點缺陷傾向于在晶界附近以團簇形式積累,因為具有較高內能和拓撲無序的晶界附近的晶格比晶粒中的晶格更易受到損傷。而點缺陷的形成與非晶化過程存在強相關性,這將下文中進行討論。這也使非晶化過程易在晶界附近發生。

圖7 模型Ⅱ、Ⅲ在輻照過程中的非晶化演化過程

以上結果說明,SiC層的輻照腫脹與非晶化存在緊密的聯系,且它們不會因晶粒結構改變而發生顯著改變。如對于模型Ⅱ,其輻照腫脹程度(圖5)和非晶結構原子比例(圖6)在輻照劑量小于0.2 dpa時迅速增加,但當輻照劑量超過0.2 dpa時增加緩慢并趨于穩定。這意味著SiC層在輻照初期輻照腫脹與非晶化程度均迅速增加,但在輻照劑量飽和后趨于穩定。輻照過程中的非晶化并非直接由晶體結構轉變為非晶結構,而是存在晶體結構轉化為中間態結構,再轉化為非晶結構的過程,且這種過程傾向于從晶界附近開始發展。

2.3 SiC輻照Frenkel對分析

在輻照過程中會產生由空位和間隙原子構成的缺陷系統,即Frenkel對。為了解SiC層在輻照過程中的Frenkel對演化情況,可用Wigner-Seitz方法[31]統計了SiC層在輻照過程中的Frenkel對數量與各類空位和間隙原子比例。鑒于模型Ⅱ和模型Ⅲ的輻照行為沒有明顯差異,此處以模型Ⅱ為例進行分析,如圖8所示。可發現Frenkel對數量隨輻照劑量的增長與圖6中的非晶結構原子類似,這說明SiC層在輻照過程中的非晶化與點缺陷的生成存在密切關系。

圖8 模型Ⅱ在輻照過程中的Frenkel對數量與各類空位和間隙原子比例

對于空位,在輻照早期C空位的比例遠大于Si空位,隨著輻照劑量的增加,二者之間的比例差距逐漸減小,在輻照劑量0.25 dpa左右時,二者的比例變得幾乎一樣。對于間隙原子,在輻照早期Si間隙原子的比例大于C間隙原子,隨著輻照劑量的增加,二者之間的比例差距先增加后減小,在輻照劑量0.25 dpa左右時,二者的比例變得幾乎一樣。這是因為C原子的離位閾值低于Si原子的離位閾值,所以在輻照早期C原子在輻照過程中更易離開原有點位,傾向于形成空位。而Si原子在輻照過程中位置的偏移量更小,所以更傾向于形成間隙原子。但隨輻照劑量的增大,C空位和Si間隙原子趨于飽和,這些傾向表現得越來越不明顯,所以C空位和Si空位、C間隙原子和Si間隙原子的比例趨于接近。

圖9示出了模型Ⅱ在輻照過程中Frenkel對的演化過程,其中圖9a～d為空位的演化過程,圖9e～h為間隙原子的演化過程。C空位用綠色標記,Si空位用紫色標記,C間隙原子用橄欖色標記,Si間隙原子用棕色標記。可看出,空位和間隙原子的演化過程與非晶化過程類似,傾向于在輻照過程中以團簇形式在晶界附近生成,然后不斷擴展并充滿晶粒內部。

綠色為C空位,紫色為Si空位,橄欖色為C間隙原子,棕色為Si間隙原子

2.4 SiC輻照反位原子分析

在輻照過程中,也會發生某種原子替代另一種原子點位的現象,即反位原子。對于SiC,如果1個Si原子占據了原本屬于C原子的點位,那么該Si原子就是Si反位原子;同理,C反位原子就是C原子占據了原本屬于Si原子的點位產生的。此處同樣以模型Ⅱ為例,采用Wigner-Seitz方法統計了SiC層在輻照過程中的反位原子數量與各類反位原子比例,如圖10所示。可發現反位原子數量曲線與圖8a中的Frenkel對數量曲線趨勢基本一致。Si反位原子與C反位原子比例的變化趨勢與圖8b中的Si空位和C空位基本一致。這也是由于C原子的離位閾值低于Si原子造成的,所以在輻照早期C原子比Si原子更易取代對方成為反位原子。但輻照劑量的增加使Si反位原子和C反位原子的比例趨近。這種趨近現象與C原子和Si原子的離位閾值差異有關[32-33]。C反位原子更易形成,其數量在輻照過程中會更易達到飽和,這主要體現在輻照早期C反位原子的比例大于Si反位原子的比例。而Si反位原子更難形成,其數量的增加相較于C反位原子有一定滯后性,但隨輻照過程的進行,Si反位原子的數量逐漸達到與C反位原子的數量接近的水平。出現圖8中現象的原因也與之類似。

圖10 模型Ⅱ在輻照過程中的反位原子數量與各類反位原子比例

圖11示出了模型Ⅱ在輻照過程中的反位原子演化過程。C反位原子用橙色標記,Si反位原子用粉色標記。可看出,反位原子的演化過程與空位和間隙原子的演化過程基本一致,這說明了輻照過程中各類點缺陷的演化之間存在強相關性。

橙色為C反位原子,粉色為Si反位原子

2.5 輻照前后的力學性能分析

為了解輻照過程中SiC層的力學性能變化,對不同輻照劑量下的模型Ⅱ和模型Ⅲ進行納米壓痕測試。計算得到的SiC層在輻照過程中的力學性能變化如圖12所示。由圖12可知,模型Ⅱ和模型Ⅲ的硬度和楊氏模量均在輻照初期隨輻照劑量的增加而迅速下降,在輻照劑量超過0.2 dpa后達到穩定,幾乎不再下降。這說明模型Ⅱ和模型Ⅲ在力學性能趨勢變化上是相似的。輻照會導致SiC層力學性能的降低,但在輻照缺陷趨于飽和后不再對SiC層的力學性能有顯著影響。

圖12 模型Ⅱ、Ⅲ在輻照過程中的力學性能

為了解輻照前后SiC層在納米壓痕過程中的力學行為,匯總模型Ⅱ和模型Ⅲ在輻照前后的載荷-深度曲線,如圖13所示。可看出,輻照前的SiC層在加載過程中的載荷隨深度的增加迅速上升,而輻照后的SiC層則上升緩慢,使前者在加載終點時的載荷遠高于后者。輻照前后的SiC層在加載過程中均發生了彈進(pop-in)現象,即加載曲線中出現了載荷的突然下降,這是彈性變形到塑性變形的轉變標志。輻照前的SiC層在加載后期出現了多次明顯的彈進現象,使其加載曲線較為粗糙;而輻照后的SiC層在加載過程中出現的彈進幅度均較小,使其加載曲線更光滑。這說明輻照會使SiC層在外力作用下的承受能力和塑性變形程度減小。

圖13 模型Ⅱ、Ⅲ輻照前后的載荷-深度曲線

為了解SiC層在納米壓痕測試過程中的應力應變分布,此處計算并給出模型Ⅰ～Ⅲ和輻照后的非晶SiC層在加載終點處沿y軸半剖面的原子Von Mises應力[34]和剪切應變,同時顯示了晶粒的分布,便于相互對照,如圖14所示。與模型Ⅰ相比,模型Ⅱ和模型Ⅲ的應力應變分布不僅局限在壓頭附近,而是與晶界的分布高度相關。模型Ⅱ的應力應變分布更傾向于沿著等軸晶的晶界橫向擴展;模型Ⅲ的應力應變分布傾向于沿著長軸晶的晶界縱向擴展;而非晶SiC層的應力應變分布則顯得沒有規律。這說明輻照會使SiC層在外力作用下的應力應變分布紊亂。

圖14 模型Ⅰ～Ⅲ和非晶SiC層的晶粒分布與應力應變分布

3 結論

本文通過分子動力學模擬研究了TRISO顆粒SiC層的輻照行為和力學性能,可得出如下結論。

1) SiC層在輻照過程中腫脹程度的理論值與實驗值基本一致,證明了Tersoff/ZBL勢和其他模擬體系參數對SiC層輻照行為研究適用。通過Oliver&Pharr方法計算的SiC單晶的楊氏模量和各種SiC多晶的硬度理論值與實驗值吻合較好,證明了Vashishta勢和Tersoff勢和其他模擬體系參數對SiC層力學性能研究適用。

2) 分子動力學模擬研究中,可通過腫脹程度、密度、原子結構類型、點缺陷演化等參量對輻照行為進行定量化分析。SiC層的晶粒結構對其輻照腫脹與非晶化影響較小。SiC層的輻照腫脹與非晶化程度在輻照初期迅速增加,但在輻照劑量飽和后趨于穩定。輻照過程中的非晶化并非直接由晶體結構轉變為非晶結構,而是存在晶體結構轉化為中間態結構,再轉化為非晶結構的過程。在點缺陷的演化過程中,C原子和Si原子的離位閾值差異導致點缺陷在早期以C空位、Si間隙原子和C反位原子為主,但隨著輻照劑量趨于飽和這些差異逐漸消失。非晶化和點缺陷傾向于從晶界附近開始發展。長軸晶和等軸晶受輻照影響規律類似。

3) 分子動力學模擬研究中,納米壓痕過程可通過載荷-深度曲線、應力應變等參量描述,更有助于分析輻照對SiC層力學性能的影響。長軸晶和等軸晶的受力主要受晶界方向影響。輻照會導致SiC層力學性能的降低,但在劑量趨于飽和后,不再對SiC層的力學性能有顯著影響。SiC層力學性能的降低與其在外力作用下的承受能力和塑性變形程度減小、應力應變分布紊亂密切相關。