基于分子动力学-格林函数法的分形粗糙表面摩擦行为研究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随着面积的增大而增大,大致呈线性趋势.