晶粒尺寸對多晶銅接觸力學特性的影響

林起崟 張瑜寒 洪 軍 王 晨 張寧靜

1.西安交通大學現代設計及轉子軸承系統教育部重點實驗室，西安，7100492.西安交通大學設計科學與基礎部件研究所，西安，710049

0 引言

機械裝備中存在大量的結合面(又稱裝配連接界面)，結合面的存在破壞了機械結構力和轉矩傳遞的連續性，因此它對整機性能具有舉足輕重的影響[1]。結合面從微觀上看是由眾多高度不同的微凸體構成的，并且結合面的表面形貌具有多尺度特征，其表征尺度可達納米尺度[2]。接觸剛度作為描述結合面接觸力學性能的重要參數之一，其大小會直接影響整機的靜動態特性[3-4]。此外，材料的表面硬度作為結合面的重要材料力學性能參數，與接觸剛度有十分密切的關系[5]。由于宏觀連續介質理論已不再適用于納米尺度接觸問題的研究，而實驗法也難以觀察到材料納米尺度下內部的原子運動情況，且材料加工費用昂貴、耗時長，因此分子動力學(molecular dynamics，MD)模擬方法[6-7]憑借其成熟的理論和高效的計算等優點成為了微納尺度下接觸研究中應用最為廣泛的數值方法[8]。

針對結合面的接觸研究，已有許多學者建立了很多分子動力學數值接觸模型。在單微凸體接觸模擬方面：黃健萌等[9]開展了兩種不同形狀的金剛石壓頭與單晶銅基體的黏著接觸分子動力學仿真，發現由于半球形壓頭底部表面積比圓錐形壓頭底部表面積更大，導致半球形壓頭與基體之間的引力更大，從而更易產生黏著接觸現象并且下壓過程中半球形壓頭所受到的法向接觸力更大；YANG等[10]在微觀尺度下用原子模擬的方法對連續介質接觸理論進行了驗證，模擬了名義光滑的球形壓頭與空間平面基體的接觸與分離過程，發現最大法向接觸力的大小與壓頭的尺寸相關，并且在壓入的過程中伴隨有不同位錯的成核、運動、湮滅；MOJUMDER等[11]開展了純鋁的納米壓痕分子動力學模擬，發現隨著壓頭尺寸的增大，壓頭與基體的接觸面積變大，壓痕力和位錯密度相應更大，并且硬度和彈性模量均增大；樊江等[12]建立了球狀微凸體剪切運動過程的分子動力學模型，分析了單晶銅微凸體在剪切運動過程中切向接觸力和法向接觸力隨微凸體半徑的變化規律，發現單晶銅在真空中發生干摩擦時滑動摩擦因數只與微凸體形狀參數有關。

在多微凸體粗糙表面的接觸模擬方面：SANTHAPURAM等[13]建立了包含有球形微凸體和圓柱形微凸體的粗糙鋁表面分子動力學接觸模型，對比了這兩種粗糙表面的接觸力學響應。黃仕平等[14]采用分子動力學-格林函數法模擬了高斯粗糙表面的接觸過程，探討了由微球體組成的粗糙面的接觸力學特性，模擬發現，在微凸體高度符合高斯分布的情況下，接觸面積與接觸力成線性關系。YASTREBOV等[15]采用分子動力學方法模擬了分形粗糙表面與彈性半空間體的無摩擦以及無黏著接觸過程，研究了Nayak參數對接觸面積演化的作用，結果表明，在給定的歸一化壓力下，接觸面積隨Nayak參數變化呈現對數遞減趨勢。

上述研究均主要著眼于粗糙表面的微觀幾何形貌結構對其接觸功能屬性的影響機理，并未針對接觸表面的材料力學性能參數(如表面硬度)對接觸性能的影響展開研究。本文基于混合勢函數(EAM和Morse)，采用Verlet算法對剛性半球形金剛石壓頭與彈塑性多晶銅基體的分子動力學黏著接觸與分離過程進行了仿真，從能量角度分析了不同晶粒尺寸對多晶銅基體的納米硬度和接觸剛度的影響規律。通過揭示多晶銅基體亞表面的變形機理，從原子尺度解釋了上述規律現象產生的原因。

1 分子動力學建模與計算方法

1.1 建模與仿真

剛性球體與彈塑性半空間體的接觸模型是接觸力學方面研究所采用的經典接觸模型之一[16]。圖1為本文構建的分子動力學接觸模擬仿真模型的XY平面示意圖。該模型由上下兩部分構成，上半部分為金剛石壓頭，鑒于典型Berkovich壓頭的尖端在納觀尺度下為半球形，因此本模擬中壓頭的形狀選取為半球形，綜合考慮壓入深度與基體尺寸，選取壓頭的半徑為3.5 nm；下半部分為多晶銅基體，其三個維度的尺寸分別為50a×50a×30a，a為晶體銅的晶格常數。由于金剛石的硬度遠遠大于銅材料的硬度，因此在接觸模擬過程中將金剛石壓頭視作剛體，即不需要用勢函數來描述C原子之間的相互作用力。多晶銅基體采用Voronoi-tesselation方法[17]生成，由許多相同晶粒尺寸的單晶銅晶粒組成，其平均晶粒尺寸G分別為12 nm、10 nm、9 nm、8 nm、7 nm和6 nm，組成多晶銅的各單晶銅晶粒取向為隨機取向。在模擬過程中多晶銅基體會發生彈塑性變形。將多晶體基體自上而下依次分為三個部分，依次為牛頓層、恒溫層和固定層，其中基體的彈塑性變形將發生在牛頓層。圖1中，hmax為壓頭在模擬過程中將要下壓的最大壓深，壓入方向豎直向下。金剛石壓頭與多晶銅基體在正式模擬前的初始間距為1 nm，此距離保證了在模擬開始前壓頭中的C原子與基體表層的Cu原子之間無相互作用力。

圖1 金剛石壓頭與多晶銅基體的分子動力學接觸模型XY平面示意圖Fig.1 Molecular dynamics contact model of diamond indenter and polycrystalline copper matrix schematic diagram of XY plane

勢函數影響著模擬的精度，決定著模擬結果的可靠性。由于物質的原子和分子構造不同，造成原子間的相互作用類型不同，因此必須采用不同類型的原子間相互作用勢。金剛石立方晶體結構的晶格常數a為0.357 nm，模擬過程中忽略C—C原子之間的相互作用。多晶銅的組成晶粒為面心立方(face center cubic，FCC)晶體結構，其晶格常數a為0.3615 nm，對Cu—Cu原子之間的相互作用力采用EAM勢函數[18]進行描述。EAM勢是一種多體勢，能夠較好地描述金屬晶體的結構特性和力學特性，在金屬體系的分子動力學模擬中應用廣泛。C—Cu原子之間的相互作用力采用Morse勢函數[19]進行描述。Morse勢是一種對偶勢，適用于表征不同原子間的相互作用力，廣泛應用于金剛石與金屬原子間相互作用力的描述。Morse勢函數表達式中包含3個可調參數，參數值的設置參考了已有相關研究報道[20-24]，分別為D=0.087 eV，α=51.4 nm-1，r0=0.205 nm，其中D為結合能系數，α為勢能曲線梯度系數，r0為原子間作用力為0時的平衡態原子間距。多晶銅基體的X、Y、Z軸的晶向分別取[100]、[001]和[010]。在模擬系統的X、Z方向上均選取周期性邊界條件以消除尺寸效應，在Y方向上選取自由邊界條件。在模擬過程中，固定層原子的位移與受力在每個時間步上均不更新；對屬于恒溫層原子的信息采用Nose-Hoover熱浴耦合算法[25]進行控制，可使得恒溫層的溫度在模擬期間保持為設定值；對牛頓層的原子采用Verlet算法[26]進行牛頓運動方程的求解并更新每時刻原子的位置與速度。模擬的時間步長取值為1 fs。

本文模擬在常溫下的接觸行為，正式模擬開始之前，先使系統經歷弛豫過程，采用Berendsen熱浴耦合算法[27]將體系加熱至模擬溫度300 K，弛豫過程中系統狀態參數的變化趨勢表明，30 ps的總弛豫時間可使得體系的溫度和總能量達到充分穩定的狀態。模擬開始時，對金剛石壓頭施加位移載荷，研究結果表明，加載的速度對仿真的結果影響不大[28]，因此使壓頭沿著Y軸負方向以100 m/s的速度做勻速運動，并逐漸壓入多晶銅基體，最大壓深為2 nm，基體在此下壓速度及壓深下可達到一定的塑性變形量。在達到最大壓深后，對金剛石壓頭施加反向且相同大小的位移載荷，使其逐漸與基體脫離接觸并勻速回到初始位置。模擬過程中記錄多晶銅基體對金剛石壓頭施加的法向載荷F。本文的模擬過程采用開源分子動力學軟件LAMMPS完成。

1.2 壓痕功法

圖2 典型加載-卸載曲線Fig.2 Typical loading-unloading curve

如圖2所示，壓痕功法[29]認為表面壓痕的形成過程是能量消耗或做功的過程。

在接觸過程中，加載曲線對X軸進行積分得到的面積是壓頭做的總功Wt，在分離過程中，卸載曲線對X軸進行積分得到的面積是彈性功We，二者之間的差值就是塑性功Wp，可表示為

Wp=Wt-We

(1)

其中，總功Wt和彈性功We可分別表示為

(2)

(3)

式中，h為壓頭下壓的位移；Fload(h)表示加載曲線；Funload(h)表示卸載曲線。

納米硬度Hp定義為

(4)

式中，Fmax為接觸過程中壓頭受到的最大法向接觸力；k為與壓頭形狀有關的常數，在本模擬中取k=0.0408。

1.3 彈性接觸剛度

由圖2可以看出，隨著壓入深度的不斷增大，壓頭受到的法向接觸力也隨之不斷增大，當壓入深度到達設定值后壓頭停止下壓。在完全脫離接觸后，基體壓深會回到一個固定值，將此時的深度稱為殘余深度hf，即壓頭在材料表面上留下的永久塑性變形。在卸載曲線頂端的起始點處作一條切線，將該切線的斜率稱為接觸剛度S。

采用如下冪函數來擬合卸載曲線頂端部分：

Funload(h)=β(h-hf)m

(5)

其中，β、m為卸載曲線的擬合參數。則接觸剛度S的計算公式為

S=βm(hmax-hf)m-1

(6)

1.4 Hertz接觸理論計算

在不考慮表面黏著力存在的情況下，Hertz接觸理論給出了壓頭與基體彈性接觸時法向載荷F、壓頭半徑R、壓入深度h之間的函數關系式：

(7)

(8)

式中，E*為復合彈性模量；Es、Ei分別為基體和壓頭的彈性模量；υs、υi分別為基體和壓頭的泊松比。

本模擬中，半球形金剛石壓頭與多晶銅基體的彈性模量和泊松比參數分別取值為：Ei=1141 GPa，Es=115 GPa；υi=0.07，υs=0.25。

2 結果與討論

2.1 法向接觸力的變化

圖3為不同晶粒尺寸G下法向接觸力隨壓頭下壓位移的變化曲線，其中黑色和藍色曲線分別為由分子動力學(MD)模擬得到的金剛石壓頭所受到的多晶銅基體在接觸和分離過程中施加的法向接觸力，紅色曲線為接觸過程中由Hertz理論擬合得到的法向接觸力。由圖3可以看出，載荷-位移曲線呈現波動的狀態，載荷的波動與基體銅的晶格變形、晶格重構還有材料的非晶相變有密切的關系。

圖3a～圖3f中的A點代表壓頭即將開始壓入多晶銅基體的臨界點。可以看出在不同晶粒尺寸條件下，壓頭未壓入基體前均出現了數值很小的負方向的力，這表明此時金剛石壓頭與多晶銅基體原子之間的吸引力起主要作用。文獻[30]指出：基體會因為引力過大而引起基體表面部分原子發生納米尺度接觸所特有的“突跳”現象。當兩表面靜止接觸或兩表面間隙處于納米量級時，由于黏著力的作用，使得兩表面黏附而不易分離。如圖4所示，不同晶粒尺寸條件下金剛石壓頭所受到的最大黏著力基本不變。這是因為黏著力的大小主要與表面能有關，相同的壓頭尺寸使得壓頭所接觸的表面積相同，所以最大黏著力的變化基本不大。

(a)G=12 nm (b)G=10 nm (c)G=9 nm

圖4 最大黏著力隨晶粒尺寸的變化曲線Fig.4 Variation curve of maximum adhesion force with grain size

隨著壓頭逐漸壓入基體，金剛石壓頭受到的多晶銅基體施加的法向接觸力逐漸增大。值得注意的是，無論何種多晶銅晶粒尺寸的條件下，在圖3a～圖3f中B點之前，由MD模擬得到的法向接觸力要比Hertz理論擬合得到的法向接觸力大。這是因為在壓深很小時，基體處于彈性變形階段，Hertz理論的適用對象是表面光滑的彈性變形體，接觸過程連續且無摩擦，不存在相對滑動，可忽略黏著力和表面力的存在，而MD模擬在原子尺度下包含了上述因素。當壓深逐漸增大且超過B點處壓深后，基體逐漸開始由彈性變形過渡到塑性變形，Hertz理論忽略了塑性變形的影響，僅適用于彈性形變，從而使得MD模擬的法向接觸力比Hertz理論擬合的法向接觸力小。

(9)

表1 不同晶粒尺寸條件下由MD模擬計算得到的復合彈性模量Tab.1 The composite elastic modulus calculated by MD simulation under different grain size conditions

圖5 最大法向接觸力隨晶粒尺寸的變化曲線Fig.5 Variation curve of maximum normal contact force with grain size

圖3a～圖3f中的C點處的法向接觸力為MD模擬得到的不同晶粒尺寸條件下半球形金剛石壓頭受到的多晶銅基體施加的最大法向接觸力。根據圖5中的曲線可以看出，當晶粒尺寸大于10 nm時，隨晶粒尺寸的減小，壓頭受到基體施加的最大法向接觸力呈現增大的趨勢，表現出細晶強化效應(即霍爾-佩奇現象)；而當晶粒尺寸小于10 nm時，最大法向接觸力隨晶粒尺寸的減小呈現減小的趨勢，顯示出明顯的反霍爾-佩奇現象。這是因為霍爾-佩奇公式是根據位錯塞積的強化作用而提出的，當晶粒尺寸為納米級時(小于10 nm)，晶粒中可存在的位錯極少，故霍爾-佩奇公式就不再適用。

當壓頭壓入基體達到最大壓深時，對壓頭施加反向且相同大小的位移載荷，壓頭與基體隨即進入分離階段。由圖3a～圖3f可以看出，不同晶粒尺寸條件下，壓頭與基體開始分離時法向接觸力迅速減小，當壓深回復到1.5 nm左右時，壓頭受到的法向接觸力均出現負值的情況，這是由于卸載過程中多晶銅基體的原子與金剛石壓頭原子之間的吸附作用導致的。隨著壓頭的上升吸附效應逐漸增強，最終壓頭原子會離開基體原子的作用范圍，即壓頭不再受力。

2.2 接觸過程的能量分析

圖6所示為根據壓痕功法計算得到的不同晶粒尺寸條件下壓頭與基體接觸以及分離過程中體系能量的變化情況。可以明顯看出，在模擬過程中，體系的總能量Wt、彈性能We以及塑性能Wp均與圖5所示的壓頭受到基體施加的最大法向接觸力的變化趨勢一致，即表現出先增大而后又減小的變化趨勢。

圖6 不同晶粒尺寸條件下總能量、彈性能、塑性能的對比Fig.6 Comparison of total energy，elastic energy，and plastic energy under different grain sizes

殘余深度hf的大小可反映壓頭在接觸過程中基體塑性變形的程度以及分離過程中彈性回復量的大小。采用Alpha-shape算法[31]對卸載后的多晶銅基體表面形貌進行擬合，得到分離接觸后的基體表面形貌，如圖7所示。

(a)G=12 nm (b)G=10 nm

通過式(4)和式(6)計算得到的不同晶粒尺寸下多晶銅的納米硬度以及接觸剛度見表 2。由表 2可知，當某一晶粒尺寸條件下的殘余深度hf相對越大，則該晶粒尺寸對應的多晶銅的納米硬度Hp相對較小，接觸剛度S相對較大。

表2 不同晶粒尺寸條件下殘余深度、納米硬度、接觸剛度計算值Tab.2 Calculated values of residual depth，nano-hardness，and contact stiffness under different grain size conditions

接觸剛度S是表征材料抵抗彈性變形能力的指標。由圖8a和圖8b可以看出，隨著晶粒尺寸的變化，接觸剛度S與彈性能We成負相關。這是由于彈性能We是基體在彈性變形中所儲存的能量，彈性能越大，表明彈性變形量越大，材料抵抗彈性變形的能力越弱，接觸剛度越低。納米硬度Hp反映了材料抵抗局部塑性變形的能力。從塑性變形的本質來說，位錯的產生及滑移是出現塑性變形的主要原因。塑性能Wp表示在外載作用下晶格畸變以及位錯滑移所產生的所有畸變能之和。然而，當晶粒尺寸為納米級(小于10 nm)時，晶粒中可存在的位錯極少，位錯塞積強化效應不再適用，因此導致納米硬度Hp與塑性能Wp成正相關，如圖8c和圖8d所示。綜合圖8b與圖8d可知，納米硬度Hp與接觸剛度S的變化總體成負相關。當晶粒尺寸為納米級時，納米硬度越小，接觸剛度越大。

(a)彈性能隨晶體尺寸變化情況 (b)接觸剛度隨晶粒尺寸變化情況

2.3 多晶銅亞表面變形機理

為研究多晶銅基體在接觸過程中的彈塑性變形過程，采用公共近鄰分析(common neighbor analysis,CNA)方法分析了多晶銅基體在與金剛石壓頭彈塑性接觸過程中位錯的形核以及擴散規律。圖9所示為不同晶粒尺寸條件下各多晶銅基體的彈塑性變化過程，從左至右依次為金剛石壓頭壓入多晶銅基體的深度h=0、0.3 nm、0.5 nm、1.0 nm、1.5 nm和2.0 nm時基體的內部原子結構形貌。圖中的深藍色原子代表晶界原子，淺藍色原子代表具有正常FCC晶體結構的銅原子，綠色原子代表位錯原子，橙色原子代表其他缺陷原子。由圖9可以看出，隨著壓頭與基體逐漸開始接觸，當壓頭壓入基體的深度h=0.3 nm時，基體上靠近壓頭正下方的部分區域開始出現位錯形核，如圖9中紅色箭頭所示；并且隨著壓頭的持續壓入，基體內萌生的位錯開始沿著晶粒的{111}滑移系滑移擴散。當壓入深度較小時，位錯的滑移運動會受到晶界的阻礙，表現為運動到晶界的位錯會被晶界吸收，如圖9中黃色箭頭所示；當壓入深度逐漸增大時，部分晶界位置會不斷地累積內應力以及原子勢能，而當累積到一定程度時，部分位錯會在晶界附近形核并穿過晶界繼續滑移擴散，從而使得塑性變形過程持續進行，如圖9中粉色箭頭所示。

(a)G=12 nm

由圖9還可以看出，較小晶粒尺寸的多晶銅擁有數量較多的晶界，由于晶界具有較大的彈性應變能，因此位錯的運動會被約束在壓頭周圍的若干晶粒中，導致較小晶粒尺寸的多晶銅的位錯分布范圍較小，產生的位錯較少。由此導致位錯塞積強化效應的不適用，使得具有較小晶粒尺寸(小于10 nm)的多晶銅材料的納米硬度相對較低、接觸剛度相對較大。

3 結論

(1)不同晶粒尺寸下，多晶銅黏著接觸過程中，最大黏著力基本保持不變。最大法向接觸力隨晶粒尺寸的減小表現出先增后減的趨勢。

(2)納米硬度和接觸剛度分別與接觸與分離過程中的塑性能和彈性能的變化相關，且納米硬度與接觸剛度變化總體成負相關。

(3)晶粒尺寸越小的多晶銅擁有數量相對更多的晶界，晶界數量會對材料的納米硬度產生影響，進而影響其接觸力學特性。