含石墨烯润滑油润滑机制的分子动力学模拟*

2019-10-23

润滑与密封 2019年10期

(1.沈阳建筑大学分析与检测技术研究中心辽宁沈阳 110168；2.沈阳建筑大学机械工程学院辽宁沈阳 110168；3.“高档石材数控加工装备与技术”国家地方联合工程实验室辽宁沈阳 110168；4.沈阳建筑大学理学院辽宁沈阳 110168)

氮化硅(Si3N4)陶瓷球轴承在高速回转机械中应用较为广泛。陶瓷球轴承的滚动体材料为Si3N4，内外圈材料为GCr15，滚动体与内外圈之间分别形成Si3N4-GCr15摩擦副。高速轴承滚道的磨损及摩擦副间的摩擦生热一直是影响轴承寿命的主要因素。石墨烯特殊的二维纳米层状结构使其具有很好的润滑特性。蒲吉斌等[1]尝试在润滑油中加入石墨烯，发现石墨烯不断地覆盖在摩擦零件表面，可有效地降低轴承摩擦因数，进而提高主轴的极限转速并降低温升。液体润滑剂中石墨烯添加剂的作用，大多通过石墨烯沉积吸附形成固体润滑膜、石墨烯片层嵌于接触粗糙峰间产生润滑作用，以及对磨粒的包裹形成滚动润滑体等宏观机制加以解析[2-3]，但从微观层面进一步解释石墨烯吸附产生的机制更具理论意义。对于纳米材料的润滑，在原子、分子尺度研究相对运动界面上的摩擦磨损与润滑行为，揭示微观摩擦磨损机制是近年来学者研究的重点[4-7]。本文作者鉴于石墨烯的纳米尺度，力图从微观机制层面揭示添加石墨烯前后，润滑剂与上下壁面间吸附作用机制。

近些年，将石墨烯作为润滑油添加剂的研究已经取得了很多显著的成果[8-16]。例如，SENATORE等[17]通过超声技术直接将石墨烯纳米片分散到矿物润滑油中，考察不同赫兹接触应力、温度和速率下，氧化石墨烯纳米片调配的矿物润滑油从边界润滑、混合润滑到弹流润滑域的摩擦磨损变化。现有的研究成果中，基于摩擦试验的石墨烯润滑油性能的研究很多，但关于石墨烯作为润滑油添加剂的润滑机制分析较少。

分子动力学模拟(Molecular Dynamics Simulation, MDS)是通过对微观离散粒子的模拟与分析，利用统计力学的方法阐释宏观连续物质的性质的技术。借助MDS技术可以用较小的实验成本和时间成本，探究现象的微观机制[18]。分子动力学模拟的方法可以帮助从分子水平上理解润滑剂与聚合物表面之间的相互作用。国内外有许多应用分子动力学研究润滑方面成果。邹鲲等人[19]应用宏观和微观特性的MDS方法，研究了超薄膜的微观结构与摩擦学特性间的关系。兰惠清和徐藏[20]采用分子动力学模拟的方法分别研究干摩擦和油润滑2种情况下不同硅含量的掺硅金刚石(Si-DLC)薄膜的摩擦过程。李义雅等[21]采用分子动力学的方法以十四烷基环为模型化合物从烃分子与Fe表面的吸附能以及烃分子之间的作用能2个方面，研究了矿物基础油分子在Fe表面的润滑作用，并比较了烃分子结构及温度对润滑作用的影响。白敏丽等[22]采用分子动力学的方法建立薄膜润滑的计算体系，以正庚烷作为润滑剂研究温度和压力对润滑油膜结构组成、滑移现象和摩擦特性的影响。曾凡林和孙毅[23]建立了作为润滑剂的正十六烷长链分子间的作用力模型，运用分子动力学模拟的方法得到了润滑剂在6种不同剪切速率下的应力、速度和温度分布，研究摩擦过程中特殊的物理现象。王慧等人[24]用分子动力学模拟研究了纳米级润滑薄膜的固液相变和界面滑移现象以及固体接触和黏着的微观机制。TANAKA等[25]采用非平衡分子动力学模拟(NEMD)建立了由2个实心壁和氟碳聚合物润滑剂组成的模型，调查薄壁润滑膜限制在墙之间的动态行为。JABBARZADEH等[26]针对C30的异构体6种不同的分子，研究了分子对润滑膜的流变性和行为的影响。GAO等[27]采用经典的分子动力学模拟和自由能计算方法对超薄油膜润滑进行计算。SIVEBAEK等[28]研究分子长度为20～1 400个碳原子的封闭烃“聚合物”固体分子的动力学摩擦计算，针对聚合物在硬衬底上滑动和聚合物在聚合物上滑动2种假设，分析不同分子长度摩擦剪切应力与剪切速度之间的关系。TARTAGLINO等[29]运用分子动力学的方法，以正丁烷和异丁烷作为润滑剂研究粗糙表面的润滑效果。上述研究表明，用分子动力学的方法研究石墨烯作为润滑油添加剂的润滑机制是可行的。

本文作者建立分子模型，研究石墨烯作为添加剂的正十六烷烃润滑油在Si3N4-GCr15摩擦副的润滑机制。模型中采用双层石墨烯(Double-layer-graphene)，即由两层以苯环结构(即六角形蜂巢结构)周期性紧密堆积的碳原子以不同堆垛方式(包括AB堆垛、AA堆垛、AA堆垛等)堆垛构成的一种二维碳材料。同时分析石墨烯含量、压力和速度对润滑油在Si3N4-GCr15摩擦副中的润滑行为的影响规律。

1 分子动力学模型及计算过程

1.1 物理模型

运用Materials Visualizer模块建立Si3N4、GCr15、C16H34和石墨烯的分子模型。首先用Amorphous Cell模块将C16H34和石墨烯混合在一起，石墨烯的质量分数根据需要设定。用Build Layers功能将Si3N4、GCr15和混合后的石墨烯润滑剂建立系统的分子模型。图1(a)是系统模型的结构图，GCr15与正十六烷烃的接触面分别是上壁面，Si3N4与正十六烷烃的接触面分别是下壁面。2个壁面的位置分别是5.645 nm和1.071 nm处。模型上下壁面原子做反向运动，即为剪切运动。对于限制剪切运动，还需在GCr15上方添加2 nm的真空层。

图1(b)所示是系统三维分子动力学模型，模型的长、宽和高分别是a、b和c。a=4.574 nm，b=4.574 nm，c=12.24 nm。

图1 系统模型的结构和三维分子动力学模型Fig 1 Model structure(a)and 3-D molecular dynamics model(b)

1.2 势能函数

分子间的作用很复杂，包括烷烃分子内部作用，润滑分子间作用，润滑分子与石墨烯的作用，润滑分子与壁面的界面作用，Si3N4和GCr15材料分子内部和分子间作用等。具体情况如图2所示。

文中在建模时，首先将石墨烯和十六烷烃作为整体，采用公式(1)计算其与Si3N4和GCr15之间的作用VLJ1和VLJ2，GCr15、Si3N4材料分子间作用为VLJ3和VLJ4。因为文中采用Universal力场，所以分子间作用只能采用LJ势。通常来说，石墨烯直接由碳原子构建，在分子动力学中采用Tersoff-Brenner势函数来描述其作用。但是在文中的模拟计算，由于石墨烯与正十六烷烃分子间具有较强的范德华力作用[30]，Tersoff-Brenner势却只考虑了碳碳化学键的作用,而忽略了范德华力作用,因此石墨烯碳原子和邻近正十六烷烃的相互作用采用LJ势函数表示[30-32]。在润滑体内部，考虑了正十六烷烃分子间作用、石墨烯分子间作用、正十六烷烃分子与石墨烯的作用VLJ5、VLJ6和VLJ7，同时也考虑了烷烃分子内部作用Eθ5和Eφ5及石墨烯分子内部的相互作用Eθ6和Eφ6，其中分子间作用采用公式(1)计算，分子内部原子的作用采用公式(2)和公式(3)计算。

公式(1)是经典的对势函数，L-J势能函数能够反映粒子的运动轨迹[33-34]。文中采用LJ/126势能函数，其表达式为

(1)

式中：VLJ为系统势能(范德华能)；ε为能量参数(r=21/6×σ的平衡位置时的能量最小值);σ为长度参数(零势能距离)；r为原子间的距离。

正十六烷烃和石墨烯内部原子之间的相互作用共价键的弯曲，采用van der Ploeg和Berendsen势[35]，表达式为

(2)

式中：Eθ为弯曲能；Kθ为系数；θ0为平衡角。

正十六烷烃的Kθ和θ0分别为62 500 K·rad2和114°，石墨烯的Kθ和θ0分别为0.9×10-18K·rad2和120°。

正十六烷烃和石墨烯内部原子之间的相互作用共价键的扭转，采用Ryckaert和Bellemans势[36]，表达式为

(3)

式中：Eφ为扭转能；φ为分子链上相邻4个“原子”确定的两平面之间的二面角；CK为系数。

正十六烷烃的C0=1 116 K，C1=1 462 K，C2=-1 578 K，C3=-368 K，C4=3 156 K，C5=-3 788 K。

模拟中，假设润滑油中靠近壁面的2个边界上的原子随着壁面做同向运动，且每一层所有原子沿着剪切方向运动速度的平均值为这一层的速度。原子的速度运用势函数计算：

(4)

式中：a为原子加速度；m为原子的质量。

(5)

(6)

式中：v为每一层的原子速度；N为每一层原子的个数。

1.3 剪切应力模型

基于维里定理，剪切应力为张量：

(7)

(8)

(9)

式中：rix、vix和fix分别表示第i个原子的位置、速度和力矢量的x方向的分量。

应力张量为三阶反对称张量，所以可以表示为

(10)

式中：SXX、SXY、SYX、SYY、SYZ、SZY、SZZ均为0。

1.4 分子动力学模型

采用Materials Studio7.0中的Forcite模块，力场为Universal，静电作用和范德华作用的非键截断均采用Atom based方法，截断距离为1.25 nm。压力和温度的控制函数分别为Berendsen和NHL(Nosé-Hoover-Langevin)方法。

模拟过程分为3部分，首先将建好的模型用Geometry Optimization进行结构优化，即采用Smart计算法找到模型的最低能量构象，以该构象作为模型进行模拟。其次，将优化后的模型进行动力学弛豫，使分子逐渐恢复到平衡态。弛豫的过程中考虑分子的热运动。势能模型选择Universal力场，同时选择恒温恒压(NPT)系综。模拟的初始温度设置为298 K。时间步长选择1 fs，总模拟时间为25 ps，模拟步数为25 000步，每5 000步输出一帧。同时设置相应的压力压缩模型。弛豫后分子动力学模型长、宽和高分别变为3.088、3.088和10.31 nm，上下壁面分别为5.67和5.271 nm。最后，进行非平衡分子动力学模拟，即剪切过程。应用Forcite模块中Confined shear功能进行剪切运动，剪切运动就是模型上下两层原子以速度v沿着x轴做反向运动，在模型中上下壁面上施加剪切速度，上层GCr15钢和下层Si3N4以相对速度向反方向运动。势能模型选择Universal力场，模拟的初始温度为298 K。设置时间步长选择0.1 fs，总模拟时间为100 ps，模拟步数为10万步，每5 000步输出一帧。考虑以上边界条件，进行分子动力学模拟并给出范德华力、类固膜厚度及剪切应力的数值分析结果。在整个润滑过程中，记录范德华能(表征吸附能力)、类固膜厚度(表征薄膜厚度)及x向的力学响应变化(表征摩擦力)。

图3(a)所示为动力学弛豫和限制剪切后的温度变化。图3(b)所示为压力102 MPa、速度141 m/s条件下，不同时刻限制剪切运动后的范德华能。

图3 弛豫及剪切后温度及能量变化Fig 3 Temperature and energy changes after relaxation and shearing(a)temperature change after relaxation and shearing；(b) van der Waals in different time

2 模拟结果讨论与分析

2.1 石墨烯含量对摩擦副润滑性能的影响

研究发现石墨烯质量分数对摩擦副的摩擦学性能有较大影响，适当的石墨烯含量可以降低摩擦因数[3,5]。设置石墨烯摩尔质量分数分别为0、2.5%、3.3%和4.2%，建立润滑油分子模型。模拟不同含量石墨烯添加到润滑油中，对润滑油薄膜润滑性能的影响。模拟时设置压力为102 MPa，剪切速度为141 m/s。

模拟中，要计算润滑剂与壁面(Si3N4+GCr15)之间范德华能。计算时，利用优化-弛豫-剪切后的模型，复制包含润滑体及上下壁面附近的原子，用Forcite模块中的energy做单点能计算；并用Forcite分析面板中View in study table选项，得到能量分项。文中建立的模型，上下分子较小，润滑剂与壁面(Si3N4+GCr15)之间范德华能与总的范德华能相差特别小。以压力102 MPa、剪切速度141 m/s条件下的范德华能为例，总的范德华能为162 374.002 kJ/mol，润滑剂与壁面(Si3N4+GCr15)之间的范德华能为163 362.162 kJ/mol，差值仅为988.16 kJ/mol，可以忽略不计。所以文中的范德华能以总的范德华能代替。用式(1)计算范德华能时，式中的ε和σ用Lorentz-Berthelot结合法则，表达式为

(11)

式中：εi为原子i的特征能量；εj为原子j的特征能量；σi为原子i的特征半径；σj为原子j的特征半径。

CH2、 CH3、石墨烯C-C 3种基团的L-J势参数如表1所示。表1中CH2、CH3数据来源于文献[37]，石墨烯C-C数据来源于文献[38]。

表1 3组基团的L-J势参数Table 1 L-J potential parameters of the three groups of groups

注：KB为玻尔兹曼常数，数值是1.381× 10-23J/K。

范德华能表征了存在于分子间的一种吸引力，范德华能越大，分子间的吸引力越大。图4表示在Si3N4-GCr15之间，不同含量石墨烯添加剂的润滑油所具有的范德华能随时间变化曲线。可以看出，Si3N4与GCr15之间的润滑油范德华能随时间变化不大。石墨烯含量会影响Si3N4与GCr15摩擦副之间范德华能的大小。其中石墨烯质量分数为2.5%时的范德华能最大，说明此时石墨烯与烷烃的混合分子与Si3N4和GCr15的表面吸附力较大；然后范德华能从大到小的排序是石墨烯质量分数为0、3.3%、4.2%时。石墨烯质量分数从0到2.5%范德华能增大，是因为不含石墨烯时只有正十六烷烃油分子的作用，石墨烯质量分数为2.5%时不只有油分子的作用，还要加上油分子与石墨烯分子之间的相互作用，所以范德华能增大。石墨烯质量分数依次为2.5%、3.3%、4.2%时建模的过程中，石墨烯的含量增加，相应的油分子含量就减少了。相对分子质量从248变化到了192和146，其余的相对分子质量不变。范德华能的不同来自于石墨烯和油分子的相互作用。根据范德华能计算公式可知，特征半径σ和特征能量ε数值越大，范德华能越大。由表1中可以看出，石墨烯C-C的特征能量ε和特征半径σ均比十六烷烃油分子中的CH2、CH3小，所以随着石墨烯质量分数变成3.3%和4.2%时，范德华能反而减小了。而且石墨烯是层状物质，比表面积较大，随着石墨烯含量的增加，石墨烯易团聚，起到的润滑作用较小。

图4 不同石墨烯含量时的范德华能随时间变化曲线Fig 4 van der Waals change with time at different graphene content

图5表示上壁面与润滑油膜之间剪切应力随石墨烯含量变化曲线。可见，添加适量石墨烯可有效降低润滑油与壁面间的剪切应力，进而降低摩擦。从图4可知，石墨烯质量分数为2.5%时，范德华能最大，说明其润滑油膜在上壁面的吸附能力最强，润滑较好，这与图5中质量分数2.5%时剪切应力最小的结果一致。其他3种质量分数下的范德华能均比质量分数为2.5%的小，吸附能力较差，剪切应力均比2.5%的大。可以看到图5中的剪切应力呈现了先降后升的现象。

图5 剪切应力随石墨烯含量变化的曲线Fig 5 The curve of shear stress with graphene content

综合以上的结果，可以看出石墨烯质量分数为2.5%时，有较好的吸附性能和较小的剪切应力。因此文中在研究压力和速度对含石墨烯润滑油润滑下Si3N4-GCr15摩擦副润滑性能的影响时，均采用质量分数为2.5%石墨烯润滑油。

2.2 压力对摩擦副润滑性能的影响

为了研究压力对润滑油润滑性能的影响，设置剪切速度为141 m/s，压力分别为25、51、76、102、127、153、178 MPa。

图6所示为压力对Si3N4与GCr15之间润滑油范德华能的影响，图6(a)和图6(b)分别为石墨烯质量分数为2.5%和不含石墨烯条件下仿真结果。可以看出，在含石墨烯和不含石墨烯的条件下，随着时间的增长，范德华能逐渐趋于稳定；且随着压力的增加，范德华能呈现增加的趋势。不同的是含石墨烯的润滑油的范德华能在相应的压力下均比不含石墨烯的范德华能要大。范德华能越大，吸附能力越好。因此，添加石墨烯有利于形成润滑膜，且润滑膜对压力更加敏感。

图6 压力对Si3N4与GCr15之间润滑油范德华能的影响Fig 6 Effect of pressure on the van der Waals energy of the lubricant between Si3N4 and GCr15(a) Graphene mass fraction 2.5%;(b) Without graphene

润滑油分子与上下壁面分子间的相互吸附作用，使得上下壁面附近润滑油分子运动速度与上下壁面速度相同，产生层状结构，出现类固性。不管是随着压力还是速度的变化，液体膜都将有一部分分子的速度与壁面的速度相同，这一部分分子被认为出现了类固性。所以，称液体膜速度与壁面速度相同部分的结构称为类固膜[21]。图7所示为不同压力下，速度沿着膜厚度方向的分布。仿真时，模型的上下壁面施加大小相等、方向相反的速度，速度值为141 m/s。图7中表明，润滑膜出现了分层结构，速度出现了不连续性的特征，由此可以断定，薄膜润滑状态存在着界面滑移，压力越大，滑移现象越明显。同时，近壁面处的分子由于吸附在摩擦副表面，润滑油分子的速度与壁面速度基本一致。这是因为在2个壁面处，具有很强烈的吸附作用，形成类固膜。图7(a)和图7(b)表明，随着压力的变化，类固膜厚度不同。压力增大，液体膜厚度降低，类固膜的厚度逐渐增大。且润滑油中含石墨烯时，类固膜厚度比不含石墨烯时要厚。同时可以看出上壁面附近润滑油类固膜较厚，说明GCr15分子对润滑油分子吸附力更大。

图7 不同压力下沿着膜厚度方向速度的分布Fig 7 Distribution of velocity along the film thickness under different pressures(a) Graphene mass fraction 2.5%；(b) Without graphene

图8所示为相同的速度下，石墨烯质量分数为2.5%和不含石墨烯时上下壁面与润滑油膜间剪切应力随压力变化曲线。模拟中上下壁面施加的速度是相反的，所以上下壁面的剪切应力大小相同，方向相反。可以看出，随着压力的增大，剪切应力是逐渐增大的。且润滑油中不含石墨烯时，润滑油与上下壁面间的剪切应力比含石墨烯时大。剪切应力越大，流动阻力越大，所以摩擦力在数值上会越大，润滑效果越差。

图8 上下壁面剪切应力随压力变化曲线Fig 8 The curves of shear stress with pressure on upper and lower wall

2.3 剪切速度对摩擦副润滑性能的影响

为了研究剪切速度对润滑油润滑性能的影响，设置压力为102 MPa，施加的剪切速度分别为47、94、141、188、236、283 m/s。

图9所示为剪切速度对Si3N4与GCr15之间润滑油范德华能的影响。

图9 剪切速度对Si3N4与GCr15之间润滑油范德华能的影响Fig 9 Effect of shear velocity on the van der Waals energy of lubricant between Si3N4 and GCr15(a) Graphene mass fraction 2.5%；(b)Without graphene

从图9可以看出，不同的条件下，范德华能都随时间增长而增长，并逐渐趋于平稳。且随着剪切速度的增加，范德华能有增长的趋势。其中剪切速度为47～141 m/s之间时，范德华能增长得较慢，当剪切速度增长到141～283 m/s之间时，范德华能增长得较快。不含石墨烯的润滑油，Si3N4与GCr15之间润滑油范德华能比含石墨烯时的范德华能小，吸附能力也较差。

图10所示为剪切速度对润滑油类固膜厚度的影响。可以看出，随着剪切速度的增加，液体膜厚度逐渐增加，类固膜厚度逐渐减小。另外，润滑油含石墨烯时，近壁面处类固膜厚度要大于不含石墨烯时的类固膜厚度。

图10 不同剪切速度下沿着膜厚度方向的速度分布曲线Fig 10 Distribution of shear velocity along the film thickness at different velocities(a) Graphene mass fraction 2.5%；(b) Without graphene

图11所示为上下壁面与润滑油膜间剪切应力随剪切速度变化曲线。可以看出，不论润滑油中是否有含石墨烯，随着剪切速度的增大，剪切应力都在逐渐增大。且在润滑油不含石墨烯时，上下壁面与润滑油膜间的剪切应力比含石墨烯时的剪切应力大。剪切应力越大，流动阻力就越大，摩擦也就越明显。