基于修正球形雙晶模型的金屬Al晶界能分子動力學(xué)計算

張明亮，楊 亮，魏承煬，李賽毅, ，張新明,

(1. 中南大學(xué) 材料科學(xué)與工程學(xué)院，長沙 410083；2. 中南大學(xué) 有色金屬材料科學(xué)與工程教育部重點實驗室，長沙 410083；3. 中南大學(xué) 有色金屬先進結(jié)構(gòu)材料與制備協(xié)同創(chuàng)新中心，長沙 410083)

作為金屬多晶材料內(nèi)相鄰晶粒間的交界區(qū)，晶界(GB)是材料顯微組織結(jié)構(gòu)的重要組成部分[1]。完整描述晶界的空間取向需要5個參數(shù)(或自由度)[2]，其中3個用來確定兩相鄰晶粒間的取向差，另兩個表示晶界面取向。晶界因存在與其空間取向相關(guān)的原子錯排而具有比晶粒更高的能量，單位面積晶界所具有的附加自由能即為晶界能[3]。晶界能作為晶界的重要特性之一，影響著材料再結(jié)晶和晶粒長大等過程中的組織演變，進而影響材料的力學(xué)及物理性能[4-10]。因此，有關(guān)晶界能的定量表征及其與晶界參數(shù)相關(guān)性的研究一直是材料科學(xué)領(lǐng)域的重要課題。通過實驗定量測定晶界能時，所得結(jié)果受到材料純度和晶界參數(shù)精度等的嚴(yán)重影響，目前僅獲得了極少數(shù)的晶界能結(jié)果[3,11-14]。READ等[15]基于晶界的位錯模型，推導(dǎo)出晶界能與取向差角的關(guān)系式，在理論上描述了晶界能受該晶界參數(shù)的影響，能夠較好地應(yīng)用于小角度晶界的能量計算。

近年來，隨著計算機模擬技術(shù)的發(fā)展，基于分子動力學(xué)的晶界能計算方法逐漸成熟，為從原子尺度探究晶界能及其各向異性提供了有效手段。該類方法的基本思想是：首先根據(jù)晶界參數(shù)構(gòu)造具有一定邊界條件的雙晶體系，然后選取合適的原子勢描述原子間的相互作用，并在絕對零度下采用經(jīng)典牛頓方程對體系進行弛豫。當(dāng)能量降低至平衡時，晶界所具有的附加能量即為晶界能。用于構(gòu)造雙晶體系的模型主要有考慮周期性邊界的塊狀(Block)模型和考慮非周期性邊界的球形(Sphere)模型。塊狀模型需保證晶界面內(nèi)兩個方向的模型尺寸與晶體原子排布的周期長度相適應(yīng)，而對于任意取向參數(shù)的晶界這點很難或無法滿足。所以，該模型主要被應(yīng)用于一些特殊取向的晶界，如重位因子較低的傾斜晶界和扭轉(zhuǎn)晶界[16-22]。為了克服塊狀模型的局限性，LEE等[23]提出了球形模型。該模型采用非周期性邊界(即自由表面)，可實現(xiàn)任意晶界的晶界能計算。但是，球形模型僅通過不同取向晶粒的剛性對接來確定晶界結(jié)構(gòu)[23]，這種對接可能使晶界處的原子相距太近甚至重疊，從而使所構(gòu)造的晶界處于非穩(wěn)態(tài)或亞穩(wěn)態(tài)，對應(yīng)的晶界能將高于理論值。對于給定參數(shù)的晶界，其結(jié)構(gòu)和能量都無法預(yù)知。因此，需要嘗試構(gòu)造大量初始晶界結(jié)構(gòu)并比較它們的能量，其中能量最低的結(jié)構(gòu)可以當(dāng)成該晶界參數(shù)所描述的晶界結(jié)構(gòu)，基于該結(jié)構(gòu)所計算的附加能量才是晶界能。塊狀模型[21-22]通過引入原子刪除操作來刪除晶界處不合理原子，能夠成功獲取初始晶界結(jié)構(gòu)。已有球形模型則尚未考慮原子刪除操作，無法排除不合理初始晶界結(jié)構(gòu)對計算結(jié)果的影響。

本文作者在基本球形雙晶建模方法的基礎(chǔ)上，引入已在塊狀模型中獲得成功應(yīng)用的原子刪除操作來修正原始球形模型，并以金屬 Al 〈011〉對稱傾斜晶界(STGB)為例，對比分析原始球形模型、修正球形模型和塊狀模型所算的晶界能，考察原子刪除操作對晶界能的影響。同時，將修正模型應(yīng)用于Al〈001〉扭轉(zhuǎn)晶界(TWGB)和〈011〉非對稱傾斜晶界(ATGB)的晶界能計算，以探討取向差角和晶界面傾斜角對晶界能的影響。最后，基于上述3類晶界的計算結(jié)果，考察晶界能與重位因子的相關(guān)性。

圖1 任意晶界的球形雙晶體系構(gòu)造過程示意圖Fig. 1 Schematic diagram showing construction of spherical bicrystal with an arbitrary GB

1 模型與方法

本研究借助 LAMMPS軟件構(gòu)造晶界能的計算模型[24]。圖1所示為球形雙晶體系構(gòu)造過程的示意圖。首先，在參考坐標(biāo)系下構(gòu)建一個直徑為D的球形初始單晶(樣品Ⅰ)，將初始單晶分別按照取向gA、gB進行旋轉(zhuǎn)分別得到單晶A和B，使它們滿足給定的取向差關(guān)系 g = gA；然后，沿給定的晶界面方位分別將A和B剖成兩個完全相同的半球，保留不同側(cè)的半球形晶粒 A′和 B′；最后，將 A′和 B′沿剖面剛性對接，組成一個球形雙晶(樣品Ⅱ)。在0 K下對單晶和雙晶進行能量最小化弛豫，弛豫后體系的總能量分別為E1和E2。由于單晶與雙晶表面原子的排布總體相同，二者所具有的表面能可認為相等[23]。所以，相對于單晶而言，雙晶所具有的額外能量為晶界的能量，即晶界能 γ0為

式中：S是晶界面積( S = π D2/4)。

由圖1可知，通過半球形晶粒A′和B′沿剖面剛性對接而構(gòu)造的雙晶體系，可能在晶界處存在相距過近或重疊的原子。從物理角度看，這些原子不能存在于實際材料的晶界結(jié)構(gòu)中。塊狀模型的應(yīng)用研究[21-22]表明，此類不合理原子的存在將導(dǎo)致晶界結(jié)構(gòu)處于非穩(wěn)態(tài)或亞穩(wěn)態(tài)，從而使晶界能量異常升高，即所得晶界能將高于其理論值。因此，需合理刪除這些原子，以獲取穩(wěn)態(tài)的晶界結(jié)構(gòu)。

本研究參照塊狀雙晶模型中的原子刪除操作[21]來刪除不合理的原子，并構(gòu)造不同的初始晶界結(jié)構(gòu)以獲取給定參數(shù)晶界的結(jié)構(gòu)和晶界能。該操作具體過程：先設(shè)定一臨界原子間距(dc)識別晶界處不合理的原子對(即間距小于 dc的相鄰兩原子)；然后刪除原子對中的任意一個來獲取晶界結(jié)構(gòu)，再通過改變dc以獲取不同的初始晶界結(jié)構(gòu)，最后選定最穩(wěn)態(tài)的晶界結(jié)構(gòu)(對應(yīng)最低晶界能量的結(jié)構(gòu))為給定參數(shù)晶界的結(jié)構(gòu)。對于具有面心立方結(jié)構(gòu)的金屬，晶體理想點陣中原子間的最近鄰距離為02 /2a (a0為晶格常數(shù))。為了避免刪除體系中間距合理的原子，同時又盡可能多地獲取不同的初始晶界結(jié)構(gòu)，本研究在 0.005a0到 0.7a0之間以0.005a0為增量依次選取dc的取值。其中，對于晶界處間距小于dc同時又分別屬于晶粒A′和B′的兩個原子，選取不同原子進行刪除將得到兩種不同的晶界結(jié)構(gòu)。因此，修正模型在原理上可構(gòu)造晶界參數(shù)相同的初始晶界結(jié)構(gòu)總數(shù)為2×(0.7-0.005)/0.005=278。顯然，進行了原子刪除操作的雙晶和初始單晶的原子總數(shù)不再相等，與雙晶具有同等原子總數(shù)單晶的能量也不再為E1。此時，對應(yīng)單晶的能量為雙晶原子總數(shù)N和平均原子勢能e1的乘積。因此，利用修正模型所構(gòu)造晶界的晶界能(dγ)計算公式為

式中：2E′為引入原子刪除操作后的雙晶的能量。

在計算某給定參數(shù)晶界的能量時要對上述278種初始晶界結(jié)構(gòu)所對應(yīng)的雙晶體系進行能量最小化弛豫，并按照式(2)計算所對應(yīng)的能量。其中，能量最低值被認為是該晶界參數(shù)下的晶界能，相應(yīng)結(jié)構(gòu)為全局最穩(wěn)態(tài)結(jié)構(gòu)，即給定晶界參數(shù)晶界的結(jié)構(gòu)。

本研究利用修正模型考察了金屬 Al中的〈011〉STGBs、〈011〉A(chǔ)TGBs和〈001〉TWGBs 3 類晶界，這些晶界對應(yīng)的取向差角(θ)、晶界面傾斜角(Ф)和重位因子(Σ)見表1。所構(gòu)造雙晶模型的直徑為D=10 nm。選定二次近鄰修正嵌入式原子勢[25-26]來描述原子的能量和相互作用力，并利用共軛梯度算法進行體系能量最小化弛豫。

表1 本研究中晶界所對應(yīng)的取向差角、晶界面傾斜角及重位因子Table 1 List of misorientation angle, inclination angle and coincidence factor of GBs considered in this work

2 結(jié)果與討論

2.1 原子刪除對晶界能的影響

圖2所示為分別采用原始球形模型、修正球形模型和塊狀模型計算得到的〈011〉對稱傾斜晶界的晶界能(分別用 γ0、γd和 γb表示)隨取向差角的變化。由圖2可以看出，3種模型所得晶界能隨取向差角的變化趨勢整體上相同，都在70.53°和129.52°處形成兩個顯著的極小值，該趨勢與實驗所測Al的晶界能結(jié)果一致[3]。70.53°和 129.52°兩晶界剛好對應(yīng)(111)和(113)孿晶界，晶界處原子錯排程度較低，所對應(yīng)的附加能量也隨之較低，即晶界能較低。圖2的結(jié)果還表明，在16.10°、38.94°和109.47°附近，原始模型和修正模型所得晶界能在絕對值大小和隨取向差角的變化趨勢上存在明顯差異。在這幾個晶界處修正模型所得晶界能呈現(xiàn)出明顯的極小值，而原始模型的結(jié)果卻呈現(xiàn)出明顯的極大值，且比前者的分別高出40%以上。修正模型所得結(jié)果與塊狀模型結(jié)果一致。

圖2 采用原始球形模型、修正球形模型和塊狀模型計算所得〈011〉對稱傾斜晶界的晶界能(分別用 γ0、γd和 γb表示)隨取向差角的變化Fig. 2 Variation of energy of 〈011〉 STGBs with misorientation angle calculated using original sphere model (γ0),modified sphere model (γd) and block model (γb)

晶界能本質(zhì)上由原子排列決定，上述原始模型和修正模型在16.10°、38.94°和109.47°所得晶界能的差別可從晶界處原子排布的不同來分析。以 θ=109.47°為例，圖3所示為不合理原子被刪除前后的晶界原子排布，分別對應(yīng)于原始模型和修正模型所獲得初始晶界結(jié)構(gòu)。如圖 3(a)所示，在原始模型構(gòu)造的晶界處存在一些相距過近的原子(圖3(a)中橢圓內(nèi)的原子)。根據(jù)Lennard-Jones兩體勢[27]可知，在一次近鄰范圍內(nèi)，原子勢能隨原子間距的減小而急劇升高。因此，晶界處間距過小的原子將具有較高的能量，直接導(dǎo)致晶界能升高。由圖3(b)可看出，在修正模型構(gòu)造的雙晶體系中，存在于原始模型構(gòu)造的雙晶體系中的不合理原子(即具有較高能量的原子)被刪除。因此，與原始模型相比，修正模型所構(gòu)造的θ=109.47° 晶界具有更穩(wěn)定的結(jié)構(gòu)和更低的晶界能。

圖2結(jié)果還表明，修正模型與塊狀模型所算晶界能在絕對值上仍然存在可見的差別。例如，修正球形模型所算結(jié)果相對于塊狀模型的平均誤差和相對平均誤差分別為14 mJ/m2和4%。這可能與兩個方面有關(guān)。首先，與塊狀模型相比，修正球形模型在計算給定取向參數(shù)晶界的能量時，未對相鄰兩晶體在晶界面內(nèi)進行剛體平移操作以獲取盡可能多的初始晶界結(jié)構(gòu)。與省略原子刪除操作類似，省略剛體平移操作也可能導(dǎo)致球形模型計算的晶界能并非全局最小值，從而產(chǎn)生一定的誤差。其次，球形模型考慮非周期性邊界條件，雙晶體系中的自由表面勢必會影響晶界結(jié)構(gòu)，導(dǎo)致晶界能的計算誤差[23]。

圖 3 基于球形模型所構(gòu)造的 109.47°〈011〉對稱傾斜晶界和初始晶界結(jié)構(gòu)對比Fig. 3 Initial configurations generated before (a) and after (b)deletion of unreasonable atoms in sphere model for 109.47°〈011〉 STGB (Misorientation axis is perpendicular to plane of paper; circles and triangles stand for atoms in (011) plane in different layers in crystals; blue and red symbols stand for crystal A′ and B′, respectively)

2.2 晶界能與取向差角、晶界面傾斜角和重位因子的相關(guān)性

圖4所示為〈001〉扭轉(zhuǎn)晶界的晶界能隨取向差角的變化情況。由圖 4可知，當(dāng)取向差角較小時(θ＜28.07°)，晶界能受取向差角的影響較大，隨取向差角的增大呈現(xiàn)出升高的趨勢，與 Read-Shockley公式所體現(xiàn)的趨勢相符；當(dāng)取向差角繼續(xù)增大時，晶界能基本保持在一個較小的范圍內(nèi)波動?；诰Ы绲奈诲e模型及相關(guān)理論可知[15]，當(dāng)取向差角較小時，晶界由離散的位錯構(gòu)成，且位錯密度(或原子錯排程度)隨著取向差角的增大而升高，從而導(dǎo)致晶界能隨取向差角的增大而升高。當(dāng)取向差角增大到一定值時，構(gòu)成晶界的位錯相互重疊，位錯密度有所上升。然而，當(dāng)取向差角進一步增大時，由位錯重疊和位錯間相互作用導(dǎo)致的位錯湮滅將在一定程度上抵消由于取向差角增大所致的位錯密度的升高，位錯密度將不再發(fā)生明顯變化，即晶界能不再升高而保持在一個較小的范圍內(nèi)波動。

圖 5所示為 Σ3〈011〉非對稱傾斜晶界的晶界能隨晶界面傾斜角Ф的變化情況。結(jié)果顯示，晶界能在Ф為0°～64.76°范圍內(nèi)隨Ф的增大而逐漸升高，在74.21°附近達到最大值，隨后緩慢降低。這種趨勢與已有金屬 Al的模擬結(jié)果[21]以及金屬 Cu的相關(guān)實驗結(jié)果[11]相符。Ф為0°和90°分別對應(yīng)于對稱共格孿晶界和對稱非共格孿晶界，晶界面兩側(cè)的原子關(guān)于晶界面成鏡面對稱排布。與相鄰晶界相比，這兩個晶界的原子錯排程度更低，對應(yīng)的晶界能也更低。當(dāng)晶界面偏離上述對稱孿晶界時，兩側(cè)晶粒在晶界上的原子排布不同(例如 Ф=70.53°時，晶界面在兩側(cè)晶粒中的取向分別為{111}和{115})，即呈現(xiàn)非對稱排布。此時，晶界處原子錯排度和對應(yīng)晶界能都將升高。

圖 4 利用修正球形模型所算〈001〉扭轉(zhuǎn)晶界的晶界能隨取向差角的變化Fig. 4 Variation of energy of 〈001〉 TWGBs with misorientation angle, calculated using modified sphere model

圖5 利用修正球形模型所算Σ3 (70.53°)〈011〉非對稱傾斜晶界的晶界能隨晶界面傾斜角的變化Fig. 5 Variation of energy of Σ3 (70.53°) 〈011〉 ATGBs with inclination angle calculated using modified sphere model

重位因子Σ是晶體點陣中重合點陣占體系總點陣比例的倒數(shù)，在一定程度上體現(xiàn)晶界結(jié)構(gòu)的特點。早期研究[3,28]認為，Σ值越小，晶界所具有的能量就越低。為了考察晶界能與重位因子的相關(guān)性，圖6所示為本研究中3類晶界的晶界能及對應(yīng)的重位因子。由圖6可以看出，高重位因子晶界的晶界能并不一定高于低重位因子晶界的，如Σ＜10晶界的能量并非都低于(甚至有時高于)Σ＞10晶界的能量。而且，同一 Σ值下(如Σ3)，不同晶界所具有的能量并不相同，甚至存在很大差異。由此看來，晶界能與Σ之間并不存在特定的相關(guān)性，無法通過Σ簡單預(yù)測或評判晶界能的相對大小。OLMSTED等[22]采用塊狀模型計算并分析了388種重位點陣晶界的能量與重位因子的關(guān)系，得到了與本研究一致的結(jié)論。

圖 6 利用修正球形模型計算的〈011〉對稱傾斜晶界、〈001〉扭轉(zhuǎn)晶界和〈011〉非對稱傾斜晶界的晶界能隨重位因子(Σ)的變化Fig. 6 Variation of energy of 〈011〉 STGBs, 〈001〉 TWGBs and 〈011〉 ATGBs with coincidence factor Σ according to predictions using modified sphere model

3 結(jié)論

1) 與原始球形模型相比，修正模型所構(gòu)造的〈011〉對稱傾斜晶界因刪除了不合理的原子而更加穩(wěn)定，所獲得的晶界能及其隨取向差角的變化趨勢與已有實驗和模擬結(jié)果相符更好。

2) 對于〈001〉扭轉(zhuǎn)晶界，在取向差角較小時晶界能隨取向差角的增大而升高，該趨勢與 Read-Shockley晶界能計算公式所預(yù)測的結(jié)果一致；當(dāng)取向差角繼續(xù)增大時，晶界能保持在較小的范圍內(nèi)波動。對于 Σ3〈011〉非對稱傾斜晶界，晶界能隨晶界面傾斜角的增大而呈現(xiàn)出先升高后減小的趨勢，在對稱孿晶界處的晶界能低于相鄰晶界的。

3) 〈001〉扭轉(zhuǎn)晶界、〈011〉對稱傾斜晶界和〈011〉非對稱傾斜晶界的晶界能結(jié)果表明，晶界能與重位因子間不存在特定的相關(guān)性，不能基于重位因子預(yù)測晶界能的相對大小。

[1] HU Hui-e, YANG Li, ZHEN Liang, SHAO Wen-zhu, ZHAN Bao-you. Relationship between boundary misorientation angle and true strain during high temperature deformation of 7050 aluminum alloy[J]. Transactions of Nonferrous Metals Society of China, 2008, 18(4): 795-798.

[2] BRANDON D G. The structure of high-angle grain boundaries[J]. Acta Metallurgica, 1966, 14(11): 1479-1484.

[3] HASSON G C, GOUX C. Interfacial energies of tilt boundaries in aluminium: Experimental and theoretical determination[J].Scripta Metallurgica, 1971, 5(10): 889-894.

[4] CHEN Da, LU Min, LIN Dong-liang. Molecular static simulation of energy features of interaction between grain boundary and dislocations in Ni3Al alloy[J]. Transactions of Nonferrous Metals Society of China, 1994, 4(3): 67-70.

[5] 盧 山, 李子然, 昝 祥. 晶界對近片層 TiAl高溫動態(tài)力學(xué)行為影響的數(shù)值模擬[J]. 中國有色金屬學(xué)報, 2012, 22(2):379-387.LU Shan, LI Zi-ran, ZAN Xiang. Numerical simulation of dynamic mechanical behavior of near lamellar TiAl at elevated temperature with influence of grain boundary[J]. The Chinese Journal of Nonferrous Metals, 2012, 22(2): 379-387.

[6] 毛衛(wèi)民, 陳 壘, 薩麗曼, 余永寧, 李云峰. 晶界對低壓電解電容器鋁箔腐蝕結(jié)構(gòu)的影響[J]. 中國有色金屬學(xué)報, 2004,14(1): 1-5.MAO Wei-min, CHEN Lei, SA Li-man, YU Yong-ning, LI Yun-feng. Influence of grain boundaries on corrosion structure of low voltage aluminum foil[J]. The Chinese Journal of Nonferrous Metals, 2004, 14(1): 1-5.

[7] YAN Wen, CHEN Jian, FAN Xin-hui. Effects of grain boundaries on electrical property of copper wires[J].Transactions of Nonferrous Metals Society of China, 2003, 13(5):1075-1079.

[8] 陳健美, 周卓夫, 唐建國, 張新明. 一種考慮晶界能各向異性模擬晶粒長大的Monte Carlo方法[J]. 材料導(dǎo)報, 2003, 17(8):77-79.CHEN Jan-mei, ZHOU Zhuo-fu, TANG Jian-guo, ZHANG Xin-ming. A monte carlo simulation of grain growth with grain boundary energy anisotropy considered[J]. Materials Review,2003, 17(8): 77-79.

[9] 李 旭, 周 慶, 陳明和, 王小芳. 基于晶界遷移率與晶界能各向異性的晶粒生長元胞自動機模擬[J]. 機械工程材料,2012, 36(7): 82-87.LI Xu, ZHOU Qing, CHEN Ming-he, WANG Xiao-fang.Cellular automata simulation for grain growth based on anisotropic grain boundary mobility and grain boundary energy[J]. Materials for Mechanical Engineering, 2012, 36(7):82-87.

[10] 魏承煬. 再結(jié)晶織構(gòu)的相場模擬與機理研究[D]. 廣州:華南理工大學(xué), 2012: 49-58.WEI Cheng-yang. Phase field modelling of texture evolution and its mechanisms during recrystallization annealing[D].Guangzhou: South China University of Technology, 2012:49-58.

[11] WOLF U, ERNST F, MUSCHIK T, FINNIS M W. The influence of grain boundary inclination on the structure and energy of Σ3 grain boundaries in copper[J]. Philosophical Magazine A, 1992, 66(6): 991-1016.

[12] GJOSTEIN N A, RHINES F N. Absolute interfacial energies of[001] tilt and twist grain boundaries in copper[J]. Acta Metallurgica, 1959, 7(5): 319-330.

[13] MORI T, MIURA H, TOKITA T, HAJI J, KATO M.Determination of the energies of [001] twist boundaries in Cu with the shape of boundary SiO2particles[J]. Philosophical Magazine Letters, 1988, 58(1): 11-15.

[14] MORI T, ISHII T, KAJIHARA M, KATO M. Determination of the energies of [001] symmetric tilt boundaries in Cu from the shape of boundary SiO2[J]. Philosophical Magazine Letters,1997, 75(6): 367-370.

[15] READ W T, SHOCKLEY W. Dislocation models of crystal grain boundaries[J]. Physical Review, 1950, 78(3): 275-289.

[16] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅰ. Boundaries on the (111) and (100) planes[J].Acta Metallurgica, 1989, 37(7): 1983-1993.

[17] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅱ. Boundaries on the (110) and (113) planes[J].Acta Metallurgica, 1989, 37(10): 2823-2833.

[18] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅲ. Symmetrical tilt boundaries[J]. Acta Metallurgica et Materialia, 1990, 38(5): 781-790.

[19] WOLF D. Structure-energy correlation for grain boundaries in f.c.c. metals - Ⅳ. Asymmetrical twist (general) boundaries[J].Acta Metallurgica et Materialia, 1990, 38(5): 791-798.

[20] RITTNER J D, SEIDMAN D N. 〈110〉 symmetric tilt grain-boundary structures in fcc metals with low stacking-fault energies[J]. Physical Review B, 1996, 54(10): 6999-7015.

[21] TSCHOPP M A, MCDOWELL D L. Structures and energies of Σ3 asymmetric tilt grain boundaries in copper and aluminium[J].Philosophical Magazine, 2007, 87(22): 3147-3173.

[22] OLMSTED D L, FOILES S M, HOLM E A. Survey of computed grain boundary properties in face-centered cubic metals: I. Grain boundary energy[J]. Acta Materialia, 2009,57(13): 3694-3703.

[23] LEE B J, CHOI S H. Computation of grain boundary energies[J].Modelling and Simulation in Materials Science and Engineering,2004, 12(4): 621-632.

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

[25] LEE B J, BASKES M I. Second nearest-neighbor modified embedded-atom-method potential[J]. Physical Review B, 2000,62(13): 8564-8567.

[26] LEE B J, SHIM J H, BASKES M I. Semiempirical atomic potentials for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, Al, and Pb based on first and second nearest-neighbor modified embedded atom method[J]. Physical Review B, 2003, 68(14): 144112.

[27] LENNARD-JONES J E. Cohesion[J]. The Proceedings of the Physical Society, 1931, 43(5): 461-482.

[28] BALLUFFI R W, MEHL R F. Grain boundary diffusion mechanisms in metals[J]. Metallurgical Transactions A, 1982,13(12): 2069-2095.