點缺陷對單晶硅力學性能影響的分子動力學研究

萬 偉， 唐昌新， 項永康， 邱 安

(南昌大學 光伏研究院， 南昌 330031)

1 引 言

硅的豐度在地球上僅次于氧，儲量豐富，易于被人類提取與利用. 隨著科技的發展，單晶硅已在大規模集成電路和太陽能電池等領域獲得了廣泛的應用，例如：單晶硅晶圓作為太陽電池等產品的原料，是由線鋸工藝從硅錠上切割而成[1]. 然而，在硅錠的生產加工過程中會產生崩邊、隱裂和碎片等不良品[2]，降低了原材料的利用率，提高了經濟成本，不利于單晶硅產業的高速發展，更會對其他以硅材料為基礎的產業產生負面影響. 因此，有必要對單晶硅的力學性能展開研究，探究產生切割不良品的原因.

從單晶硅力學性能面臨的問題入手，Liu等人[3]采用壓痕法等實驗手段對單晶硅加工過程中產生的變形與表面損傷等問題進行研究，并對比結合分子動力學等手段所給出的模擬值來驗證晶體亞結構的變形與破壞理論. 由于分子動力學模擬方法可以提供原子微觀結構演變的詳細信息，為理解材料微觀性能提供了一個有力的研究工具. 在單晶硅產業發展的數十年歷程中，通過多學科綜合研究等一系列方法，我們對氮[4]、氧[5]等雜質原子以及位錯[6]等因素對硅力學性能的強化效應及其影響機理有了充分的認識. 另外在硅片的生產過程中，由于熱應力、快速結晶等因素，單晶硅中會形成點缺陷、點缺陷團簇甚至原子規模的小孔洞，這些因素都會對產品的質量造成不可忽略的影響. 因此，本文將點缺陷納入研究范圍，對單晶硅力學性能的影響因素加以討論和完善. 目前在點缺陷對金屬[7, 8]、化合物[9-11]和碳基材料[12, 13]的力學性能研究中，已經通過理論分析、實驗測試和分子動力學模擬等方式開展了較多高水平的工作. 作為單晶硅主要競爭對象的鑄造多晶硅，點缺陷對其力學性能的影響亦有過系統性研究[14]. 但點缺陷對單晶硅的力學性能的影響如何，未見相關報道.

另一方面，單晶硅因其具有的高理論比容量和清潔環保等優點，被廣泛應用在以硅基負極為首的新型電極材料上. 而硅在作為負極材料發生的嵌鋰過程中會發生巨大的體積膨脹現象，造成電池首次庫侖效率低和容量衰減等問題[15]. 因此，抑制硅負極電池充放電過程中的體積膨脹效應，探究單晶硅在循環應力作用下的力學性能變化情況，已經成為硅負極研究的關鍵問題之一. 為解決困擾硅負極電池的體積膨脹效應，Darbaniyan等人[16]研究了應變率和鋰質量分數對硅晶體的力學性能影響. Kim等人[17]以預鋰化的方式成功抑制了體積的變化，還有一些研究者以負極結構設計為著力點，提出了多種改善硅負極力學性能的優化方式. 例如，天津大學楊全紅教授研究組在硅負極的構建中提出了一種通過可控收縮結合碳籠網絡穩定碳/硅界面的力學緩沖增強策略[18]，有效地改善了嵌鋰過程中的體積膨脹效應. Xie等人[19]將納米硅顆粒包覆到富氮多孔炭基體中. 所制備的蜂窩狀多孔復合材料具有顯著的循環穩定性和優異的倍率性能，并且能夠很好地適應嵌鋰過程中的膨脹效應. 但尚未有研究者從原子級別的微觀尺度，以分子動力學手段討論在較大原子尺度的孔洞或點缺陷團簇作用下的單晶硅晶體力學性能的變化情況.

本文運用分子動力學方法，選擇硅的典型勢函數，在恒定應變率作用下模擬單晶硅的拉伸破壞行為，通過在晶體中構造不同原子大小的點缺陷，以此來分析單個點缺陷，點缺陷團簇乃至大型孔洞對單晶硅力學性能的影響. 同時利用數值模擬給出的數據驗證了米塞斯強度理論，討論了單晶硅中的應力集中效應與點缺陷之間的聯系，闡明了點缺陷對單晶硅力學性能的影響機理. 期望能為單晶硅力學性能的同類研究問題提供參考.

2 模型與方法

2.1 模型構造

單晶硅具有金剛石型晶體構型，晶格常數為0.543 nm. 本文采用LAMMPS軟件構建理想的單晶硅超晶胞模型，坐標軸X、Y、Z方向分別對應單晶硅的[100]、[010]、[001]晶向，計算區域取40×40×40個晶胞，包含512000個硅原子，模擬盒子幾何大小為217×217×217 nm3. 模型在三個方向上均運用周期性邊界條件，因為周期性邊界條件可以消除表面效應，使模擬更加符合實際情況[20]. 構建完成后的理想晶體如圖1所示. 體系建立后再在理想晶體的幾何中心位置劃分出不同半徑的球形區域，以刪除該區域內晶硅原子的方式構建不同大小的點缺陷. 點缺陷的大小則使用其生成時刪除的晶硅原子數度量.

圖1 單晶硅晶胞模型和6個含點缺陷的樣品

2.2 勢函數

分子動力學模擬是基于勢函數描繪的相互作用進行的，因此選擇一個能較好地描繪硅晶體性質的勢函數對于模擬的準確性而言至關重要. 硅的原子間勢函數發展至今已有三十多個版本[21]，經過一系列比較與參考后，本文決定使用Tersoff勢的第二個版本進行分子動力學模擬. Tersoff 勢是一種鍵級勢, 首次報道于1986年, 后分別于1988年和1989 年先后經過兩次修改, 形成了 T1[22], T2[23]和 T3[24]三個版本, 該勢函數除了能較好地描述金剛石結構外, 對硅的非正四面體構型也能描述, 如團簇、晶向和液態等結構. 其系統的總勢能為：

(1)

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

(2)

式中：fR與fA分別表示排斥作用函數和吸引作用函數，rij為原子i，j之間的距離，fc為光滑截斷函數，作用是使兩原子在距離較遠時的相互作用光滑地趨向于零，aij和bij分別表示排斥和吸引作用的系數. Zhou等人[25]對比了Tersoff，SW和MEAM等勢函數對硅晶體熔化特性的描述，結果表明Tersoff 勢相對更適合描述硅的熔化和凝固過程.

2.3 模擬方法

單晶硅拉伸過程采用恒定應變率加載方式，應變率加載方向為X軸上的[100]晶向，應變率為1×108/S-1. 模擬過程分為弛豫與加載兩個階段. 由于初始構造的原子體系是不平衡的，為使體系達到平衡狀態，將體系在NPT系綜(等溫等壓恒定原子數目)弛豫300 ps，使其溫度達到300 K. 弛豫完畢后開始加載應變：該過程采用NVE系綜，利用Berendsen[26]恒溫器控制溫度，Berendsen恒壓器控制Y軸與Z軸上的壓力趨于零以消除邊界壓力對模擬結果的影響. 本模型時間步長設置為1 fs；使用Verlet算法計算原子的運動軌跡，結果的可視化分析則通過Ovito[27]軟件完成.

表1 模擬參數設置

3 結果與討論

3.1 點缺陷對力學性能的影響現象

材料的力學性能受晶體缺陷的影響最為嚴重[28]，而硅晶體又易受加工時產生的熱應力等因素的影響而產生常見的點缺陷. 為分析點缺陷對單晶硅力學性能的影響現象，將所有單晶硅樣本沿X軸方向施以恒定應變率加載的應力—應變曲線繪制為圖2.

圖2 所有樣品的應力—應變曲線

理想晶體的應力-應變曲線表明，單晶硅在恒定應變率的拉伸過程包括彈性階段和屈服階段. 在彈性階段中應變不斷增加，應力也呈線性增加的趨勢，此階段中的變形為可恢復的彈性變形. 當應力超過屈服強度之后，晶體發生破壞，應力迅速下降至下屈服點，轉入屈服階段. 對比不同大小點缺陷的應力-應變曲線則表明，相較于不含點缺陷的理想晶體，點缺陷使單晶硅的屈服強度、極限應變顯著降低，下屈服點升高.

進一步討論點缺陷對屈服強度的影響時發現，點缺陷大小與屈服強度之間存在著函數關系. 圖3的擬合結果表明，點缺陷大小與屈服強度σs的變化服從公式(3)所述的指數函數關系.

圖3 屈服強度與點缺陷大小的指數擬合曲線

σs(c)=σ0+A×Exp(R×c)

(3)

式中，c為點缺陷的大小，σ0為屈服強度下限值，約為12.05±0.46 GPa，A=5.41±0.64，R=-0.0051±0.0015.A與R均是與點缺陷自身性質有關的參數. 至于上述兩個參數的影響因素，還有待更深入的研究和討論.

3.2 從應力變化的角度討論點缺陷的影響機理

Bullegas等人[29]研究了應力集中在復合材料破壞時的作用，從中歸納出內應力的聚集與釋放是材料破壞的主要原因. 而材料的破壞通常離不開應力的影響. 所以為更深入地討論單晶硅的破壞機理，就有必要先觀察和統計其在破壞前一刻的應力分布情況.

圖4和圖5利用0.5 ps內的應力張量數據將σx在晶體中的分布情況繪制為可視化圖像. 從中可見，在單個點缺陷和點缺陷團簇中，σx的分布大致是相同的. 圖5(a)說明應力的集中現象總是分布在垂直于外加應力方向的危險截面上. 在圖5(b)，(c)和(d)中，隨著應變的增加，與[100]方向成一定夾角的(111)晶面上的應力水平也逐步增加. 圖4和圖5已經從原子尺度上展示了σx的集中現象和點缺陷之間的聯系，但僅憑應力云圖還無法獲取達成破壞的所需的極限內應力范圍，故還需要對應力進行一些統計工作.

圖6與圖7定量地給出了所有單晶硅樣本中破壞前1 ps內的最大應力與點缺陷大小的變化關系，圖6中σx的變化情況表明，晶體內部的最大正應力是始終大于應力—應變圖像上的屈服強度的，兩者的變化趨勢相同，說明單晶硅的破壞可能取決于該點的σx，但不能說是絕對. 另外σy和σz雖然被控制為接近0，但隨著點缺陷大小的增大，卻有著上升的趨勢. 圖7中的剪應力τxy，τyz和τxz則隨著樣品中點缺陷大小的增大而增大. 這一變化情況則表明，含有點缺陷的單晶硅晶體在單向拉伸時其內部存在著復雜應力狀態，僅單純考慮σx的變化情況與屈服強度的關系過于片面，晶體的破壞應該是各向應力共同作用的結果. 為討論正應力和剪應力變化情況之間的聯系以及此類復雜應力狀態的決定因素，還需從強度理論著手；使用分子動力學模擬得到的應力張量數據，對晶體內每個原子的受力情況按各強度理論的公式進行校核，確定一種能夠較好地描述單晶硅屈服強度降低現象的強度理論，以供后續的分析與討論.

圖6 破壞前1 ps內所有晶體樣本內最大正應力與點缺陷大小的關系

圖7 破壞前1 ps內所有晶體樣本內最大剪應力與點缺陷大小的關系

圖8 按四類強度理論給出的理論值與單晶硅實際屈服強度的校核

(4)

式中，σs為米塞斯強度理論給出的理論屈服強度，σx，σy和σz分別指某點X，Y和Z軸上的正應力，τxy，τyz和τxz分別對應某點XY，YZ和XZ平面上的切應力. 當材料內部某一點的正應力和切應力滿足米塞斯等效應力公式給出的理論屈服強度時，該點即進入塑性狀態. 因為已經發現米塞斯強度理論可以很好地描述單晶硅的屈服強度的變化，所以為進一步討論點缺陷對單晶硅屈服強度的影響機理，就有必要對晶體內部的米塞斯應力水平進行分析.

圖9 樣本3中面在發生破壞時的裂紋和米塞斯應力的演化情況，(a) 時間步2042500的裂紋快照；(b) 時間步2044500的裂紋快照；(c) 時間步2046000的裂紋快照；(d) 時間步2047500的裂紋快照；(a1) 時間步2042500的米塞斯應力快照；(b1) 時間步2044500的米塞斯應力快照；(c1) 時間步2046000的米塞斯應力快照；(d1) 時間步2047500的米塞斯應力快照

在了解了點缺陷對單晶硅力學性能影響的一系列現象后，將點缺陷對屈服強度的影響機理概括為以下結論：點缺陷會改變晶體內部的應力分布，在其周圍區域產生應力集中效應，使這些區域的破壞條件符合米塞斯屈服準則. 因為應力集中區域的米塞斯應力總是最容易達到理論屈服強度，所以實際屈服強度會按米塞斯強度理論的預計趨勢變化. 在討論單晶硅的理論屈服極限時，可按米塞斯強度理論所給出的參考值來預測其力學性能.

3.3 破壞后的單晶硅微結構

3.2節中提到，小規模的破壞是在<111>晶面族上開始的，經過圖9所示的延展后在晶體中產生宏觀裂紋. 通過Ovito軟件觀察破壞后的結構時，發現晶體發生破壞的結構又表現出一定的規律性. 本文通過去除晶硅原子的方式，將破壞后的微結構幾何構造繪制成圖10.

圖10(a)與圖10(c)中的破壞微結構說明，于[100]方向施加恒定應變率時，裂紋在(100)面與(110)面上表現出明顯的周期性重復. 而從圖10(b)中的[110]晶向視角來看，發生破壞后的結構類似于一種非正交二維網格. 進一步觀察后發現該網格是沿著解理面分布的. 另一方面，圖10(d)中所給出的單晶胞破壞結構表明，二維網格在單個樣品中的基本組成單元是兩個夾角為70.53 °，且同屬于<111>晶面族的晶面. 上述現象說明破壞不完全由拉力引起，剪力也參與了破壞，側面驗證了米塞斯強度理論對單晶硅力學性能描述的準確性.

圖10 單晶硅晶體受[100]晶向拉力破壞后的微結構，(a) [100]晶向視角；(b) [110]晶向視角；(c) [010]晶向視角；(d) 單樣品中的網格結構

表2 破壞后結構的面分布情況

4 結 論

本文通過分子動力學模擬研究了點缺陷對單晶硅力學性能的影響和機理，結果歸納為以下幾點：

(1) 點缺陷對單晶硅彈性模量的影響很小，但會降低其屈服強度. 屈服強度σs與點缺陷大小c符合公式(3)中所述的指數函數關系.

(2) 通過對硅晶體中每個原子的應力進行計算和可視化處理，結果發現在點缺陷附近區域存在著應力集中現象，該現象說明點缺陷能夠改變其附近區域的應力狀態. 接著使用原子的六個應力張量對四大強度理論進行了校核，給出了理論屈服強度的變化情況. 其中，第四強度理論所給出的理論屈服強度最符合單晶硅實際屈服強度的變化情況，其描述的米塞斯應力的演化情況也能較好地解釋單晶硅裂紋的擴展現象. 所以本文認為點缺陷對單晶硅屈服強度的影響機理是：點缺陷通過改變晶體內部的應力分布，使這些區域的發生破壞的應力條件符合米塞斯屈服準則. 因為米塞斯應力的最大值總是在應力的集中區域出現，所以實際屈服強度會按第四強度理論的理論預計趨勢變化.

(3)通過直觀地觀察發生破壞后的微結構，發現破壞后的單晶硅結構類似于一種二維網格. 產生的微結構沿<111>解理晶面族分布，但在<111>晶面族中的4個可能發生破壞的晶面中，每次僅會有兩個夾角為70.53°的晶面發生破壞，通過統計各晶面上的應力推測該現象的產生原因是由于破壞總傾向于在應力滿足了破壞條件的晶面上產生.