含有缺陷3C-SiC陶瓷拉伸性能的分子動(dòng)力學(xué)模擬

馬小強(qiáng)， 徐喻瓊， 蘇華山， 杜曉超， 袁顯寶， 周建軍

( 三峽大學(xué)機(jī)械與動(dòng)力學(xué)院 水電機(jī)械設(shè)備設(shè)計(jì)與維護(hù)湖北省重點(diǎn)實(shí)驗(yàn)室， 宜昌 443002)

1 引 言

隨著核能發(fā)展，新一代核能系統(tǒng)對(duì)反應(yīng)堆用材料提出了更高要求，比如，要有較好的抗輻照性能、抗蠕變以及耐腐蝕等. 碳化硅(SiC)具有耐高溫、常規(guī)化學(xué)惰性、低密度、低熱膨脹系數(shù)以及中子俘獲截面小等優(yōu)良性能[1]，是新一代核能和反應(yīng)堆結(jié)構(gòu)材料的重要候選材料之一[2]. SiC是以共價(jià)鍵為主的化合物，因缺少能促使晶體變形的滑移系統(tǒng)而呈現(xiàn)出明顯的脆性. SiC纖維增韌SiC基復(fù)合材料(SiCf/SiC)克服了陶瓷材料的脆性[3]，保持了高溫強(qiáng)度、抗蠕變性能、耐腐蝕等優(yōu)點(diǎn)，作為燃料包殼材料或反應(yīng)堆內(nèi)構(gòu)件材料，在輕水反應(yīng)堆[4]、高溫氣冷堆[5]、聚變堆[6]等有良好的應(yīng)用前景.

在反應(yīng)堆環(huán)境中，由于中子、裂變碎片等對(duì)SiCf/SiC復(fù)合材料輻照會(huì)產(chǎn)生缺陷，致使材料中出現(xiàn)微空洞、元素偏析、位錯(cuò)、非晶等微觀結(jié)構(gòu)變化[7-10]. 另外，由于材料制備工藝的限制，在制備SiC材料時(shí)不可避免地會(huì)在材料中產(chǎn)生空位、微空洞以及非化學(xué)計(jì)量比等缺陷[11]. 所有這些缺陷都可能使材料的力學(xué)性能改變而影響材料的正常使用，甚至可能導(dǎo)致發(fā)生事故[12].

SiCf/SiC是由SiC纖維和SiC陶瓷晶體經(jīng)過(guò)材料復(fù)合技術(shù)制備的具有優(yōu)越性能的復(fù)合材料. SiC晶體是由四個(gè)碳原子和一個(gè)硅原子(SiC4，硅原子位于正四面體中心，四個(gè)碳原子分別位于正四面體的四個(gè)頂點(diǎn))或者四個(gè)硅原子和一個(gè)碳原子(CSi4，C原子位于正四面體中心，四個(gè)硅原子分別位于正四面體的四個(gè)頂點(diǎn))交替以sp3鍵結(jié)合成的正四面體結(jié)構(gòu). 在不同物理、化學(xué)環(huán)境下，SiC形成的晶型不同，除了具有立方結(jié)構(gòu)的β-SiC(即3C-SiC)，其他的為α-SiC. 3C-SiC晶體由周期為3層的SiC正四面體結(jié)構(gòu)密排堆積而成的立方晶體結(jié)構(gòu)[12]，只有3C-SiC能夠滿足核能需要[2].

分子動(dòng)力學(xué)模擬是研究材料力學(xué)性能的重要手段的重要手段[13, 14]，對(duì)于3C-SiC力學(xué)性能的分子動(dòng)力學(xué)研究已經(jīng)有些成果[15-17]，這些成果主要研究的是無(wú)缺陷的單晶體體系、孿晶納米線和薄膜，而對(duì)一定缺陷(空位、微空洞以及反位替代)的晶體拉伸性能研究鮮有報(bào)道.

本文應(yīng)用分子動(dòng)力學(xué)程序LAMMPS[18]對(duì)含有空位、微空洞、以及反位替代三種不同缺陷體系的3C-SiC的拉伸性能進(jìn)行了分子動(dòng)力學(xué)模擬. 通過(guò)模擬計(jì)算，得到了體系能量和晶體晶格常數(shù)以及應(yīng)力-應(yīng)變曲線隨缺陷“濃度”的變化關(guān)系. 通過(guò)分析應(yīng)力-應(yīng)變曲線，得到了不同缺陷體系的楊氏模量、斷裂應(yīng)變、拉伸強(qiáng)度隨缺陷“濃度”的變化關(guān)系，最后分析了3C-SiC拉伸斷裂機(jī)理. 數(shù)據(jù)分析過(guò)程中應(yīng)用分子可視化軟件Atomeye[19]和Ovito[20]直觀的觀察了拉伸過(guò)程中3C-SiC結(jié)構(gòu)的變化情況.

2 計(jì)算模型與方法

3C-SiC原子間相互作用使用Tersoff勢(shì)函數(shù)描述[21]，該勢(shì)函數(shù)已成功應(yīng)用于3C-SiC納米線的單軸壓縮和拉伸等[22]的力學(xué)性能模擬，Tersoff勢(shì)函數(shù)具體討論見(jiàn)第2節(jié).

以晶胞的[100]、[010]、[001]方向?yàn)槟M體系的x，y，z坐標(biāo)軸方向建立10a×10a×10a(a=0.4318 nm，為晶格常數(shù))的3C-SiC單晶模型，該模型包含8000個(gè)原子，結(jié)構(gòu)如圖1所示. 為了模擬不同的缺陷，分別在無(wú)缺陷的單晶模型中隨機(jī)的刪除不同個(gè)數(shù)的原子形成不同“濃度”的空位型缺陷，或者在無(wú)缺陷晶體中刪除含有8個(gè)原子的團(tuán)簇，產(chǎn)生不同數(shù)量的微空洞. 由于制備工藝和方法的限制，制備的3C-SiC晶體和纖維的碳硅原子的比例不是嚴(yán)格的1:1，碳的所占比例大于硅[1]，因此，本工作以不同比例碳原子隨機(jī)的取代硅原子位置，形成反位替代缺陷. 三種不同的缺陷類型以及缺陷的“濃度”列于表1.

圖 1 3C-SiC晶體結(jié)構(gòu)的3D模型，藍(lán)色表示Si原子，黃色表示C原子Fig. 1 The 3D model of 3C-SiC, The blue balls and yellow balls denote the Si and C atoms, respectively

在拉伸模擬之前，對(duì)三種不同類型缺陷的所有模擬體系在NPT系綜下進(jìn)行1 ns的弛豫(弛豫時(shí)間步長(zhǎng)為1 fs)，使體系溫度達(dá)到300 K并形成穩(wěn)定結(jié)構(gòu). 以時(shí)間步長(zhǎng)為0.1 fs，拉伸應(yīng)變率為10-4/ps進(jìn)行拉伸模擬，拉伸時(shí)每100個(gè)積分步對(duì)原子速度重新調(diào)節(jié)以 控制系統(tǒng)溫度不變，而后再次加載應(yīng)力拉伸. 每一次加載在x軸([100]方向)方向的應(yīng)力會(huì)使體系在x軸方向產(chǎn)生一個(gè)微小的應(yīng)變?chǔ)牛搼?yīng)變使模擬體系的x軸方向大小變?yōu)?+ε，在拉伸變形的同時(shí)體系中的原子重新匹配盒子大小. 拉伸過(guò)程中計(jì)算每個(gè)原子應(yīng)力的x分量，體系的拉伸應(yīng)力是對(duì)模擬體系中所有原子應(yīng)力的x分量的平均值.

表1 不同缺陷體系參數(shù)

3 原子間勢(shì)函數(shù)

SiC原子間相互作用不僅取決于原子之間的距離，還和原子鍵的方向有關(guān). Tersoff多體勢(shì)能夠很好地描述SiC體系中的C-Si、C-C、Si-Si間的相互作用關(guān)系. Tersoff勢(shì)描述的系統(tǒng)勢(shì)能為任意兩個(gè)原子之間的勢(shì)能總和：

(1)

其中Vij(rij)是由吸引勢(shì)和排斥勢(shì)構(gòu)成的第i個(gè)原子和第j個(gè)原子之間的相互作用勢(shì)能：

Vij(rij)=fc(rij)[fR(rij)+bijfA(rij)]

(2)

其中，rij為第i個(gè)原子和第j個(gè)原子之間的距離；bij描述了除第j個(gè)原子之外所有原子對(duì)第i個(gè)原子的作用；fc(rij)為原子間相互作用的細(xì)化截?cái)嗪瘮?shù)，與截止半徑有關(guān)；fR(rij)是第i個(gè)原子和第j個(gè)原子之間的吸引對(duì)勢(shì)，其主要是指兩體相互作用；fA(rij)第i個(gè)原子和第j個(gè)原子之間的排斥對(duì)勢(shì)，其主要是指三體相互作用. 各個(gè)函數(shù)的具體形式表示為：

fR(rij)=-Ae-λijrij

(3)

fA(rij)=Be-μijrij

(4)

(5)

式中，A和B為吸引對(duì)勢(shì)和排斥對(duì)勢(shì)的結(jié)合能；R為截?cái)喟霃剑籖為fC(rij)的參數(shù)；λij和μij為吸引對(duì)勢(shì)與排斥對(duì)勢(shì)的勢(shì)能曲線梯度系數(shù).

上述各量的具體表達(dá)式為：

(6)

(7)

(8)

式中θijk為鍵ij和鍵ik之間的鍵角；β為鍵角系數(shù)；ζij為角勢(shì)能；c,d,h1為彈性常數(shù). 文中所用的參數(shù)的取值列于表2.

表2 Tersoff勢(shì)函數(shù)參數(shù)

4 結(jié)果與討論

4.1 體系能量和晶格常數(shù)

拉伸模擬之前，對(duì)構(gòu)建的不同“濃度”的缺陷體系在NPT系綜下弛豫1 ns，使體系的溫度達(dá)到300 K并使結(jié)構(gòu)穩(wěn)定. 對(duì)含有不同缺陷類型的穩(wěn)定結(jié)構(gòu)的勢(shì)能和晶格常數(shù)計(jì)算，計(jì)算結(jié)果如圖2(a)和圖3所示. 從圖2(a)可以看出，缺陷類型為空位和微空洞時(shí)，體系勢(shì)能隨著缺陷“濃度”增大而呈線性增加，反位替代型缺陷的體系勢(shì)能和沒(méi)有缺陷的體系勢(shì)能基本一致，沒(méi)有發(fā)生明顯的變化. 分析空位缺陷和微空洞缺陷的體系中原子間的成鍵情況可以發(fā)現(xiàn)，在相同缺陷“濃度”時(shí)，空位型缺陷體系中沒(méi)有完全成鍵的原子數(shù)要比微空洞體系的多，相同缺陷“濃度”時(shí)，空位缺陷體系勢(shì)能要高于的微空洞缺陷體系的勢(shì)能，相應(yīng)地，隨著缺陷“濃度”增加，空位缺陷體系的勢(shì)能增加要更快一些. 相同缺陷“濃度”時(shí)反位替代缺陷的能量低于其他兩種類型，這說(shuō)明反位替代型缺陷相對(duì)要穩(wěn)定一些.

圖 2 體系能量曲線 (a)穩(wěn)定結(jié)構(gòu)的原子勢(shì)能曲線； (b)體系能量隨應(yīng)變變化曲線(8%缺陷“濃度”為例)Fig. 2 energy Curve. (a) atomic potential energy with defect concentration； (b) total energy of system as a function of strain

圖 3 晶格常數(shù)隨缺陷“濃度”的變化Fig. 3 Lattice constant with defect concentration

從圖3可以看出，缺陷類型不同，晶格常數(shù)隨缺陷類型變化不同. 和沒(méi)有缺陷的晶格常數(shù)相比，微空洞體系的晶格常數(shù)變化不明顯，空位體系的晶格常數(shù)隨著缺陷“濃度”的增加略微減小(缺陷“濃度”為1.6%時(shí)減小了0.5%). 反位替代缺陷弛豫后形成穩(wěn)定結(jié)構(gòu)晶格常數(shù)隨著被替代原子數(shù)增加而減小，發(fā)生了較大的變化，在這種缺陷中，當(dāng)SiC中部分Si原子被C原子替代后，晶體中的部分Si-C鍵也相應(yīng)地被替換為C-C鍵，由于Si-C鍵(鍵長(zhǎng)為0.189 nm)變成C-C鍵(鍵長(zhǎng)為0.137 nm)[23]，鍵長(zhǎng)相應(yīng)減小，晶格常數(shù)也變小.

拉伸過(guò)程中，外力要對(duì)體系做功而使體系的能量發(fā)生變化，如圖2(b)所示為沒(méi)有缺陷的體系和缺陷“濃度”為8%時(shí)體系總能量(所有原子動(dòng)能和勢(shì)能的和)隨應(yīng)變變化關(guān)系. 由于模擬過(guò)程中保持溫度(300 K)不變，所以體系所有原子的總動(dòng)能不變，體系總能量的變化是由原子之間的勢(shì)能變化引起的. 從圖中可以看出，體系總能量隨著應(yīng)變的增加分為三個(gè)階段. 第一階段，應(yīng)變較小時(shí)(ε<0.07)原子勢(shì)能較低，此階段晶體形變屬于彈性階段，體系能量隨著應(yīng)變?cè)黾佣徛?，這個(gè)過(guò)程中系統(tǒng)的能量增加是由于原子之間的鍵長(zhǎng)發(fā)生了變化引起的. 第二階段(應(yīng)變?chǔ)?0.07到體系勢(shì)能最大)，這個(gè)階段體系能量隨著應(yīng)變?cè)黾映示€性增加，對(duì)此階段的原子構(gòu)型分析可以發(fā)現(xiàn)，該階段晶體體系中原子之間不只是鍵長(zhǎng)發(fā)生變化，鍵角和鍵結(jié)構(gòu)也會(huì)相應(yīng)發(fā)生變化(圖7)，從圖7還可以看出隨著應(yīng)變?cè)黾硬糠衷渔I斷裂，相應(yīng)的使體系能量進(jìn)一步增加而形成不穩(wěn)定結(jié)構(gòu). 第三階段，體系能量突然減小階段，這個(gè)突變是由于晶體發(fā)生脆性斷裂，大量原子之間失去相互作用，局部應(yīng)力卸載而釋放能量產(chǎn)生的. 由于已經(jīng)發(fā)生了斷裂，即使在此階段卸載應(yīng)力，體系不可能再恢復(fù)到初始狀態(tài).

4.2 應(yīng)力-應(yīng)變曲線

應(yīng)力—應(yīng)變曲線反應(yīng)了材料的基本力學(xué)性能，圖4是三種不同缺陷體系(a為空位缺陷，b為微空洞缺陷，c為反位替代缺陷)沿x方向拉伸的應(yīng)力-應(yīng)變曲線. 從該組曲線可以看出不同類型缺陷的應(yīng)力-應(yīng)變具有相同的變化趨勢(shì)，尤其是SiC在拉伸過(guò)程中會(huì)出現(xiàn)的脆性斷裂. 體系應(yīng)變較小時(shí)，應(yīng)力-應(yīng)變曲線滿足胡克定律，隨著應(yīng)變?cè)黾?，三種缺陷體系都存在應(yīng)力最大值，這個(gè)最大應(yīng)力就是 拉伸強(qiáng)度，應(yīng)變進(jìn)一步增加應(yīng)力會(huì)急劇減小，體系發(fā)生脆性斷裂，對(duì)應(yīng)的應(yīng)變?yōu)閿嗔褢?yīng)變.

圖 4 應(yīng)力-應(yīng)變曲線(a)空位缺陷體系； (b)微空洞缺陷體系；(c)反位替代缺陷體系 Fig. 4 Stress-strain curves for (a) vacancies； (b) micro-voids and； (c) antisite with defect concentration

楊氏模量是表征固體材料抵抗形變能力的物理量，取決于材料本身的性質(zhì)，其值為應(yīng)變?chǔ)?0.03是應(yīng)力-應(yīng)變曲線的斜率[15]，對(duì)模擬得到的應(yīng)力-應(yīng)變曲線中應(yīng)變小于0.03的部分?jǐn)M合求斜率，即可得到不同缺陷體系的楊氏模量. 例如，沒(méi)有缺陷的晶體材料，應(yīng)力隨著應(yīng)變可以增加到77.4 GPa(拉伸強(qiáng)度)后突然減小，此時(shí)應(yīng)變?yōu)?.311(斷裂應(yīng)變). 擬合應(yīng)變-應(yīng)力曲線，得到?jīng)]有缺陷的3C-SiC模擬體系的楊氏模量為485.9 GPa，而溫度為300 K時(shí)實(shí)驗(yàn)值為392-694 GPa[24]，計(jì)算結(jié)果和實(shí)驗(yàn)符合的很好.

4.3 力學(xué)性能參數(shù)

4.3.1楊氏模量

圖5(a)為不同缺陷體系的楊氏模量隨缺陷“濃度”變化的關(guān)系曲線. 從變化曲線可以看出，空位、和微空洞體系的楊氏模量和反位替代體系的楊氏模量隨缺陷“濃度”變化是不同的，前兩種情況的楊氏模量隨著缺陷“濃度”增加而呈現(xiàn)線性關(guān)系減小，缺陷“濃度”為1.6%時(shí)的楊氏模量分別減小了6.58%和5.35%，分析拉伸過(guò)程的體系原子構(gòu)型不難發(fā)現(xiàn)，空位和微空洞缺陷體系中有部分原子并沒(méi)有形成四個(gè)化學(xué)鍵而成為穩(wěn)定的正四面體結(jié)構(gòu)，這些沒(méi)有完全成鍵的原子不但使體系穩(wěn)定性變差(體系勢(shì)能較高)，而且影響了材料的力學(xué)性能. 反位替代缺陷的楊氏模量隨替代“濃度”增大呈現(xiàn)線性增加，缺陷“濃度”為1.6%的楊氏模量和沒(méi)有缺陷體系比較增加了7.1%，通過(guò)分析反位替代體系原子間成鍵情況可以發(fā)現(xiàn)，楊氏模量增大是因?yàn)楫?dāng)碳原子替代硅原子數(shù)增加時(shí)Si-C鍵數(shù)量減少，C-C鍵數(shù)量增加，在形變較小時(shí)C-C 的sp3鍵的強(qiáng)度Si-C 的sp3鍵強(qiáng)度大的緣故，因此隨著反位替“濃度”的增加楊氏模量也相應(yīng)地增加.

4.3.2拉伸強(qiáng)度和斷裂應(yīng)變

圖5(b)為不同缺陷體系的拉伸強(qiáng)度隨缺陷“濃度”變化的關(guān)系曲線，圖5(c)對(duì)應(yīng)的斷裂應(yīng)變. 從變化曲線可以看出，具有各種缺陷體系的拉伸強(qiáng)度隨著缺陷“濃度”增加而減小，并且減小速率有變慢. 對(duì)于空位和微空洞缺陷，在缺陷“濃度”小于0.8%時(shí)，這兩種缺陷的拉伸強(qiáng)度基本相同，缺陷“濃度”大于0.8%時(shí)，空位體系的拉伸強(qiáng)度比微 空洞體系的拉伸強(qiáng)度小，隨缺陷“濃度”的增加，這個(gè)差值變大，這是由于空位缺陷和微空洞缺陷相比，不完全成鍵的原子數(shù)量較多，在外力作用下，這些不完全成鍵原子更容易移動(dòng)，拉伸強(qiáng)度要降低的更快一些. 對(duì)于反位替代缺陷體系，當(dāng)應(yīng)變?cè)黾拥?.1時(shí)，C-C 的sp3鍵變?yōu)镃-C 的sp2弱鍵，隨著應(yīng)變進(jìn)一步增加C-C的sp2鍵處會(huì)形成空洞，以至于斷裂[13]，所以當(dāng)反位替代原子“濃度”增加時(shí)，拉伸強(qiáng)度相應(yīng)減小，并且減小的更快.

圖5 楊氏模量(a)，拉伸強(qiáng)度(b)和斷裂應(yīng)變(c)隨缺陷“濃度”的變化曲線Fig. 5 Young’s modulus (a) ，tensile strength(b)and(c)elongation percentage with defect concentration

4.4 斷裂機(jī)理

圖6為拉伸過(guò)程中包含182個(gè)原子的局部體系圖像，圖中紅色箭頭表示此位置原子之間沒(méi)有成鍵. 從圖6中可以看出在應(yīng)變小于0.075時(shí)，大部分鍵只是拉伸，沒(méi)有出現(xiàn)斷裂的情況，此過(guò)程彈性變形過(guò)程，只改變?cè)娱g鍵長(zhǎng)，如果此階段卸載外力，結(jié)構(gòu)還可以恢復(fù)到完美的穩(wěn)定結(jié)構(gòu). 隨著應(yīng)變進(jìn)一步增加(但小于0.25時(shí))，原子間所成的鍵有拉伸變形外，有部分原子鍵已經(jīng)斷裂，但有部分結(jié)構(gòu)還是晶體結(jié)構(gòu)，此過(guò)程是一個(gè)塑性變形過(guò)程，體系不完全成鍵原子數(shù)增多，甚至在局部出現(xiàn)空洞，這時(shí)如果卸載外力，結(jié)構(gòu)已經(jīng)不能回到?jīng)]有缺陷的完美狀態(tài). 當(dāng)應(yīng)變超過(guò)斷裂應(yīng)變式，發(fā)生脆性斷裂. 在室溫的納米壓痕實(shí)驗(yàn)中也觀察到了與上述過(guò)程相一致的結(jié)果[25]，在一次加載—卸載實(shí)驗(yàn)中，即在發(fā)生脆性斷裂以前，卸載曲線和加載曲線能夠完全重合.

單晶拉伸是基于剪切應(yīng)力在一個(gè)滑移面的一種變形機(jī)制，3C-SiC是金剛石結(jié)構(gòu)，滑移面為{111},模擬時(shí)的加載方向?yàn)閇100]，在(010)面上的剪切應(yīng)力不為零，所以，斷裂面不是完全的{111}方向，而是和{111}方向之間有一定夾角. 圖7是沒(méi)有缺陷體系脆性斷裂前后(ε=0.25為斷裂前，ε=0.31斷裂應(yīng)變，ε=0.33為斷裂后)體系構(gòu)型變化，從圖7可以看出斷裂面的方向不是嚴(yán)格的{1 1 1}.

5 結(jié) 論

本文采用LAMMPS程序?qū)锌瘴弧⑽⒖斩春头次惶娲N缺陷體系的3C-SiC中沿[100]方向的拉伸變形過(guò)程進(jìn)行了分子動(dòng)力學(xué)模擬，通過(guò)對(duì)拉伸過(guò)程中體系能量、應(yīng)力-應(yīng)變曲線、楊氏模量、拉伸強(qiáng)度以及斷裂機(jī)制的分析，得到了以下有用的結(jié)論：

(1)不同缺陷體系穩(wěn)定結(jié)構(gòu)的能量隨缺陷“濃度”增加變化趨勢(shì)不同，空位和微空洞體系的勢(shì)能隨缺陷“濃度”增加成線性增大；反位替代缺陷的體系能量和無(wú)缺陷體系能量基本相同，隨替代“濃度”基本沒(méi)有發(fā)生變化. 在拉伸過(guò)程中模擬的所有體系的勢(shì)能變化都經(jīng)歷了緩慢增大，近線性增大和突然減小的三個(gè)階段.

圖 6 拉伸過(guò)程中原子構(gòu)形(包含182個(gè)原子)Fig. 6 Atomic configurations of 182 atoms of the bulk 3C-SiC with the tensile

圖 7 脆性斷裂前后原子構(gòu)型Fig. 7 Atomic configurations before and after brittle fracture for 3C-SiC

(2) 空位和微空洞缺陷體系的楊氏模量隨著缺陷“濃度”增加呈線性減小，反位替代缺陷體系的楊氏模量隨缺陷“濃度”增加呈線性增加.

(3) 三種缺陷體系的拉伸強(qiáng)度和斷裂應(yīng)變都是隨著缺陷“濃度”增加而減小，空位和微空洞缺陷體系減小的更快一些.

(4)單軸拉伸斷裂面不是嚴(yán)格的在{111}滑移面.