王云天 曾祥國 楊鑫

(四川大學(xué)建筑與環(huán)境學(xué)院, 深地科學(xué)與工程教育部重點實驗室, 成都 610065)

采用嵌入原子勢的分子動力學(xué)模擬方法, 研究了5 × 109 s—1應(yīng)變率下, 溫度效應(yīng)對單晶鐵中孔洞成核與生長的影響, 并對 NAG (nucleation and growth)模型在單晶鐵中的適用性進行了探討.結(jié)果表明：隨著溫度的升高, 單晶鐵的抗拉強度峰值降低, 1100 K溫度下單晶鐵抗拉強度峰值比100 K溫度下降低了35.9%.在100—700 K溫度下, 拉應(yīng)力時程曲線表現(xiàn)出雙峰值特點, 分析表明, 第一峰值是由于拉應(yīng)力升高引起內(nèi)部結(jié)構(gòu)發(fā)生相變而產(chǎn)生, 第二峰值則是因發(fā)生孔洞成核與生長而產(chǎn)生; 900—1100 K溫度下, 拉應(yīng)力時程曲線表現(xiàn)為單峰值, 孔洞成核與生長是拉應(yīng)力下降的主要原因.分析發(fā)現(xiàn), 孔洞在高溫下更容易成核, 高應(yīng)變率下單晶鐵中孔洞成核與生長和NAG模型有較好的符合度, 單晶鐵中孔洞成核閾值與生長閾值都遠高于低碳鋼,并且孔洞成核閾值與生長閾值隨著溫度的升高而逐漸降低.研究結(jié)果可為建立高應(yīng)變率下金屬材料動態(tài)損傷演化模型提供借鑒.

1 引 言

金屬材料在高速沖擊下的動態(tài)損傷問題是力學(xué)中最困難、最復(fù)雜的問題之一, 且在國防與民用領(lǐng)域中都有非常重要的應(yīng)用背景, 因此該問題廣受關(guān)注[1].在高速沖擊下, 材料內(nèi)部加載稀疏波和自由面反射稀疏波相互作用在內(nèi)部產(chǎn)生拉應(yīng)力區(qū)域,引起材料內(nèi)部微孔洞成核、生長和貫穿, 最終由于損傷的累積導(dǎo)致材料斷裂[2-6].

金屬材料動態(tài)損傷受多種因素影響, 例如材料的微結(jié)構(gòu)、初始缺陷、溫度、加載應(yīng)變率等, 其中加載應(yīng)變率的影響非常重要, 金屬材料的斷裂強度隨加載應(yīng)變率的增加而增大, 在極高的應(yīng)變率下材料的斷裂強度與理論強度十分接近[7].在實驗研究中,常用的研究金屬材料動態(tài)損傷的方法主要有霍普金森桿、輕氣炮、激光加載等, 加載應(yīng)變率在102—1010s—1之間, 利用上述加載手段, 已經(jīng)對不同應(yīng)變率下金屬材料的動態(tài)損傷進行了大量的研究[8-15].Zaretsky和Kanel[14]采用輕氣炮加載方法, 研究了應(yīng)變率達到 105s—1量級時多晶鐵樣品在不同初始溫度下的動態(tài)響應(yīng), 分析了溫度對層裂強度的影響以及沖擊引起的相變行為.Remington等[2]利用激光加載方法, 研究了不同晶粒尺寸與加載應(yīng)變率 (106—107s—1)下鉭的動態(tài)響應(yīng), 結(jié)果表明晶粒尺寸對層裂強度有顯著影響, 高應(yīng)變率下孔洞在晶界處更容易成核.

在實驗研究的基礎(chǔ)上, 建立了描述金屬材料動態(tài)損傷的理論模型, 代表性的理論模型包括：美國洛斯阿拉莫斯國家實驗室提出的Gurson模型[16]和Johnson模型[17], 以及Curran等[3]提出的NAG(nucleation and growth)模型, 這些理論模型描述了宏觀尺度下金屬材料的動態(tài)損傷行為.金屬材料在高速沖擊下的動態(tài)損傷是一個非常復(fù)雜的問題,需要考慮的因素眾多, 包括孔洞的成核與生長、相變、熱耗散等, 并且時間尺度從皮秒到毫秒, 空間尺度從納米到毫米[18], 因此, 難以在一個理論模型中將所有因素考慮在內(nèi).

目前, 分子動力學(xué)模擬已成為實驗研究的重要補充手段, 為高應(yīng)變率下金屬材料動態(tài)損傷研究提供了新的方法.分子動力學(xué)模擬可以在微觀尺度揭示金屬材料的損傷過程, 并且可以直接測量相關(guān)物理量的變化, 例如溫度、拉應(yīng)力、相變和孔洞演化等數(shù)據(jù)[18-24], 有助于從根本上了解金屬材料在高速沖擊下的動態(tài)損傷機制.利用分子動力學(xué)模擬,Mayer[25]對單晶鐵在應(yīng)變率為 5 × 109s—1、單軸拉伸加載條件下的動態(tài)剪切與拉伸強度進行了研究,結(jié)果表明位錯的均勻成核不會導(dǎo)致剪切應(yīng)力的松弛, 同時發(fā)現(xiàn)拉應(yīng)力時程曲線中有兩個峰值出現(xiàn),Mayer認為拉應(yīng)力的第一個峰值是由于位錯成核導(dǎo)致的塑性變形而產(chǎn)生的, 第二個峰值是由孔洞成核與生長而產(chǎn)生的.Rawat和Raole[19]分析了應(yīng)變率為 1.5 × 108—1.5 × 1010s—1、三軸拉伸加載條件下單晶鐵中的孔洞演化情況, 結(jié)果表明在高應(yīng)變率下孔洞更容易成核, 但孔洞在較低的應(yīng)變率下有更快的生長速度.Shao等[26]研究了應(yīng)變率為2 ×109s—1、壓縮加載條件下溫度和預(yù)置空洞對單晶鐵中結(jié)構(gòu)變化的影響, 結(jié)果表明在相變過程中密排六方結(jié)構(gòu) (hexagonal close packed crystal structure,HCP)相的剪切應(yīng)力總是低于體心立方結(jié)構(gòu)(bodycentered cubic crystal structure, BCC)相的剪切應(yīng)力.馬文等[27]對沖擊壓縮條件下納米多晶鐵的相變過程進行了研究, 結(jié)果表明晶粒尺寸對晶界變形及相變原子比例有顯著影響.分子動力學(xué)模擬的結(jié)果也可以為理論模型提供參數(shù), Sugandhi等[28]通過粒子群優(yōu)化算法從分子動力學(xué)模擬的結(jié)果中獲得NAG模型的最佳擬合參數(shù).

Rawat等[29]采用分子動力學(xué)方法研究了不同溫度下單晶銅在 5 × 109s—1應(yīng)變率下的動態(tài)損傷情況, 結(jié)果表明溫度對單晶銅中孔洞成核與生長有顯著影響, 同時對NAG模型在微觀尺度的適用性進行了分析.目前已經(jīng)進行的關(guān)于鐵的動態(tài)損傷研究中, 已有關(guān)于不同應(yīng)變率、相變、晶粒尺寸等因素的研究, 但在高應(yīng)變率下溫度的影響尚少有涉及, 并且NAG模型在微觀上能否描述單晶鐵的動態(tài)損傷過程的研究也比較少.因此, 本文使用分子動力學(xué)模擬方法, 建立了單晶鐵的三軸拉伸模型,該模型廣泛應(yīng)用于研究金屬材料在高應(yīng)變率下的動態(tài)損傷中[30-32].本文分析了溫度對單晶鐵在高應(yīng)變率下孔洞成核與生長的影響, 并探討了NAG模型在單晶鐵中的適用性.研究結(jié)果可為建立高應(yīng)變率下金屬材料動態(tài)損傷演化模型提供參考.

2 分子動力學(xué)模型

本文采用Mendelev等[33]提出的嵌入原子法(embedded atom method, EAM)描述單晶鐵原子間相互作用, 該勢函數(shù)廣泛應(yīng)用于高應(yīng)變率下單晶鐵的分子動力學(xué)模擬研究中[34,35].分子動力學(xué)模擬采用Lammps軟件[36]進行, 首先建立單晶鐵的分子動力學(xué)模型, 模型尺寸為 28.55 nm × 28.55 nm ×28.55 nm, 模型x,y和z坐標分別沿 [100], [010],[001]方向, 模型共包含 2 × 106個原子, 如圖1 所示.在3個方向施加周期性邊界條件, 在加載之前采用NPT系綜對模擬系統(tǒng)進行充分弛豫(30 ps),以確保模擬系統(tǒng)在加載前達到指定溫度的平衡狀態(tài), 其控制溫度分別為 100, 300, 500, 700, 900,1100 K.弛豫結(jié)束后, 使用速度標定法控制溫度,采用NVE系綜并在x,y和z坐標軸方向施加相等的應(yīng)變率使系統(tǒng)處于三軸拉伸狀態(tài).

圖1 單晶鐵三軸拉伸模型Fig.1.Model of single crystal iron under triaxial tension(atoms are colored by CNA).

使用Ovito軟件[37]對分子動力學(xué)模擬結(jié)果進行可視化分析, 采用Python語言編寫的程序進行孔洞體積分析, 其基本思想是將整個模擬區(qū)域劃分為小的立體方格并計算空方格的個數(shù), 然后統(tǒng)計在任意時刻的空方格數(shù)量來計算孔洞體積, 這種孔洞體積統(tǒng)計方法已被用于計算BCC晶體與面心立方結(jié)構(gòu) (face-centered cubic crystal structure, FCC)晶體的孔洞體積分析中[29,32,38]; 模型內(nèi)部結(jié)構(gòu)變化使用位錯分析 (dislocation analysis, DXA)[39]、徑向分布函數(shù) (radial distribution function, RDF)分析與共鄰分析 (common neighbor analysis, CNA)[40]等方法來進行; 使用基于C語言編寫的改進遺傳算法程序確定NAG模型最佳擬合參數(shù).

3 計算結(jié)果與討論

本節(jié)討論了在應(yīng)變率為 5 × 109s—1、三軸拉伸加載下單晶鐵在100—1100 K溫度下的動態(tài)響應(yīng).分析了溫度對拉應(yīng)力、孔洞體積分數(shù)、系統(tǒng)勢能的影響, 對拉應(yīng)力在不同溫度下的變化特點進行了探討, 同時還討論了NAG模型在高應(yīng)變率下單晶鐵孔洞成核與生長的適用性問題, 確定了NAG模型的最佳擬合參數(shù), 并分析了溫度對NAG模型參數(shù)的影響.

3.1 溫度對拉應(yīng)力的影響

通過分子動力學(xué)模擬計算了單晶鐵在100—1100 K溫度下拉應(yīng)力隨時間變化情況, 拉應(yīng)力時程曲線見圖2.分子動力學(xué)模擬得到單晶鐵在300 K 溫度下最大拉應(yīng)力約為 18.3 GPa, Mayer[25]研究單晶鐵在應(yīng)變率為 5 × 109s—1、單軸拉伸下得到的最大拉應(yīng)力約為 17.6 GPa, Ashitkov 等[13]在鐵(樣品厚度為50 mm)的層裂實驗中測得的最大拉應(yīng)力為 17.8—20.3 GPa.從圖2 可以看出, 在100—700 K溫度下, 拉應(yīng)力時程曲線有兩個峰值出現(xiàn), 而在900—1100 K溫度下拉應(yīng)力時程曲線只有單峰值出現(xiàn).

100—1100 K溫度下拉應(yīng)力峰值如圖3所示.從圖3可以看出, 隨著溫度的升高, 拉應(yīng)力峰值逐漸降低, 最大拉應(yīng)力為 100 K 時的 21.4 GPa, 最小拉應(yīng)力為 1100 K 時的 13.7 GPa, 降低幅度為35.9%, 表明溫度對單晶鐵的拉應(yīng)力峰值有明顯的影響, 相對而言, 應(yīng)變率的變化對單晶鐵的拉應(yīng)力峰值的影響并不明顯[19].

圖2 100—1100 K 溫度下拉應(yīng)力隨時間的變化Fig.2.Tensile stress as a function of time at 100—1100 K.

圖3 拉應(yīng)力峰值隨溫度的變化Fig.3.Relationship of peak pressure and temperature.

雙拉應(yīng)力峰值時程曲線以100 K溫度為例,從圖2可以看出, 拉應(yīng)力隨時間持續(xù)增長, 在15.7 ps時出現(xiàn)第一個峰值 (約為 21.4 GPa)后開始下降, 在 16.5 ps 時, 拉應(yīng)力停止降低, 隨后有所增長 (約增長 0.03 GPa), 19.1 ps時達到第二個峰值 (約為 16.9 GPa), 19.1 ps之后, 由于孔洞成核與生長, 拉應(yīng)力快速下降.在 5 × 109s—1、單軸拉伸下對單晶鐵的研究中, Mayer[25]也觀察到拉應(yīng)力時程曲線的雙峰值現(xiàn)象, 認為第一個拉應(yīng)力峰值內(nèi)結(jié)構(gòu)發(fā)生塑性變形, 第二個拉應(yīng)力峰值內(nèi)孔洞成核.雙拉應(yīng)力峰值時程曲線不僅在BCC結(jié)構(gòu)的單晶鐵中出現(xiàn), 在FCC結(jié)構(gòu)的單晶銅中也觀察到類似現(xiàn)象, Rawat等[29]在單晶銅的三軸拉伸分子動力學(xué)模擬中觀察到單晶銅在1250 K溫度下(接近銅的熔點)拉應(yīng)力時程曲線有兩個峰值出現(xiàn),Rawat認為拉應(yīng)力的第一個峰值是由于發(fā)生了相變, 第二個峰值是由于孔洞成核與生長, 這說明Mayer與Rawat等對拉應(yīng)力時程曲線的雙峰研究結(jié)果一致.不同的是, FCC結(jié)構(gòu)的單晶銅在1250 K溫度下出現(xiàn)拉應(yīng)力雙峰值現(xiàn)象, 而BCC結(jié)構(gòu)的單晶鐵是在100—700 K溫度下出現(xiàn), 這表明在三軸拉伸下溫度對BCC和FCC晶體結(jié)構(gòu)的微觀損傷演化的影響是不同的.

為探究拉應(yīng)力時程曲線出現(xiàn)雙峰值的原因, 對模擬結(jié)果進行了如下分析：1)分析拉應(yīng)力時程曲線與孔洞體積分數(shù)關(guān)系, 即確認孔洞的成核閾值與成核時間; 2)通過DXA, 辨別塑性變形發(fā)生時間;3)借助晶體結(jié)構(gòu)分析, 利用RDF分析和CNA方法來分析結(jié)構(gòu)變化.為便于對比分析, 在圖2中以100 K溫度為代表標示出了拉應(yīng)力雙峰值之間的四個特征點, 分別為：第一個峰值點(A點)、第一個峰值后拉應(yīng)力快速下降點(B點), 第一個峰值點后拉應(yīng)力再次增長點(C點)、第二個峰值點(D點), 100—700 K溫度下四個特征點時間如表1所列.

表1 100—700 K 溫度下四個特征點時間Table 1.Four characteristic points time at 100 —700 K.

從表1 可以看出, 在 100—500 K 溫度下, 第一個拉應(yīng)力峰值出現(xiàn)的時間逐漸提前, 500 K溫度下第一個拉應(yīng)力峰值出現(xiàn)的時間比100 K溫度下提前了 2.9 ps, 700 K 溫度下第一個拉應(yīng)力峰值出現(xiàn)時間較500 K溫度下有所延后, 但相差不大, 僅為 0.5 ps; 100—700 K 溫度下拉應(yīng)力再次增長的時間以及第二個拉應(yīng)力峰值出現(xiàn)的時間隨溫度的升高逐漸提前, 700 K溫度下拉應(yīng)力第二個峰值較100 K溫度下提前了1.5 ps.

3.2 溫度對孔洞體積分數(shù)的影響

圖4給出了100—1100 K溫度下孔洞體積分數(shù)隨時間的變化.從圖4可以看出, 盡管100—700 K與900—1100 K溫度下拉應(yīng)力時程曲線有不同的特點, 但是孔洞體積分數(shù)隨時間的變化趨勢基本相同, 隨著溫度的升高, 孔洞出現(xiàn)的時間提前,孔洞體積分數(shù)增長速率有所加快, 100—1100 K溫度下孔洞體積分數(shù)在 30 ps時基本相同, 約為29%.100—1100 K溫度下孔洞成核時間如圖5所示, 從圖5 可以看出, 1100 K 溫度下孔洞成核時間是 12.8 ps, 100 K 溫度下孔洞成核時間為 17.4 ps,1100 K溫度下孔洞成核時間比100 K溫度下提前了26.4%, 這表明溫度對孔洞成核時間有顯著影響, 孔洞成核時間隨溫度升高而逐漸提前.圖6給出了體積分數(shù)為0.1時, 100—1100 K溫度下系統(tǒng)內(nèi)部孔洞分布情況.從圖6可以看出, 系統(tǒng)內(nèi)部的孔洞總體積是多個孔洞的成核、生長和貫穿累積而來, 1100 K溫度下孔洞數(shù)量明顯高于100 K溫度,這表明隨著溫度的升高系統(tǒng)內(nèi)部孔洞成核數(shù)量更多, 導(dǎo)致孔洞體積分數(shù)增長速率隨溫度升高而加快.

圖4 100—1100 K 溫度下孔洞體積分數(shù)隨時間的變化Fig.4.Void volume fraction as a function of time at 100—1100 K.

圖5 100—1100 K 溫度下孔洞成核時間Fig.5.Void nucleation time at 100—1100 K.

圖6 100—1100 K 溫度下內(nèi)部孔洞分布 (孔洞體積分數(shù)為 0.1 時)Fig.6.Distribution of voids at 100—1100 K (void volume fraction is 0.1).

圖7給出了孔洞體積分數(shù)隨時間的變化情況,并與拉應(yīng)力時程曲線進行了對比, 圖中虛線對應(yīng)孔洞成核的時間(以tn表示), 孔洞成核后體積分數(shù)開始增長.從圖7可以看出, 拉應(yīng)力時程曲線在特征點C之前孔洞體積分數(shù)一直保持為0; 在特征點C附近, 孔洞成核, 孔洞體積分數(shù)開始增長, 但增長速度較慢, 直到特征點D之后, 增長速率才明顯加快.圖8給出了100—700 K溫度下拉應(yīng)力時程曲線特征點D現(xiàn), 即孔洞沒有成核; 在第一個峰值之后, 拉應(yīng)力下降到特征點C附近時的孔洞分布情況, 可以看出系統(tǒng)內(nèi)部有明顯的孔洞產(chǎn)生.因此, 100—700 K溫度下拉應(yīng)力達到第一個峰值時沒有孔洞出時, 孔洞開始成核, 之后孔洞體積分數(shù)持續(xù)增長.

100—1100 K溫度下孔洞體積分數(shù)的發(fā)展可以分為3個階段, 以900和1100 K為例進行說明,如圖7所示.從圖7可以看出, 孔洞體積分數(shù)增長可以分為3個階段.

1)初始增長階段 (S1), 與 100—700 K 溫度下相同, 在拉應(yīng)力達到峰值之前已有孔洞成核, 此時孔洞體積增長主要是孔洞成核帶來, 初始增長階段孔洞體積分數(shù)增長緩慢, 表明孔洞還未開始大規(guī)模生長, 直到拉應(yīng)力時程曲線峰值拐點之前.從圖5可以看出, 隨著溫度的升高S1階段逐漸提前,1100 K比100 K溫度下孔洞出現(xiàn)時間提前了約4.6 ps, 這表明孔洞在高溫下更容易成核.

2)快速增長階段(S2), 當孔洞體積分數(shù)達到0.01后, 成核孔洞開始快速生長, 是孔洞體積分數(shù)迅速升高的主要原因.S2階段孔洞體積分數(shù)增長速率最快, 孔洞體積分數(shù)變化趨勢與指數(shù)函數(shù)特點相似.通過與拉應(yīng)力時程曲線的對比, 可以發(fā)現(xiàn)孔洞體積分數(shù)快速增長與拉應(yīng)力時程曲線的快速下降對應(yīng), 表明孔洞體積分數(shù)的快速增長導(dǎo)致了拉應(yīng)力的快速下降.1100 K溫度下S2階段持續(xù)時間約為 3 ps, 100 K 溫度下 S2 階段約為 7.5 ps, 二者相差2.5倍, 這表明高溫下孔洞生長的速度更快, 這也導(dǎo)致了高溫下拉應(yīng)力下降速率更快.

3)線性增長階段(S3), 此時孔洞生長減慢, 孔洞貫穿帶來體積增長, 增長速率適中, 在 30 ps時不同溫度下孔洞體積分數(shù)數(shù)值相差不大, 約為29%.這表明單晶鐵在 5 × 109s—1應(yīng)變率三軸拉伸下, 溫度對最終孔洞體積分數(shù)的影響并不明顯, 這與應(yīng)變率對孔洞體積分數(shù)的影響不同[19].

3.3 位錯分析

圖7 100-1100 K 溫度下孔洞體積分數(shù)與拉應(yīng)力的關(guān)系 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.7.Void volume fraction as a function of time at 100-700 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 K.

孔洞體積分數(shù)分析表明拉應(yīng)力時程曲線的特征點C附近孔洞開始成核與生長, 拉應(yīng)力時程曲線第一峰值內(nèi)沒有孔洞出現(xiàn).使用Ovito軟件對分子動力學(xué)模擬結(jié)果進行了DXA, 圖9給出了100—1100 K溫度下位錯密度隨時間的變化情況, 圖中時間軸上的標記點(紅色圓點)是拉應(yīng)力時程曲線峰值出現(xiàn)的時間, 在100—700 K溫度下是第二個峰值出現(xiàn)的時間.

DXA 表明, 單晶鐵在應(yīng)變率為 5 × 109s—1、三軸拉伸加載下有三種類型位錯出現(xiàn)：1/2位錯、位錯、位錯, 這些位錯類型在單晶鐵及其合金的實驗中都有觀測到[41].除了這三種類型的位錯之外, 還有一部分無法識別的位錯類型,定義為其他(Others).從圖9可以看出, 在100—1100 K溫度下, 單晶鐵中 1/2位錯是最主要的位錯類型, 占比最大,位錯次之, 占比略高于其余兩種位錯類型.100 K溫度下位錯密度在21 ps開始增長, 300—900 K 溫度下位錯密度在20 ps開始增長, 1100 K 溫度下位錯密度在18 ps開始增長.1100 K 溫度下位錯出現(xiàn)時間比100 K 溫度提前了 3 ps左右, 這表明溫度的升高使位錯出現(xiàn)的時間有所提前.通過與拉應(yīng)力時程曲線的對比分析, 可以發(fā)現(xiàn)位錯密度增長時間比拉應(yīng)力時程曲線峰值(100—700 K拉應(yīng)力第二個峰值)出現(xiàn)時間滯后 2—3 ps左右.對于 100—700 K溫度下拉應(yīng)力時程曲線, 位錯密度開始增長的時間位于拉應(yīng)力時程曲線第二個峰值拐點之后, 此時拉應(yīng)力開始快速下降, 這表明拉應(yīng)力時程曲線中第一個峰值內(nèi)沒有出現(xiàn)因位錯而引起的塑性變形, 這與Mayer[25]的結(jié)論不同.Rawat等[29]在單晶銅的三軸拉伸分子動力學(xué)模擬中認為拉應(yīng)力時程曲線第一個峰值內(nèi)發(fā)生了相變, 為了探討拉應(yīng)力時程曲線第一個峰值內(nèi)是否發(fā)生了相變, 對分子動力學(xué)模擬結(jié)果進行了RDF分析和CNA.

圖8 100-700 K 溫度下拉應(yīng)力時程曲線第二峰值點孔洞分布Fig.8.Void distribution on the second-peak of tensile stress at 100-700 K.

3.4 徑向分布函數(shù)分析

徑向分布函數(shù)通過給定某個粒子的坐標, 計算其他粒子在空間的分布概率隨距離變化情況, 常用于研究物質(zhì)內(nèi)部的有序性.晶體具有規(guī)則的周期性結(jié)構(gòu), 原子在其晶格位置附近波動, 由于晶體有序的結(jié)構(gòu), 徑向分布函數(shù)有長程的峰, 在晶體結(jié)構(gòu)轉(zhuǎn)變或被破壞后, 長程的峰消失, 因此可以用徑向分布函數(shù)來反應(yīng)模擬體系內(nèi)部晶體結(jié)構(gòu)變化.使用Ovito軟件計算了100—1100 K溫度下整個系統(tǒng)RDF 變化情況, 100—700 K 溫度下選取了 0 ps、拉應(yīng)力時程曲線的四個特征點對應(yīng)時間, 900—1100 K溫度下選取了0 ps、孔洞成核時間、拉應(yīng)力峰值時間以及拉應(yīng)力峰值后1 ps, 如圖10所示.可以看出, 在 100—700 K 溫度下 0 ps時系統(tǒng)內(nèi)部長程有序, 遠距離處依然有較明顯的曲線波動, 在拉應(yīng)力時程曲線第一個峰值時遠距離處的曲線波動開始減弱; 特征點B過后, 遠距離的曲線波動基本消失, 只有近距離內(nèi)出現(xiàn)明顯的曲線波動, 系統(tǒng)內(nèi)部長程有序消失.在 900—1100 K 溫度下, 0 ps時系統(tǒng)內(nèi)部長程有序, 孔洞成核使系統(tǒng)內(nèi)部長程有序減弱, 在拉應(yīng)力達到峰值后, 內(nèi)部長程有序完全消失.

RDF分析表明, 100—700 K溫度下拉應(yīng)力達到第一個峰值后系統(tǒng)內(nèi)部長程有序開始消失, 表明在拉應(yīng)力達到第一個峰值后系統(tǒng)內(nèi)部發(fā)生了結(jié)構(gòu)改變; 在 900—1100 K 溫度下, 孔洞成核引起系統(tǒng)內(nèi)部長程有序消失.

3.5 共鄰分析

RDF分析表明, 拉應(yīng)力時程曲線的第一個峰值內(nèi)系統(tǒng)內(nèi)部可能發(fā)生了相變, 使用Ovito軟件中的自適應(yīng)共鄰分析方法, 對系統(tǒng)的晶體結(jié)構(gòu)變化進行了分析.圖11給出了100—1100 K溫度下系統(tǒng)內(nèi)部晶體結(jié)構(gòu)變化情況, 圖中FCC代表面心立方結(jié)構(gòu)、BCC代表體心立方結(jié)構(gòu)、HCP代表密排六方結(jié)構(gòu)、Other代表未知結(jié)構(gòu).100—700 K 溫度下選取拉應(yīng)力時程曲線的四個特征點對應(yīng)的時間,900—1100 K溫度下選取孔洞成核時間與拉應(yīng)力峰值時間, 將模型中間進行切片觀察, 切片厚度為5 ?, 圖中使用黑色圓圈將內(nèi)部孔洞標示出來.

圖9 100-1100 K 溫度下位錯密度變化情況 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.9.Evolution of dislocation density as a function of time at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K;(f) 1100 K.

從圖11可以看出：在 100—700 K溫度下拉應(yīng)力達到第一個峰值時, 系統(tǒng)內(nèi)部出現(xiàn)FCC,HCP和未知結(jié)構(gòu), 隨著拉應(yīng)力的降低, FCC,HCP和未知結(jié)構(gòu)數(shù)量逐漸增多; 在拉應(yīng)力達到特征點C時, 大部分BCC結(jié)構(gòu)已經(jīng)轉(zhuǎn)變?yōu)镕CC,HCP和未知結(jié)構(gòu), 只有少量BCC結(jié)構(gòu)留存; 在拉應(yīng)力達到第二個峰值時, 可以觀察到多個孔洞出現(xiàn), 孔洞被未知結(jié)構(gòu)包圍.RDF 分析表明, 拉應(yīng)力時程曲線的第一個峰值內(nèi)系統(tǒng)內(nèi)部發(fā)生了相變,BCC結(jié)構(gòu)轉(zhuǎn)變?yōu)镕CC, HCP和未知結(jié)構(gòu), 沒有孔洞出現(xiàn); 在第二個峰值時, 只有少量BCC結(jié)構(gòu)留存, 孔洞在未知結(jié)構(gòu)區(qū)域內(nèi)成核.在 900—1100 K溫度下系統(tǒng)內(nèi)部未知結(jié)構(gòu)數(shù)量增多, FCC結(jié)構(gòu)數(shù)量有所減少, 特別是在 1100 K 溫度下, 系統(tǒng)內(nèi)部出現(xiàn)大量的未知結(jié)構(gòu), 在拉應(yīng)力達到峰值時可以看到明顯的孔洞出現(xiàn), 孔洞被未知結(jié)構(gòu)包圍.在鐵納米管和納米晶體鐵的單軸拉伸分子動力學(xué)模擬中也出現(xiàn)了BCC結(jié)構(gòu)到HCP結(jié)構(gòu)的相變出現(xiàn)[42,43].本文中出現(xiàn)的BCC結(jié)構(gòu)到HCP結(jié)構(gòu)的轉(zhuǎn)變可能是α-ε相變, 在鐵的沖擊加載實驗和分子動力學(xué)模擬中都有觀察到[44-47].

圖10 100-1100 K 溫度下徑向分布函數(shù)變化情況 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.10.Radial distribution function of the system at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K;(f) 1100 K.

圖11定性地給出了系統(tǒng)內(nèi)部結(jié)構(gòu)變化情況,為定量地描述系統(tǒng)的結(jié)構(gòu)變化, 計算了100—1100 K溫度下FCC, BCC, HCP和未知結(jié)構(gòu)隨時間的變化情況, 如圖12 所示.從圖12 可以看出, BCC 結(jié)構(gòu)首先轉(zhuǎn)變?yōu)槲粗Y(jié)構(gòu), 在拉應(yīng)力時程曲線達到第一個峰值后, FCC, HCP 結(jié)構(gòu)開始出現(xiàn); FCC,HCP結(jié)構(gòu)在拉應(yīng)力時程曲線達到第二個峰值后占比達到最大.在 100—300 K 溫度下, 二者占比約為30%, 隨著溫度的升高, FCC結(jié)構(gòu)占比下降,500 K 時占比約為 25%, 700 K 時占比約為 20%;BCC結(jié)構(gòu)占比最低點出現(xiàn)在FCC, HCP占比峰值之后, 100 K 時約為 10%, 并且隨著溫度的升高BCC結(jié)構(gòu)占比逐漸降低, 900 K溫度下約為2%.從圖7可以看出孔洞體積分數(shù)在第二個拉應(yīng)力峰值后開始快速增長, 對比圖12結(jié)構(gòu)占比變化, 可以發(fā)現(xiàn)孔洞體積的快速增長導(dǎo)致FCC, HCP結(jié)構(gòu)占比降低, BCC結(jié)構(gòu)占比再次升高, 引起這種變化的機理將在以后的文章中進行進一步的分析.

對100—700 K溫度下拉應(yīng)力時程曲線的分析表明系統(tǒng)內(nèi)部首先發(fā)生了相變, 隨后孔洞成核與生長, 但是隨著溫度的升高, 在 900—1100 K 溫度下拉應(yīng)力時程曲線只有一個峰值出現(xiàn).從圖12可以看出, 在孔洞成核之前系統(tǒng)內(nèi)部已經(jīng)發(fā)生相變,FCC, HCP結(jié)構(gòu)占比已經(jīng)開始增長, 但與 100—700 K 溫度下不同的是, BCC 到 FCC, HCP 的相變并未導(dǎo)致拉應(yīng)力的降低, 原因可能是FCC,HCP相變發(fā)生的時間與孔洞成核的時間非常接近,孔洞成核時FCC, HCP相變的數(shù)量非常少, 因此由FCC, HCP相變引起的拉應(yīng)力降低幅度很小,拉應(yīng)力只出現(xiàn)一個峰值.結(jié)合圖7孔洞體積分數(shù)變化情況, 在拉應(yīng)力達到峰值后孔洞體積分數(shù)開始快速增長, 這表明孔洞成核與生長是導(dǎo)致900—1100 K溫度下拉應(yīng)力降低的主要原因[30,48].因此, 在100—700 K溫度下, 拉應(yīng)力下降的原因首先是系統(tǒng)內(nèi)部發(fā)生相變, BCC結(jié)構(gòu)向FCC, HCP和未知結(jié)構(gòu)轉(zhuǎn)變, 然后是孔洞成核與生長; 在 900—1100 K 溫度下, 孔洞成核與生長是導(dǎo)致拉應(yīng)力下降的主要原因.

圖11 100-1100 K 溫度下內(nèi)部結(jié)構(gòu)變化 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.11.Snapshots for the structural changes at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 K.

為了對比初始BCC結(jié)構(gòu)與重生成的BCC結(jié)構(gòu)之間的區(qū)別, 選取模型中部進行切片觀察, 為了能夠更加明顯地辨識孿晶結(jié)構(gòu), 選擇切片厚度為2 ?.圖13 給出了 300 K 溫度下系統(tǒng)內(nèi)部結(jié)構(gòu)變化情況, 紅色箭頭標示出孿晶出現(xiàn)的區(qū)域.從圖13可以看出, FCC和HCP結(jié)構(gòu)在拉應(yīng)力第二個峰值之后轉(zhuǎn)變?yōu)锽CC結(jié)構(gòu), 孔洞被未知結(jié)構(gòu)包圍, 當孔洞周圍的未知結(jié)構(gòu)原子隨著孔洞生長而移動時, 導(dǎo)致原本屬于FCC和HCP結(jié)構(gòu)的原子重新排列, 再次產(chǎn)生 BCC 結(jié)構(gòu).對比 24, 30 和 0 ps時的BCC結(jié)構(gòu), 可以發(fā)現(xiàn)重新生成的BCC結(jié)構(gòu)與初始BCC結(jié)構(gòu)并不相同, 重新出現(xiàn)的BCC結(jié)構(gòu)可能是變形孿晶的結(jié)果[49], 與BCC結(jié)構(gòu)Mo和BCC結(jié)構(gòu)Ta的孔洞生長分子動力學(xué)模擬中出現(xiàn)的變形孿晶情形相似[50].

圖12 100-1100 K 溫度下內(nèi)部晶體結(jié)構(gòu)占比 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.12.Crystal structure fraction as a function of time at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K;(f) 1100 K.

通過孔洞體積分數(shù), DXA, RDF, CNA 等分析方法對分子動力學(xué)模擬結(jié)果進行了分析, 結(jié)果表明拉應(yīng)力時程曲線的第一個峰值內(nèi)系統(tǒng)內(nèi)部發(fā)生了相變, 單晶鐵原有的BCC結(jié)構(gòu)轉(zhuǎn)變?yōu)镕CC,HCP和未知結(jié)構(gòu), 沒有出現(xiàn)位錯引起的塑性變形;在拉應(yīng)力時程曲線的第二個峰值內(nèi)孔洞開始成核與生長, 孔洞體積的快速增長使FCC和HCP結(jié)構(gòu)占比降低, 孔洞的生長導(dǎo)致BCC結(jié)構(gòu)重新產(chǎn)生,但與初始的BCC結(jié)構(gòu)不同, 新出現(xiàn)的BCC結(jié)構(gòu)可能是變形孿晶的結(jié)果.

3.6 溫度對原子鍵能及系統(tǒng)勢能的影響

圖14是鐵原子鍵能在不同體系溫度下的變化情況, 圖中ε是摩爾鍵能, 單位 kJ/mol, 摩爾鍵能的計算參照 Eberhart和Horner[51]的文章提出的方法, 計算時間點為0 ps, 即模擬體系充分弛豫后的時刻, 計算所需數(shù)據(jù)來自美國國家標準與技術(shù)研究院(NIST)[52].從圖14可以看出, 鐵原子鍵能在100 K 時為 104.68 kJ/mol, 隨著模擬體系溫度的升高, 鐵原子鍵能逐漸降低, 在 1100 K 時鐵原子鍵能為 93.99 kJ/mol, 降低幅度為 10.21%.這表明溫度的升高顯著影響了鐵原子鍵能, 在高溫下鐵原子鍵能明顯降低.

圖13 300 K 溫度下系統(tǒng)內(nèi)部結(jié)構(gòu)變化Fig.13.Structural changes at different time at 300 K.

圖14 300-1100 K 溫度下鐵原子鍵能Fig.14.Bond energy of iron at 300-1100 K.

圖15 100-1100 K 溫度下系統(tǒng)勢能Fig.15.Potential energy of the system at 100-1100 K.

圖15是不同溫度下系統(tǒng)勢能隨時間變化情況.從圖15可以看出, 系統(tǒng)總勢能隨著溫度的升高而增加.在 100—700 K 溫度下, 拉應(yīng)力時程曲線第一個峰值(圖中紅色圓圈標示)導(dǎo)致系統(tǒng)勢能減小, 但幅度不大, 并且隨著溫度升高這種影響越來越低, 此時系統(tǒng)內(nèi)部并未有孔洞成核, 是相變引起的系統(tǒng)勢能變化; 系統(tǒng)勢能峰值出現(xiàn)在拉應(yīng)力時程曲線第二個峰值拐點之后, 此時孔洞體積分數(shù)快速增長, 表明孔洞成核與生長引起勢能降低.在900—1100 K溫度下, 系統(tǒng)勢能隨時間持續(xù)增長,勢能曲線的峰值出現(xiàn)在拉應(yīng)力時程曲線拐點之后,隨后快速降低, 表明在高溫下孔洞成核與生長是導(dǎo)致勢能降低的主要原因.當有孔洞成核時, 原子間的相互作用會變得更強, 同時鍵長縮短, 導(dǎo)致勢能的快速下降.

3.7 NAG模型

Curran等[3]和Seaman等[53]通過對鋁和銅等金屬材料的研究, 揭示出其斷裂是大量準球形空洞成核和生長的結(jié)果.通過研究發(fā)現(xiàn)單位體積的孔洞數(shù)量以及孔洞體積與拉應(yīng)力呈指數(shù)關(guān)系, 從而提出一種考慮微孔洞成核和生長效應(yīng)的損傷模型—NAG模型.NAG模型是一種微細觀物理模型, 主要用于描述金屬材料中孔洞的成核和生長而發(fā)生的損傷過程.通過3.2部分對孔洞體積的分析, 發(fā)現(xiàn)單晶鐵在 5 × 109s—1應(yīng)變率 100—1100 K 溫度下孔洞體積非線性增長階段與NAG模型相符合.因此計算了100—1100 K溫度下NAG模型參數(shù)并且分析了溫度對NAG模型參數(shù)的影響.

3.7.1 NAG 模型介紹

孔洞的總體積由兩部分組成：孔洞成核體積和孔洞生長體積.

1)孔洞成核體積

當材料中的拉應(yīng)力Ps超過材料成核閾值Pn0時, 新的孔洞生成, 成核率按下式計算：

在Δt時間間隔內(nèi)孔洞成核體積由下式計算：

式中Rn為形核尺寸參數(shù).在NAG模型的計算中,形核尺寸參數(shù)Rn的值為3.1 ?, 即1.01倍鐵的晶格常數(shù), 也是孔洞體積分析中立方體方格的尺寸.

2)孔洞生長體積

當拉應(yīng)力Ps超過孔洞生長閾值Pg0時, 成核的孔洞就會生長, 這時新的孔洞體積由下式計算：

式中,Vv0為在每個時間間隔初始時的孔洞體積,η為材料黏度, 在Δt時間間隔內(nèi)孔洞成核與生長的總體積為

每個數(shù)據(jù)點的相對誤差φ定義如下：

式中,VMD為分子動力學(xué) (molecular dynaimics,MD)模擬中得到的孔洞體積分數(shù),VNAG為NAG模型得到的孔洞體積分數(shù).

相對平方誤差Σ定義如下：

式中M為數(shù)據(jù)點總數(shù).

3.7.2 NAG 參數(shù)分析

對于任意給定的一組NAG參數(shù), 在已知拉應(yīng)力時程曲線的前提下, 可以得到孔洞總體積隨時間變化的函數(shù)關(guān)系.參數(shù)擬合采用基于改進遺傳算法編寫的C語言程序進行, 得到了100—1100 K溫度下 NAG 模型擬合參數(shù)使用擬合得到的NAG模型參數(shù)計算了孔洞的體積分數(shù)變化情況, 并與MD模擬的結(jié)果進行了對比,如圖16所示.從圖16可以看出, 不同溫度下單晶鐵的孔洞體積分數(shù)變化與NAG模型符合較好.

圖16 100-1100 K 溫度下 NAG 與 MD 的孔洞體積分數(shù)結(jié)果的對比 (a) 100 K; (b) 300 K; (c) 500 K; (d) 700 K; (e) 900 K; (f) 1100 KFig.16.Comparison of void volume fraction between the NAG model and MD at 100-1100 K：(a) 100 K; (b) 300 K; (c) 500 K;(d) 700 K; (e) 900 K; (f) 1100 K.

表2 100-1100 K 溫度下 NAG 模型最佳擬合參數(shù)Table 2.Best-fit NAG parameters at 100-1100 K.

除了與MD結(jié)果進行對比之外, 還將100—1100 K 溫度下 NAG 模型擬合參數(shù) (Pn0,P1,,Pg0,η)與實驗中得到的低碳鋼 (mild steel) NAG模型參數(shù)[54]進行了對比, 詳見表2.

從表2可以看出：

1)單晶鐵的孔洞成核閾值(Pn0)遠高于低碳鋼中的成核閾值, 差距在10倍以上.這是由于低碳鋼是多晶結(jié)構(gòu), 存在大量的缺陷, 孔洞在晶界處更容易成核[2], 而在單晶鐵中不存在這些缺陷, 所以孔洞成核閾值更高.

2)單晶鐵的孔洞成核拉應(yīng)力靈敏度系數(shù)(P1)較低, 約為低碳鋼的20%, 這表明單晶鐵中的成核速率對拉應(yīng)力的敏感性比對多晶材料的敏感性高得多.

3)與低碳鋼相比, 單晶鐵的材料黏度η非常低, 相差約1000倍;

4)隨著溫度的升高, 孔洞成核閾值Pn0和孔洞生長閾值Pg0都逐漸減小.高溫下, 原子運動更加劇烈, 加速了原子間相互分離的趨勢, 導(dǎo)致孔洞成核與生長閾值降低;

5)成核閾值與生長閾值的比值在單晶鐵與低碳鋼中相差不大, 單晶鐵的成核閾值Pn0約為生長閾值Pg0的6倍, 低碳鋼中成核閾值是生長閾值的5倍[54].

4 結(jié) 論

本文使用分子動力學(xué)模擬方法, 研究了100—1100 K 溫度下單晶鐵在 5 × 109s—1應(yīng)變率三軸拉伸加載的動態(tài)響應(yīng), 通過拉應(yīng)力變化、孔洞體積分數(shù)以及微結(jié)構(gòu)分析研究了單晶鐵的動態(tài)響應(yīng)以及溫度效應(yīng)的影響, 并且對NAG模型在高應(yīng)變率下單晶鐵中的適用性進行了探討, 得到以下主要結(jié)論.

1)在 100—700 K溫度下, 拉應(yīng)力時程曲線出現(xiàn)兩個峰值, 采用 DXA, RDF, CNA 和孔洞體積分數(shù)等方法對結(jié)果進行了分析.分析表明, 100—700 K溫度下拉應(yīng)力時程曲線的第一個峰值內(nèi)系統(tǒng)內(nèi)部發(fā)生了相變, 原有的BCC結(jié)構(gòu)轉(zhuǎn)變?yōu)镕CC, HCP和未知結(jié)構(gòu), 在第一個峰值之后孔洞開始成核與生長.在 900—1100 K 溫度下, 拉應(yīng)力時程曲線只有一個峰值出現(xiàn), 孔洞成核與生長是導(dǎo)致拉應(yīng)力下降的主要原因.

2)隨著溫度的升高, 系統(tǒng)內(nèi)部拉應(yīng)力峰值降低, 最大拉應(yīng)力為 100 K 時的 21.4 GPa, 最小拉應(yīng)力為 1100 K 時的 13.7 GPa, 降低幅度為 35.9%;系統(tǒng)總勢能隨溫度升高而增加, BCC結(jié)構(gòu)轉(zhuǎn)變?yōu)镕CC, HCP和未知結(jié)構(gòu)會導(dǎo)致系統(tǒng)勢能降低, 但降幅不大, 孔洞快速生長是導(dǎo)致勢能快速降低的主要原因.

3)孔洞成核在發(fā)生HCP, FCC相變之后與拉應(yīng)力達到峰值之前(100—700 K是拉應(yīng)力時程曲線第二個峰值之前)發(fā)生, 孔洞在未知結(jié)構(gòu)區(qū)域內(nèi)成核, 孔洞的生長推動周邊原子的移動, 使FCC和HCP結(jié)構(gòu)排列轉(zhuǎn)變?yōu)樾碌腂CC結(jié)構(gòu), 新產(chǎn)生的BCC結(jié)構(gòu)不同于初始的BCC結(jié)構(gòu), 可能是變形孿晶的結(jié)果.

4)孔洞體積分數(shù)的增長可以分為三個階段,第一階段由孔洞成核主導(dǎo), 此時孔洞體積分數(shù)增長較慢, 第二階段由孔洞生長主導(dǎo), 此時孔洞體積分數(shù)增長最快, 呈指數(shù)型增長, 第三階段由孔洞貫穿主導(dǎo), 此時孔洞體積分數(shù)呈線性增長, 速度適中.孔洞在高溫下更容易成核與生長, 100 K溫度下孔洞體積分數(shù)快速增長階段S2的時間是1100 K溫度下的2.5倍.孔洞的快速生長階段與拉應(yīng)力的峰值相對應(yīng), 因此高溫下拉應(yīng)力下降速率隨孔洞體積分數(shù)快速增長階段時間縮短而加快.

5) NAG 模型可以用來描述單晶鐵在 5 × 109s—1應(yīng)變率100—1100 K溫度下孔洞體積非線性增長,在100—1100 K溫度下擬合結(jié)果與MD結(jié)果都符合較好, 誤差在 20% 以內(nèi).與低碳鋼相比, 單晶鐵的拉應(yīng)力敏感系數(shù)與黏度更低.隨著溫度的升高,單晶鐵的孔洞成核和生長閾值逐漸降低.

本文的分析結(jié)果有助于更全面地了解單晶鐵在高速沖擊下的動態(tài)響應(yīng), 研究結(jié)果可以為建立考慮溫度影響的單晶鐵在高應(yīng)變率下的損傷模型提供參考.本文研究只考慮了溫度的影響, 然而金屬材料在高應(yīng)變率下的動態(tài)響應(yīng)受多個因素影響, 溫度只是重要因素之一, 因此后續(xù)研究將考慮多種因素共同作用, 例如晶粒尺寸、預(yù)置孔洞、應(yīng)變率等因素與溫度共同作用對單晶鐵的動態(tài)響應(yīng)的影響.