基于分子動(dòng)力學(xué)
--格林函數(shù)法的微凸體接觸數(shù)值分析1)

黃仕平吳杰胡俊亮鄭恒斌王衛(wèi)鋒,2)

?(華南理工大學(xué)土木與交通學(xué)院，廣州510640)?(清華大學(xué)摩擦學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100084)??(華南農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，廣州510642)

生物、工程及交叉力學(xué)

基于分子動(dòng)力學(xué)
--格林函數(shù)法的微凸體接觸數(shù)值分析1)

黃仕平?,?吳杰?胡俊亮?鄭恒斌??王衛(wèi)鋒?,2)

?(華南理工大學(xué)土木與交通學(xué)院，廣州510640)?(清華大學(xué)摩擦學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100084)??(華南農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，廣州510642)

表面接觸是摩擦的先決條件，其真實(shí)接觸面積、壓應(yīng)力大小、空間分布等一直是接觸力學(xué)關(guān)注的核心問(wèn)題.采用分子動(dòng)力學(xué)--格林函數(shù)法(GFMD)模擬粗糙面的接觸過(guò)程，驗(yàn)證了其在大規(guī)模接觸分析中的高效及準(zhǔn)確性，同時(shí)探討了由微球體組成的粗糙面的接觸力學(xué)特性，并分析了分子尺度下的結(jié)果和傳統(tǒng)力學(xué)模型計(jì)算結(jié)果的差異.結(jié)果表明，單個(gè)微凸體接觸結(jié)果和分子動(dòng)力學(xué)--格林函數(shù)法模擬所得非常接近，誤差在5%以?xún)?nèi).數(shù)值模擬發(fā)現(xiàn)，在微凸體高度符合高斯分布的情況下，接觸面積和接觸力成線性關(guān)系；在相同接觸面積下，微凸體模型得出的接觸力偏高，是上限值.微凸體模型沒(méi)有考慮微凸體間的相互影響，實(shí)際是高估了彈性體的剛度；實(shí)際接觸過(guò)程中微凸體相互影響，微凸體對(duì)臨域形變影響尤其大，使接觸區(qū)域更加離散.GFMD模型可以準(zhǔn)確計(jì)算數(shù)十億量級(jí)別分子、原子接觸過(guò)程中真實(shí)接觸面積及分布，為后續(xù)摩擦、滑移等分析提供可靠的參考.

微凸體模型,粗糙面接觸,分子動(dòng)力學(xué),格林函數(shù)法,接觸力學(xué)

引言

接觸摩擦現(xiàn)象無(wú)時(shí)不在，無(wú)處不有，其理論廣泛應(yīng)用于機(jī)械控制、樁基擠壓、熱傳導(dǎo)、潤(rùn)滑與密封、光干涉等領(lǐng)域.據(jù)統(tǒng)計(jì)，世界上能源的1/3～1/2消耗于摩擦與磨損，約80%的機(jī)器零件失效是由摩擦與磨損引起的[1]；經(jīng)調(diào)查，合理運(yùn)用摩擦學(xué)原理可以節(jié)約1%～1.4%的國(guó)民生產(chǎn)總值[2].表面間的接觸是摩擦產(chǎn)生的先決條件，其接觸面積、應(yīng)力分布是研究摩擦的基礎(chǔ)，因此一直是摩擦學(xué)關(guān)注的課題.沒(méi)有完全光滑的表面，所有表面都是粗糙面，因此表面接觸實(shí)際是粗糙面間的接觸[34].表面粗糙度對(duì)黏附力、接觸面積、摩擦性能等均有顯著的影響[58].粗糙面可以看成由無(wú)數(shù)微凸體組成，其內(nèi)部接觸摩擦機(jī)理遠(yuǎn)比光滑表面復(fù)雜.微凸體的形狀、接觸方向、微凸體間的相互作用都是未知的，而且接觸、摩擦過(guò)程可能產(chǎn)生彈塑性變形、應(yīng)力集中、斷裂等問(wèn)題，因此表面接觸、摩擦成為力學(xué)領(lǐng)域中極具挑戰(zhàn)性的難題[910].

近半世紀(jì)以來(lái)，國(guó)內(nèi)外眾多學(xué)者對(duì)表面接觸理論進(jìn)行了研究，提出了各種理論模型.這些模型大致可以分為4類(lèi)：(1)微凸體模型(aspeirty model)；(2)半解析模型；(3)降維模型；(4)數(shù)值計(jì)算模型.粗糙面接觸最有影響力的模型是微凸體模型，又稱(chēng)統(tǒng)計(jì)學(xué)模型[1112].微凸體模型假設(shè)粗糙面由許多半徑相同但高度不同的球狀微凸體組成，微凸體高度服從高斯分布.微凸體模型以微凸體接觸力學(xué)為基礎(chǔ)，對(duì)壓力和剪力進(jìn)行統(tǒng)計(jì)意義上的累加.之后，把隨機(jī)過(guò)程理論應(yīng)用于描述表面形貌[1315]，即對(duì)微凸體的密度、高度和曲率分布及其相關(guān)性進(jìn)行理論推導(dǎo)，其結(jié)果被后續(xù)研究者廣泛采用.微凸體模型被很多研究者不斷改進(jìn)，如將球狀微凸體替換成橢球狀微凸體的Bush,Gibson和Thomas的模型[16]，在微凸體模型基礎(chǔ)上考慮了微凸體曲率半徑等分形特征發(fā)展起來(lái)的分形模型[1718]，考慮了微凸體接觸對(duì)接觸方向的M isra-Huang模型[3].

半解析模型中比較著名的是德國(guó)學(xué)者Persson等提出的Persson模型[1920]，其認(rèn)為在彈性范圍內(nèi)，接觸體的彈性能與外力成正比，然后根據(jù)能量微分方程，得出了接觸體的豎向和水平向剛度.降維模型由學(xué)者Popov等提出[21]，他們認(rèn)為3D模型可簡(jiǎn)化成1D彈簧模型，然后在1D邊界上劃分單元進(jìn)行數(shù)值計(jì)算，可以快速對(duì)粗糙度、移動(dòng)速度、壓力等參數(shù)進(jìn)行數(shù)值計(jì)算，得到和3D模型相似的結(jié)果.其中，Persson和Popov兩個(gè)研究團(tuán)隊(duì)曾多次公開(kāi)發(fā)表學(xué)術(shù)論文質(zhì)疑對(duì)方模型的合理性[22].

隨著計(jì)算機(jī)計(jì)算能力的提升，近年來(lái)數(shù)值計(jì)算方法也成為研究熱點(diǎn).早期的數(shù)值計(jì)算方法主要是有限元法[23]、邊界元法[24].近年來(lái)，為更深入地探索接觸理論，分子動(dòng)力學(xué)--格林函數(shù)法(GFMD)在接觸摩擦中得以應(yīng)用[2527].分子動(dòng)力學(xué)從原子、分子力場(chǎng)出發(fā)分析接觸、摩擦機(jī)理，獲得了超潤(rùn)滑、多尺度效應(yīng)等用連續(xù)介質(zhì)力學(xué)理論難以模擬的力學(xué)行為[8,26,28].

本文以分子動(dòng)力學(xué)--格林函數(shù)法為工具，驗(yàn)證了其作為大規(guī)模分子、原子級(jí)計(jì)算方法的高效及準(zhǔn)確性，同時(shí)探討了由微球體組成的粗糙面的豎向接觸力學(xué)特性，分析了分子尺度下的結(jié)果和傳統(tǒng)力學(xué)模型計(jì)算結(jié)果的差異，并研究了造成這種差異的主要原因.

1 理論部分

1.1 分子動(dòng)力學(xué)--格林函數(shù)法

從分子、原子力場(chǎng)出發(fā)模擬表面接觸行為，有助于從根源上發(fā)現(xiàn)問(wèn)題，近年來(lái)得到廣泛認(rèn)可.分子動(dòng)力學(xué)與格林函數(shù)法[25,27,29]的大致思路是：在接觸層采用分子動(dòng)力學(xué)模擬表面力學(xué)效應(yīng)；接觸層以外，采用格林函數(shù)法模擬其彈性響應(yīng)，如圖1所示.該方法的優(yōu)點(diǎn)是既沒(méi)有忽略表面接觸層復(fù)雜的、分子級(jí)別的力學(xué)效應(yīng)，同時(shí)又利用格林函數(shù)的降維特性極大地提高了計(jì)算效率，從而使模擬大規(guī)模的表面接觸行為成為可能.應(yīng)注意到，兩個(gè)粗糙面的接觸可以簡(jiǎn)化成一個(gè)剛性粗糙面和一個(gè)完全光滑的彈性面間的接觸[20].

圖1 分子動(dòng)力學(xué)--格林函數(shù)法示意圖Fig.1 Schematic diagram ofmoleculardynamics-Green’s functionmethod

表面層原子間的勢(shì)能與他們之間的距離相關(guān)，比較通用的非鍵結(jié)勢(shì)能形式是Lennard-Jones(LJ)模型[30]，其數(shù)學(xué)表達(dá)式為

其中，r為原子間的距離，ε,σ為勢(shì)能參數(shù)，計(jì)算中分別取其為勢(shì)能、距離的基本單位.若無(wú)特別說(shuō)明，本文的力、距離的單位分別為ε/σ,σ.由式(1)可知當(dāng)r=21/6σ≈1.12σ時(shí)，原子間相互作用力為0.當(dāng)r＞1.12σ時(shí)，原子間作用力為吸引力；當(dāng)r＜1.12σ時(shí)，原子間作用力為排斥力.因此，文中表面原子的最初間距取為1.12σ.此外，本文不考慮黏附作用，LJ勢(shì)能的截?cái)嗑嚯x(cuto ff distance)亦取為1.12σ.若僅考慮分子間的非鍵結(jié)勢(shì)能，則其總勢(shì)能函數(shù)為

式中，n為原子個(gè)數(shù).由于各原子位置的改變，引起下層原子受力，如圖1(c)所示，其表達(dá)式為

下層原子在受力作用下將產(chǎn)生彈性變形，利用格林函數(shù)法[31]，可以寫(xiě)出其表達(dá)式

其中，u?為彈性力學(xué)基本解，P為原子的位置，Ω為積分邊界，即接觸面.基本解的表達(dá)式為

同樣，基本解對(duì)應(yīng)的彈性力t?ij表達(dá)式為

隨著原子位置的不斷迭代，原子力(分子力)和彈性力將達(dá)到平衡，即滿足收斂準(zhǔn)則.

1.2 微凸體模型

微凸體模型是Greenwood等提出的模型[12]，該模型假設(shè)粗糙面由無(wú)數(shù)半徑相同但高度不同的球體組成，單個(gè)球體的接觸應(yīng)力滿足Hertz解[32]，并假定球體間不存在相互影響.兩個(gè)粗糙面的接觸可以簡(jiǎn)化成一個(gè)剛性粗糙面和一個(gè)完全光滑的彈性面間的接觸[20]，此時(shí)僅需要把彈性體的彈性模量換算成等效彈性模量，其表達(dá)式為

式中，E1,E2分別是上、下兩彈性體的彈性模量，ν1,ν2分別為上、下兩彈性體的泊松比.據(jù)Hertz理論，有

式中，f為微球體的接觸力，a為圓形接觸面積半徑，R為球體的半徑，d為球體壓入深度.由于整個(gè)表面由無(wú)數(shù)微球體組成，根據(jù)中心極限定理，表面球體高度必然滿足高斯分布，假定其概率密度函數(shù)為p(z)，則表面的接觸力F的表達(dá)式為

式中，N表示微凸體個(gè)數(shù)，z為微凸體高度，D為兩表面剛接觸時(shí)光滑面到粗糙面參考面間的距離[12].微凸體模型邏輯簡(jiǎn)單，計(jì)算方便，特別是引入隨機(jī)過(guò)程后，其各個(gè)統(tǒng)計(jì)參數(shù)的計(jì)算在數(shù)學(xué)上是精確的，因此獲得廣泛認(rèn)可.

2 本文模型的建立

本文的原子模型采用面心立方，由于其對(duì)稱(chēng)性，僅取晶胞的第一層原子構(gòu)成表面形貌.表面形貌由大量半徑相同但凸出高度不同的球體組成，如圖2(a)所示；同時(shí)球體高度符合高斯分布，如圖2(b)所示.采用的計(jì)算程序是美國(guó)Sandia國(guó)家實(shí)驗(yàn)室開(kāi)發(fā)的開(kāi)源分子動(dòng)力學(xué)程序Lammps[33],把格林函數(shù)當(dāng)成了分子的一種力場(chǎng)整合進(jìn)該程序.計(jì)算前先生成剛性粗糙面和彈性光滑面，然后導(dǎo)入主程序，利用位移加載，使粗糙接觸面逐步向彈性光滑面靠近.每次位移加載步長(zhǎng)為0.01σ，然后進(jìn)行迭代分析，直至彈性層分子(原子)受力平衡.在四周邊界處，原子可能溢出，因此設(shè)置了周期性邊界；在每次進(jìn)行迭代分析時(shí)，設(shè)置最大計(jì)算步為5×104步，收斂準(zhǔn)則為原子力的1-范數(shù)等于0.01ε/σ.由于勢(shì)能函數(shù)是長(zhǎng)程力，因此我們?cè)O(shè)置了截?cái)嗑嚯x為1.12σ，即當(dāng)他們的距離大于1.12σ時(shí)，兩原子受力為0，當(dāng)兩原子受力小于1.12σ時(shí)，處于受力(接觸)狀態(tài).為簡(jiǎn)單起見(jiàn)，文中的算例中，彈性模量E1=∞，E2=3ε/σ3，泊松比v1=0.5，v2=0.5；球體半徑R=50σ，高斯分布方差為5σ.

圖2 表面形貌及概率分布Fig.2 Surfacemorphology and probability distribution

3 結(jié)果與討論

3.1 單個(gè)微凸體壓入數(shù)值模擬實(shí)驗(yàn)

利用球體方程生成如圖3(a)所示的剛性球體，球體半徑為50σ，壓入512σ×512σ×1024σ的彈性體中.這里需要指出的是，在彈性體中，僅在第1層建立原子，第1層以下利用格林函數(shù)法進(jìn)行模擬.這樣的優(yōu)點(diǎn)是只需要布置512×512=262144個(gè)原子，而常規(guī)的分子動(dòng)力學(xué)計(jì)算需要布置512×512×1024=368435456個(gè)原子，這種計(jì)算量即使是超級(jí)計(jì)算機(jī)都難以完成.每次加載完之后，輸出原子的位置和力.接觸面積與力關(guān)系的數(shù)值計(jì)算結(jié)果和Hertz理論結(jié)果如圖3(b)所示，兩者最大誤差為5%.造成兩者差距的主要原因在于Hertz理論是基于半球體和無(wú)限半空間彈性體的接觸，而本次數(shù)值模擬的球體半徑比較小，因此兩者誤差是合理的.本算例利用4核普通計(jì)算機(jī)并行計(jì)算，30min左右完成計(jì)算，可見(jiàn)GFMD計(jì)算效率很高.

圖3 單個(gè)微凸體接觸行為Fig.3 Singleasperity contactbehavior

3.2 粗糙面接觸數(shù)值模擬實(shí)驗(yàn)

這個(gè)實(shí)驗(yàn)是將一個(gè)1024σ×1024σ(1048576個(gè)原子組成)的剛性粗糙面壓入1024σ×1024σ×2048σ的彈性體中.同樣，首先生成符合高斯分布的球體高度數(shù)據(jù)1024個(gè)(32×32)，然后將1024σ×1024σ的剛性平面分成1024個(gè)網(wǎng)格，每個(gè)網(wǎng)格大小亦為32σ×32σ，之后在每個(gè)網(wǎng)格隨機(jī)放入之前建立的球體，如圖2(a)和圖4(a)所示.雖然只有1024個(gè)微凸體，但接觸面積--力仍然成線性關(guān)系，如圖4(b)所示，這和大量數(shù)值模擬及試驗(yàn)結(jié)果吻合[11,17].可見(jiàn)，只要微凸體高度符合高斯分布，接觸面積與力基本成線性比例[3436].在相同接觸面積下，微凸體模型接觸力計(jì)算值偏大，最大的差值為20%.其偏大的原因?yàn)椋何⑼贵w模型忽略了微凸體間的相互影響，這樣實(shí)際是高估了彈性體的剛度.在接觸面積較小的時(shí)候，微凸體間距離較大，他們的相互影響可以忽略不計(jì)，因此，低接觸面積下，兩者的接觸力幾乎一致.

圖4 粗糙表面接觸行為Fig.4 Rough surface contactbehavior

3.3 微凸體空間分布對(duì)結(jié)果的影響

為研究微凸體空間分布對(duì)接觸面積與力關(guān)系曲線的影響，研究團(tuán)隊(duì)做了一個(gè)對(duì)比試驗(yàn).分別在512σ×512σ和1024σ×1024σ大小的表面布置16×16=256個(gè)球體，兩者球體高度一致，但是間距不一樣，后者的間距是前者的兩倍.接觸面積與力關(guān)系如圖5所示，在相同接觸面積下，微凸體間距大的模型則接觸力較小，間距小的模型則接觸力較大.在接觸面積為7%左右時(shí)，512σ×512σ算例比1024σ×1024σ算例的接觸力大8%左右；在接觸面積較小時(shí),兩個(gè)算例結(jié)果幾乎一致，因?yàn)榇藭r(shí)微凸體之間的影響較小.圖6分別顯示了微凸體模型和GFMD模型在5%接觸面積下計(jì)算的接觸點(diǎn)的空間分布情況(模型參數(shù)見(jiàn)3.2算例).在圖6(a)中，由于假設(shè)微凸體間完全沒(méi)有相互影響，其接觸面積可以直接從剛性粗糙面截?。粓D6(b)是考慮微凸體間的

圖5 不同微凸體間距下粗糙表面接觸行為Fig.5 Rough contactbehaviorw ith di ff erentasperity distance

圖6 在5%接觸面積下的接觸點(diǎn)分布情況Fig.6 Contactspotsdistribution at5%contactarea

影響，計(jì)算結(jié)果來(lái)自于GFMD.顯而易見(jiàn)，考慮微凸體間相互影響時(shí)的接觸區(qū)域更加離散，接觸點(diǎn)更多.圖6(a)沒(méi)有考慮微凸體間的相互影響，接觸點(diǎn)容易連成一片，其接觸點(diǎn)僅為97個(gè)；而圖6(b)考慮微凸體間的相互影響，接觸點(diǎn)更加離散，其接觸點(diǎn)為152個(gè).可見(jiàn)，微凸體模型沒(méi)有考慮微凸體的相互影響，其實(shí)際接觸區(qū)域并不準(zhǔn)確.準(zhǔn)確的接觸面積分布、壓力大小等對(duì)后續(xù)的摩擦、滑移等分析非常重要[8,26]，因此GFMD是接觸、摩擦分析的可靠計(jì)算方法.

4 結(jié)論

本文利用分子動(dòng)力學(xué)--格林函數(shù)法對(duì)大規(guī)模粗糙表面接觸進(jìn)行數(shù)值分析，得出以下結(jié)論：

(1)分子動(dòng)力學(xué)--格林函數(shù)法計(jì)算效率高，僅在接觸層進(jìn)行分子、原子間的接觸，且可利用并行計(jì)算加速；分子層以下采用格林函數(shù)法，可以起到降維效果.該方法從分子、原子力場(chǎng)出發(fā)，對(duì)數(shù)十億級(jí)的大規(guī)模原子尺度接觸分析效果良好.

(2)在微凸體高度符合高斯分布的情況下，接觸面積和接觸力成線性關(guān)系；在相同接觸面積下，微凸體模型得出的接觸力偏高，是上限值.

(3)微凸體模型沒(méi)有考慮微凸體間的相互影響，實(shí)際是高估了彈性體的剛度；微凸體之間的距離決定了他們之間的相互影響程度，距離大時(shí)影響小，距離小時(shí)影響較大.實(shí)際接觸過(guò)程中微凸體相互影響，微凸體對(duì)臨域形變影響尤其大，使接觸區(qū)域更加離散.

(4)微凸體模型可以快速地預(yù)測(cè)粗糙面豎向接觸中的基本接觸特性，對(duì)于精確度不高的分析仍是一種簡(jiǎn)明的計(jì)算方法.然而，對(duì)于后續(xù)的摩擦、滑移等分析，能計(jì)算準(zhǔn)確接觸面積、壓力分布等的GFMD模型，是更可靠的計(jì)算方法.

1楊艷峰,鄭堅(jiān),狄長(zhǎng)春等.基于微觀分析的火炮擋彈裝置磨損失效機(jī)理研究.摩擦學(xué)學(xué)報(bào)，2014，34(3)：304-310(Yang Yanfeng,Zheng Jian,DiChangchun,etal.Wearmechanism of gun cartridge stop device based onm icroanalysis.Tribology,2014,34(3)：304-310(in Chinese))

2 Bhushan B,Ko PL.Introduction to tribology.Applied Mechanics Reviews,2013,56(1)：136

3 Huang HP,M isra A.M icro-Macro-Shear-Displacementbehavior of contacting rough solids.Tribology Letters,2013,51：431-436

4 M isra A,Huang SP.M icromechanical stress-displacementmodel for rough interfaces：E ff ectofasperity contactorientation on closure and shear behavior.International JournalofSolids and Structures,2012,49(1)：1111-1120

5 Peng ZL,Wang C,Chen SH.Them icrostructure morphology on ant footpadsand itse ff ecton antadhesion.Acta Mechanica,2016,227(7)：2025-2037

6 Peng ZL,Chen SH.E ff ectsof surface roughnessand fil thickness on theadhesion ofabioinspired nanofilm PhysicalReview E,2011,83(5)：051915

7 Ben-David O,Rubinstein SM,Fineberg J.Slip-stick and the evolution of frictionalstrength.Nature,2010,463(7277)：76-79

8 Mo YF,Turner KT,Szlufarska I.Friction lawsat thenanoscale.Nature,2009,457：1116-1119

9瓦倫丁L,波波夫.接觸力學(xué)與摩擦學(xué)的原理及其應(yīng)用.北京：清華大學(xué)出版社,2011(Valentin L.Popov.ContactMechanics and Friction：Physical Principles and Application.Beijing：Tsinghua University Press,2011(in Chinese))

10黃平,孟永鋼,徐華.摩擦學(xué)教程.北京：高等教育出版社,2007(Huang Ping,Meng Yonggang,Xu Hua.Tribology Course.Beijing：Higher Education Press,2007(in Chinese))

11 Greenwood JA,Tripp JH.The contact of two nom inally fla rough surfaces//Proceedings of the Institution of Mechanical Engineers,1970,185：625-634

12 Greenwood JA,Williamson JB.Contact of nom inally fla surfaces//Proceedings of the Royal Society of London Series AMathematicaland PhysicalSciences,1966,295：300-319

13 Nayak PR.Random processmodelof rough surfaces in plastic contact.Wear,1973,26(3)：305-333

14 Nayak PR.Random processmodel of rough surfaces.Journal of Lubrication Technology,1971,93(3)：398-407

15 Longuet-Higgins MS.The statistical analysis of a random,moving surface.Philosophical Transactions ofthe Royal Society ofLondon SeriesAMathematicaland PhysicalSciences,1957,249(966)：321-387

16 Bush AW,Gibson RD.Elastic contact of a rough surface.Wear,1975,35(1)：87-111

17 Xie HP,Wang JN,XieWH.Fractal e ff ectsof surface roughnesson themechanicalbehaviorof rock joints.Chaos,Solitons&Fractals,1997,8(2)：221-252

18 Majumdar A,Bhushan B.Fractalmodel of elastic-plastic contact between rough surfaces.Journal of Tribology-Transactions of the ASME,1991,113(1)：1-11

19 Campana C,Persson BNJ,M¨uer MH.Transverse and normal interfacial sti ff ness of solidswith random ly rough surfaces.Journal of Physics:Condensed Matter,2011,23(8)：085001

20 Persson BNJ.Relation between interfacial separation and load：A generaltheory of contactmechanics.PhysicalReview Letters,2007,99(12)：125502

21 Li Q,Popov M,Dimaki A,et al.Friction between a viscoelastic body and a rigid surfacewithrandom self-a ffi ne roughness.PhysicalReview Letters,2013,111(3)：034301

22 Persson BN.Contactmechanics for random ly rough surfaces：on the validity of themethod of reduction of dimensionality.Tribology Letters,2015,58(11)：1-4

23王勖成.有限單元法.北京：清華大學(xué)出版社,2003(Wang Xucheng.Finite Element Method.Beijing：Tsinghua University Press,2003(in Chinese))

24姚振漢,王海濤.邊界元法.北京：高等教育出版社，2010(Yao Zhenhan,Wang Haitao.Boundary Element Methods.Beijing：Higher Education Press,2010(in Chinese))

25 Pastewka L,Sharp TA,Robbins MO.Seam less elastic boundaries for atom istic calculations.PhysicalReview B,2012,86：075459

26 Luan B,RobbinsMO.Thebreakdown of continuum models formechanical contacts.Nature,2005,435：929-932

27 Campa?n′a C.Using Green’s function molecular dynam ics to rationalize the success of asperity models when describing the contact between self-a ffi ne surfaces.PhysicalReview E,2008,78：026110 28 LiSZ,LiQY,Carpick RW,etal.Theevolving quality of fractional contactw ith graphene.Nature,2016,539：541-545

29 Luan BQ,RobbinsMO.Hybrid atom istic/continuum study of contactand friction between rough solids.Tribology Letters,2009,36：1-16

30陳正隆,徐為人,湯立達(dá).分子模擬的理論與實(shí)踐.北京：化學(xué)工業(yè)出版社,2007(Chen Zhenglong,XuWeiren,Tang Lida.The Theory and Practice of Molecular Simulation.Beijing：Chem ical Industry Press,2007(in Chinese))

31王元淳.邊界元法基礎(chǔ).上海：上海交通大學(xué)出版社，1988(Wang Yuanchun.Foundation of Boundary Element Method.Shanghai：Shanghai Jiao Tong University Press,1988(in Chinese))

32 Hertz H.On the contactof elastic solids.Journal fur de reine und angewandteMathematik,1881,92：156-171

33 Plimpton S,Crozier P,Thompson A.LAMMPS-large-scale atom ic/molecular massively parallel simulator.Sandia National Laboratories,2007,18

34 Zavarise G,PaggiM.Reliability ofm icromechanical contactmodels：a still open issue//W riggers P,Laursen TA,eds.CISM International Centre for Mechanical Sciences,Springer Vienna,CISM,Udine,2007,498：39-82

35 Hyun S,PeiL,MolinariJ-F,etal.Finite-elementanalysisof contact between elastic self-a ffi ne surfaces.Physical Review E,2004,70：026117

36 Buzio R,Boragno C,Biscarini F,etal.The contactmechanics of fractalsurfaces.NatureMaterials,2003,2(4)：233-236

NUMERICAL ANALYSISOFASPERITY CONTACTMODEL BASED ONMOLECULAR DYNAM ICS-GREEN’SFUNCTIONMETHOD1)

Huang Shiping?,?Wu Jie?Hu Junliang?Zheng Hengbin??WangWeifeng?,2)?

(SchoolofCivilEngineering and Transportation，South China University ofTechnology，Guangzhou 510640，China)?(State Key Laboratory ofTribology atTsinghua University，Beijing 100084，China)??(College ofWaterConservancy and CivilEngineering，South China AgriculturalUniversity，Guangzhou 510642，China)

Rough contact is a prerequisite for surface friction.The rough contact behaviour such as the contact area,the pressure distribution and spatial distributionshasbeen one of the core issues in contactmechanicsand tribology.In this paper,themolecular dynamics-Green’s functionmethod(GFMD)isused to simulate the contactmechanism of the rough surface,where the asperitymodel is used for the rough surface,i.e.,the surface is composed of numerous spherical asperities.Starting w ith the atom ic ormolecular force file to consider the rough contactbehaviour,themolecular dynamics-Green’s function method is able to capture themechanisms such as super-lubrication and multi-scale e ff ectbehaviour,which are not found in traditional continuum mechanics.Themolecular dynam ics-Green’s functionmethod demonstrates its high e ffi ciency in large scalemolecular dynam ics simulations and is able to simulate the system composed of billions of atoms.The results of single asperity contactbased on Hertz contact theory are very close to those simulated by themolecular dynam ics-Green’s functionmethod,and the di ff erence is less than 5%.It is found by numerical simulation that the contactarea is linearly related to the contact force if the asperity heights follow the Gaussian distribution,and the contact forceobtained by theasperitymodel is theupper lim itgiven the same contactarea.A lthough Asperitymodel is fast,itoverestimates the sti ff nessof theelastomer due to the neglection of the interaction between the asperities.In real contact process,asperities have considerable e ff ects on each other,especially on the deformation of the adjacentarea,whichmakes the contact spotsmore discrete.The information of the real contactarea and its spatial distributions,isof importance for the follow ing simulation on surface friction.

asperitymodel,rough contact,molecular dynam ics,Green’s function,contactmechanics

O343.3

A

10.6052/0459-1879-17-084

2017-03-04收稿，2017-04-10錄用,2017-04-10網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金(11202080,11672108)，清華大學(xué)摩擦學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(SKLTKF15B05)和交通運(yùn)輸部建設(shè)科技項(xiàng)目基金(2014318363230)資助項(xiàng)目.

2)王衛(wèi)鋒,教授,主要研究方向：實(shí)驗(yàn)力學(xué)、結(jié)構(gòu)健康監(jiān)測(cè).E-mail：ctw fwang@scut.edu.cn

黃仕平,吳杰,胡俊亮,鄭恒斌,王衛(wèi)鋒.基于分子動(dòng)力學(xué)--格林函數(shù)法的微凸體接觸數(shù)值分析.力學(xué)學(xué)報(bào),2017,49(4)：961-967

Huang Shiping,Wu Jie,Hu Junliang,Zheng Hengbin,WangWeifeng.Numericalanalysisof asperity contactmodelbased onmolecular dynam ics-Green’s functionmethod.Chinese JournalofTheoreticaland Applied Mechanics,2017,49(4)：961-967