宓思恩 劉小明 魏悅廣，2)

* (中國科學(xué)院力學(xué)研究所非線性力學(xué)國家重點實驗室，北京 100190)

? (中國科學(xué)院大學(xué)工程科學(xué)學(xué)院，北京 100049)

** (北京大學(xué)工學(xué)院力學(xué)與工程科學(xué)系，北京 100871)

引言

先進(jìn)結(jié)構(gòu)材料基于微納米結(jié)構(gòu)設(shè)計，其宏觀力學(xué)性能依賴于微納米結(jié)構(gòu)，具有跨尺度力學(xué)特征[1-4].對其跨尺度力學(xué)性能的準(zhǔn)確模擬，需要采用從微觀的離散MD (molecular dynamics)模擬到宏觀的連續(xù)介質(zhì)表征.然而從離散體系到連續(xù)介質(zhì)體系的過渡是面臨的挑戰(zhàn)性科學(xué)和技術(shù)難題.原子模擬理論和方法的發(fā)展使得材料在納米尺度下的力學(xué)行為得到了充分的研究[5-9].如何將原子模擬方法與基于連續(xù)介質(zhì)力學(xué)理論的表征相關(guān)聯(lián)，把握先進(jìn)材料的跨尺度力學(xué)行為，有效地建立先進(jìn)材料的跨尺度力學(xué)理論，從而有效表征其在不同尺度下的力學(xué)性能成為一個重要的研究方向[10-13].其中，如何在離散原子體系中計算針對連續(xù)介質(zhì)定義的力學(xué)變量也是亟需解決的一個重要的科學(xué)技術(shù)問題.由于在原子模擬中(如MD)，模型由離散的原子構(gòu)成，體系中對于應(yīng)力和應(yīng)變等連續(xù)介質(zhì)力學(xué)理論框架下的力學(xué)變量并沒有統(tǒng)一的定義.許多研究都試圖在原子模擬中定義與傳統(tǒng)連續(xù)介質(zhì)力學(xué)等效的應(yīng)力和應(yīng)變，由于原子結(jié)構(gòu)及原子運(yùn)動的多樣性和離散特征，基于該過程給出的應(yīng)力和應(yīng)變的定義是難以直接與連續(xù)介質(zhì)理論的應(yīng)力-應(yīng)變定義形成有效過渡的，在過渡區(qū)不可避免地引起非平衡“鬼力”.為此，本文試圖建立一種從離散模擬(MD) 到連續(xù)介質(zhì)彈性FEA (finite element analysis)的過渡方法.在介紹這一過渡方法之前，下面首先簡要回顧一下若干MD 模擬研究中相關(guān)應(yīng)力和應(yīng)變概念的定義和原理.

在離散體系MD 方法中，人們也曾對應(yīng)力和應(yīng)變的概念進(jìn)行過定義和相關(guān)研究，例如，Irving 和Kirkwood[14]在流體力學(xué)方程的推導(dǎo)中，根據(jù)統(tǒng)計力學(xué)原理定義了位力應(yīng)力(virial stress)的概念，其中連續(xù)流體中的每一微小部分都被視為由N個微粒所組成的系統(tǒng)，應(yīng)力則被描述為該系統(tǒng)所受的面力，其中包含動能項與微粒間的作用力項.這種方法在計算過程中引入了適當(dāng)?shù)臅r間和空間量的平均，并且在計算微粒間的相互作用力時，為了將統(tǒng)計力學(xué)與流體力學(xué)運(yùn)動方程相聯(lián)系，只考慮通過微粒中心的作用力項，因此在計算微粒間作用力項時，只計算了微粒的正應(yīng)力.Hardy 等[15-16]注意到在Irving 和Kirkwood[14]的工作中，離散系統(tǒng)中守恒定律的準(zhǔn)確性依賴于整體的平均，研究過程中引入了無窮級數(shù)和狄拉克δ函數(shù)，需要通過近似求解確定.Hardy 使用了一種局域函數(shù)來代替狄拉克δ函數(shù)，以克服Irving 和Kirkwood 的工作在數(shù)學(xué)上的復(fù)雜性.這種局域函數(shù)代表對一個空間原子的運(yùn)動產(chǎn)生顯著影響的所有原子所占有體積的特征尺度，Hardy 證明了這種局域函數(shù)的取值范圍可以任意選擇，通過取極限，獲得穩(wěn)定的位力應(yīng)力值.Zimmerman 等[17]將Hardy 提出的計算方法與局部平均的位力應(yīng)力作對比，顯示Hardy 的方法在絕大多數(shù)情況下的精確性更優(yōu)于局部平均的位力應(yīng)力.在這兩種表達(dá)式中，選取不同特征體積時計算得到的應(yīng)力值不同，但其波動的幅度都可以通過適當(dāng)擴(kuò)大參與計算的特征體積，或引入適當(dāng)?shù)臅r間平均降到最低，以此得到一個穩(wěn)定值.Tsai[18]提出了另一種計算應(yīng)力的方法，該方法考慮了通過特定邊界的力和動量通量，來計算邊界所包含部分的應(yīng)力.通過邊界的合力的法向分量計算體內(nèi)的正應(yīng)力，邊界面內(nèi)的切向分量計算剪應(yīng)力.Lutsko[19]進(jìn)一步發(fā)展了Irving 和Kirkwood[14]的研究工作，通過局部的動量平衡方程和體積平均推導(dǎo)出應(yīng)力張量.Zhou[20]在其理論中指出，傳統(tǒng)的Cauchy 應(yīng)力應(yīng)該與原子的速度無關(guān)，只有分子間作用力部分才能真正對應(yīng)于傳統(tǒng)Cauchy 應(yīng)力.位力應(yīng)力的動能項依賴于原子的質(zhì)量和速度，會違背動量守恒，且不具備明確的物理意義.Liu 和Qu[21]則從柯西應(yīng)力的原始定義出發(fā)，論證了在經(jīng)過合適的時間平均后，基本的拉格朗日原子應(yīng)力與柯西應(yīng)力等效，進(jìn)一步提出了不含原子速度項的拉格朗日應(yīng)力和歐拉應(yīng)力，并證明經(jīng)典的位力應(yīng)力實際上對應(yīng)歐拉位力應(yīng)力.Falk 和Langer[22]在研究材料不可逆的剪切變形時，將局部的應(yīng)變定義為一定范圍內(nèi)，鄰近原子相對于中心原子的實際位移與發(fā)生均勻變形時的位移的均方差的最小值.Shimizu 等[23]在測量局部彈性應(yīng)變時引用了其中均勻應(yīng)變的計算.這是通過局部變形矩陣計算的拉格朗日應(yīng)變，本質(zhì)上反應(yīng)原子與其鄰近原子間的相對運(yùn)動.雖然應(yīng)變計算的形式簡單，但是Gullett 等[24]指出，其中剛體的轉(zhuǎn)動部分可能產(chǎn)生虛構(gòu)的剪應(yīng)變，導(dǎo)致計算錯誤.他們同時在文中介紹了變形梯度的計算，并與傳統(tǒng)的連續(xù)變形梯度進(jìn)行對比.同時，為了區(qū)分最鄰近原子與稍遠(yuǎn)的原子的不同貢獻(xiàn)，在變形梯度計算時引入了權(quán)函數(shù).為了建立材料的原子運(yùn)動離散體系與連續(xù)介質(zhì)體系的應(yīng)力和應(yīng)變概念的關(guān)聯(lián)，Mott 等[25]通過考慮原子結(jié)構(gòu)特征給出了一種定義變形梯度的幾何方法，以此獲得對應(yīng)變的定義.

由于原子的位置坐標(biāo)本質(zhì)上是離散的，計算所得的原子尺度的應(yīng)力與應(yīng)變值也是離散的，同時，以上所有的計算方法都需要選取一個合適范圍來確定計算所包含的鄰近原子(選取截斷半徑)，該范圍的確定在一定程度上是任意的，其明顯地對計算結(jié)果產(chǎn)生影響.在應(yīng)力計算時采用的局域函數(shù)，以及計算變形梯度時的權(quán)函數(shù)也是通過不精確地擬合原子相互作用與原子間距離的關(guān)系所得，都會使計算結(jié)果受到人為因素的影響.因此，通過傳統(tǒng)離散體系MD 模擬計算的應(yīng)力與應(yīng)變與連續(xù)介質(zhì)力學(xué)框架下定義應(yīng)力與應(yīng)變之間存在一定的差異.

本文針對晶體材料的微觀周期性結(jié)構(gòu)特征，提出一種協(xié)同考慮材料原子結(jié)構(gòu)、運(yùn)動、模擬以及連續(xù)介質(zhì)有限元分析的耦合計算方法，即MD-FEA 方法，將原子模擬結(jié)果直接轉(zhuǎn)變?yōu)橐阎邢拊?jié)點位移的應(yīng)力應(yīng)變場的后處理問題，避免了在后處理過程中仍需選取截斷半徑的問題.該方法的顯著優(yōu)點體現(xiàn)在:既直接將原子運(yùn)動離散體系轉(zhuǎn)化為連續(xù)介質(zhì)體系的力學(xué)場計算，又避免了傳統(tǒng)原子模擬方法在后處理過程中由于選擇截斷半徑而引起的誤差.對于許多具有周期性晶格排布的材料(如FCC，BCC和HCP 等)，將基本原子結(jié)構(gòu)選為連續(xù)介質(zhì)有限單元十分容易，同時以原子位移作為單元結(jié)點位移，計算有限單元內(nèi)的力學(xué)場也十分方便.因此，本文建立了一種結(jié)合了分子動力學(xué)與有限元分析的計算方法(MD-FEA)，針對典型FCC 晶體材料結(jié)構(gòu)的力學(xué)問題(如，Al/Ni 雙材料納米桿的拉伸與納米壓痕)，開展新方法與MD 方法的比對研究，以驗證新方法的有效性.

1 MD-FEA 方法

由于晶體材料中原子的排布具有周期性(如FCC，BCC 和HCP 等)，有利于將整個晶體材料的原子模型分割為有規(guī)律的單元體，進(jìn)而整個模型的力學(xué)行為便可以運(yùn)用有限元分析，通過各個單元的變形進(jìn)行描述.以圖1 中的面心立方晶體為例，分子動力學(xué)模型中的原子可以分割成若干有限元中的六面體單元，每個MD 模型中原子都可以視為單元中的節(jié)點.因此，每個FEA 六面體單元包含14 個節(jié)點，8 個位于六面體的頂點，6 個位于六面體的面心.每個頂點的原子同屬于8 個單元，面心的原子同屬于2 個單元.因此，通過MD 模擬得到每個原子的位移即為對應(yīng)單元的節(jié)點位移，六面體單元的應(yīng)變場可以通過屬于該單元的14 個節(jié)點的位移進(jìn)行插值計算.與傳統(tǒng)FEA 方法不同的是，這種方法不需要計算模型的剛度矩陣，每個節(jié)點的位移都是通過MD 模擬獲得的.

圖1 14 節(jié)點單元示意圖Fig.1 Schematic diagram of a 14-node element

1.1 MD 模擬的相關(guān)基本關(guān)系

求解原子的牛頓運(yùn)動基本方程為

其中ri，mi，Pi分別為原子i的位置、質(zhì)量和其受到的作用力，U為n個原子所組成的系統(tǒng)的總勢能，由系統(tǒng)中每個原子的作用勢求和得到.對晶體結(jié)構(gòu)金屬材料，常用的兩種原子作用勢為:Lennard-Jones(L-J)作用勢和鑲嵌原子勢(EAM 勢).

對于Lennard-Jones 作用勢，單個原子的作用勢Ei為

其中rij為原子i與j之間的距離，ε，σ和r0分別為兩個原子的結(jié)合能，平衡距離和截斷距離.

對于EAM 原子作用勢[17-18]

其中Ei是原子i的總勢能，α和β分別為原子i和j的類別，F(xiàn)代表原子的鑲嵌能，是該原子電子密度ρ的函數(shù)，rij是原子i和j之間距離，φ為對勢.

1.2 FEA 相關(guān)基本關(guān)系

對于一個由面心立方晶體構(gòu)成的14 個原子的結(jié)構(gòu)，轉(zhuǎn)換成有限元的單元結(jié)構(gòu)即為由14 節(jié)點構(gòu)成的空間六面體單元，其結(jié)構(gòu)如圖1 所示.常采用的等參元的位移場可以描述為

其中，- 1 ≤ξ，η，ζ ≤1 為等參元的坐標(biāo)分量.

由此可以確定單元的形函數(shù)N為

因此單元內(nèi)一點處的應(yīng)變分量為

其中δe是由原子位移確定的節(jié)點的位移矢量，應(yīng)變-位移矩陣B可以分割為

雅可比矩陣J描述了實際單元與等參元的坐標(biāo)轉(zhuǎn)換關(guān)系.

如果考慮材料為線性彈性情況，單元的應(yīng)力場與應(yīng)變能密度分別為

可應(yīng)用式(4)、式(6)和式(9)通過節(jié)點位移插值直接計算單元內(nèi)高斯點處的位移、應(yīng)變和應(yīng)力.若考慮3 × 3 × 3 的高斯積分點方案，則單元內(nèi)有27 個高斯積分點(通常簡稱高斯點).

2 MD-FEA 方法的應(yīng)用

下面通過兩個典型物理問題的研究，通過比較本文MD-FEA 方法的結(jié)果與MD 模擬結(jié)果，以檢驗本文方法的有效性.

2.1 軟硬兩種金屬(Al/Ni)等截面納米桿串聯(lián)結(jié)構(gòu)的拉伸問題

研究的問題如圖2 所示.

首先采用分子動力學(xué)軟件LAMMPS[26]建立由Al 和Ni 構(gòu)成的雙材料納米桿拉伸模型，原子間的相互作用采用EAM 勢函數(shù)[26-27](見式(3)).整個系統(tǒng)由233 568 個原子(93 296 個Al 原子和140 272 個Ni 原子) 構(gòu)成，模型在x方向上長為44 nm (Al 和Ni 段的長度各為22 nm)截面為邊長8 nm 的正方形.模型的初始構(gòu)型如圖2 所示.Al 和Ni 晶體的[1 0 0]晶向與x軸平行，模型在x方向的兩端設(shè)為剛體，y，z方向的4 個面為自由表面.整個系統(tǒng)在NVT 系綜下保持溫度為0.01 K.

圖2 Al/Ni 雙材料納米桿拉伸分子動力學(xué)模型.界面及端部標(biāo)注為黃色Fig.2 Molecular Dynamics model of the Al/Ni composite nano cylinder system.The interface layers and the boundaries are marked yellow

需要指出的是，許多研究[28-30]表明，實際情況下的Al/Ni 界面十分復(fù)雜，難以進(jìn)行準(zhǔn)確的原子模型建模.為了方便研究，Gurtin 與Murdoch[31]建立了簡化的材料表/界面模型:假定材料的表/界面厚度為0 且具有一定的力學(xué)性質(zhì)，并在連續(xù)介質(zhì)力學(xué)框架下根據(jù)這個模型建立了材料的表/界面彈性理論.根據(jù)這個假設(shè)，同時也為了簡化計算，在本模型中，Al 和Ni 在界面處各取一層原子層，并在x方向上進(jìn)行固定，使得界面原子在運(yùn)動中整體保持在同一平面內(nèi).界面處原子的相互作用采用L-J 勢式(2)，以改變不同工況下的界面性能.另外，在納米桿拉伸模型中，通常采用桿整體的應(yīng)變變化作為外部載荷，由于計算規(guī)模的限制，應(yīng)變率通常在106～ 108s-1之間.而在本模型中，為了保證界面在變形過程中保持平面，限制了界面原子在x方向的位移，Al 和NI 之間的變形無法通過界面進(jìn)行有效地傳遞.因此，本模型采用力加載的方式，在模型x方向兩端施加p=160.2 N/s 的外部載荷取代模型整體的應(yīng)變加載，系統(tǒng)在整個加載過程中保持平衡.系統(tǒng)的時間步長設(shè)為0.001 ps.由圖3 的應(yīng)變-時間關(guān)系可知，Al，Ni 兩相的應(yīng)變率分別為4×107s-1和1.7×107s-1，并且在整個模擬過程中基本保持一致.

圖3 拉伸載荷下Al 和Ni 的應(yīng)變-時間曲線Fig.3 The strain-time relation of Al and Ni under uniform force loading

模型中Al 和Ni 的應(yīng)力應(yīng)變曲線如圖4 所示.在達(dá)到極限應(yīng)力σxx=7.56 GPa 之前，Al 和Ni 整體上都符合彈性變形.

達(dá)到極限應(yīng)力后，Al 在界面附近發(fā)生破壞，如圖5 所示.

圖5 超過極限應(yīng)力σxx=7.56 GPa 后Al-Ni 系統(tǒng)的變形Fig.5 Deformed Al-Ni system after the critical tensile stress σxx=7.56 GPa

根據(jù)模型的應(yīng)力-應(yīng)變關(guān)系，可知Al 和Ni 的楊氏模量分別為72.2 GPa 與134.6 GPa，這與其他研究[27，32]中，根據(jù)EAM 勢的下的Al 和Ni 體積模量分別為79 GPa 和181 GPa，泊松比分別為0.35 與0.385計算所得的彈性模量分別為71 GPa 和125.3 GPa 相近.整個模型在彈性階段的楊氏模量為92.1 GPa.

通過MD-FEA 方法計算模型每個單元高斯點處的拉伸應(yīng)變εxx代表整個模型的連續(xù)應(yīng)變場.在極限應(yīng)力σxx=7.56 GPa 下，模型在z=4 nm 處的x-y截面的拉伸應(yīng)變εxx如圖6 所示，界面位置為x=0 nm，邊界位置分別為x=-22 nm 和x=22 nm.材料Al內(nèi)部的拉伸應(yīng)變?yōu)?.12，材料Ni 內(nèi)部拉伸應(yīng)變?yōu)?.05，與圖4 的結(jié)論一致，而在材料的界面和端部，由于邊界條件的影響，拉伸應(yīng)變與材料內(nèi)部有所區(qū)別.

圖6 拉伸載荷σxx=7.56 GPa 下，Al 和Ni 在中截面處的應(yīng)變εxx，界面位置為x=0 nmFig.6 εxx at the middle section of Al and Ni under σxx=7.56 GPa.The interface at x=0 nm

圖6 為材料在z=4 nm 處的中截面的拉伸應(yīng)變εxx.在截面上，材料Al 和Ni 從模型底部(y=0 nm)到模型中部(y=4 nm)的應(yīng)變分布如圖7 所示.應(yīng)變在剛體邊界附近趨向于零，并朝著材料內(nèi)部急劇增加，并經(jīng)過少許的降低后，在材料Al 和Ni 的體內(nèi)分別穩(wěn)定在0.122 和0.053.在材料的界面附近，材料Al 處的應(yīng)變隨著與界面距離的減少整體上有少許的降低，而材料Ni 處的應(yīng)變有少許增加，且與體內(nèi)應(yīng)變相比具有較大的分散性.界面附近應(yīng)變的變化集中在距離界面4 nm 內(nèi)的區(qū)域中.

圖7 材料中截面內(nèi)不同位置的拉伸應(yīng)變.(a) Al，(b) NiFig.7 Tensile strain in (a) Al and (b) Ni at the different location of the middle section

在相同的拉伸應(yīng)力下，界面附近的拉伸應(yīng)變εxx在材料Al 處的降低和材料Ni 處的升高主要是由于材料性能的不同引起的.根據(jù)EAM 勢函數(shù)的結(jié)論與應(yīng)力應(yīng)變關(guān)系曲線，可知Al 和Ni 在[1 0 0]晶向上的楊氏模量分別為72.2 GPa 和134.6 GPa，泊松比分別為0.35 和0.385.因此，Al 和Ni 在y方向平行于界面的正應(yīng)變分別為

設(shè)Al 和Ni 在界面處平行于界面的正應(yīng)變的差值為失配應(yīng)變

根據(jù)界面處的無滑移邊界條件，需要在界面處引入額外的平行于界面的正應(yīng)力，以保證在x=0 的界面處，滿足

假設(shè)界面處完全是彈性變形，且界面兩邊材料中額外引入的應(yīng)力使得界面整體保持平衡.因此

在z方向上同理.由此可得Al 和Ni 在x方向上的拉伸應(yīng)變

由式(18)計算可得，模型中Al 在x 方向上的拉伸應(yīng)變在界面處減小6.5%，Ni 的拉伸應(yīng)變在界面處增大7.2%，這與圖7 的結(jié)論一致.

在分子動力學(xué)模擬中可以得到模型的原子應(yīng)變.根據(jù)文獻(xiàn)[16-17]的研究，原子的應(yīng)變可以通過變形梯度計算得到

變形梯度Fi通過求下式的最小值得到

其中dji和分別為即時構(gòu)型和初始構(gòu)型下原子j到原子i的矢量.計算原子應(yīng)變的截斷半徑設(shè)為0.7 nm，且不設(shè)置加權(quán)函數(shù).拉伸載荷下，在模型的整體拉伸應(yīng)變從0.02 到0.08 的變形過程中，通過兩種方法計算得到的模型中心(y=z=4 nm)處沿x方向的拉伸應(yīng)變εxx如圖8 所示.使用MD-FEA 方法計算得到的材料內(nèi)部的拉伸應(yīng)變(MF)與通過變形梯度計算的原子應(yīng)變(AS)在模型整體應(yīng)變達(dá)到0.12之前基本相同.在之后的變形過程中，原子應(yīng)變計算的拉伸應(yīng)變大于使用MD-FEA 方法得到的拉伸應(yīng)變.這主要是由于MD-FEA 計算了在彈性變形中常用的小應(yīng)變張量，而原子應(yīng)變基于變形梯度計算的是格林-拉格朗日應(yīng)變張量.這兩種應(yīng)變在變形較小時的差異可忽略不計，但對于計算應(yīng)變達(dá)到0.12 的材料Al，其結(jié)果有明顯的差異.

圖8 模型整體應(yīng)變0.02～0.08 的變形過程中MD-FEA 方法計算的應(yīng)變與原子應(yīng)變的對比Fig.8 Strain computed by MD-FEA method and deformation gradient with the system strain from 0.02 to 0.08

同時，本文計算了7.56 GPa 拉伸載荷下，模型中心在截斷半徑0.4～0.7 nm 下原子應(yīng)變計算的拉伸應(yīng)變值，如圖9 所示.不同截斷半徑下，計算的原子拉伸應(yīng)變在模型端部和體內(nèi)基本保持一致，但在距離界面2 原子層內(nèi)產(chǎn)生明顯的差異.考慮到圖7 所示的兩種材料的拉伸應(yīng)變在界面附近的橫截面內(nèi)有較大的分散性，可以認(rèn)為這種復(fù)雜的應(yīng)變分布是導(dǎo)致不同截斷半徑下計算得到的界面處的拉伸應(yīng)變有較大差異的主要原因.

圖9 不同截斷半徑下計算的原子應(yīng)變Fig.9 Atom strain computed in different cutoff radius

在模型兩端均勻的拉伸載荷下，材料Al 和Ni 在界面附近的拉伸應(yīng)變在同一橫截面內(nèi)表現(xiàn)出一定的周期性變化，如圖10 所示.除了由于應(yīng)力集中出現(xiàn)在4 個角落的最大應(yīng)變外，在材料Al 界面附近第2 層原子處，中心附近的拉伸應(yīng)變大于四周邊界附近的拉伸應(yīng)變，邊界附近的部分區(qū)域有著較大的拉伸應(yīng)變，并在靠近材料內(nèi)部時急劇減小.第2 層Ni 原子處，四角附近的拉伸應(yīng)變大于材料中心區(qū)域，并且有規(guī)律地，在某些位置的拉伸應(yīng)變明顯減小.在兩種材料界面附近的第4 層原子處，大致能觀察到與各自第2 層原子處相反的應(yīng)變分布，且應(yīng)變在橫截面內(nèi)的變化幅度明顯減小.這表明界面附近的拉伸應(yīng)變場在平行與垂直界面方向上都有波動，且隨著與界面距離的增加，波動幅度逐漸降低.在距離界面第6 層原子處，應(yīng)變波動基本消失，除邊界附近應(yīng)變與內(nèi)部稍有變化外，材料內(nèi)部的應(yīng)變與遠(yuǎn)離界面的應(yīng)變基本相同.

圖10 (a) Al 和(b) Ni 在距界面第2，4，6 層截面原子處的拉伸應(yīng)變Fig.10 Tensile strain of the 2nd，4th and 6th layers of the cross section near the interface for (a) Al and (b) Ni

圖10 (a) Al 和(b) Ni 在距界面第2，4，6 層截面原子處的拉伸應(yīng)變(續(xù))Fig.10 Tensile strain of the 2nd，4th and 6th layers of the cross section near the interface，for (a) Al and for (b) Ni (continued)

界面原子在弛豫過程中發(fā)生重構(gòu)時，在對其應(yīng)力應(yīng)變的研究中也觀察到了類似的波動現(xiàn)象[33].但是，本文選擇了弛豫后的構(gòu)型為參考構(gòu)形，消除了弛豫過程中殘余應(yīng)變的影響，這表明界面處的應(yīng)變波動也會由外部施加的拉伸載荷引起.材料Al 在橫截面上有20 × 20 個晶格，材料Ni 有23 × 23 個晶格.界面處的晶格失配與應(yīng)變波動的周期相吻合.因此可以認(rèn)為界面的原子失配不僅僅引起界面處的殘余應(yīng)變，同時也改變了界面處的材料性能，使得界面截面的不同位置在相同的拉伸載荷下產(chǎn)生不同的拉伸應(yīng)變.

拉伸應(yīng)變在界面附近處的劇烈變化是原子應(yīng)變在不同的截斷半徑下得到的結(jié)果產(chǎn)生巨大差異的主要原因之一.圖11 展示了Al 在距界面第2 層原子處的原子應(yīng)變在不同截斷半徑下的分布，作為與MD-FEA 方法得到的結(jié)果(圖10(a1))的對比.截斷半徑rc設(shè)置在0.35～0.38 nm 之間.結(jié)果顯示截斷半徑rc=0.36 nm 和rc=0.37 nm 下的原子應(yīng)變分布與MD-FEA 方法的結(jié)果有相似之處，而rc=0.35 nm和rc=0.38 nm 的結(jié)果卻截然不同.同時，相比于MD-FEA 方法，原子應(yīng)變的分布規(guī)律顯得較為粗略，忽略了許多細(xì)節(jié).對于FCC 晶體，每個晶體單元擁有4 個原子，因此每個單元根據(jù)原子位移得到4 個原子應(yīng)變值，而MD-FEA 方法準(zhǔn)確計算了每個單元在27 個高斯點處的應(yīng)變，因此能夠提供更多的應(yīng)變細(xì)節(jié)，得到的應(yīng)變場也更加精確.

圖11 截斷半徑為0.35～0.38 nm 時，距界面第2 層Al 原子截面處的拉伸應(yīng)變分布Fig.11 Atom strain at 2th layer in phase Al with the cutoff radius rc from 0.35 to 0.38

圖11 截斷半徑為0.35～0.38 nm 時，距界面第2 層Al 原子截面處的拉伸應(yīng)變分布(續(xù))Fig.11 Atom strain at 2th layer in phase Al with the cutoff radius rc from 0.35 to 0.38 (continued)

當(dāng)計算一種材料的原子應(yīng)變時，截斷半徑內(nèi)可能同時包含其他材料的原子，這也會明顯影響計算的結(jié)果.如圖12 所示，當(dāng)計算界面附近Al 原子的應(yīng)變時，截斷半徑內(nèi)參與計算的不僅有附近的Al 原子，同時也包含了部分Ni 原子.隨著截斷半徑的增大，使得界面附近的Al 原子計算結(jié)果誤差也隨之增大.在相同的拉伸載荷下，材料Al 的應(yīng)變明顯大于Ni 的應(yīng)變，這使得界面附近Al 原子的應(yīng)變隨著截斷半徑的增大而減小，Ni 的原子應(yīng)變隨截斷半徑的增大而增大.

圖12 rc=0.8 nm 下，界面原子附近參與原子應(yīng)變計算的Al 原子(黃色和橙色)和Ni 原子(綠色和藍(lán)色)數(shù)量示意圖Fig.12 A schematic diagram of the number of atom Al (yellow and orange) and Ni (blue and green) that participated in the calculation of the atom strain of the selected atom near the interface with rc=0.8 nm

以距界面第2，4，6 層原子層截面處中線上(z=4 nm)的拉伸應(yīng)變?yōu)槔?.56 GPa 的拉伸載荷下，由MD-FEA 方法計算的結(jié)果與不同截斷半徑下原子應(yīng)變的結(jié)果如圖13 所示.

圖13 不同截斷半徑下，距離界面第2，4，6 層原子層截面中線上的(a1)-(a3) Al 原子和(b1)-(b3) Ni 原子的原子應(yīng)變分布以及同一位置下MD-FEA 方法計算的應(yīng)變Fig.13 MD-FEA strain and atom strain of the 2nd，4th and 6th layers of (a1)-(a3) Al and (b1)-(b3) Ni in the middle line of interface section computed by different cutoff radius

圖13 不同截斷半徑下，距離界面第2，4，6 層原子層截面中線上的(a1)-(a3) Al 原子和(b1)-(b3) Ni 原子的原子應(yīng)變分布以及同一位置下MD-FEA 方法計算的應(yīng)變(續(xù))Fig.13 MD-FEA strain and atom strain of the 2nd，4th and 6th layers of (a1)-(a3) Al and (b1)-(b3) Ni in the middle line of interface section computed by different cutoff radius (continued)

在截斷半徑為0.4 nm 時，距離界面第2 層原子層的截面內(nèi)計算的Al 和Ni 原子應(yīng)變表現(xiàn)出明顯的波動，并且隨著截斷半徑的增大，這種波動幅度明顯減小.同時，隨著截斷半徑的增大，Al 的平均原子應(yīng)變呈現(xiàn)逐漸減小的趨勢，而Ni 的平均原子應(yīng)變逐漸增大.在此截面上，通過MD-FEA 方法計算的應(yīng)變大于原子應(yīng)變的計算值，且在截面中線上，應(yīng)變的變化趨勢也與某一特定截斷半徑下的原子應(yīng)變類似(rc=0.6 nm 的Al 與rc=0.4 nm 的Ni).而在距離界面第4 層原子層處，不同截斷半徑下計算的原子應(yīng)變較為一致，隨著截斷半徑的增大，Al 和Ni 原子的原子應(yīng)變計算值均稍有增大，在距離界面第6 層原子層處，不同截斷半徑下計算的原子應(yīng)變基本相同.在第4 層與第6 層原子層處，原子應(yīng)變的計算值均大于MD-FEA 方法的計算值，而這兩種方法計算的應(yīng)變沿著截面中線的變化趨勢較為一致，且原子應(yīng)變的變化幅度明顯小于MD-FEA 方法計算的應(yīng)變.隨著與界面距離的增加，不同方法計算的Al 原子的應(yīng)變稍有增大，Ni 原子的應(yīng)變稍有減小，這與圖7 和圖8的結(jié)論相吻合.

在MD-FEA 方法中，應(yīng)變是以參考構(gòu)型中的一個單元為基本單元計算的，這種方法消除了人為選取的截斷半徑的影響.單元的體積由材料的晶格常數(shù)確定，每個單元的節(jié)點(原子)數(shù)量相同，且單元的體積足夠小，從而能得局部的精確結(jié)果.在MDFEA 方法中，每個單元都只包含了同一種原子，界面附近的Al 原子只能通過相互作用改變界面另一邊Ni 原子的構(gòu)型來間接反應(yīng)對Ni 原子的影響，而不能Al 原子構(gòu)型的變化直接參與Ni 原子處應(yīng)變的計算，反之亦然.用這種方法計算界面附近的應(yīng)變場更加合理.

圖14 展示了模型在1.26 GPa 的拉伸載荷下中截面的應(yīng)力分布.位力應(yīng)力由下式計算得到

其中mi，vi分別為原子i的質(zhì)量與速度，rij，fij分別為原子i與原子j之間的距離與作用勢，V為原子體積.

根據(jù)圖14(a)和圖14(b)可知，MD-FEA 方法計算的應(yīng)力明顯大于原子的位力應(yīng)力.在加載的拉伸應(yīng)力為1.26 GPa 時，MD-FEA 方法計算的Al 和Ni 在模型體內(nèi)的應(yīng)力為1.32 GPa 見圖14(a)，而位力應(yīng)力計算的應(yīng)力值明顯較小，且在兩種材料中并不一致見圖14(b).這主要是由于在MD-FEA 方法中，選擇弛豫后的構(gòu)型為參考構(gòu)形，而在參考構(gòu)形中計算的位力應(yīng)力有一定殘余的壓應(yīng)力存在見圖14(c).

圖14 拉伸應(yīng)力為1.26 GPa 時，不同方法計算的模型中截面應(yīng)力分布Fig.14 Distribution of the stress on the cross section computed by different methods under the applied tensile stress 1.26 GPa

圖14(d)展示了MD-FEA 方法與消除了殘余應(yīng)力的位力應(yīng)力在模型軸心處的對比，可見兩種方法計算的結(jié)果較為一致.在使用MD-FEA 方法計算應(yīng)力時，假設(shè)模型的楊氏模量和泊松比在變形過程中保持不變，而通過模型的應(yīng)力-應(yīng)變關(guān)系可知，材料的楊氏模量在變形較大時有輕微的變化.因此，鑒于選取的材料參數(shù)為小變形下的參數(shù)，可以認(rèn)為用MD-FEA 方法計算的應(yīng)力在模型的變形較小時較為精確.同時，MD-FEA 方法計算的應(yīng)力在模型的表面與內(nèi)部較為一致，而位力應(yīng)力的結(jié)果顯示了由于鄰近的原子構(gòu)型的差異，在同一拉伸載荷下，模型表面原子和內(nèi)部原子計算的應(yīng)力有顯著的差異.因此，MD-FEA 方法更能準(zhǔn)確地計算模型的應(yīng)力分布.

2.2 金剛石壓頭/Al 基體的納米壓痕實驗分析

使用分子動力學(xué)軟件LAMMPS 建立納米壓痕模型，基體長和寬分別為40 個晶格，高為15 個晶格的Al，壓頭為由半徑為10 nm 的金剛石球形壓頭.初始構(gòu)形中壓頭尖端與基體表面距離為1.156 nm，并在每時間步中對壓頭施加0.05 晶格的位移(約0.02 nm) 后，對整個模型弛豫.模型中的Al 基體采用EAM 勢[27]，金剛石壓頭采用tersoff 勢[34]，Al 和金剛石之間的相互作用使用L-J 勢描述，其中，結(jié)合能ε=31.5×10-3eV，平衡距離為2.976 ?，截斷半徑為10 ?[35].基體的底部與四周為固定邊界，整個系統(tǒng)采用NVE 系綜.溫度設(shè)定為0.01 K.模型的初始構(gòu)型見圖15(a).圖15(b)展示了載荷與壓入深度D之間的關(guān)系.

圖15 (a)納米壓痕的分子動力學(xué)模型和(b)載荷-壓入深度曲線Fig.15 (a) Molecular dynamics model of the indentation and (b) the force-depth relation during the simulation

從D=-0.55 nm 開始，基體和壓頭之間表現(xiàn)出相互吸引，壓頭受到基體向下的拉力，并在D=-0.4 nm 時達(dá)到最大值.直到D=-0.1 nm 時，拉力轉(zhuǎn)變?yōu)閴毫?，且在D=0.15 nm 之前不斷增大，在D=0.15 nm時壓力出現(xiàn)第一次波動，基體開始出現(xiàn)位錯，此時壓力為150 nN.

基體的應(yīng)變分量εxx和γxz在臨界壓入深度下的分布如圖16 所示.在這兩種計算方法中，εxx的影響區(qū)域較為接近，而MD-FEA 方法計算的反對稱的剪應(yīng)變γxz的大小和影響范圍都明顯大于原子應(yīng)變的結(jié)果.

圖16 MD-FEA 方法計算的應(yīng)變分量(a1)-(a2) εxx 與(b1)-(b2) γxz 和對應(yīng)原子應(yīng)變的對比Fig.16 Comparison of (a1)-(a2) εxx and (b1)-(b2) γxz computed by the atom strain and MD-FEA method

圖17 展示了用MD-FEA 方法計算的基體中截面在壓痕過程中，位錯出現(xiàn)前的各個階段(圖15(a)-圖15(e))的應(yīng)變場，并與原子應(yīng)變進(jìn)行了對比.計算原子應(yīng)變的截斷半徑取為rc=0.45 nm.

圖17 壓痕過程中(a1)-(e1)原子應(yīng)變與(a2)-(e2) MD-FEA 應(yīng)變分量εzz 的對比Fig.17 (a1)-(e1) Atom strain and (a2)-(e2) the strain field computed by MD-FEA method of εzz during the indentation

圖17 壓痕過程中(a1)-(e1)原子應(yīng)變與(a2)-(e2) MD-FEA 應(yīng)變分量εzz 的對比(續(xù))Fig.17 (a1)-(e1) Atom strain and (a2)-(e2) the strain field computed by MD-FEA method of εzz during the indentation (continued)

由MD-FEA 方法計算得到的基體內(nèi)部的應(yīng)變場εzz和原子應(yīng)變的結(jié)果基本相同.基體的應(yīng)變主要集中在壓頭下方，并隨著壓入深度的增加不斷擴(kuò)大.壓坑周圍出現(xiàn)一圈拉伸應(yīng)變區(qū)域，并隨著壓入深度的增加不斷往外擴(kuò)張，這與其他研究的結(jié)論類似[36].在壓入深度為0.1 nm 位錯出現(xiàn)之前，壓縮應(yīng)變主要分布在壓痕表面下方約3.5 nm 的區(qū)域內(nèi)，壓縮應(yīng)變的最大值并不在壓痕與基體的接觸表面，而是出現(xiàn)在壓痕表面下方約1 nm 的區(qū)域中.在出現(xiàn)位錯時，位錯位置計算得到的應(yīng)變值顯著增大，因此能明確判別出位錯產(chǎn)生的位置.

由兩種方法計算得到的基體中心的應(yīng)變分量εzz和εxx的在壓痕過程中的對比如圖18 所示.在應(yīng)變較小時，MD-FEA 方法計算得到應(yīng)變值和原子應(yīng)變完全一致，在應(yīng)變較大時，MD-FEA 方法計算的壓應(yīng)變εzz較原子應(yīng)變稍大.應(yīng)變分量εzz和εxx的最大值都出現(xiàn)在壓頭下方的基體內(nèi)部而不是壓痕表面，且這兩者出現(xiàn)的位置并不一致.εzz的最大值出現(xiàn)在壓痕表面約1 nm 的位置，而εxx的最大值出現(xiàn)在壓痕表面約1.5 nm

圖18 由MD-FEA 應(yīng)變和原子應(yīng)變在基體中心處的應(yīng)變分量εzz 和εxx 的對比Fig.18 Comparison of εzz and εxx at the center of the substrate between MD-FEA method (MF) and atom strain (AS) during the indentation

截斷半徑在計算壓痕過程中基體原子應(yīng)變的影響較小.在截斷半徑rc=0.4 nm 和rc=0.8 nm 下計算的基體中截面的應(yīng)變分布如圖19 所示，應(yīng)變分量εzz和εxx在不同截斷半徑下的結(jié)果在基體內(nèi)部基本一致，僅在壓痕表面的原子處有所差異.rc=0.8 nm 下計算的應(yīng)變分量εzz在壓痕表面的值較rc=0.4 nm 時稍大，而分量εxx在壓痕表面的值則較rc=0.4 nm 時稍小.這種差異在分析壓痕過程中基體的整體變形時可忽略不計，但在關(guān)注壓痕表面的變化時必須加以考慮.

圖19 截斷半徑rc=0.4 nm 和rc=0.8 nm 下中截面上的原子應(yīng)變分量(a1)-(a2) εzz 與(b1)-(b2) εxx 分布圖Fig.19 Atom strain of (a1)-(a2) εzz and (b1)-(b2) εxx with the cutoff radius rc=0.4 nm and rc=0.8 nm

圖19 截斷半徑rc=0.4 nm 和rc=0.8 nm 下中截面上的原子應(yīng)變分量(a1)-(a2) εzz 與(b1)-(b2) εxx 分布圖(續(xù))Fig.19 Atom strain of (a1)-(a2) εzz and (b1)-(b2) εxx with the cutoff radius rc=0.4 nm and rc=0.8 nm (continued)

圖20 展示了由兩種方法計算的基體中心不同深度的應(yīng)力分量σzz和σxx在壓痕各階段的對比.在位力應(yīng)力的計算中，基體的楊氏模量設(shè)為72.2 GPa，泊松比為0.35.σzz在基體不同深度均有體現(xiàn)，σxx則主要影響了距離壓痕表面2 nm 以內(nèi)的區(qū)域.除了基體表面位力應(yīng)力計算不準(zhǔn)確的部分，兩種方法計算的基體內(nèi)部的應(yīng)力在變形較小時基本一致，在變形接近產(chǎn)生位錯的臨界狀態(tài)時，計算的位力應(yīng)力值略大于MD-FEA 方法計算的結(jié)果.

圖20 壓痕過程中由MD-FEA 方法(MF)計算的應(yīng)力分量與位力應(yīng)力(AS)分量σzz 和σxx 的對比Fig.20 Comparison of σzz and σxx between MD-FEA method (MF) and atom virial stress (AS) during indentation

圖21 展示了由MD-FEA 方法計算的圖15(b)中從a～e 狀態(tài)下的正應(yīng)力σzz在中截面處的分布，并與對應(yīng)位置處原子的位力應(yīng)力進(jìn)行了對比.由于原子構(gòu)型的差異，基體表/界面處原子的位力應(yīng)力結(jié)果與體內(nèi)原子有明顯差異，因此明顯影響了壓痕表面的應(yīng)力結(jié)果.忽略起始階段由于壓頭與基體相互吸引產(chǎn)生的拉應(yīng)力，基體最大壓應(yīng)力出現(xiàn)在壓痕表面下方，且比同一狀態(tài)下最大壓應(yīng)變εzz出現(xiàn)的位置更接近壓痕表面.在產(chǎn)生位錯的臨界狀態(tài)下，壓應(yīng)力σzz的最大值達(dá)到10 GPa.

圖21 壓痕過程中基體截面處(a1)-(e1)原子位力應(yīng)力與(a2)-(e2) MD-FEA 應(yīng)力場的對比Fig.21 (a1)-(e1) Virial stress and (a2)-(e2) the stress field computed by MD-FEA method on the middle section of the substrate during indentation

圖21 壓痕過程中基體截面處(a1)-(e1)原子位力應(yīng)力與(a2)-(e2) MD-FEA 應(yīng)力場的對比(續(xù))Fig.21 (a1)-(e1) Virial stress and (a2)-(e2) the stress field computed by MD-FEA method on the middle section of the substrate during indentation (continued)

3 結(jié)論及討論

本文提出了一種通過結(jié)合分子動力學(xué)(MD)和有限元分析(FEA)來計算原子模型中的應(yīng)力和應(yīng)變的新方法(MD-FEA).這種方法通過離散的原子坐標(biāo)和位移數(shù)據(jù)，計算連續(xù)介質(zhì)力學(xué)框架下的計算連續(xù)的應(yīng)力和應(yīng)變場.本文通過建立Al/Ni 雙材料納米桿拉伸模型和金剛石壓頭/鋁基體的納米壓痕模型，將MD-FEA 方法計算得到的應(yīng)力/應(yīng)變場和在傳統(tǒng)的原子模型中廣泛使用的位力應(yīng)力，以及通過離散的變形梯度計算得到的原子應(yīng)變進(jìn)行對比.結(jié)果顯示，在模型體內(nèi)，MD-FEA 方法計算的應(yīng)變分量和原子應(yīng)變相一致，但在原子位移急劇變化的模型表/界面上有明顯的差異.這主要是由于原子應(yīng)變在計算過程中平均了鄰近原子的應(yīng)變，且包含了截斷半徑內(nèi)不同類型，材料參數(shù)有差異的原子.除了模型表/界面處由于原子構(gòu)型的差異，位力應(yīng)力計算不準(zhǔn)確外，MD-FEA 方法計算的應(yīng)力與位力應(yīng)力在模型體內(nèi)變形較小時完全一致.由于MD-FEA 方法通過楊氏模量和泊松比將原子模型中的應(yīng)力和應(yīng)變相聯(lián)系，在大變形下，由于材料參數(shù)可能發(fā)生變化，計算的應(yīng)力不精確.因此通過這種方法計算應(yīng)力時應(yīng)限定在小變形下.

結(jié)果顯示:在應(yīng)變較為均勻的區(qū)域，MD-FEA 方法計算的應(yīng)變場與基于MD 變形梯度計算的原子應(yīng)變相一致，而在界面附近應(yīng)變變化劇烈的區(qū)域，MDFEA 方法能夠提供精確的應(yīng)變場，而原子應(yīng)變只能得到較大范圍內(nèi)的平均值.用MD-FEA 方法計算的應(yīng)力場與位力應(yīng)力在模型體內(nèi)相一致，但在模型界面與表面區(qū)域，由于鄰近原子的構(gòu)型發(fā)生變化，位力應(yīng)力的計算有很大的誤差，而MD-FEA 方法仍能得到準(zhǔn)確的應(yīng)力場.MD-FEA 方法避免了用應(yīng)變梯度算法計算原子應(yīng)變時截斷半徑的選取，以及計算位力應(yīng)力時原子體積的確定，也不需要引入加權(quán)函數(shù)，消除了人為因素對應(yīng)力應(yīng)變結(jié)果的影響.同時，MDFEA 方法得到的是原子模型中連續(xù)的應(yīng)力應(yīng)變場，將離散的原子模型和傳統(tǒng)連續(xù)介質(zhì)力學(xué)理論框架之間建立了聯(lián)系.

需要指出的是:按本文提出的計算模型，原子的位置完全是按照MD 計算給出的，所以本文方法與MD 計算結(jié)果對于原子的位移來說完全相同，而應(yīng)變和應(yīng)力則不同;本文將一個晶格處理為一個有限元單元，事實上這里的有限單元可以擴(kuò)大至一個單元代表多個晶格的情況，只要初始劃分的單元的節(jié)點能夠與原子的初始位置對應(yīng)即可，但此時由兩種模型給出的應(yīng)力和應(yīng)變結(jié)果的差異將會更大一些，這個將在后續(xù)研究中作進(jìn)一步研究和討論.