基于分子动力学--格林函数法的微凸体接触数值分析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