張學陽 胡望宇 戴雄英
1) (湖南工程學院計算科學與電子學院,湘潭 411104)
2) (湖南大學材料科學與工程學院,長沙 410082)
由于鐵在人類社會中,特別是在國防和工業(yè)領域扮演著重要的角色,因此,鐵在動態(tài)高壓下的行為一直是研究的熱點問題之一,其中α → ε 馬氏體相變尤其受到人們的關注,首次實驗研究早在1956 年就被發(fā)表了[1].之后,人們開展了一系列關于鐵的沖擊加載的研究,特別是動態(tài)加載下鐵的塑性活動和相變的實驗研究[2–4],引發(fā)了關于相變和塑性相互影響的諸多討論[5–11].但是,由于動態(tài)高壓實驗的瞬時性,研究人員應很難得到動態(tài)細節(jié),計算機模擬則補充了實驗研究所欠缺的細節(jié),逐漸成為了主流的研究方法之一.其中,分子動力學(MD)模擬基于原子間相互作用勢直接模擬材料的結構演化,不需要構建復雜的物理或力學模型,可以提供沖擊下微觀結構的動態(tài)變化細節(jié).但是,長久以來沒有合適的勢能夠在模擬中重現(xiàn)出鐵在沖擊加載下的塑性過程,因此動態(tài)加載下鐵的塑性和相變的耦合一直是MD 研究的盲點.直到最近幾年,一些新開發(fā)的勢克服了這一缺陷,特別是Wang等[12]開發(fā)的MAEAM 勢不僅成功模擬相變前的塑性,還克服了相變產(chǎn)物中有很多面心立方相變后的(FCC)原子的缺點(實驗中已經(jīng)證明FCC 結構不會大量出現(xiàn)).這些為研究沖擊下鐵的相變與塑性耦合的微觀細節(jié)奠定了基礎.
很早就有研究發(fā)現(xiàn),沖擊加載下的相變與塑性都喜歡在缺陷附近形核[13],之后人們陸續(xù)開展了很多關于缺陷對相變與塑性影響的研究[14–17].晶界作為自然界中最常見的缺陷,其對相變與塑性的影響自然也受到了人們極大的關注.例如Zong 等[18]研究了鈦的晶界在沖擊下對相變的影響,并且提出共格孿晶界對鈦的α → ω 相變是有利的.對于鐵的MD 研究來說,前期人們因為缺乏合適的勢函數(shù),只能把研究的重點都集中在相變上,后來Gunkelmann 等[19–21]和Wang 等[22–24]在分別開發(fā)了適合動態(tài)高壓下鐵的勢函數(shù)之后,都在第一時間研究了沖擊加載下鐵的相變與塑性的耦合問題.Gunkelmann 等[25]指出在晶界附近位錯先于相變形核,且可以在模擬中觀察到明顯的三波結構,同時測得的相變壓力閾值高于靜態(tài)下.而Wang 等[24]發(fā)現(xiàn)晶界能夠通過影響相變與塑性的耦合從而改變相變的機制,并提出新的機制與實驗中觀察到的缺陷誘導機制的特點相符,人們對多晶內相變與塑性耦合的認識有了很大的進步.但是,模擬中多晶的響應會受限于晶粒的尺寸,使得對晶界附近相變與塑性的認識還不夠清晰,而雙晶模型可以突破這個限制.因此,我們在早些時候就采用了雙晶模型來研究晶界如何影響相變以及塑性,總結了典型晶界對相變形核的重要影響[26,27].但是晶界的種類十分豐富,除了經(jīng)常被研究的典型晶界,還存在很多的研究幾乎沒有涉及到的普通晶界,而這些晶界附近的沖擊響應是否會符合已有的認知目前是不清楚的,因此關于晶界對相變與塑性影響的認知還不完善,有待于進一步研究.
在以往的關于材料的沖擊加載MD 研究中,活塞法等產(chǎn)生的沖擊波普遍被近似為一維沖擊波,沖擊方向以外的應力分量對沖擊響應的影響被完全忽視.此外,MD 研究中討論的常見沖擊晶向都是中心對稱的,例如〈001〉,〈110〉和〈111〉對于體心立方結構(BCC)和面心立方結構(FCC),〈0001〉,〈1010〉和〈1210〉對于六方密排結構(HCP)[28–30],而沿非中心對稱晶向沖擊的研究是罕見的.那么,非常見晶向在動態(tài)沖擊加載中的表現(xiàn)是否一直都是可被忽略的呢? Zhakhovsky 等[29]的研究指出沖擊波的結構與晶體的晶向有密切的聯(lián)系,因此關注沿非中心對稱的晶向沖擊的過程其實也是有必要的.對此我們研究了沿非中心對稱晶向沖擊單晶鐵的過程[31],發(fā)現(xiàn)選取非中心對稱晶向為沖擊方向時,沖擊波不能再被簡單視為一維的,晶格的各向異性在沖擊過程中有明顯的影響,甚至會影響相變路徑的選擇.考慮到多晶中晶向的隨機性和多樣性,顯然晶向的各向異性在多晶的沖擊加載中將會產(chǎn)生更重大的影響,只不過沒有引起人們足夠的重視而已.此外,在設計雙晶模型研究不同晶界對相變的影響的時候,發(fā)現(xiàn)某些模型中的晶界附近的行為與我們預期的完全不同,聯(lián)系到上述我們關于沿非中心對稱晶向沖擊的研究,晶格的各向異性很有可能是導致其反常的重要因素.本文的重點就是基于上述現(xiàn)象討論沖擊加載下鐵的各向異性如何影響雙晶模型中的晶界附近的相變.
綜合上述,雖然已經(jīng)對鐵的沖擊響應過程做出了一定的探索,但是還有很多問題有待于人們去揭示,例如晶格的各向異性在沖擊加載中的影響就沒有引起足夠的重視,關于沖擊加載下晶向的各向異性的研究還有待進一步開展.基于此,本文主要研究了沖擊下包含對稱晶界的雙晶模型中非常規(guī)晶向對晶界附近相變的影響.由于多晶取向的隨機性和合金中各種原子對中心對稱性的破壞,本文能夠為多晶和合金的沖擊實驗提供一定的理論支持.
使用Lammps 軟件[32]模擬了包含對稱晶界的雙晶鐵模型在沖擊加載下的非平衡分子動力學過程,之前的研究已經(jīng)證明MAEAM 勢[22]在研究沖擊下的鐵的相變與塑性的問題上有很大的優(yōu)勢.所構建的模型如圖1 所示,模型1 沿X,Y和Z方向的尺寸分別為19.74 nm,19.90 nm 和162.81 nm,晶體取向分別為〈100〉,〈031〉和〈013〉.模型2 沿X,Y和Z方向的尺寸分別為20.23 nm,20.13 nm 和151.45 nm,晶體取向分別為〈110〉,〈332〉和〈113〉.這兩個模型所包含的原子數(shù)都超過了500 萬,且測試結果表明模型的尺寸是足夠的,不會有邊界效應的影響.模型弛豫的方法如下: 首先,使用共軛梯度法使模型的能量處于最小化;然后,將樣本的初始條件設置為等溫等壓系綜(NPT),樣品的溫度設置為10 K,較低的初始溫度可以消除噪聲的影響,使得圖像更加清晰,而且經(jīng)過測試發(fā)現(xiàn)樣品初始溫度為室溫的模擬結果與10 K 的差異很小,表明設置較低的溫度不會影響本文的結論,對于X,Y和Z的邊界采用周期性邊界條件,將模型置于上述條件下弛豫20 ps;最后,將弛豫好的模型Z軸方向的邊界條件改為非周期性邊界條件,并且將系綜改變?yōu)槲⒄齽t系綜(NVE)后,進入沖擊加載階段.弛豫階段的時間步長為1 fs,沖擊階段的時間步長為0.5 ps.沖擊波由活塞法產(chǎn)生,并沿Z軸從左向右掃過模型,還對比模擬了沖擊波在同樣的模型中沿Z軸從右向左掃過的過程,排除了兩側微觀結構不嚴格對稱造成的影響.本文使用的可視化工具為Ovito[33],依據(jù)維里應力計算方法計算應力[34].通過將模型沿Z軸切片來統(tǒng)計模型速度和應力張量等[35].本文中的c軸為垂直于HCP 晶胞的正六邊形底面的向量,模型中大量相變產(chǎn)物的c軸需要通過我們自主開發(fā)的微結構分析程序來確認,其簡要思路是先通過Ovito 軟件篩選出屬于HCP 結構的原子并保留,循環(huán)選取其中的每個原子作為中心原子,通過原子間的距離找出中心原子的12 個最近鄰原子,找出12 個最近鄰原子中與中心原子處于同一個平面的6 個原子,并求解該平面的法線向量,即為c軸.

圖1 模型示意圖 (a) 包含sigma5 晶界的模型1;(b) 包含sigma11 晶界的模型2Fig.1.Model schematic figures: (a) Model 1 containing sigma5 grain boundaries;(b) model 2 containing sigma11 grain boundaries.
在關于沖擊加載下晶界附近相變的研究中[26,27],對稱晶界兩側的響應幾乎都是空間對稱的,例如在沿垂直于晶界的方向沖擊時,對稱sigma5〈100〉扭轉晶界、sigma〈110〉傾側晶界、sigma〈110〉扭轉晶界以及共格孿晶晶界兩側的相變產(chǎn)物和相變路徑都是相同的,盡管由于反射波的影響,晶界左側的響應稍微“弱”于右側.上述沖擊響應滿足空間對稱性與晶界兩側微觀結構的對稱保持一致(如圖1所示),然而發(fā)現(xiàn)有些對稱晶界的沖擊響應是不對稱的,如圖2 所示.圖2(a)是包含sigma5 晶界的模型1 在沖擊加載下某個時刻的微觀結構圖,代表HCP 結構的紅色原子是相變發(fā)生后的主要產(chǎn)物,該晶界兩側在受到?jīng)_擊后的響應有著明顯的差異,即右側有較多的HCP 原子,而左側觀察不到HCP原子的存在.而且隨著模擬時間的增加,左側也不會發(fā)生相變,我們還通過模擬更長模型的沖擊來確認了這一點.該結果表明,晶界兩側的相變閾值存在一定差異.晶界左側因為有反射波的影響,因此左側的沖擊響應比右側要“弱”一些是可以預期的,但是反射波引起的差異應該是十分有限的,模型1中晶界兩側相變閾值的差異不能用反射波來解釋.圖2 還展示了模型1 同時刻的剪切應力和溫度的分布,可以看出,沖擊波到達的所有區(qū)域中,波前左側的區(qū)域1 是剪切應力最高的,對應這一區(qū)域的結構正在發(fā)生彈性形變.區(qū)域1 左側是區(qū)域2,與區(qū)域1 相比剪切應力較低,對應該區(qū)域的結構發(fā)生了塑性形變,塑性過程會釋放剪切應力.區(qū)域3 對應著相變區(qū)域,發(fā)生了相變的區(qū)域剪切應力達到最低,因為相變過程中有大量原子發(fā)生滑移,原子間的滑移能夠進一步釋放剪切應力.根據(jù)上述剪切應力與形變過程的對應關系可以判斷,區(qū)域4 是沒有發(fā)生相變的,這與圖2(a)呈現(xiàn)的微觀結構是相符的.根據(jù)圖2(c)溫度分布可知,塑性行為是引起溫度上升的主要原因,溫度最高的是上述各區(qū)域的邊界以及晶界,對應沖擊加載下這一部分原子擁有最高的活躍度.

圖2 沖擊加載下含sigma5 晶界的模型1 在24 ps 時刻的(a)微觀結構圖,紅色表示HCP 結構的原子,藍色表示BCC 結構的原子,綠色表示FCC 結構的原子,白色代表無序原子,晶界的位置可以通過觀察白色原子的位置得到;(b)剪切應力分布圖;(c) 溫度分布圖Fig.2.(a) Microstructure figures of the model 1 containing sigma5 grain boundaries under shock loading at 24 ps,where red represents HCP structure atoms,blue represents BCC structure atoms,green represents FCC structure atoms,and white represents disordered atoms,the position of the grain boundaries can be obtained by observing the position of white atoms.(b) The shear stress distribution;(c) the temperature distribution.
為了研究影響晶界兩側的響應差異的因素,還模擬了含sigma11 晶界的模型2 的沖擊響應過程,如圖3 所示.圖3(a)表明模型2 中晶界兩側的響應也是不對稱的,右側的相變產(chǎn)物多于左側.經(jīng)過對微觀結構的分析發(fā)現(xiàn),晶界兩側相變產(chǎn)物的類型也是不一樣的,需要進一步分析導致晶界兩側相變產(chǎn)物存在差異的原因.據(jù)目前已有的研究,不同的相變模式下相變路徑、相變過程以及相變的產(chǎn)物都有很大的差異,而目前最常見的兩種相變模式為應力援助模式[26]和應變誘導模式[27].應力援助模式的特點: 壓縮過程中能量累積到能克服相變的勢壘時開始發(fā)生相變,相變發(fā)生較慢.應變誘導模式的特點: 多發(fā)生在缺陷附近,缺陷附近非常規(guī)的局域結構可作為相變的核,因此相變閾值低,相變開始較快.根據(jù)圖3(a)可知,晶界左側的相變符合應力援助模式的特點,而右側的相變符合應變誘導模式的特點.圖3 還展示了模型2 沖擊加載下的剪切應力和溫度的分布情況,與圖2(b)不同的是晶界右側沒有塑性形變區(qū)(圖2(b)中的區(qū)域2),即剪切應力減小的相變區(qū)與彈性形變區(qū)直接相鄰,這表明晶界右側的晶格不需要被壓縮太多就可以開啟相變.而且,模型2 晶界右側的相變閾值要明顯低于左側,右側相變的這些特點也是符合缺陷誘導模式的.晶界左側由于有反射波的影響其剪切應力的分布較為復雜,但是整體比晶界右側要高,這對應于左側的相變閾值高和相變產(chǎn)物少.

圖3 沖擊加載下含sigma11 晶界的模型2 在20 ps 時刻的(a) 微觀結構圖,顏色標示同圖2(a);(b) 剪切應力分布圖;(c) 溫度分布圖Fig.3.(a) Microstructure of the model 2 containing a sigma11 grain boundary under shock loading at 20 ps,and the color marking is the same as that of Figure 2(a);(b) the shear stress distribution of the model 2;(c) the temperature distribution of the model 2.
為了進一步確認晶界兩側的相變模式,分析了晶界兩側相變產(chǎn)物的c軸和相變動力學.圖4 是BCC 結構轉變至HCP 結構的示意圖,左側是BCC結構的一組原子,右側是HCP 結構的一組原子,選取的原子相變后剛好組成HCP 晶格的一個晶胞.BCC 晶格要轉變至HCP 晶格需要兩個過程,一是晶格沿〈001〉方向壓縮,對應原子2 和3 的距離減小,同時原子1 和2 之間的距離小幅度地增加.二是相鄰層的原子之間發(fā)生相對滑移,在圖4中滑移即黃色原子相對紫色原子整體移動.相變需要這兩個過程結合起來,但是這兩個過程的時序卻不是固定的,比如說,在缺陷誘導模式相變過程中,由于有缺陷幫助形核,相變閾值較低,壓縮與滑移過程幾乎同時進行,而應力援助模式相變過程中,由于需要靠壓縮積累相變的能量來克服勢壘,因此滑移落后于壓縮開始的時刻,模型2 中晶界兩側相變的差異將會對照上述兩種相變模式的特點來討論.圖5(a)和圖5(b)分別展示了模型2 中所選取的晶界兩側附近的HCP 晶胞,通過c軸的定向來確定相變前后的晶向關系.經(jīng)分析圖5(a)中晶界左側晶胞的c軸與模型中X,Y和Z軸的夾角分別約為90o,25o和125o,與模型2 晶界左側晶體中的晶向幾乎一致;而圖5(b)中晶界右側晶胞的c軸與X軸一致,即與晶界右側晶體中的[110]晶向一致,由此可知,相變后晶界兩側HCP 晶胞的c軸是不對稱的.通過觀察圖5 中兩個晶胞的形態(tài)也可以得到相變后晶界兩邊晶胞的c軸是不對稱的,即兩邊經(jīng)歷了不同的轉變路徑.

圖5 模型2 中sigma11 晶界兩側發(fā)生相變后分別選取的一個HCP 晶胞的原子形貌圖 (a) 左側;(b) 右側Fig.5.Atomic morphology of an HCP cell selected after phase transition on both sides of grain boundary in model 2:(a) On the left;(b) on the right.
為了了解晶界兩側相變動力學的差異,追蹤了晶界兩側原子轉變過程中的位置.圖6 中的晶格參數(shù)a,b和r的定義對應于圖4:a是原子2 和3 之間的距離,a減小可表征轉變過程中晶格的壓縮過程,b是原子1 和2 之間的距離,r是原子1 和4 之間的距離比上原子2 和4 之間的距離,r的變化則可以表征原子4 等黃色原子相對于紫色原子的滑移過程.圖6(a)和圖6(b)分別展示了模型2 中晶界兩側的相變過程中a,b和r隨時間的變化,圖中p點表示晶格受到?jīng)_擊后壓縮至最大值附近的時刻,s點表示相鄰層的原子滑移過程開啟的時刻.從圖6 可知,晶界左側的p點與s點之間存在一個時間差,即相變過程中壓縮開始一段時間后才開啟滑移過程,而右側的p點與s點幾乎同時出現(xiàn),表明兩個過程幾乎同時發(fā)生.上述結果表明右側滑移發(fā)生的時刻要明顯早于左側,這對應于前面討論過的右側有更低的相變閾值.聯(lián)系前面兩種不同的相變模式對應壓縮和滑移過程不同的時序,左側的規(guī)律符合應力援助模式,右側的則符合缺陷誘導模式.至此,從多個角度證實了模型2 中晶界兩側的相變存在著很大的差異,結合模型1 的討論可以得出,對稱晶界兩側在受到?jīng)_擊時確實存在著不對稱的響應.

圖6 展示了沖擊過程中晶界兩側的晶格參數(shù)(lattice parameter)隨時間的變化 (a) 對應晶界左側;(b) 對應晶界右側.圖中a 和b 已經(jīng)被歸一化處理了Fig.6.Variation of lattice parameters on the (a) left and(b) right sides of grain boundaries over time during the shock process,where a and b have been scaled.
為了探究影響對稱晶界兩側不對稱響應的原因,關注了晶界兩側在沖擊波到達后未發(fā)生相變的時刻.選取一組沿Z軸均勻分布的原子,記錄其沖擊過程中在X,Y軸方向的坐標與初始坐標的偏差,模型1 中所選取的原子在t=12 ps 時的位置偏差如圖7(a)所示.從圖7(a)可知,原子在X軸方向相對于初始位置沒有移動,而在Y軸方向沖擊波掃過區(qū)域的原子都有了不同的偏移.活塞法產(chǎn)生的沖擊波普遍被近似為一維沖擊波,因此很少關注沖擊波在垂直于沖擊方向的X和Y軸方向產(chǎn)生的影響,圖7(a)表明“一維沖擊波”在晶體內傳播時對垂直傳播方向的X和Y軸方向也有不可忽視的影響.在圖7(a)中,沖擊過程中原子只在Y軸上有移動,模型1 的X軸為〈100〉晶向,Y軸為〈031〉晶向,根據(jù)BCC 晶格的特點,晶格關于〈100〉晶向是對稱的,而關于〈031〉晶向是非對稱的,這表明原子在受到?jīng)_擊時在X和Y軸方向不同的表現(xiàn)與其是否是對稱晶向有關.

圖7 活塞速度為0.5 km/s 的加載下,模型1 在12 ps 時內部原子的(a)位置偏差和(b)剪切應力隨時間的變化Fig.7.Time dependent variation of (a) position deviation and (b) shear stress of atoms in model 1 under loading with a piston speed of 0.5 km/s at 12 ps.
圖7(a)中晶界兩側的原子在Y軸方向的移動幅度并不是相等的,這對應晶界兩側的剪切應力有較大的差異.圖7(b)展示了模型1 在12 ps 時的剪切應力的分布,這個時刻沖擊波已到達晶界,但是晶界處還沒有發(fā)生相變,因此研究此時晶界兩側的差異對于理解兩邊不同的響應有重要的意義.從圖7(b)可知,晶界兩側附近的剪切應力都低于離晶界較遠的區(qū)域,但是晶界右側的剪切應力明顯低于左側,這對應模型1 中晶界兩側只有右側發(fā)生了相變.在模型2 中,因為晶界右側發(fā)生的是缺陷誘導模式相變,相變速度較快,因此捕捉不到?jīng)_擊波到達但相變未開始的圖像,所以將活塞速度降低至0.2 km/s,該沖擊強度下晶界兩側都不會發(fā)生相變,這樣可以排除相變對晶界兩側剪切應力的影響.圖8 對應模型2,與圖7 類似,所追蹤的原子在沖擊下也只在Y軸方向有移動,Y軸對應的〈332〉晶向是非對稱的,X軸的〈110〉晶向是對稱晶向.而且圖8 中晶界兩側移動偏差和剪切應力的差異與模型2 在沖擊加載下晶界兩側不一樣的相變過程的聯(lián)系也是類似于圖7,這表明移動偏差的差異與該方向的晶向是否對稱有密切的聯(lián)系,即沖擊的晶向是影響晶界兩側發(fā)生不對稱相變的重要原因.

圖8 活塞速度為0.2 km/s 的加載下,模型2 在17.5 ps 時原子的(a)位置偏差和(b)剪切應力隨時間的變化Fig.8.Time dependent variation of (a) position deviation and (b) shear stress of atoms in model 2 at 17.5 ps with a piston speed of 0.2 km/s.
之前研究了沿非中心對稱的晶向沖擊時單晶鐵的響應,發(fā)現(xiàn)X軸和Y軸方向對應的晶向是否對稱對相變路徑的選擇產(chǎn)生了重要的影響[31],結合本文的研究可以得到鐵的各向異性對其沖擊響應有重大的影響.特別是,之前的研究發(fā)現(xiàn)晶格的各向異性的影響并不只存在于鐵的模擬中,而是廣泛存在于BCC 晶格金屬的沖擊加載中,因此晶格的各向異性對沖擊加載下相變的影響應該得到進一步的重視.
本文利用分子動力學方法模擬了沖擊加載下對稱晶界兩側不對稱的響應過程,結論如下:
1)有對稱微觀結構的晶界兩側存在著不對稱的沖擊響應,本文中的具體表現(xiàn)為模型1 中sigma5晶界兩側的相變閾值存在較大的差異,而模型2中sigma11 晶界兩側的相變閾值和相變路徑都存在著較大的差異,特別是該晶界兩側不同的相變動力學表明晶界兩側觸發(fā)了不同的相變模式,即應力援助模式與缺陷誘導模式;
2)沖擊下兩個模型中的原子都會沿Y軸方向偏移,這表明沿非中心對稱晶向沖擊時,活塞法產(chǎn)生的沖擊波不應該再被簡單視為一維的;還發(fā)現(xiàn)沖擊下原子在垂直于沖擊方向的X和Y軸上是否偏移與晶向的對稱性密切相關,并且在有偏移時晶界兩側的偏移是不一樣的,這對應晶界兩側的剪切應力分布也是不對稱的,剪切應力的差異最終導致晶界兩側的相變閾值和相變模式出現(xiàn)差異.綜上所述,對稱晶界兩側的沖擊響應與晶向的對稱性密不可分,因此晶格的各向異性是造成晶界兩側不對稱響應的重要因素.