基于分子動力學模擬的單晶硅沖擊壓縮相變研究*

劉夢婷，李旺輝，奉蘭西，張曉晴，姚小虎

（華南理工大學土木與交通學院, 廣東 廣州 510641）

晶體硅具有純度高、晶體取向良好等特點，在電子等領域應用廣泛，同時由于其復雜的結構相變機制，在相圖研究中也被廣泛關注。

在準靜態(tài)荷載下，單晶硅表現(xiàn)出顯著的多態(tài)性，目前已觀察到13 種相結構。Mujica 等指出，單晶硅在10～11 GPa 的載荷下發(fā)生了從立方金剛石（cubic diamond, CD）結構到β-白錫（β-tin）結構的相變，并通過第一性原理計算出存在中間暫穩(wěn)態(tài) Imma 相。而增加壓力時，晶體硅繼續(xù)發(fā)生向簡單六角形（simple hexagonal，SH）結構的轉變；再進一步增加壓力，變?yōu)镃mca 空間群結構，而在41 GPa 的高壓下，晶體硅又從Cmca 結構轉變?yōu)榱矫芏逊e結構（hexagonal close pack，HCP）。在動態(tài)荷載實驗中，Gilev 等研究了在沖擊壓縮下單晶硅的金屬化現(xiàn)象，發(fā)現(xiàn)在沖擊壓力達到23 GPa 時，晶體結構的缺陷導致電導率發(fā)生了偏離。Loveridge-Smith 等用X 射線衍射技術測量了沖擊加載下單晶硅的正交平面晶格參數(shù)，發(fā)現(xiàn)沿壓縮方向單晶硅的晶格間距減少了近11%，而與壓縮方向正交平面的晶格間距沒有顯著變化。Turneaure 等研究了單晶硅[100]、 [111]晶向在15.9、21.7 GPa 沖擊壓力下的動力學響應，發(fā)現(xiàn)其彈性波和非彈性波晶向依賴性很大，相變波晶向依賴性很小。Zhao 等對單晶硅[001]晶向進行了激光沖擊壓縮和回收實驗，觀察到兩個塑性響應：靠近沖擊表面的大塊非晶層和沿{111}平面產(chǎn)生的滑移帶，該結果揭示了單晶硅在一維應變條件下的塑性響應。Smith 等對單晶硅[001]、[110]和[111]晶向進行了沖擊實驗，揭示了單晶硅材料沖擊響應的各向異性，發(fā)現(xiàn)了其在沖擊壓力高于13 GPa 時轉變?yōu)楹唵瘟切谓Y構。Liu 等研究了單晶硅沿[100]晶向沖擊壓縮的塑性變形，發(fā)現(xiàn)其塑性變形能力隨著溫度增加而增加。Kishimura 等對在38 GPa 沖擊壓力下的單晶硅進行回收并進行了X 射線分析和拉曼光譜分析。XRD 分析顯示，單晶硅在11 GPa 沖擊壓力下有少量亞穩(wěn)態(tài)相產(chǎn)生，在38 GPa 沖擊壓力下有CuSi 這一物質生成，而在其余加載壓力下并沒有明顯的高壓相和亞穩(wěn)相存在。拉曼光譜的分析顯示，其峰值的整體偏移是由晶體尺寸減少引起的，該結果說明了單晶硅沖擊壓縮下的晶粒化過程。Turneaure 等用實時XRD 檢測了單晶硅在沖擊壓力低于54 GPa 下的晶體結構，發(fā)現(xiàn)單晶硅在沖擊壓力為31～33 GPa 時產(chǎn)生沖擊熔化，并從熔融邊界再結晶到密堆積六角形結構。由于硅和鍺屬于同族元素，具有相似的物理力學特性，因此Renganathan 等用原位XRD 檢測了鍺在[100]晶向的相變，發(fā)現(xiàn)鍺在沖擊壓力為15.7 GPa 時轉變?yōu)榘族a結構，在沖擊壓力為31.5 GPa 時轉變?yōu)槿廴趹B(tài)；應力卸載后，又轉變?yōu)榱⒎浇饎偸Y構，該結果表明鍺從白錫結構到立方金剛石結構的相變是可逆的。McBride 等用實時X 射線檢測技術揭示了多晶硅的沖擊相變行為，觀察到不同壓力下的多種固-固相變以及固-液相變。Paul 等則利用密度泛函理論、進化算法和晶格動力學等多種方法，構建了4 TPa 和6 000 K 下硅的相圖，為極高壓相圖的研究提供了參考。

盡管以上實驗與計算研究結果豐富，但實驗之間也存在較顯著的差異。一方面，實驗直接獲取的物理量過少，難以直接表征結果；另一方面，各個實驗設備和實驗過程也存在差異，變量如應變率、加載時間等不一致，而這些因素都會對單晶硅的相變有顯著影響。同時，單晶硅作為一種脆性材料，不僅加工困難，在沖擊壓縮時樣品也很難回收，而且單晶硅的相圖比較復雜，存在穩(wěn)態(tài)相和非穩(wěn)態(tài)中間相等多個相，因此它在沖擊壓縮下的物理和力學特性尚無完全統(tǒng)一的結論。此外，單晶硅沖擊實驗中加載時間較短，很難實時測量其相變過程中的晶格結構。所以，盡管檢測微觀結構變化的原位實時診斷技術取得了很大的進展，但對于可逆相變往往還是求助于把宏觀力學量的實驗測試結果與靜高壓數(shù)據(jù)進行比較，以確認新相的結構、相變的類型及性質，這也會導致研究結果的差異。

由于分子動力學具有超高的時間和空間分辨率，在微納米尺度研究材料的物理和力學特性時具有獨特的優(yōu)勢。而在當前的超算能力下，分子動力學模擬規(guī)模可以達到千萬甚至數(shù)十億原子，并能同時獲取各個原子運動的軌跡，進而從原子尺度到微納米尺度揭示材料的變形與力學參量的關系，使得其能夠對理論計算和實驗研究提供重要的補充。在分子動力學模擬單晶硅方面也有諸多發(fā)現(xiàn)。Swift 等用第一性原理計算得到單晶硅的物態(tài)方程，這個物態(tài)方程可以反應單晶硅從金剛石結構到體心四方結構的轉變。Demkowicz 等用SW 勢函數(shù)模擬非晶硅在施加恒定壓力下的變形情況，發(fā)現(xiàn)壓縮過程中存在類晶體和類液體部分，后者更容易發(fā)生塑性流動。Kumagai 等改進了Tersoff 勢函數(shù)，使之既能模擬彈性常數(shù)又能準確估測熔點。Higginbotham等通過分子動力學模擬，表明壓力誘導相變是彈性波異常的原因。Mogni 等則基于Tersoff-like 勢函數(shù)模擬[001]晶向單晶硅在激光沖擊下的相變，發(fā)現(xiàn)了新的Imma 相。Zhao 等研究了基于分子動力學模擬的沖擊壓縮下單晶硅的力學響應，發(fā)現(xiàn)硅在沖擊壓力低于10 GPa 時為彈性響應，在稍高壓力下，由于位錯和缺陷沿著特定的晶體學方向擴展，會進一步產(chǎn)生新的缺陷，此外在更高壓力下會產(chǎn)生交叉的非晶帶。

盡管基于分子動力學的單晶硅變形與相變的模擬研究成果豐富，但仍有若干關鍵問題尚未明確，包括單晶硅的動態(tài)壓縮各向異性，以及單晶硅復雜相變機制的物理描述和不同因素對其相變的影響規(guī)律和機制等。深入研究這些關鍵問題有利于更加全面地認清單晶硅在動態(tài)壓縮下的變形和相變的物理和力學特性。本文中將基于大規(guī)模分子動力學方法模擬研究單晶硅的沖擊壓縮響應，從微納米尺度認識晶體硅的沖擊壓縮特性，揭示其在動加載下的各向異性和結構相變等特征。

1 研究方法

1.1 模型的構建與沖擊壓縮模擬方法

基于廣泛使用的LAMMPS平臺，采用分子動力學方法模擬研究單晶硅的沖擊壓縮響應，主要對其非彈性變形和相變行為展開分析。晶體硅通常情況下具有立方金剛石結構，在300 K 溫度下晶格常數(shù)為0.543 nm。模型的建立是由單個晶胞沿、、晶向分別復制30、30 和390 個晶胞得到的超胞結構。因此，模型尺寸約為16.3 nm×16.3 nm×212 nm，見表1。

表1 單晶硅計算模型詳細參數(shù)Table 1 Parameters of single crystal Si sample for MD simulation

晶體硅中Si-Si 原子間的相互作用采用改進的半經(jīng)驗Tersoff 勢函數(shù)描述，并使用Erhart 等優(yōu)化的參數(shù)。與改進前的勢函數(shù)相比，該勢函數(shù)可以成功模擬晶體硅的多種性質，包括彈性常數(shù)和一般熱力學性質，以及在靜水壓縮和沖擊壓縮條件下立方金剛石相到β-白錫相的結構特征。

在施加沖擊壓縮前，先對構建的初始模型進行弛豫，弛豫過程中保持三維方向均為周期性邊界條件，使整個晶體模型在NPT 系綜下弛豫了20 ps。從能量和殘余應力的曲線上確認了材料模型為能量趨于穩(wěn)定且殘余應力在誤差允許范圍內趨于零的狀態(tài)，此時獲得了合理穩(wěn)定的初始結構。在沖擊壓縮加載階段，為了消除邊界對應力波的影響，在沿軸施加自由邊界條件的同時，沿、方向施加周期性邊界條件。模擬平板沖擊采用的是活塞法，活塞對自由表面產(chǎn)生沖擊波。單晶硅沖擊加載的示意圖如圖1所示。沖擊壓縮模擬過程中系統(tǒng)處于NVE 系綜，并且不施加溫度控制，時間步長為1 fs。這樣的加載方式和效果，與典型的平面沖擊實驗相似，均為一維應變問題。為揭示單晶硅沖擊響應的各向異性，尤其是沖擊結構的演化，分別對[001]、[110]和[111]晶向進行系列沖擊壓縮研究，沖擊粒子速度從0.3 km/s逐漸增加至3.2 km/s。

圖1 分子動力學模擬中單晶硅的沖擊壓縮示意圖Fig. 1 Schematic of shock compression of single crystal Si in molecular dynamic simulations

1.2 主要物理量的計算和可視化分析方法

為了深入分析單晶硅的沖擊壓縮響應，利用OVITO進行可視化分析，并采用配位數(shù)、中心對稱參數(shù)、鍵角分析等多種分析技術，對晶體硅的變形和結構變化進行識別和分析。

2 結果與討論

2.1 沖擊響應與各向異性

分子動力學模擬計算結果與不同文獻中實驗所得的晶體硅沖擊Hugoniot 應力和體積變化的關系如圖2 所示。從圖2 中可以看出，在體積壓縮較小時，所有結果高度吻合；壓縮較大時，實驗的應力-體積關系呈現(xiàn)兩大分歧：一是以Goto 等和Gust 等的結果為代表，這兩個獨立研究團隊的結果比較接近，二是以Turmeaure 等和McBride 等的實驗結果為代表。因此，對晶體硅不同沖擊實驗研究的差異值得關注，這可能與實驗手段、技術和材料樣品差異有關。一方面，Gust 等早期采用爆炸驅動的沖擊波壓縮實驗，而Goto 等和Turmeaure 等采用輕氣炮驅動的平板沖擊裝置對不同晶向的單晶硅進行沖擊實驗，McBride 等則采用激光驅動的應力脈沖加載對多晶硅進行沖擊壓縮實驗；另一方面，盡管所有實驗都是運用Rankine-Hugoniot 關系推算試樣的內部壓力，但壓力推算中所需的粒子速度或相變波速度測定的差異顯著，例如Gust 等提及其結果與相關歷史數(shù)據(jù)的對比時，對應于HEL 和相變點的臨界體積（或密度）是很吻合的，但壓力數(shù)值卻存在差異，他們將此歸因于粒子速度的測定差異。值得注意的是，兩大存在分歧的壓力-體積關系均有相互獨立的研究團隊的結果支撐，可能預示著深層次的原因在于速度測定的技術差異與硅材料響應有關。Goto 等的輕氣炮驅動的平面沖擊實驗中采用了與Gust 等相似的測量技術（斜鏡與條紋相機）對硅試件的自由面測速，其Hugoniot 壓力推算值也與Gust 等非常吻合；然而，Turneaure 等的輕氣炮驅動平面沖擊實驗與McBride 等的激光沖擊實驗則采用了VISAR 測試系統(tǒng)對硅/LiF 窗口或者硅自由面速度測定，其結果相接近。此外，Smith 等指出硅自由表面端存在的應力松弛和初始的速度濺射下的自由表面測量結果可能會影響應力和體積的關系測定。因此，這些結果在壓力數(shù)值上的差異深刻反映了實驗技術差異和數(shù)據(jù)處理方法對實驗結果的重要影響。

圖2 不同實驗[13,20,32-33]與模擬中晶體硅的沖擊壓縮Hugoniot 應力和體積變化的關系Fig. 2 Shock Hugoniot stress as a function of volume change in various experiments[13,20,32-33] and simulations of Si crystals

與實驗結果相比，本文計算結果與Goto 等和Gust 等的結果更接近，但也在一定程度上高于實驗值。盡管如此，基于該原子間勢函數(shù)模擬計算的應力-體積整體變化趨勢與實驗結果相似。例如，隨著體積壓縮的增加，可以明顯看出在體積比低于0.8 時單晶硅已經(jīng)處于非彈性階段；而進一步壓縮至體積比為約0.67 時，材料壓力急劇增加顯示出典型的高壓相變特征。因此，盡管該原子間勢函數(shù)對單晶硅的Hugoniot 壓力存在高估的不足，但其可以成功反映晶體硅在沖擊壓縮下的若干重要變形響應特征，對本文研究動高壓下硅的變形和相變仍有價值。

在不同沖擊壓縮下，單晶硅的剪切應力隨粒子速度變化如圖3 所示。在[001]、[110]和[111]晶向中，其剪切應力在一定范圍內均隨著粒子速度的增加而逐漸增加，達到一定臨界粒子速度后開始急劇下降。但是，該臨界粒子速度在不同晶向差異顯著，在[001]晶向最小，為1.4 km/s，在[110]晶向最大，為1.8 km/s，而在[111] 晶向時介于兩者之間，為1.5 km/s。在這3 個晶向中，最大剪切應力分別為5.0、7.8 和6.2 GPa。剪切應力的大幅度下降意味著單晶硅在這些沖擊加載下發(fā)生了重要的變形（如塑性或相變等），使得剪切應力釋放。可以觀察到 [001]和[111]晶向中剪切應力經(jīng)歷了由正值轉為負值的階段：分別在1.8、2.4 km/s 達到最大負剪切應力，而后逐漸減小。負剪切應力現(xiàn)象存在于單晶碳化硅和單晶金屬鐵中，學者們普遍認為這是由于材料過度弛豫釋放的結果，并且具有晶向依賴性，多發(fā)生在[001]晶向。從剪切應力的計算方法角度分析，負剪應力意味著在沖擊加載中材料內部的橫向應力高于縱向應力。最終，[001]、[110]和[111]晶向的剪切應力分別在2.0、1.9、2.6 km/s 之后趨近零，因此材料失去剪切強度。此外，結合圖1 中應力和體積壓縮的關系，單晶硅在[110]晶向的奇異表現(xiàn)還體現(xiàn)在彈塑性或相變轉變的臨界狀態(tài)，相比于[001]和[111]晶向的在約0.85 體積壓縮比時的“彈塑性”轉變，[110]晶向則表現(xiàn)出奇異的彈性特征。這種沖擊壓縮下奇異的彈性響應現(xiàn)象在金剛石等材料的分子動力學模擬研究也有過發(fā)現(xiàn)，這些材料呈現(xiàn)出非單調的剪應力-應變關系。由于剪切應力是塑性等變形的驅動力，因此，單晶硅沖擊壓縮下奇異的彈性響應應該與其剪切應力的非單調性相關。

圖3 沖擊Hugoniot 狀態(tài)的剪切應力(τ)與粒子速度(up) 曲線Fig. 3 Shock Hugoniot shear stress as a function of particle velocity

為了進一步探究單晶硅沖擊加載下的材料特性及其各向異性，圖4 中給出了沿[001]、[110] 和[111]晶向加載時單晶硅的縱向應力、剪切應力和密度波剖面。根據(jù)單晶硅的沖擊Hugoniot 曲線，單晶硅3 個晶向的彈塑性或相變閾值是不同的，[110]是3 個晶向相變閾值最大的，閾值粒子速度為1.8 km/s，其次是[111]晶向，閾值粒子速度為1.5 km/s，[001]晶向最小，為1.4 km/s。因此，波剖面分析中選取了各晶向中具有代表性的3 個粒子速度狀態(tài)。如圖4(a)所示，在[001]晶向中，粒子速度分別為1.3、1.4、1.5 km/s的波剖面體現(xiàn)了單晶硅從單一彈性波到彈性-塑性/相變波的沖擊波結構轉變過程。在縱向應力波剖面中，當粒子速度為1.3 km/s 時，波剖面平整，且波陣面十分狹窄陡峭。當粒子速度增加至1.4 km/s 時，彈性前驅波幅值增加，但同時材料出現(xiàn)軟化現(xiàn)象，導致彈性波之后出現(xiàn)應力幅值下降的現(xiàn)象，同時由于塑性或相變的激發(fā)，在靠近活塞端的波剖面區(qū)域呈現(xiàn)比軟化區(qū)域更高幅值的塑性或相變波。當粒子速度為1.5 km/s 時，類似的現(xiàn)象依然存在，但由于沖擊強度增加，彈性前驅波幅值增加，這可能是由于單晶硅沖擊Hugoniot 彈性極限對加載速率敏感性很強導致的；同時，更高的沖擊強度使單晶硅的軟化越發(fā)明顯，使得彈性前驅波之后的波幅值下降更快；同樣，塑性或相變波部分幅值增加，對應的波速度也在增加。從[001]晶向的剪切應力波剖面可以進一步揭示該現(xiàn)象的變化。剪切應力對材料的響應非常靈敏，單晶硅在沖擊壓縮下呈現(xiàn)復雜的材料響應，包括彈性波及其軟化和塑性或相變波。剪切應力隨著彈性波的出現(xiàn)急劇上升，然后隨著軟化作用而略有降低，并在塑性或相變波區(qū)域迅速下降，更強的塑性或相變波引起更大幅值的剪切應力釋放。密度波剖面的變化對應了包括剪切應力波剖面在內的各物理量波剖面情況，即密度隨著彈性波抵達而快速上升，同時由于軟化作用略有回落，并在相變波區(qū)域由于材料的進一步壓縮而顯著上升。本文中研究獲得的單晶硅的波剖面特征在Higginbotham 等和Stubley等的研究工作也得到了體現(xiàn)。

圖4 單晶硅中分別沿[001]、[110]和[111]晶向的沖擊波剖面Fig. 4 Shock profiles in single crystal Si along [001], [110] and [111] crystal orientation

在[110]和[111]晶向中，復雜的多波結構依舊存在，但軟化作用更加微弱。同時，與[001]晶向不同的是，在[110]和[111]晶向中波剖面彈性波幅值幾乎保持不變，顯示出與加載速率無關，塑性或相變波更加復雜。

基于波剖面的分析，采用可視化分析技術進一步探究沖擊引起的單晶硅非彈性變形和相變等材料響應。圖5 中分別給出了沿著[001]、[110]和[111]晶向沖擊壓縮時的可視化結果，粒子速度分別為1.4、1.8、1.5 km/s。結合上文分析可知，這些粒子速度沖擊壓縮下單晶硅3 個晶向中的剪切應力均開始產(chǎn)生急劇的下降，意味著單晶硅在這些加載區(qū)域發(fā)生了重要變形響應。中心對稱參數(shù)和配位數(shù)是反映晶體結構特征的重要指標，可以體現(xiàn)單晶硅的發(fā)生與結構相關的變化。圖5(a)中單晶硅根據(jù)中心對稱參數(shù)分析著色，而圖5(b)則根據(jù)原子最小配位數(shù)著色。在中心對稱參數(shù)分析中，3 個晶向的單晶硅均發(fā)生了顯著的中心對稱參數(shù)變化，但呈現(xiàn)出不同的帶狀紋理。未干擾區(qū)域具有比壓縮區(qū)域更高的中心對稱參數(shù)，而在帶狀條紋區(qū)域該參數(shù)則更低，并且在交叉區(qū)域該參數(shù)呈現(xiàn)出趨于零的變化。與此同時，在原子最小配位數(shù)分析中也可以發(fā)現(xiàn)3 個晶向的差異顯著。[001]單晶硅存在配位數(shù)為5 和6 的帶狀和交叉區(qū)域；而[110]單晶硅中在靠近加載端一段區(qū)域幾乎全部變成具有更大配位數(shù)的原子；在[111]單晶硅中則僅有極少數(shù)原子發(fā)生最小配位數(shù)的變化。這說明在3 個晶向中單晶硅的晶體結構均發(fā)生了結構性的變化，但是呈現(xiàn)明顯的晶體各向異性。從不同的配位數(shù)和中心對稱參數(shù)看，單晶硅沖擊壓縮下發(fā)生了多種復雜結構相變，這體現(xiàn)在某些結構相變可引起中心參數(shù)的變化但不會引起配位數(shù)的變化。

圖5 單晶硅沿[001]、[110]、[111]晶向分別在1.4、1.8、1.5 km/s 沖擊壓縮下的可視化Fig. 5 Visualizations of single crystal Si in [001], [110] and [111] crystal orientations with shock particle velocities of 1.4, 1.8, 1.5 km/s

2.2 沿[001]晶向沖擊的變形與相變

為進一步探究單晶硅的沖擊結構相變等特征，重點分析[001]晶向單晶硅的情況。圖6 中給出了不同粒子速度下的單晶硅沖擊響應可視化分析，通過晶格原子的偏轉角可以顯著地區(qū)分出沖擊壓縮下單晶硅變形和相變的局部變化，藍色區(qū)域原子的偏轉角非常小，黃色區(qū)域原子較初始結構偏轉約54°，綠色區(qū)域原子較初始結構偏轉45°。圖中[001]晶向的相變帶狀區(qū)域延伸方向與加載方向呈57°夾角，與偏轉角接近，這樣的少許差異與材料處于壓縮狀態(tài)相關。這表明相變波的相變帶遵循特定的晶體學取向，是沿著靠近{111}晶面向前傳播，其中{111}晶面與{001}晶面夾角為54°。隨著粒子速度的增加，相變帶狀區(qū)域寬度有略微增加的趨勢。當粒子速度為2.0 km/s 時，靠近沖擊加載端的區(qū)域出現(xiàn)無序的狀態(tài)，這個無序區(qū)域隨著粒子速度的增加而擴展。當粒子速度為2.4 km/s 時，一個顯著的無序狀態(tài)與晶體狀態(tài)同時存在于相同加載方向位置。考慮到強沖擊可能引起劇烈溫升和更高的壓力狀態(tài)，因此可能引起單晶硅的熔化，這意味著在此沖擊條件下單晶硅處于固-液共存的區(qū)域。值得注意的是，本文中所發(fā)現(xiàn)的沖擊波陣面后方固-固相變和固-液相變共存的現(xiàn)象也出現(xiàn)在其他兩個晶向，該發(fā)現(xiàn)與McBride 等的激光沖擊實驗診斷結果一致，本文中沖擊壓縮下固-液共存甚至完全熔化所對應的臨界體積壓縮比，也與McBride 等的實驗高度吻合。如圖6 所示，當粒子速度高達2.8 km/s 時熔化區(qū)域已經(jīng)擴展很大，并在3.2 km/s 時波陣面后幾乎完全熔化。

圖6 不同沖擊粒子速度下[001]單晶硅的結構演化Fig. 6 Structural evolutions of single crystal Si in [001]orientation at various shock particle velocities

基于圖6 的代表性局部區(qū)域，進一步通過鍵角分析和晶向分布函數(shù)分析沖擊壓縮下單晶硅的結構相變。圖7(a)為幾個代表性的局部區(qū)域的鍵角分析。區(qū)域 Ⅰ代表沖擊波陣面處于帶狀區(qū)域之前的部分，區(qū)域Ⅱ代表黃色帶狀區(qū)域，區(qū)域Ⅲ代表綠色的交叉區(qū)域，區(qū)域Ⅳ代表中間合圍區(qū)域，區(qū)域Ⅴ代表無定型區(qū)域。初始結構中，鍵角呈現(xiàn)109°，區(qū)域Ⅰ主要鍵角峰值在105°，區(qū)域Ⅱ和Ⅲ則依次主要呈現(xiàn)96.5°和接近90°，而區(qū)域Ⅳ則與初始結構非常相似。這說明除了區(qū)域Ⅳ，其余區(qū)域均出現(xiàn)了與結構變化相關的鍵角變化。圖7(b)的徑向分布函數(shù)進一步揭示了可能的結構變化。從圖7 可知，與初始結構相比，相變區(qū)域Ⅱ～Ⅲ、Ⅴ變化明顯。其中區(qū)域Ⅱ～Ⅲ中的峰值位置和數(shù)量與初始結構差異顯著，體現(xiàn)了晶體結構的固-固轉變，而區(qū)域Ⅴ中峰值之后趨于平滑的變化則說明硅材料向無定型狀態(tài)的變化。

圖7 對代表性局部區(qū)域的鍵角分析和徑向分布函數(shù)分析Fig. 7 Bond-angle and radial distribution functions analysis for representative local regions

2.3 沿[110]和[111]晶向沖擊的變形與相變

利用相似的后處理技術，對單晶硅沿[110]和[111]晶向沖擊壓縮的變形和相變進行分析。圖8～9分別為單晶硅[110]晶向中不同粒子速度下的結構演化和代表性局部的鍵角統(tǒng)計和徑向分布函數(shù)分析。當粒子速度為1.8 km/s 時，在靠近左端一定范圍內出現(xiàn)了兩種變化：垂直于加載方向呈帶狀條紋區(qū)域的晶格原子有很小的晶向偏轉，緊隨其后的是一定程度的原子無序狀態(tài)區(qū)域。帶狀區(qū)域中區(qū)域Ⅰ的鍵角分布顯示了兩個主峰值，分別對應約90°和120°的鍵角，而區(qū)域Ⅱ的主峰值仍為109°。區(qū)域Ⅲ的徑向分布函數(shù)在第1 個峰值后趨于平緩且鍵角分布沒有獨立峰值，該區(qū)域為無序狀態(tài)，這可能是沖擊引起的熔化。對比圖4 的波剖面可知，區(qū)域Ⅰ和區(qū)域Ⅲ這兩種變化與波陣面后方的相變波對應，而且是分離的兩列波，對應著材料內部不同程度的剪應力釋放和密度上升。隨著沖擊粒子速度的增加，兩列波波速均在增加。并且當粒子速度高于2.4 km/s 時，固-液相變區(qū)域逐漸占據(jù)主導。

圖8 不同沖擊粒子速度下[110]單晶硅的結構演化Fig. 8 Structural evolutions of single crystal Si in [110]orientation at various shock particle velocities

圖9 對圖8 代表性局部區(qū)域的鍵角分析和徑向分布函數(shù)分析Fig. 9 Bond-angle and radial distribution functions analysis for representative local regions

圖10～11 為單晶硅[111]晶向中不同沖擊粒子速度下的結構演化和代表性局部的鍵角統(tǒng)計與徑向分布函數(shù)分析。當粒子速度為1.6 km/s 時，在波陣面后一定區(qū)域內也出現(xiàn)了帶狀條紋以及靠近端部的無序狀態(tài)區(qū)域。與[001]和[110]晶向不同的是，[111]晶向中帶狀條紋出現(xiàn)后，隨著沖擊粒子速度的增加，條紋寬度逐漸減小，并在粒子速度為2.4 km/s 時逐漸消失。無序狀態(tài)區(qū)域隨著粒子速度的增加也不斷擴展，并在3.2 km/s 時擴展速度急劇增加。這一現(xiàn)象可能是沖擊加載下溫升越發(fā)顯著引起的軟化作用和壓力引起的結構變化共同作用的結果。在相對低的沖擊粒子速度下軟化作用相對較弱，而壓力誘發(fā)的晶格原子偏轉更顯著，當沖擊強度增加時，溫升會更顯著，軟化作用影響增加，晶格偏轉受到一定的削弱。類似于壓縮下的變形孿晶，高溫會抑制變形孿晶的形成。由于硅晶體材料的變形和相變對溫度、壓力、應變率等因素比較敏感，因此其動態(tài)物理和力學行為會受到這些競爭性作用因素的影響。當粒子速度高達3.2 km/s 時，固-液相變波波速急劇增大，這是軟化區(qū)域大面積激活轉化為熔化狀態(tài)導致的。

圖10 不同沖擊粒子速度下[111]單晶硅的結構演化Fig. 10 Structural evolutions of single crystal Si in [111]orientation at various shock particle velocities

圖11 對代表性局部區(qū)域的鍵角分析和徑向分布函數(shù)分析Fig. 11 Bond-angle and radial distribution functions analysis for representative local regions

2.4 討論

從早期實驗與模擬研究以及本文的模擬工作中可以看到晶體硅的復雜相變。在單軸加載條件下，沿著[001]、[110]、[111]晶向的加載過程都會發(fā)生相變，但是相變特征和機制具有各向異性。實驗中已探明晶體硅在壓力逐漸增加作用下發(fā)生一些列相變，如cd 相→β-tin 相→Imma 相→sh 相等。其中cd 相是四配位結構，sh 相為六配位結構。β-tin 相是介于四配位和六配位之間的一種體心四角結構，而Imma 相又是介于β-tin 相和sh 相之間的正交相。因此，本文模擬中涉及的固-固相變可能包括β-tin 相和Imma 相。Mogni 等基于相同的原子間勢函數(shù)的模擬分析指出，[001]晶向單晶硅的沖擊壓縮固-固相變結構是Imma 相，而且他們認為僅[001]晶向中帶狀區(qū)域為Imma 相。通過本文中對[001]單晶硅沖擊壓縮下的交叉區(qū)域的偏轉角和配位數(shù)以及鍵角的分析，可以看出交叉區(qū)域的結構應該與Imma 相有所不同。該區(qū)域的配位數(shù)顯示可達5 甚至6，與sh 相具有相同配位，但本文中尚無足夠完整的結構信息用于確定具體結構。從圖5 中可以發(fā)現(xiàn)，單晶硅[001]、[110]和[111]晶向具有明顯各向異性，如沿[001]晶向加載時配位數(shù)發(fā)生變化，但在[110]和[111]晶向時配位數(shù)沒有變化，該現(xiàn)象是由于不同加載方式的原子堆積方式不同導致的。此外，不同晶向的變形與相變受到滑移系的影響，位錯阻力小的滑移系更容易滑移。對于單晶硅這種金剛石結構，密排面為(111)晶面，密排方向為<110>方向，所以位錯更容易沿著滑移系{111}<110>產(chǎn)生。所以相較于[001]晶向，[110]晶向偏轉角非常小，發(fā)生結構變化的平面幾乎平行于(110)平面。

晶體硅的沖擊壓縮研究揭示了不同的物理和力學行為，但不同實驗之間的確存在較大差異。在已經(jīng)報道出的晶體硅變形中除了熔化還包括位錯滑移、無定型化等。McBride 等指出非原位探測很難識別出單晶硅沖擊壓縮階段的熔化，一些中間相變也很難通過回收試樣的方法準確觀察到，因為晶體硅在高壓卸載時可能產(chǎn)生其他相結構。對分子動力學模擬工作，盡管本文定性地描述了單晶硅的沖擊壓縮變形和相變，觀測到彈性、非彈性、固-固相變以及固-液相變和固-液共存等物理特征，對比了晶體的各向異性，但是也應指出，本文所采用的勢函數(shù)仍然不能準確包含單晶硅的全部變形和相變特征。因此，應對不同原子間勢函數(shù)進行系統(tǒng)評估和對相關勢函數(shù)進行修改。

總之，本文已經(jīng)初步明確單晶硅沖擊壓縮下可依次呈現(xiàn)彈性、非彈性到固-固相變和固-液相變的變化趨勢及其各向異性特征，對進一步理解單晶硅的沖擊壓縮變形和相變提供了一定的支持。針對動態(tài)加載下單晶硅相變的復雜性，后續(xù)仍需對單晶硅在動態(tài)條件下的變形和相變研究提供更深入全面的納觀尺度信息。

3 結 論

采用大規(guī)模分子動力學方法，模擬研究了單晶硅的沖擊壓縮響應，獲得了其沖擊壓力和剪切應力的變化，分析了其沖擊波傳播特性各向異性，并重點研究了不同沖擊強度下單晶硅變形的非彈性轉變和相變行為，得到主要研究結論如下：

（1）單晶硅在沖擊壓縮下呈現(xiàn)彈性到相變的變化趨勢，其剪切應力隨著沖擊粒子速度的增加先增加，達到臨界點后由于結構的變化導致應力急劇下降，其剪應力和相應的臨界點存在明顯的各向異性；

（2）基于偏轉角和鍵角以及徑向分布函數(shù)的分析發(fā)現(xiàn)，沿[001]晶向沖擊壓縮下可發(fā)生多種固-固相變以及固-液相變，并觀察到與McBride 等的實驗高度一致的固-液共存的重要現(xiàn)象。

本文研究工作可為動加載下晶體硅的相變研究提供新的納米尺度的結果支撐，為相關實驗提供新的補充。