鋁多孔洞演變行為的分子動(dòng)力模擬

方 洲，劉曉波，徐慶軍

（南昌航空大學(xué) 航空制造工程學(xué)院，江西 南昌 330063）

1 引言

近五十年來(lái)，鋁已成為世界上最為廣泛應(yīng)用的金屬之一。但是，鋁很容易受到溫度，壓力等因素的影響；同時(shí)在制造過(guò)程中，存在著不可避免的各種缺陷，如縮孔、砂眼、氣孔、夾渣、微觀裂紋、雜質(zhì)、孔洞、刻痕、切口等，這些缺陷容易降低鋁材的質(zhì)量，并在宏觀上影響成品力學(xué)性能，甚至導(dǎo)致材料斷裂失效[1-2]。斷裂失效造成的事故往往使生命和財(cái)產(chǎn)蒙受巨大的損失。我們常?？吹降牧熊嚦鲕?、飛機(jī)失事、橋梁坍塌以及各種大小車禍等均是由材料缺陷演變的斷裂失效所致[3]。因此，掌握金屬內(nèi)部的裂紋、孔洞等缺陷的演變規(guī)律，對(duì)于預(yù)防類似事故的發(fā)生具有重要的指導(dǎo)作用。

隨著科學(xué)技術(shù)的快速發(fā)展，特別是計(jì)算機(jī)技術(shù)的發(fā)展，使需要大量計(jì)算的分子動(dòng)力學(xué)方法得以興起和廣泛應(yīng)用[4-8]。文獻(xiàn)[9]對(duì)鋁單晶中心貫穿微裂紋的愈合過(guò)程進(jìn)行了分子動(dòng)力學(xué)模擬。結(jié)果表明，當(dāng)加熱溫度超過(guò)臨界溫度，或外加壓應(yīng)力Kl超過(guò)臨界值時(shí)微裂紋將完全愈合。文獻(xiàn)[10]利用嵌入原子模型方法在型加載情況下，研究單晶Ni在有缺陷和沒(méi)有缺陷下的裂紋擴(kuò)展情況。獲得了裂紋尖端的臨界負(fù)載和應(yīng)變能的分布。分子動(dòng)力學(xué)模擬根據(jù)原子間勢(shì)函數(shù)不同的獲取方法，又可以分為經(jīng)典分子動(dòng)力學(xué)和現(xiàn)代分子動(dòng)力學(xué)。經(jīng)典分子動(dòng)力學(xué)是以牛頓力學(xué)為基礎(chǔ)來(lái)計(jì)算原子間的相互作用，而現(xiàn)代分子動(dòng)力學(xué)里的原子相互作用的計(jì)算是基于量子力學(xué)考慮了電子的影響。針對(duì)不同的研究對(duì)象和體系，兩種方法各有優(yōu)劣?？偟膩?lái)說(shuō)，面對(duì)簡(jiǎn)單體系及金屬材料時(shí)經(jīng)典分子動(dòng)力學(xué)模擬計(jì)算量更小同時(shí)與實(shí)驗(yàn)結(jié)果十分吻合，而面對(duì)復(fù)雜體系時(shí)，由于勢(shì)函數(shù)獲取更為準(zhǔn)確故常采用現(xiàn)代分子動(dòng)力學(xué)方法。分子動(dòng)力學(xué)方法是對(duì)某一體系內(nèi)某一時(shí)刻下所有已知位置和速度的粒子通過(guò)運(yùn)動(dòng)方程求解出粒子下一時(shí)刻的位置和速度的一種“確定性方法”，并由此推斷出在當(dāng)前條件下溫度、壓力、動(dòng)能、勢(shì)能等熱力學(xué)性質(zhì)。借助分子動(dòng)力學(xué)模擬鋁制品內(nèi)部的裂紋、孔洞等的演變可以使我們深入的了解復(fù)雜的機(jī)制，對(duì)于預(yù)防類似事故的發(fā)生具有積極的指導(dǎo)作用。便是在原子尺度下使用經(jīng)典分子動(dòng)力學(xué)方法在分子動(dòng)力學(xué)模擬軟件LAMMPS（Large-scale Atomic/Molecular Massively Parallel Simulator）中，對(duì)含有預(yù)制微觀孔洞的鋁在受速度加載拉伸作用時(shí)的變形機(jī)理進(jìn)行了研究。

2 計(jì)算模型

單晶鋁的結(jié)構(gòu)及物理特性，如表1所示。

表1 單晶Al結(jié)構(gòu)及物理特性Tab.1 The Structure and Physical Properties of Al

圖1 預(yù)制等間距孔洞分子動(dòng)力學(xué)初始模型Fig.1 Molecular Dynamics Initial Model with Prefabricate Equally Spaced Voids

在LAMMPS中根據(jù)表1數(shù)據(jù)建立模型，如圖1所示。x、y、z方向的尺寸分別為50a0×50a0×0.25a0，a0為晶格系數(shù)，將原子模型分為三個(gè)區(qū)域：上層、中間層、下層，在XY平面上，上層和下層對(duì)應(yīng)的區(qū)域大小為50a0×2a0，中間層的中間水平位置上預(yù)制4個(gè)相同大小的孔洞，孔洞為邊長(zhǎng)為2a0的正方形；同時(shí)4個(gè)孔洞的孔中心間距相同為L(zhǎng)=10a0，對(duì)應(yīng)的孔洞的坐標(biāo)，如表2所示。將系統(tǒng)初始溫度設(shè)為1K，低溫可以避免原子的熱激活效應(yīng)。本模型采用嵌入原子法進(jìn)行分子動(dòng)力學(xué)模擬，X和Y方向采用非周期性邊界條件，Z方向采用周期性邊界條件。以時(shí)間步長(zhǎng)為1fs對(duì)整個(gè)體系充分弛豫，達(dá)到平衡狀態(tài)。根據(jù)本模型的大小，將中間層原子和上層原子在Y方向上的初始速度設(shè)定為0.35A°/ps，X和Z方向的初始速度設(shè)置為0，下層原子設(shè)置為固定不動(dòng)。模擬的時(shí)間步長(zhǎng)設(shè)定為0.001ps，整個(gè)體系采用微正則系綜（NVE），整個(gè)程序運(yùn)行60000步，每間隔一千步記錄一次原子的坐標(biāo)位置、動(dòng)能、勢(shì)能和總能量。

表2 圖1模型中孔洞中心在x-y坐標(biāo)Tab.2 The X-Y Coordinates of Voids in the Fig.1

3 模擬的勢(shì)函數(shù)和算法

分子動(dòng)力學(xué)的對(duì)象是一個(gè)粒子系統(tǒng)，系統(tǒng)中的原子間的相互作用用勢(shì)函數(shù)來(lái)描述，因此采用EMA[11]勢(shì)函數(shù)進(jìn)行模擬，計(jì)算鋁原子之間的相互作用，系統(tǒng)的總勢(shì)能表示為：

式中：Fi—嵌入能函數(shù)；ρi—除了第i個(gè)原子以外所有原子在i處產(chǎn)生的電子云密度之和；φij—第i個(gè)原子與第j個(gè)原子之間的對(duì)勢(shì)作用函數(shù)；rij—第i原子與第j個(gè)原子之間的距離。運(yùn)用分子動(dòng)力學(xué)方法計(jì)算有N個(gè)粒子系綜的運(yùn)動(dòng)方程時(shí)，根據(jù)牛頓力學(xué)建立微分方程如下：

式中：mi—第i個(gè)原子的質(zhì)量；ri—第i個(gè)原子的原子坐標(biāo)；Fi—作用在第i個(gè)原子上的力；ai—第i個(gè)原子的加速度。

原子的初始位置以及初始速度為已知條件，利用微分方式（2），即可以求出任意時(shí)刻的原子位置和速度。再利用統(tǒng)計(jì)力學(xué)便可求解模型體系的平均動(dòng)能和勢(shì)能，而總能量為兩者之和。

模擬使用Velocity-Verlet算法進(jìn)行計(jì)算。其基本算法形式為：

Velocity-Verlet是通過(guò)原子在t時(shí)刻的位置R（t），加速度a（t）及 t-Δt時(shí)刻的位置計(jì)算出t+Δt時(shí)刻的位置R（t+Δt）。

4 等間距的多孔洞體系演變行為的結(jié)果與分析

首先，對(duì)模型進(jìn)行運(yùn)算，運(yùn)算完后對(duì)模擬結(jié)果采用可視化操作，以獲得等間距孔洞的整個(gè)原子運(yùn)動(dòng)的軌跡圖。提取了幾個(gè)特殊時(shí)刻下體系等間距孔洞的原子圖，如圖2所示。

體系運(yùn)行3000步時(shí)的原子圖，如圖2（a）所示。此時(shí)可以看到體系中編號(hào)2、3孔洞之間已經(jīng)產(chǎn)生微裂紋缺陷，編號(hào)2、3孔洞內(nèi)邊緣已經(jīng)開(kāi)始產(chǎn)生微小的滑移以及微裂紋，但是孔洞的大小與初始模型相比沒(méi)有明顯的變化，其它區(qū)域的原子沒(méi)有非常明顯的變化。

圖2 特殊時(shí)刻下的裂紋擴(kuò)展原子圖Fig.2 Under the Special Moment of Crack Extension Atomic Figure

12000步時(shí)的原子圖，如圖2（b）所示。圖中可以很明顯的看出，編號(hào)2、3孔洞的內(nèi)邊緣產(chǎn)生的裂紋的尖端原子出現(xiàn)了大面積的混亂排列，原子晶格位置發(fā)生重組；此時(shí)編號(hào)2、3孔洞之間的裂紋不斷產(chǎn)生子裂紋，子裂紋繼續(xù)擴(kuò)展又產(chǎn)生新的子裂紋，由于編號(hào)2、3孔洞內(nèi)邊緣裂紋擴(kuò)展產(chǎn)生的子裂紋在不斷的接近，導(dǎo)致編號(hào)2、3孔洞之間區(qū)域被裂紋不斷擴(kuò)展，此時(shí)編號(hào)2、3孔洞已經(jīng)快被裂紋貫穿。另外，編號(hào)1、2孔洞的內(nèi)邊緣開(kāi)始發(fā)射位錯(cuò)、產(chǎn)生微小的滑移及微裂紋。

21000步時(shí)的原子圖，如圖2（c）所示。此時(shí)非常明顯的變化是編號(hào)2、3孔洞已經(jīng)被其內(nèi)邊緣裂紋擴(kuò)展完全貫穿。同時(shí)隨著編號(hào)2、3孔洞內(nèi)邊緣產(chǎn)生的裂紋進(jìn)一步擴(kuò)展，裂紋尖端變得更加鈍化。上層與中間層的邊緣產(chǎn)生了撕裂，編號(hào)1、2孔洞內(nèi)邊緣產(chǎn)生裂紋不斷相互擴(kuò)展使其間距越來(lái)越??；編號(hào)3、4孔洞也相同。其他區(qū)域沒(méi)有明顯的變化。

41000步時(shí)的原子圖，如圖2（d）所示。此時(shí)產(chǎn)生明顯的變化是編號(hào)1與2、3與4孔洞由于內(nèi)邊緣裂紋擴(kuò)展已經(jīng)完全貫穿，此時(shí)編號(hào)1、2、3、4孔洞之間完全貫穿。另外體系的左右中間邊緣已經(jīng)開(kāi)始產(chǎn)生滑移并進(jìn)一步向內(nèi)擴(kuò)展，同時(shí)編號(hào)1、4孔洞內(nèi)邊緣產(chǎn)生的裂紋也進(jìn)一步向外擴(kuò)展，體系經(jīng)過(guò)一系列的循環(huán)過(guò)程：母裂紋擴(kuò)展—裂紋尖端鈍化—母裂紋產(chǎn)生子裂紋—子裂紋尖端鈍化—產(chǎn)生次子裂紋，最終加載體系被完全貫穿。

5 孔洞大小及孔間距對(duì)含多孔洞的鋁演變行為的影響

5.1 孔洞大小對(duì)含多孔洞的鋁演變行為的影響

為了探求孔洞的大小對(duì)體系的影響，在初始模型的基礎(chǔ)上改變預(yù)制孔洞的大小，分別設(shè)置了孔洞邊長(zhǎng)均為0.5a0、1.5a0、2a0、2.5a0四種情況下體系運(yùn)行的情況，將體系運(yùn)行的結(jié)果繪制為時(shí)間步數(shù)—總能量曲線圖，不同孔洞大小的體系總能量演變圖，如圖3所示。

圖3 不同孔洞大小的體系總能量演變圖Fig.3 System Energy Evolution with Different Voids

間距相等、孔洞邊長(zhǎng)均為0.5a0的能量曲線，如圖3（a）所示。能量曲線經(jīng)歷快速上升至最高點(diǎn)然后下降至一個(gè)平穩(wěn)的位置保持了大致10000步的時(shí)間步長(zhǎng)，然后能量又開(kāi)始下降。直到大概44000步后，能量曲線保持在一個(gè)穩(wěn)定水平至最后。

間距相等、孔洞邊長(zhǎng)均為1.5a0的能量曲線，如圖3（b）所示。能量曲線大致經(jīng)過(guò)12000步快速上升至最高點(diǎn)。然后能量曲線下降過(guò)程中又經(jīng)歷了兩次輕微下降并且后兩次的能量的下降都是先保持一定的時(shí)間步后才繼續(xù)下降。最后隨著體系繼續(xù)下降，在40000步左右時(shí)體系被完全拉開(kāi)，能量曲線保持在一個(gè)穩(wěn)定水平至最后。

間距相等、孔洞邊長(zhǎng)均為2a0的能量曲線，如圖3（c）所示。能量曲線先快速上升至最高點(diǎn)，此時(shí)時(shí)間步大致為12000步，然后隨著體系進(jìn)一步加載，能量曲線開(kāi)始下降。對(duì)孔洞邊長(zhǎng)為2a0的能量曲線進(jìn)行觀察，我們發(fā)現(xiàn)曲線在下降區(qū)域內(nèi)變化較為平緩，沒(méi)有出現(xiàn)如孔邊長(zhǎng)為0.5a0、1.5a0曲線那樣的平穩(wěn)一段時(shí)間后再下降的情況。在時(shí)間步為37000步左右時(shí)體系被完全拉開(kāi)，能量曲線保持在一個(gè)穩(wěn)定水平至最后。

間距相等、孔洞邊長(zhǎng)均為2.5a0的能量曲線，如圖3（d）所示。能量曲線同樣是先快速上升至最高點(diǎn)，此時(shí)時(shí)間步約為12000步。然后隨著體系進(jìn)一步加載能量曲線開(kāi)始下降，對(duì)比其他三個(gè)孔的曲線，發(fā)現(xiàn)孔邊長(zhǎng)為2.5a0的曲線的下降過(guò)程最為平緩。在體系加載到時(shí)間步為35000時(shí)體系被完全拉開(kāi)，能能量曲線保持在一個(gè)穩(wěn)定水平至最后。

綜合分析圖3不同孔洞大小的體系總能量演變圖，不同的孔洞大小對(duì)應(yīng)的能量曲線都是在12000步左右時(shí)體系能量達(dá)到相應(yīng)的最高點(diǎn)。但是孔洞的邊長(zhǎng)越大，對(duì)應(yīng)的孔洞的整體能量水平越低，能量曲線在下降的過(guò)程表現(xiàn)為越來(lái)越平穩(wěn)沒(méi)有波動(dòng)?？梢酝茢嚯S著模型中4個(gè)等間距、等大小的孔洞邊長(zhǎng)變大，孔洞邊緣發(fā)射的位錯(cuò)與相鄰孔洞以及模型左右兩側(cè)邊界發(fā)射的位錯(cuò)匯合的概率更大，所需的能量也越低，孔洞的演變路徑在Y方向也更加規(guī)則，因此從孔洞到達(dá)最高點(diǎn)開(kāi)始至體系完全被拉開(kāi)的區(qū)間內(nèi)能量曲線趨于平穩(wěn)下降。最后，隨著孔洞邊長(zhǎng)的增加，體系被拉開(kāi)所用的時(shí)間步更少。

5.2 孔間距對(duì)含多孔洞的鋁演變行為的影響

為了探究不同的孔間距對(duì)含預(yù)制孔洞的體系演變的影響，本節(jié)分別設(shè)置初始模型中的L為6a0、8a0、10a0對(duì)比在這三種孔洞等間距情況下體系的演變規(guī)律，其他條件與初始模型一致。得到能量曲線圖，如圖4所示。

圖4 不同孔間距下的能量演變圖Fig.4 Energy Evolution Under Different Void Spacing

對(duì)于圖4中曲線a：其先經(jīng)歷一個(gè)上升階段在12000步左右時(shí)達(dá)到能量最高點(diǎn)，然后開(kāi)始下降，直到40000步左右時(shí)體系被完全拉開(kāi)，此后能量曲線處于一個(gè)穩(wěn)定位置。對(duì)比b、c曲線，a曲線的最高點(diǎn)最高，下降的階段是十分平穩(wěn)沒(méi)有明顯的波動(dòng)。

對(duì)于圖4中曲線b：其先經(jīng)歷一個(gè)上升階段在12000步左右時(shí)達(dá)到能量最高點(diǎn)，對(duì)比c曲線最高點(diǎn)，b曲線的最高點(diǎn)大于c曲線的最高點(diǎn)。隨后曲線開(kāi)始下降，直到40000步左右時(shí)體系被完全拉開(kāi)，此后能量開(kāi)始處于一個(gè)穩(wěn)定的位置。但在它的下降階段，曲線出現(xiàn)了一個(gè)波動(dòng)后又繼續(xù)平穩(wěn)下降。

對(duì)于圖4中曲線c：其先經(jīng)歷一個(gè)上升階段在12000步左右時(shí)達(dá)到能量最高點(diǎn)，然后開(kāi)始下降，直到40000步左右時(shí)體系被完全拉開(kāi)，此后能量曲線處于一個(gè)穩(wěn)定的位置。但在它的下降階段，曲線出現(xiàn)了多個(gè)較大的能量波動(dòng)并有短暫的穩(wěn)定階段。

綜合分析曲線a、b、c；第一，它們都是在12000步左右時(shí)被加載到能量最高值，此刻模型中4個(gè)孔洞剛開(kāi)始演變；但隨著孔洞間距L的減小，a、b、c曲線的最高值依次減小。這是因?yàn)殡S著4個(gè)孔洞之間的三處間隔區(qū)域直線距離L減小，使得孔洞之間的區(qū)域整體原子結(jié)合力減小從而導(dǎo)致體系內(nèi)孔洞（因?yàn)轶w系內(nèi)含4個(gè)孔洞，各個(gè)孔洞開(kāi)始演變的時(shí)間有間隔）開(kāi)始演變所需要吸收的能量減小，即間距越大越難被拉開(kāi)。第二，隨著孔洞間距的減小，能量曲線在下降階段的波動(dòng)數(shù)量越來(lái)越多，幅度越來(lái)越大，變化趨于復(fù)雜。這是因?yàn)榭锥磧?nèi)邊緣原子間結(jié)合力是非常弱的，所以導(dǎo)致孔洞的內(nèi)邊緣在孔洞剛開(kāi)始演變時(shí)會(huì)出現(xiàn)很多微小位錯(cuò)以及微小裂紋，孔洞之間的間隔越近就越容易出現(xiàn)彼此內(nèi)邊緣發(fā)射的位錯(cuò)及產(chǎn)生的微裂紋相互匯集從而導(dǎo)致體系能量出現(xiàn)多個(gè)波動(dòng)。

6 結(jié)論

建立含多孔洞的鋁的分子動(dòng)力學(xué)模擬模型，通過(guò)觀察特殊時(shí)刻下的原子運(yùn)動(dòng)軌跡圖，詳細(xì)的說(shuō)明了特殊時(shí)刻體系的裂紋擴(kuò)展規(guī)律和孔洞的生成及匯集等情況，對(duì)于等間距的多孔洞鋁的研究表明，隨著孔洞直徑的增大，其越容易與相鄰的孔洞匯合，體系會(huì)更快、更容易被拉開(kāi)。在孔洞尺寸相同的情況下，孔洞間距越大即體系內(nèi)孔洞分布越均勻，孔洞開(kāi)始演變所需要吸收的能量也越多，也就是含有均勻分布的孔洞的模型更難被拉開(kāi)。這也表示集中分布的孔缺陷比均勻分布的孔缺陷危害更大。同時(shí)，相鄰孔洞內(nèi)邊緣產(chǎn)生的微裂紋、發(fā)射的位錯(cuò)更容易因?yàn)殚g隔較大互相引導(dǎo)發(fā)展成為一條大的、規(guī)則裂紋，從而使體系被拉開(kāi)時(shí)的演變十分平穩(wěn)。

［1］盧功輝.金屬型鋁合金鑄件裂紋產(chǎn)生原因及對(duì)策［J］.特種鑄造及有色合金，2006，26（5）：308-310.（Lu Gong-hui.Metal type aluminum alloy casting crack causes and countermeasures［J］.Special Casting and Nonferrous Alloys，2006，26（5）：308-310.）

［2］陳明，李革，張文飛.納米尺度孔洞周圍應(yīng)力集中現(xiàn)象分析［J］.材料導(dǎo)報(bào) B：研究篇，2011，25（3）：131-134.（Chen Ming，Li Ge，Zhang Wen-fei.Discussion on the phenomenon of stress concentration around hole in nano-scale［J］.Material Review B：Research，2011，25（3）：131-134.）

［3］張俊善.材料強(qiáng)度學(xué)［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2004：107-109.（Zhang Jun-shan.Strength of Engineering Materials［M］.Harbin：Harbin Institute of Technology Press，2004：107-109.）

［4］M.M aki-Jaskari，K.Kaski，A.Kuronen.Simulations of crack initiation in silicon［J］.Computational Materials Science，2000（17）：336-342.

［5］W.P.Wu，Z.Z.Yao.Molecular dynamics simulation of stress distribution andmicrostructureevolutionaheadofagrowingcrackinsinglecrystalnickel［J］.Theoretical and Applied Fracture Mechanics，2013（62）：67-75.

［6］Aude Simar，Hyon-Jee Lee Voigt，Brian D.Wirth.Molecular dynamics simulations of dislocation interaction with voids in nickel［J］.Computational Materials Science，2011（50）：1811-1817.

［7］劉曉波，徐慶軍，劉劍.鋁裂紋擴(kuò)展行為的分子動(dòng)力學(xué)模擬［J］.中國(guó)有色金屬學(xué)報(bào)，2014，24（6）：1408-1413.（Liu Xiao-bo，Xu Qing-jun，Liu Jian.Molecular dynamics simulation of the aluminum crack propagation behavior［J］.Transactions of Nonferrous Metals Society of China，2014，24（6）：1408-1413.）

［8］余震，陳定方，朱宏輝.基于分子動(dòng)力學(xué)的微納米切削模擬研究進(jìn)展［J］.機(jī)械設(shè)計(jì)與制造，2009（11）：217-219.（Yu Zhen，Chen Ding-fang，Zhu Hong-hui.Development of research in micro-nanometric cutting simulation based on molecular dynamics［J］.Machinery Design&Manufacture，2009（11）：217-219.）

［9］周國(guó)輝，高克瑋，萬(wàn)發(fā)榮.微裂紋愈合過(guò)程的分子動(dòng)力學(xué)模擬［J］.自然科學(xué)進(jìn)展，2001（3）：300-305.（Zhou Guo-hui，Gao Ke-wei，Wan Fa-rong.Molecular dynamics simulation of micro crack healing process［J］.Progress in Natural Science，2001（3）：300-305）

［10］Majid Karimi，Tom Roarty，Theodore Kaplan.Molecular dynamics simulations of crack propagation in Ni with defects ［J］.Modelling Simul.Mater.Sci.Eng.，2006（14）：1409-1420.

［11］Ackland G J，Vitek V Many-body potentials and atom-scale relaxations innoble-metalalloys［J］.PhysicalReviewB，1990，41（15）：10324-10333.