基于分子動力學-格林函數法的分形粗糙表面摩擦行為研究1)

黃仕平 陳 梟 蕭明強

(華南理工大學土木與交通學院,廣州 510640)

引言

接觸摩擦現象無處不在,其廣泛發生在機械、能源和電子等領域,合理利用摩擦學原理,能有效地減少能源消耗和提升作業效率[1].宏觀上光滑和平整的表面,實際上是由許多不規則的微凸體組成的[2-3].當兩個粗糙固體表面接觸時,實際接觸只發生在表觀面積的很小一部分,真實的接觸面積由大量離散的接觸團簇組成,接觸團簇的形狀、大小、分布和相互作用是揭示接觸摩擦行為的關鍵,其內部接觸摩擦機理相比光滑表面更加復雜,因此表面接觸和摩擦問題成為力學領域極具挑戰性的課題[4].

近年來,學者們在接觸摩擦研究領域也已經取得了許多重要的成果.由Greenwood 等[5]提出的GW 接觸模型是統計接觸模型的典型代表,其假設粗糙面由許多半徑相同但高度不同的球狀微凸體組成,以赫茲理論為基礎,對接觸力和摩擦力進行統計意義上的累加.Majumdar 等[6]基于W-M 分形函數提出的MB 接觸模型是分形接觸模型的典型代表,其結果更接近粗糙表面的非穩定隨機特性.Chang等[7]考慮微凸體的塑性變形,并給出單個微凸體發生塑性變形時的接觸面積和接觸載荷計算公式.Pérez-Ràfols 等[8]研究具有非高斯高度分布的自仿射分形表面的接觸剛度.Popov 等[9]使用降維方法模擬分形粗糙表面與彈性體之間的摩擦力,發現摩擦系數與表面輪廓的均方梯度成正比.李玲等[10]采用二次函數建立微凸體接觸半徑與接觸變形的解析關系,并重新推導了微凸體發生彈性、彈塑性和塑性變形的接觸表達式.馮燕等[11]考慮熱應力的影響,建立了粗糙表面的熱彈塑性接觸模型,用于表征粗糙表面的熱力學特性.陳少華等[12-13]研究剛性球體與彈性梯度半空間的接觸和摩擦問題.占旺龍等[14]基于Iwan 模型的非線性力學性質,推導出臨界摩擦力分布函數.周華等[15]考慮接觸過程中微凸體的相互作用,建立分形接觸模型,并推導了粗糙表面接觸的切向位移和切向載荷的關系.

隨著科學計算方法的發展,采用數值方法更準確模擬粗糙表面的接觸摩擦行為已成為研究熱點.Kogut 等[16]通過有限元模擬了球體與平面的彈塑性接觸問題,并構建了KE 彈塑性模型,將赫茲理論拓展到完全塑性接觸.Sofuoglu 等[17]運用有限元法研究了二維剛性粗糙體與彈塑性平面的熱滑動接觸問題.Sellgren 等[18]利用有限元軟件建立了考慮工程表面特性的接觸模型,發現微凸體的高度分布對接觸剛度有顯著影響,但其曲率對接觸剛度影響較小.Nadimi 等[19]利用邊界元法研究均方根高度、赫斯特指數和分形維數等參數對實際接觸面積的影響.阮曉光等[20]使用Abaqus 軟件研究粗糙表面的接觸壓力和摩擦力分布.王權等[21]基于分子動力學模擬探究摩擦深度對銅鎳合金摩擦行為的影響.范立峰等[22]采用彈塑性有限元方法分析了界面分形參數對法向接觸剛度的影響.宋劍鋒等[23]使用Abaqus軟件模擬兩個粗糙表面的接觸過程,分析了表面位移和表面粗糙度對接觸載荷和接觸面積的影響.

雖然已經有許多學者對粗糙表面的接觸摩擦行為進行了研究,但大部分研究都忽略了摩擦過程中接觸團簇之間的相互作用和影響,并且大規模的數值模擬計算開銷大.為了減少計算開銷,文獻[25-27]提出分子動力學-格林函數法(GFMD),利用格林函數起到了降維的效果,并已成功運用于表面接觸和摩擦分析.分子動力學從分子和原子尺度出發分析接觸摩擦機理,能獲得很多連續介質力學難以模擬的力學行為[28-29].

本工作以分子動力學-格林函數法為工具,首先利用二維傅里葉變換方法生成分形粗糙表面,建立分子動力學模型;其次進行分形粗糙表面和彈性光滑平面的接觸模擬,利用廣度優先搜索算法識別接觸團簇;最后模擬分形粗糙表面的滑動摩擦,從不同尺度計算最大摩擦系數和摩擦力,利用影響矩陣法分析接觸團簇之間的相互作用,并結合彈性力學理論模型進一步分析該相互作用的影響因素,研究接觸團簇之間的距離和接觸團簇面積對相互作用的影響.本文結果可為粗糙表面的界面分析和優化提供一定的理論依據.

1 分子動力學-格林函數法理論

分子動力學-格林函數法的核心是在表面邊界層采用分子動力學方法模擬其表面力學效應,表面邊界層以外則采用格林函數法直接計算其彈性響應.該方法的優點是既沒有忽略表面邊界層復雜及分子級別的力學效應,同時又利用格林函數法提高了計算效率,從而使得大規模表面接觸摩擦行為的研究成為可能.

兩個粗糙表面的接觸可以簡化為一個彈性光滑面與一個剛性粗糙面之間的接觸[30],如圖1 所示,此時彈性光滑面的彈性模量根據以下公式換算

圖1 剛性粗糙面與彈性光滑面的接觸Fig.1 Rigid rough surface in contact with elastic smooth surface

式中,E1和E2分別是上下表面的彈性模量,v1和v2分別是上下表面的泊松比.

邊界層和基體層的動力學方程為

式中,m表示原子的質量矩陣,uiα表示第 α層的第i個原子的位移向量,Diαjβ表示力常數矩陣,fi0表示原子所受的力矢量.

通過傅里葉變換簡化動力學方程,結果如下

如果基底的變形很小,此時假設整個基底的位移是施加在邊界層上的力的線性函數

式中,G表示格林函數矩陣.只要提前計算出G00,即可用來表示邊界層原子的位移,格林函數矩陣G的具體計算可參考文獻[26].

2 計算模型及方法

2.1 分子動力學模型

本文采用Lammps 軟件模擬計算分子動力學模型,對Lammps 軟件進行二次開發,增加了格林函數法的計算模塊[31].整體模型由一個剛性分形粗糙表面和一個彈性光滑表面組成,如圖2 所示.采用的原子類型為gfmd 類型,由二次開發所得,相鄰原子間距d=1.12σ,模型尺寸為 512d×512d.彈性光滑表面采用各向同性連續彈性體的格林函數解[32-33],等效彈性模量E*=4ε/σ3,泊松比為0.3,四周設置周期性邊界條件來模擬無限大半空間.分形粗糙表面通過二維傅里葉變換法生成,表面特征由4 個參數控制,分別是空間頻率上限值 ωH、空間頻率下限值 ωL、赫斯特指數H、幅度的標準差P.本文的分形粗糙表面 模型取 ωH=1/48,ωL=1/128,H=0.5,P=0.5,其高度的均方根為0.99,斜率的均方根為0.05,具有很好的普適性.

圖2 分子動力學模型Fig.2 Molecular dynamics model

2.2 模擬方法

本文采用位移加載控制模型的接觸和摩擦過程,模擬時間步長為0.005 ps,位移加載步長為 0.01σ,接觸過程的時間步為1000 步,滑動摩擦過程的時間步為500 步,沿x方向滑動.每一位移步都要進行迭代分析,直到彈性表面原子受力平衡,每次迭代的最大計算步為 5.0×104步,迭代收斂準則為原子的1-范數小于等于 0.01ε/σ.

采用Lennard-Jones (LJ)勢來模擬原子間的相互作用[34],其表達式如下所示

式中,r表示原子間的距離,ε 和σ為勢能參數,采用國際單位制會導致物理量的數值非常小,可能會引起較大的舍入誤差,故本文模擬采用LJ 約化單位制.在LJ 約化單位制中,ε=1,σ=1,計算中分別取其為勢能、距離的基本單位,力的單位為 ε/σ.LJ 勢只考慮原子間的吸引力和排斥力,不考慮電荷作用,當r=21/6σ ≈1.12σ時,原子間相互作用力為0,當r＜1.12σ時,原子間作用力為排斥力,當r＞1.12σ時,原子間作用力為吸引力.本文模型不考慮原子的黏附作用,因此LJ 勢能的截斷距離取 1.12σ.

摩擦過程中的變化規律和特征可以通過摩擦系數、摩擦力等參數表征.在模擬過程中,本文從3 個尺度來分析界面的摩擦行為: (1)原子尺度,反應接觸原子的個體效應;(2)接觸團簇尺度,反應接觸原子的群體效應;(3)界面尺度,反應接觸團簇的群體效應.

定義原子i受到的法向接觸力Fi為其z方向的分力,受到的摩擦力fi為其x方向的分力;界面受到的法向接觸力F和摩擦力f分別為界面所有原子的法向接觸力和摩擦力之和;接觸團簇受到的的法向接觸力Fcc和摩擦力fcc分別為該接觸團簇所有原子的法向接觸力和摩擦力之和,表達式如下

式中,n表示界面總原子數,ncc表示該接觸團簇總原子數.

引入摩擦系數,其值定義為摩擦力與法向接觸力的比值,表達式如下

式中,μi表示原子i的摩擦系數,μcc表示接觸團簇的摩擦系數,μ表示界面的摩擦系數.

2.3 接觸團簇識別算法

接觸團簇的大小、形狀、分布是研究界面接觸摩擦行為的重要特性.定義z方向的受力fz＞0的原子為接觸原子,為了研究摩擦過程中接觸團簇之間的相互作用,需要分別識別組成各個接觸團簇的接觸原子,對其進行單獨的編號.本文采用廣度優先搜索算法,避免了遞歸搜索的棧溢出問題.定義相鄰的接觸原子屬于同一個接觸團簇,搜索步驟如下.

(1) 遍歷整個模型原子,當遇到接觸原子時,從此原子開始做廣度優先搜索.創建一個輔助隊列,將該原子加入隊列,此時接觸團簇數量加1,且在廣度優先搜索中刪除此原子;

(2) 廣度優先搜索時,判斷隊列首部原子是否是接觸原子,若是則刪除該原子,并將該原子的鄰接原子加入隊列,若不是則跳過該原子;

(3) 循環出隊列的首原子,直到整個隊列為空,則該接觸團簇遍歷完畢.

3 結果與討論

3.1 接觸團簇分布

隨著壓入深度不斷增加,接觸面積也不斷增大,當接觸面積達到5%時,采用廣度優先搜索算法識別到17 個接觸團簇,分別對其進行編號,其結果如圖3 所示.

圖3 接觸團簇分布Fig.3 Contact cluster distribution

3.2 摩擦行為

3.2.1 摩擦系數

采集每一位移步下所有原子的法向接觸力和摩擦力,以此求得當前位移步下,原子、接觸團簇和界面的摩擦系數.之后再根據每一位移步的結果,統計得到在摩擦過程中原子、接觸團簇和界面的最大摩擦系數,其結果分別如圖4 所示.

圖4 摩擦系數Fig.4 Friction coefficient

對原子尺度和接觸團簇尺度的最大摩擦系數分別取平均值,最終得到原子尺度的平均摩擦系數為0.59,接觸團簇尺度的平均摩擦系數為0.54,而界面尺度的摩擦系數為0.35.結果表明原子摩擦系數 ＞接觸團簇摩擦系數 ＞ 界面摩擦系數,摩擦系數從小尺度到大尺度逐漸減小.

3.2.2 摩擦力-位移曲線

分別統計3 個尺度下的摩擦力值,得到摩擦力-位移曲線如圖5 所示.從圖5 可以看到,在滑動過程中摩擦力呈現周期變化的趨勢,波動周期約等于原子間距 1.12σ.界面尺度下,如圖5(a)所示,當位移達到 0.92σ時,摩擦力達到最大,為 314.52ε/σ;接觸團簇尺度下,如圖5(b)所示,用圓點標記了部分接觸團簇的最大摩擦力值,接觸團簇1 在位移達到 0.92σ時達到最大摩擦力,為 62.99ε/σ,而接觸團簇5 在位移達到 0.70σ時達到最大摩擦力,為14.39ε/σ,接觸團簇并非同時達到最大摩擦力,而是發生局部滑移;原子尺度下,如圖5(c)所示,原子也是分批達到最大摩擦力,從而發生滑移.

圖5 摩擦力-位移曲線Fig.5 Friction-displacement curve

整體滑移模型認為,微凸體整體滑移,摩擦力是同時到達峰值的.但在本文的分子動力學模擬中,原子和接觸團簇都是分批滑移的,這在原子和接觸團簇層面解釋了局部滑移現象.在摩擦過程中,接觸團簇并非同時達到最大摩擦力,而是發生局部滑移現象.發生局部滑移的接觸團簇會對其他接觸團簇產生相互作用,加速其他接觸團簇滑移,因此整體滑移模型預測的摩擦力,實際上是分子模擬結果的上限值,而摩擦力的下限值,則取決于接觸團簇的大小、分布等,這可能是界面設計的關鍵因素.

3.3 影響矩陣法求相互作用

3.3.1 影響矩陣法

為了進一步探究摩擦過程中接觸團簇之間的相互作用,本文使用影響矩陣法來建立各個接觸團簇之間的相互作用機制.影響矩陣法引入了目標向量、受調向量、施調向量和影響矩陣的概念,本文目標向量即接觸團簇的摩擦力F,受調向量即當前接觸團簇的摩擦力F0,施調向量即接觸團簇的位移向量D.

定義影響矩陣K,矩陣中的元素kij表示第j個接觸團簇滑動單位位移時對第i個接觸團簇摩擦力的影響.其計算方法如下: 給第i個接觸團簇施加單位位移,同時固定其他接觸團簇,計算在單位位移下各接觸團簇的摩擦力,記為Ki={k1i,k2i,···,kni}T,同理得到K1～Kn,即可得到接觸團簇之間的影響矩陣為

當模型滿足線形疊加原理時,則目標向量等于受調向量加上影響矩陣與施調向量之積,即接觸團簇的摩擦力F等于當前接觸團簇的摩擦力F0加上影響矩陣K與位移向量D之積,表達式如下

3.3.2 影響矩陣計算

以接觸團簇1 為例,對其施加單位位移,同時固定其他接觸團簇,計算在單位位移下各接觸團簇的摩擦力.為了更好的體現接觸團簇之間相互作用的大小,對影響矩陣進行歸一化得到K1={1,0.097,0.031,0.047,0.040,0.026,0.036,0.043,0.012,0.014,0.010,0.084,0.026,0.011,0.041,0.028,0.005}T,其結果如圖6 所示,同理計算得到K2～K17.

圖6 接觸團簇相互作用Fig.6 Interaction of contact clusters

3.3.3 摩擦力計算

摩擦力計算結果如圖7 所示.根據影響矩陣法計算得到的最大摩擦力為 331.20ε/σ,與GFMD 模擬結果較為接近,說明在一定程度上,可以用影響矩陣法來模擬接觸團簇之間的相互作用.Ki′j表示不考慮局部滑移的影響,根據各個接觸團簇的滑移時刻和對應的影響矩陣系數計算得到的摩擦力,最大摩擦力為 379.11ε/σ,比GFMD 模擬結果大了20%.由此可見在摩擦過程中,接觸團簇之間的相互作用對滑動摩擦力的影響很大,由于局部滑移的存在,先滑移的接觸團簇會對其他接觸團簇產生力的作用,并加速其達到局部最大摩擦力,從而使得系統整體的摩擦力下降,下降的程度受到接觸團簇分布和大小等的影響.

圖7 摩擦力計算Fig.7 Calculation of friction force

3.3.4 影響因素分析

kij的計算結果受到接觸團簇形狀、分布和距離等的影響,目前很難從理論上明確計算.本文建立兩個簡單模型,一個是有兩個球狀接觸團簇的GFMD數值模型,另一個是基于彈性力學中切向集中力作用下彈性半空間的變形解的理論模型,嘗試分析接觸團簇之間的距離和面積對kij的影響.

(1) 數值模型

數值模型的正視圖如圖8 所示.原子初始間距d=1.12σ,彈性光滑平面的尺寸為 512d×512d,兩個球體的半徑均為 25d,間距為150d.根據赫茲理論,一個剛性球體與彈性半空間體接觸時會生成圓形的接觸團簇,其壓力分布形式為

圖8 GFMD 模型Fig.8 GFMD model

(2) 理論模型

考慮有一個彈性半空間體,在半徑為a的圓上作用切向非均布力,其中q0表示圓心處的切向力,r表示離圓心的距離.則任意點B點的位移可由彈性力學推導得到.

當B點在圓外時,B點的切向位移為

當B點在圓內時,B點的切向位移為

(3)kij計算

由于微觀模型和宏觀模型的差異,本文不研究kij的絕對值,轉為研究kij/kii的相對關系,kij/kii體現了第j個接觸團簇滑移對第i個接觸團簇的加速效果.數值模型的kij計算同上述情況,理論模型中kij/kii計算如下式

(4) 距離的影響

分析接觸團簇之間的距離對kij的影響.根據上述方法,保持接觸團簇的半徑不變,改變兩個接觸團簇之間的距離,分別計算GFMD 數值模型和彈性力學理論模型的k12/k11大小,結果如圖9(a)所示.從圖中可以看出,kij隨著距離的增加而減小,大致呈反比趨勢.

圖9 kij 變化Fig.9 Changes of kij

(5) 面積的影響

分析接觸團簇面積對kij的影響.保持接觸團簇之間的距離不變,改變第二個接觸團簇面積大小,分別計算GFMD 數值模型和彈性力學理論模型的k12/k11大小,結果如圖9(b)所示.從圖中可以看出,kij隨著面積的增加而增大,大致呈線性趨勢.

4 結論

本文基于分子動力學-格林函數法建立了分形粗糙表面的微觀模型,分析摩擦過程中的力學特性,并利用影響矩陣法建立接觸團簇的相互作用機制,在接觸團簇層面解釋了局部滑移現象及其對界面摩擦行為的影響,主要結論如下.

(1) 分別從原子、接觸團簇和界面3 個尺度分析摩擦行為,結果表明原子摩擦系數 ＞ 接觸團簇摩擦系數 ＞ 界面摩擦系數,摩擦系數從小尺度到大尺度逐漸減小.并且在摩擦過程中,摩擦力呈周期性波動,接觸團簇并非同時達到最大摩擦力,而是發生局部滑移,先滑移的接觸團簇會加速其它接觸團簇的滑移.因此,整體滑移模型預測的摩擦力實際上是分子模擬結果的上限值.

(2) 利用影響矩陣法建立了接觸團簇的相互作用機制.根據影響矩陣計算得到的最大摩擦力與GFMD 模型的結果較為一致,說明影響矩陣法的正確性,而不考慮局部滑移影響計算得到的最大摩擦力比GFMD 模型結果大20%,可見接觸團簇之間的相互作用影響之大,不可忽略.

(3) 利用球狀數值模型和彈性力學方法,探究接觸團簇之間的距離和接觸團簇面積對影響矩陣元素kij的影響.結果表明kij隨著距離的增大而減小,大致呈反比趨勢;kij隨著面積的增大而增大,大致呈線性趨勢.