馬小強， 徐喻瓊， 蘇華山， 杜曉超， 袁顯寶， 周建軍

( 三峽大學機械與動力學院 水電機械設備設計與維護湖北省重點實驗室， 宜昌 443002)

1 引 言

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

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

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

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

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

2 計算模型與方法

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

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

圖 1 3C-SiC晶體結構的3D模型，藍色表示Si原子，黃色表示C原子Fig. 1 The 3D model of 3C-SiC, The blue balls and yellow balls denote the Si and C atoms, respectively

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

表1 不同缺陷體系參數

3 原子間勢函數

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

(1)

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

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

(2)

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

fR(rij)=-Ae-λijrij

(3)

fA(rij)=Be-μijrij

(4)

(5)

式中，A和B為吸引對勢和排斥對勢的結合能；R為截斷半徑；R為fC(rij)的參數；λij和μij為吸引對勢與排斥對勢的勢能曲線梯度系數.

上述各量的具體表達式為：

(6)

(7)

(8)

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

表2 Tersoff勢函數參數

4 結果與討論

4.1 體系能量和晶格常數

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

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

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

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

拉伸過程中，外力要對體系做功而使體系的能量發(fā)生變化，如圖2(b)所示為沒有缺陷的體系和缺陷“濃度”為8%時體系總能量(所有原子動能和勢能的和)隨應變變化關系. 由于模擬過程中保持溫度(300 K)不變，所以體系所有原子的總動能不變，體系總能量的變化是由原子之間的勢能變化引起的. 從圖中可以看出，體系總能量隨著應變的增加分為三個階段. 第一階段，應變較小時(ε<0.07)原子勢能較低，此階段晶體形變屬于彈性階段，體系能量隨著應變增加而緩慢升高，這個過程中系統(tǒng)的能量增加是由于原子之間的鍵長發(fā)生了變化引起的. 第二階段(應變ε>0.07到體系勢能最大)，這個階段體系能量隨著應變增加呈線性增加，對此階段的原子構型分析可以發(fā)現，該階段晶體體系中原子之間不只是鍵長發(fā)生變化，鍵角和鍵結構也會相應發(fā)生變化(圖7)，從圖7還可以看出隨著應變增加部分原子鍵斷裂，相應的使體系能量進一步增加而形成不穩(wěn)定結構. 第三階段，體系能量突然減小階段，這個突變是由于晶體發(fā)生脆性斷裂，大量原子之間失去相互作用，局部應力卸載而釋放能量產生的. 由于已經發(fā)生了斷裂，即使在此階段卸載應力，體系不可能再恢復到初始狀態(tài).

4.2 應力-應變曲線

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

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

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

4.3 力學性能參數

4.3.1楊氏模量

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

4.3.2拉伸強度和斷裂應變

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

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

4.4 斷裂機理

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

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

5 結 論

本文采用LAMMPS程序對含有空位、微空洞和反位替代三種缺陷體系的3C-SiC中沿[100]方向的拉伸變形過程進行了分子動力學模擬，通過對拉伸過程中體系能量、應力-應變曲線、楊氏模量、拉伸強度以及斷裂機制的分析，得到了以下有用的結論：

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

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

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

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

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

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