極端應(yīng)變率下鎢孔洞演化及其力學(xué)響應(yīng)的分子動力學(xué)模擬

百曉丹，陳遜，馬運柱，劉文勝，梁超平

極端應(yīng)變率下鎢孔洞演化及其力學(xué)響應(yīng)的分子動力學(xué)模擬

百曉丹，陳遜，馬運柱，劉文勝，梁超平

(中南大學(xué) 輕質(zhì)高強(qiáng)國家重點實驗室，長沙 410083)

極端應(yīng)變率；分子動力學(xué)模擬；孔洞；塑性變形；位錯

鎢及鎢合金具有密度大、熔點高和高溫強(qiáng)度高等優(yōu)異性能，是國防軍事領(lǐng)域的關(guān)鍵材料之一，也是目前對付裝甲目標(biāo)等的主要動能武器彈芯材料[1]。隨著防御工事和武器裝備的發(fā)展，高超聲速飛行器技術(shù)及其軍事應(yīng)用成為現(xiàn)今研究的熱點[2]。因此，應(yīng)用于未來高超聲速領(lǐng)域的鎢及鎢合金材料將承受更加嚴(yán)苛的作業(yè)環(huán)境，如極端應(yīng)變率等?，F(xiàn)有高應(yīng)變率下的鎢及鎢合金變形行為研究主要通過Hopkinson拉/壓桿[3-4]、輕氣炮[5]和爆炸加載[1]等技術(shù)實現(xiàn)，應(yīng)變率范圍在102~ 106s?1。近年來，隨著激光能量的提高和美國等國家強(qiáng)激光裝置對學(xué)術(shù)界的開放，科學(xué)家得以通過激光驅(qū)動撞擊實驗，揭示材料在更高應(yīng)變率106~1010s?1下的塑性行為[6]。在激光驅(qū)動撞擊的極高應(yīng)變率下，材料出現(xiàn)與低應(yīng)變率下截然不同的塑性變形機(jī)制，如隨應(yīng)變率升高，塑性變形由位錯滑移為主轉(zhuǎn)變?yōu)閷\晶變形為主；隨應(yīng)變率進(jìn)一步升高，材料又發(fā)生結(jié)構(gòu)相變或非晶化轉(zhuǎn)變[7]。然而，高應(yīng)變率使得直接對材料塑性變形行為的觀察變得極為困難，使用第三代窄脈沖同步輻射光源能原位觀察材料在激光驅(qū)動撞擊下的晶體結(jié)構(gòu)信息，但無法獲得具體的位錯運動、孿晶變形以及相變過程[8]。分子動力學(xué)模擬可準(zhǔn)確紀(jì)錄原子在飛秒、皮秒尺度下的運動軌跡，已被證明是研究材料在極端應(yīng)變率下變形行為的最有效手段[6]，如ZEPEDA- RUIZ[9]和WEHRENBERG等[10]通過分子動力學(xué)預(yù)測了金屬鉭在極端應(yīng)變率下位錯到孿晶的轉(zhuǎn)變過程，并通過實驗驗證了模擬結(jié)果。鎢及鎢合金通常采用粉末冶金方法制備，存在大量孔洞，而孔洞演化又是材料在極端應(yīng)變率下碎裂失效的主要原因[6]。為此，本文作者采用分子動力學(xué)模擬來研究極端應(yīng)變率(108~1011s?1)下單晶鎢中孔洞的演化規(guī)律，獲得孔洞尺寸和應(yīng)變率對鎢動態(tài)力學(xué)性能的影響，并從原子尺度揭示孔洞演化的塑性機(jī)理，為鎢合金在極端應(yīng)變率下的應(yīng)用提供理論指導(dǎo)。

1 模擬方法

采用LAMMPS軟件進(jìn)行分子動力學(xué)模擬[11]，通過嵌入原子方法描述鎢原子間的相互作用[12?15]?？锥茨Ｐ褪窃诔叽鐬?8×48×48(為鎢的晶格常數(shù)，=0.316 5 nm)的體心立方(BCC)鎢單晶塊的中心掏半徑分別為12，15，17和19的圓孔(為鎢的原子半徑，=0.137 0 nm)，具體的模型細(xì)節(jié)列于表1。對所建立的模型，首先通過共軛梯度法靜態(tài)優(yōu)化其結(jié)構(gòu)和原子位置，并進(jìn)一步通過等溫等壓系綜(isoth- ermal-isobaric ensemble, NPT)獲得優(yōu)化后的結(jié)構(gòu)在300 K下的平衡結(jié)構(gòu)。其中時間步長設(shè)置為0.001 fs，運行50 000步后的孔洞結(jié)構(gòu)如圖1所示。最終，通過在正則系綜(canonical ensemble, NVT)下對模型的X軸(對應(yīng)[100]晶向)以恒定應(yīng)變率(108~1011s?1)施加體積變化來獲得體系的應(yīng)力狀態(tài)，模擬沿X軸壓縮到初始長度的1/4時終止。在模擬過程中，通過Langevin控溫器來保持體系溫度為300 K。為了獲得變形過程中的原子信息，根據(jù)不同應(yīng)變率來選取模擬的時間步長(2~0.5 fs)，以保證相同應(yīng)變間隔下的原子信息一致。變形后的缺陷通過OVITO軟件中的位錯提取算法(dislocation analysis, DXA)和公共近鄰分析(common neighbor analysis, CNA)進(jìn)行表征[16]。

表1 不同孔徑模型的結(jié)構(gòu)參數(shù)

Note:is the radius of tungsten atom,=0.137 0 nm; ais the lattice constant of tungsten, a=0.316 5 nm.

圖1 等溫等壓系綜(NPT)300 K弛豫后的孔洞結(jié)構(gòu)

2 結(jié)果與討論

2.1 孔洞尺寸效應(yīng)

對不同孔徑的孔洞模型進(jìn)行計算，結(jié)果表明孔洞的大小不會影響鎢的彈性和塑性變形行為。圖2所示為無孔洞的理想單晶模型和具有不同尺寸孔洞的模型在應(yīng)變率為108s?1時的應(yīng)力?應(yīng)變曲線。從圖中看出，在無孔洞的理想狀態(tài)下，鎢的屈服強(qiáng)度為58.5 GPa，孔洞的存在使得屈服強(qiáng)度急劇降低，孔半徑=12時屈服強(qiáng)度下降至29.0 GPa。隨孔徑增大，屈服強(qiáng)度進(jìn)一步降低，在=19時屈服強(qiáng)度降至26.0 GPa。產(chǎn)生這一現(xiàn)象的主要原因是孔洞處產(chǎn)生的高應(yīng)力集中和低的位錯啟動力。封閉孔的應(yīng)力集中計算公式為[17]：

式中：σmax為孔洞處的最大應(yīng)力；σ￥為遠(yuǎn)離孔洞處的應(yīng)力；R為孔的半徑。由式(1)可知孔洞處的應(yīng)力集中隨孔徑增大而增加，R=19r的孔洞處的應(yīng)力集中程度較R=12r孔洞處的增加25%，這是導(dǎo)致鎢的屈服強(qiáng)度下降的主要原因?？锥刺幐叩膽?yīng)力導(dǎo)致鎢的塑性變形提前發(fā)生，無孔洞的完美晶體在應(yīng)變ε=0.068時發(fā)生屈服，而在含有R=12r的孔洞時，在ε=0.039就已經(jīng)發(fā)生塑性變形。

從應(yīng)力?應(yīng)變曲線的彈性階段看出，孔洞的存在也導(dǎo)致鎢的彈性模量減小，并且隨孔洞尺寸增大，彈性模量進(jìn)一步降低。從原理上看，彈性模量是材料的本征屬性，反映了原子、離子或分子之間的鍵合強(qiáng)度。中心對稱參數(shù)(centrosymmetry parameter, CSP)的計算公式為[11]：

式中：為成鍵數(shù)；e為原子平衡間距；為受力面積；d/d為原子間勢能的一階導(dǎo)數(shù)。從式(3)可知彈性模量與成鍵數(shù)呈線性關(guān)系。表2所列為模型弛豫前斷1~4根鍵時的CSP值和總斷鍵數(shù)，式(3)中的成鍵數(shù)等于無孔洞時的總成鍵數(shù)total減去有孔洞時的斷鍵數(shù)broken，即=total?broken。因此彈性模量與斷鍵數(shù)呈線性關(guān)系，斜率為負(fù)值，即隨斷鍵數(shù)增加，彈性模量降低。圖3所示為彈性模量與總斷鍵數(shù)的對應(yīng)關(guān)系，圖中的彈性模量由圖2中不同孔徑的模型彈性變形階段線性擬合得到。從圖中3看出，彈性模量與斷鍵數(shù)呈很好的線性關(guān)系，隨孔洞半徑增大，斷鍵數(shù)增多，彈性模量由理想致密晶體(斷鍵數(shù)為0)的387.05 GPa降低到=19(斷鍵數(shù)為21 356)的370.83 GPa。從以上分析可知，內(nèi)封閉孔洞處的斷鍵導(dǎo)致單晶鎢的彈性模量降低。

表2 不同孔徑的孔洞模型斷1~4根鍵時的中心對稱參數(shù)和總斷鍵數(shù)

圖3 含孔洞單晶鎢的彈性模量與斷裂金屬鍵數(shù)量及孔洞半徑的關(guān)系

2.2 應(yīng)變率的影響

圖4 應(yīng)變率與孔洞半徑對鎢壓縮性能的影響

(a) Stress-true strain curve of tungsten under different strain rates (void radius=12); (b) Yield strength of tungsten with different void radius as a function of strain rate

2.3 孔洞塑性變形機(jī)理

為了研究孔洞在應(yīng)力集中下的變形機(jī)理，采用通過位錯提取算法和公共近鄰分析來研究孔洞形貌變化和缺陷。與實驗不同，分子動力學(xué)模擬得到變形過程中的原子結(jié)構(gòu)信息，圖6所示為典型應(yīng)力?應(yīng)變曲線上各點的原子結(jié)構(gòu)信息。從圖6看出，在屈服點a時，孔洞沿[100]發(fā)生一定的壓縮，且位錯開始在某些面上形核；壓縮進(jìn)行到b和c點時，一個位錯環(huán)開始向外發(fā)射；到達(dá)d點時，第一個位錯環(huán)完成一次剪切回到孔洞，而其他方向的位錯環(huán)繼續(xù)向前發(fā)展；不同方向的位錯環(huán)不斷發(fā)射和運動，導(dǎo)致孔洞發(fā)生不同方向的剪切，最終孔洞坍塌，呈現(xiàn)e點處不規(guī)則形狀的孔洞，如圖6(e)所示。

圖5 單晶鎢在不同應(yīng)變率下屈服時的原子結(jié)構(gòu)快照和徑向分布函數(shù)

通過對應(yīng)力矩陣進(jìn)行晶體坐標(biāo)系到滑移系轉(zhuǎn)換，得

圖6 含孔洞R=12r的單晶鎢在應(yīng)變率108 s?1下壓縮變形過程中的原子快照

圖7 R=12r孔洞在應(yīng)變率108 s?1下壓縮變形時的位錯環(huán)結(jié)構(gòu)(對應(yīng)圖6中的c點)

不同{112}面上發(fā)射的1/2[111]位錯在一定條件下可發(fā)生如下位錯分解：

從而在初始滑移面上留下的不全位錯，所以當(dāng)1/2[111]螺位錯發(fā)生交滑移到相鄰的{112}面時，體系中留下一段孿晶。通過位錯提取算法分析發(fā)現(xiàn)，孔洞處發(fā)射位錯的伯氏矢量均為1/2[111]，當(dāng)位錯接觸到另一{112}滑移面時，發(fā)生交滑移，從而留下一段孿晶。圖8所示為圖6變形過程中d點的變形孿晶原子結(jié)構(gòu)圖。由圖可知圖7中的發(fā)生分解，形成了孿晶片段。由此判斷孔洞在受到[100]方向極高應(yīng)變率壓縮時，通過位錯和孿晶協(xié)同發(fā)生變形。

表3 R=12r孔洞模型在應(yīng)變率108 s?1下壓縮屈服點處不同滑移系的剪應(yīng)力

*The star indicates that this system is active at yielding.

3 結(jié)論

1) 對單晶鎢進(jìn)行壓縮變形，隨孔徑增大，孔洞處的應(yīng)力集中增大，導(dǎo)致鎢的壓縮屈服強(qiáng)度降低，同時由于孔洞處斷鍵數(shù)增加，導(dǎo)致材料的彈性模量減小。

2) 含有孔洞的單晶鎢的屈服強(qiáng)度隨應(yīng)變率增加而增大，呈現(xiàn)出很強(qiáng)的應(yīng)變率硬化行為，而彈性模量不隨應(yīng)變率變化而變化。

4) 在非晶化之前，孔洞的塑性變形是由位錯滑移和孿晶變形的共同作用造成的。

[1] 王星, 李樹奎, 王迎春, 等. 不同加載應(yīng)變率下鎢合金變形及破壞機(jī)理研究[J]. 北京理工大學(xué)學(xué)報, 2010, 30(11): 112?116. WANG Xing, LI Shukui, WANG Yingchun, et al. Study on deformation and failure mechanism of tungsten alloy under different loading strain rates[J]. Journal of Beijing Institute of Technology, 2010, 30(11): 112?116.

[2] 蔡國飆, 徐大軍. 高超聲速飛行器技術(shù)[M]. 北京: 科學(xué)出版社, 2012: 243?272. CAI Guobiao, XU Dajun. Hypersonic Vehicle Technology[M]. Beijing: Science Press, 2012: 243?272.

[3] 郭文啟, 劉金旭, 呂翠翠, 等. 90W-7Ni-3Fe合金拉伸力學(xué)行為的應(yīng)變率效應(yīng)[J]. 稀有金屬材料與工程, 2013, 42(4): 793?796. GUO Wenqi, LIU Jinxu, Lü Cuicui, et al. Strain rate effect of tensile mechanical behavior of 90W-7Ni-3Fe alloy[J]. Rare Metal Materials and Engineering, 2013, 42(4): 793?796.

[4] MA Y, ZHANG J, LIU W, et al. Microstructure and dynamic mechanical properties of tungsten-based alloys in the form of extruded rods via microwave heating[J]. International Journal of Refractory Metals & Hard Materials, 2014, 42: 71?76.

[5] 劉海燕, 宋衛(wèi)東, 寧建國. 不同晶粒度鎢合金動態(tài)力學(xué)性能研究[J]. 材料工程, 2007(6): 4?7. LIU Haiyan, SONG Weidong, NING Jianguo. Research on dynamic mechanical properties of tungsten alloys with different grain sizes[J]. Materials Engineering, 2007(6): 4?7.

[6] REMINGTON T P, REMINGTON B A, HAHN E N, et al. Deformation and failure in extreme regimes by high-energy pulsed lasers: A review[J]. Materials Science & Engineering A, 2017, 688(Complete): 429?458.

[7] BOURNE N K. Atomistic views of deformation[J]. Nature, 2017, 550(7677): 461?463.

[8] MILATHIANAKI D, BOUTET S, WILLIAMS G J, et al. Femtosecond visualization of lattice dynamics in shock- compressed matter[J]. Science, 2013, 342(6155): 220?223.

[9] ZEPEDA-RUIZ L A, STUKOWSKI A, OPPELSTRUP T, et al. Probing the limits of metal plasticity with molecular dynamics simulations[J]. Nature, 2017, 550(7677): 492?495.

[10] WEHRENBERG C E, MCGONEGLE D, BOLME C, et al. In situ X-ray diffraction measurement of shock-wave-driven twinning and lattice dynamics[J]. Nature, 2017, 550(7677): 496?499.

[11] PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of Computational Physics, 1995, 117(1): 1?19.

[12] ZHOU X, JOHNSON R, WADLEY H. Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers[J]. Physical Review B, 2004,69(14): 144113.

[13] REN Q Q, FAN J L, GONG H R. Work function and cohesion properties of W-Fe interfaces[J]. Materials Letters, 2015, 145: 205?208.

[14] MI S T, GONG H R, FAN J L. Structural stability and mechanical property of Fe-W solid solutions from a constructed Fe-W potential[J]. Journal of Applied Physics, 2019,126(11): 115102.

[15] WEI W, CHEN L, GONG H, et al. Strain-stress relationship and dislocation evolution of W-Cu bilayers from a constructed n-body W-Cu potential[J]. Journal of Physics Condensed Matter, 2019, 31(30): 305002.

[16] STUKOWSKI A, ALBE K. Extracting dislocations and non-dislocation crystal defects from atomistic simulation data[J]. Modelling & Simulation in Materials Science & Engineering, 2010, 18(8): 2131?2145.

[17] JR W D C. Materials Science and Engineering: An Introduction [M]. United States of America: John Wiley & Sons, 1994: 150? 196.

[18] DUPONT V, GERMANN T C. Strain rate and orientation dependencies of the strength of single crystalline copper under compression[J]. Physical Review B Condensed Matter, 2012, 86(13): 71?75.

[19] TANG Y, BRINGA E M, REMINGTON B A, et al. Growth and collapse of nanovoids in tantalum monocrystals[J]. Acta Materialia, 2011, 59(4): 1354?1372.

[20] TANG Y, BRINGA E M, MEYERS M A. Ductile tensile failure in metals through initiation and growth of nanosized voids[J]. Acta Materialia, 2012, 60(12): 4856?4865.

[21] REMINGTON T P, RUESTES C J, BRINGA E M, et al. Plastic deformation in nanoindentation of tantalum: A new mechanism for prismatic loop formation[J]. Acta Materialia, 2014, 78: 378? 393.

[22] BRINGA E M, TRAIVIRATANA S, MEYERS M A. Void initiation in fcc metals: Effect of loading orientation and nanocrystalline effects[J]. Acta Materialia, 2010, 58(13): 4458? 4477.

[23] DUPONT V, GERMANN T C. Strain rate and orientation dependencies of the strength of single crystalline copper under compression[J]. Physical Review B, 2012, 86(13): 134111.

Molecular dynamics simulation of void evolution and its mechanical response in single-crystal tungsten under extreme strain rate

BAI Xiaodan, CHEN Xun, MA Yunzhu, LIU Wensheng, LIANG Chaoping

(National Key Laboratory of Science and Technology for National Defense on High–Strength Structural Materials, Central South University, Changsha 410083, China)

extreme high strain rate; molecular dynamics simulation; void; plastic deformation; dislocation

TG146.4

A

1673-0224(2020)04-288-08

國家自然科學(xué)基金青年科學(xué)基金項目(51801236)；國家自然科學(xué)基金重點項目(51931012)

2020?06?01；

2020?06?17

梁超平，副教授，博士。電話：13357312107；E-mail: cpliang@csu.edu.cn

(編輯 湯金芝)