基于分子動力學
--格林函數法的微凸體接觸數值分析1)

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

?(華南理工大學土木與交通學院，廣州510640)?(清華大學摩擦學國家重點實驗室，北京100084)??(華南農業大學水利與土木工程學院，廣州510642)

生物、工程及交叉力學

基于分子動力學
--格林函數法的微凸體接觸數值分析1)

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

?(華南理工大學土木與交通學院，廣州510640)?(清華大學摩擦學國家重點實驗室，北京100084)??(華南農業大學水利與土木工程學院，廣州510642)

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

微凸體模型,粗糙面接觸,分子動力學,格林函數法,接觸力學

引言

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

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

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

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

本文以分子動力學--格林函數法為工具，驗證了其作為大規模分子、原子級計算方法的高效及準確性，同時探討了由微球體組成的粗糙面的豎向接觸力學特性，分析了分子尺度下的結果和傳統力學模型計算結果的差異，并研究了造成這種差異的主要原因.

1 理論部分

1.1 分子動力學--格林函數法

從分子、原子力場出發模擬表面接觸行為，有助于從根源上發現問題，近年來得到廣泛認可.分子動力學與格林函數法[25,27,29]的大致思路是：在接觸層采用分子動力學模擬表面力學效應；接觸層以外，采用格林函數法模擬其彈性響應，如圖1所示.該方法的優點是既沒有忽略表面接觸層復雜的、分子級別的力學效應，同時又利用格林函數的降維特性極大地提高了計算效率，從而使模擬大規模的表面接觸行為成為可能.應注意到，兩個粗糙面的接觸可以簡化成一個剛性粗糙面和一個完全光滑的彈性面間的接觸[20].

圖1 分子動力學--格林函數法示意圖Fig.1 Schematic diagram ofmoleculardynamics-Green’s functionmethod

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

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

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

下層原子在受力作用下將產生彈性變形，利用格林函數法[31]，可以寫出其表達式

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

同樣，基本解對應的彈性力t?ij表達式為

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

1.2 微凸體模型

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

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

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

式中，N表示微凸體個數，z為微凸體高度，D為兩表面剛接觸時光滑面到粗糙面參考面間的距離[12].微凸體模型邏輯簡單，計算方便，特別是引入隨機過程后，其各個統計參數的計算在數學上是精確的，因此獲得廣泛認可.

2 本文模型的建立

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

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

3 結果與討論

3.1 單個微凸體壓入數值模擬實驗

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

圖3 單個微凸體接觸行為Fig.3 Singleasperity contactbehavior

3.2 粗糙面接觸數值模擬實驗

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

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

3.3 微凸體空間分布對結果的影響

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

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

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

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

4 結論

本文利用分子動力學--格林函數法對大規模粗糙表面接觸進行數值分析，得出以下結論：

(1)分子動力學--格林函數法計算效率高，僅在接觸層進行分子、原子間的接觸，且可利用并行計算加速；分子層以下采用格林函數法，可以起到降維效果.該方法從分子、原子力場出發，對數十億級的大規模原子尺度接觸分析效果良好.

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

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

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

1楊艷峰,鄭堅,狄長春等.基于微觀分析的火炮擋彈裝置磨損失效機理研究.摩擦學學報，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,波波夫.接觸力學與摩擦學的原理及其應用.北京：清華大學出版社,2011(Valentin L.Popov.ContactMechanics and Friction：Physical Principles and Application.Beijing：Tsinghua University Press,2011(in Chinese))

10黃平,孟永鋼,徐華.摩擦學教程.北京：高等教育出版社,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王勖成.有限單元法.北京：清華大學出版社,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陳正隆,徐為人,湯立達.分子模擬的理論與實踐.北京：化學工業出版社,2007(Chen Zhenglong,XuWeiren,Tang Lida.The Theory and Practice of Molecular Simulation.Beijing：Chem ical Industry Press,2007(in Chinese))

31王元淳.邊界元法基礎.上海：上海交通大學出版社，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網絡版發表.

1)國家自然科學基金(11202080,11672108)，清華大學摩擦學國家重點實驗室開放基金(SKLTKF15B05)和交通運輸部建設科技項目基金(2014318363230)資助項目.

2)王衛鋒,教授,主要研究方向：實驗力學、結構健康監測.E-mail：ctw fwang@scut.edu.cn

黃仕平,吳杰,胡俊亮,鄭恒斌,王衛鋒.基于分子動力學--格林函數法的微凸體接觸數值分析.力學學報,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