分子模拟技术

目前, 关于层状矿物的研究已有诸多报道。Zeng等[6]利用分子动力学 (Molecular dynamics, MD) 方法研究了有机蒙脱石的层间结构, 结果表明烷基链在单链插层的蒙脱石层间呈现单层、双层、假三层分布, 而在双链插层的蒙脱石层间能够出现假四层分布。Hackett等[7]采用MD方法研究了烷基铵改性层状硅酸盐的结构, 发现层间距随烷基铵浓度的增加而增大, 且烷基铵在硅酸盐层间分层排列。Boek等[8,9]用蒙特卡罗 (Monte Carlo, MC) 法对蒙脱石的膨胀特性进行了研究, 结果表明蒙脱石的层间距与水含量成一定的函数关系, 且理论分析结果和实验结果吻合。Pablo等[10]研究了钠蒙脱石的水化情况, 结果显示水分子在钠蒙脱石层间域中分为两层。在有机物/氧化石墨插层复合物的研究中, 林舜嘉等[11]研究了CnTAB链长、浓度对CnTAB/GO插层复合物层间距的影响, 结果表明, 随着CnTAB链长和浓度的增加, 插层复合物的层间距逐渐增大。Matsuo等[4]通过聚苯胺与氧化石墨中十六烷基铵之间的交换反应制备了聚苯胺/氧化石墨插层复合物, 结果表明复合物的层间距为1.14~1.68nm, 且聚苯胺在氧化石墨层间呈两层平行排列。Bissessur等[12]利用剥离重组法成功制备了聚吡咯/氧化石墨插层复合物。然而, 关于有机物/氧化石墨插层复合物的研究主要集中在实验方面, 在分子模拟[13]方面的研究尚未见文献报道。对CnTAB插层GO进行分子力学模拟的研究, 有利于进一步揭示CnTAB/GO插层复合物的结构及作用机理。

本工作在课题组前期实验研究基础上[11], 运用分子力学 (Molecular mechanics, MM) 方法在实验值的范围之内对CnTAB (n=12, 14, 16, 18) /GO插层复合物的层间距、CnTAB分子在GO层间的分布、插层复合物的价键能、非价键能进行模拟分析, 旨在从分子力学的角度对CnTAB分子在GO层间的分布及插层作用机理进行研究。

1 模型建立与模拟方法

1.1 模型建立

实验研究表明GO结构中含氧官能团以羟基 (-OH) 和环氧基 (C-O-C) 为主, 仅在结构层边缘接有羧基 (-COOH) 或羰基 (-C=O) [14,15]。根据前人关于氧化石墨烯结构的理论研究[16], 所建立的GO结构模型 (图1) 主要包含-OH、C-O-C两种官能团, 且官能团为异侧邻位分布, n (C) /n (O) 为2/1。GO模型采用12×12×1的超晶胞结构, 在XOY面上的大小为29.52×29.52, c值为GO的实验值, 由X射线衍射分析为8.45[11]。

模拟盒子应用三维周期性边界条件, 在空间3个方向无限重复来模拟宏观体系。初始时, 把1个CnTAB分子 (图2) 以平躺方式置入GO层间, 对其进行结构优化, 然后逐个加入CnTAB分子, 直到插层复合物优化后的层间距达到实验中插层饱和时的值为止, 研究其相关变化, 插层复合物的结构模型如图3所示。

1.2 模拟方法

模拟软件采用Materials Studio[17], 结构优化运用UFF (Universal Force Field) 力场。UFF力场几乎包含了周期表中所有的元素, 它对预测有机分子、主族无机物分子等几何构象具有很好的精度。UFF力场的函数形式[18]为:

式中:Etotal为体系总能量;Evalence为价键作用能;Enonbond为非价键作用能;Ebond为键伸缩能;Eangle为键角弯曲能;Etorsion为二面角扭曲能;Einversion为反转能;EvdW为范德华能;Ecoulomb为静电能。

能量优化采用Forcite模块中的Smart Minimizer法。优化的参数设置为:超软 (Ultra-fine) 收敛, 其中能量偏差为2.0×10-5kcal/mol, 原子均方根力为0.001kcal/ (mol·) , 晶胞的均方根应力为0.001GPa, 均方根位移为1.0×10-5, 迭代步数为10万步, 长程静电作用采用Ewald求和法, 范德华力的截断半径为0.9nm[19]。

1.3 模拟结果

分子力学模拟获得了CnTAB/GO插层复合物优化后的层间距 (图4) 、CnTAB分子在GO层间的分布 (图5) 及插层复合物优化前后的价键能、非价键能 (表1) 。

2 结果与讨论

2.1 CnTAB分子数和碳原子数对CnTAB/GO插层复合物层间距的影响

利用分子力学方法对构建的模拟体系进行结构优化, 根据布拉格公式2dsinθ=nλ, 采用Reflex模块得到模拟体系的层间距 (图4) 。

图4是CnTAB/GO插层复合物优化后的层间距与CnTAB分子数和碳原子数的关系图。由图4可以看出:随着CnTAB分子数的增加, CnTAB/GO插层复合物优化后的层间距 (d (001) ) 呈阶梯状逐步增大, 达到实验插层饱和值时[11], 插层的CnTAB分子数分别达到24、26、28和30, 形成的阶梯数分别为4、5、5和6。当CnTAB分子数分别为1~8、1~7、1~7和1~7时, d (001) 分别达到15.6、15.7、15.6和15.7 (处于图4中第一个阶梯) , 碳原子数对d (001) 的影响极小;当CnTAB分子数分别为9~14、8~12、8~12和8~12时, d (001) 分别达到19.6、19.7、19.8和20.0 (第二个阶梯) , 碳原子数对d (001) 有较小影响;当CnTAB分子数分别为15~18、15~18、15~17和14~16时, d (001) 分别达到22.4、23.1、23.3和24.6 (第三个阶梯) , 碳原子数对d (001) 的影响比较明显, d (001) 随着碳原子数的增加而增加;当CnTAB分子数分别为20~24、20~22、19~22和18~21时, d (001) 分别达到26.4、27.2、27.6和27.8 (第四个阶梯) , d (001) 随着碳原子数的增加而增加;当CnTAB (n=14、16、18) 分子数分别为23~26、26~28和23~25时, d (001) 分别达到30.0、34.0和31.6 (第五个阶梯) , 碳原子数对d (001) 有明显影响;当CnTAB (n=18) 分子数为28~30时, d (001) 达到36.5 (第六个阶梯) 。

分析认为, d (001) 随着CnTAB分子数的增加呈非线性的阶梯状逐步增大, 这与Boek等[20]研究的蒙脱石层间距呈现阶梯状增加相类似。这是因为GO结构层之间是通过范德华力连接的, 当插层的CnTAB分子数较少时, CnTAB分子在GO层间域中呈单层分布, 相互之间的作用力较弱, 只能有限地撑开GO的结构层 (图4中第一个阶梯) , 当继续增加CnTAB分子数时, CnTAB分子之间的相互作用力也随之增加, 在单层分布饱和时开始形成第二层, 使GO的层间域进一步变大, 进而使d (001) 增大 (第二个阶梯) , 再继续增加CnTAB分子数, 当双层分布达到饱和时开始形成第三层 (第三个阶梯) , 接着再继续不断地增加CnTAB分子数, CnTAB分子将形成四、五和六层分布 (第四、五和六个阶梯) , 从而出现插层复合物的层间距随着CnTAB分子数的增加而呈阶梯状增加的变化趋势。因此, CnTAB/GO插层复合物优化后的层间距与CnTAB分子数和碳原子数有关。

2.2 CnTAB分子在GO层间的分布

结合图4的分析, 以C18TAB/GO插层复合物为例, 在每个阶梯阶段随机选取 (C18TAB分子数分别为6、12、16、22、26、30) 一个C18TAB/GO插层复合物优化后的结构 (图5) 。

图5为C18TAB/GO插层复合物在结构优化后所形成的能量相对较低、结构相对稳定的构型。由图5可以看出:插层复合物优化后的结构相对于初始结构 (图3) 有较大变化。首先, 初始结构层面较为平整, 且初始超晶格为规则的六方格子, 而优化后的结构层面出现不同程度的弯曲及倾斜。其次, 最明显的变化是, 初始时CnTAB分子以单层平躺方式分布在GO层间 (图3) , 优化后CnTAB分子在GO层间出现分层, 这与Zeng等[6]研究的烷基链在蒙脱石层间出现单层、双层、假三层, 甚至假四层分布很相似。随着C18TAB分子数由少到多, C18TAB分子在GO层间呈现1~6层分布, 而CnTAB (n=12, 14, 16) /GO插层复合物具有相似的分层情况, CnTAB分子在GO层间分别呈现1~4、1~5和1~5层分布。

分析认为, CnTAB分子之间以及与GO的结构层之间存在多种相互作用力 (排斥力、吸引力、范德华力等) , 在结构优化的过程中, 为了使插层复合物的结构达到最稳定及能量最低化, CnTAB分子在多种力的作用下在GO层间开始运动, 寻找能量最低的位置, 从而导致优化后的结构发生不同程度的弯曲和倾斜;当结构稳定时, 随着CnTAB分子数的增加, CnTAB分子在GO层间呈现不同层数的分布, 说明此时的结构能量最低, CnTAB分子在GO层间最佳的位置是分层分布, 与层间距的分析一致, 反映了结构变化对层间距的影响。

2.3 CnTAB/GO插层复合物优化前后价键能、非价键能的变化

为了阐明插层复合物的作用机理, 以C18TAB/GO插层复合物为例, 对插层复合物优化前后的价键能、非价键能进行分析 (表1) 。

表1为C18TAB/GO插层复合物的能量优化参数表。由表1可以看出, 当C18TAB/GO插层复合物结构稳定时, 它们的总能量较优化前有很大降低, 模拟结果符合能量最低原理。插层复合物中的能量, 优化前以非价键作用的范德华能为主, 其次为价键作用的键角弯曲能, 优化后以价键作用的键角弯曲能为主, 其次为非价键作用的范德华能, 这说明C18TAB/GO插层复合物结构优化前以非价键作用的范德华作用为主, 优化后以价键作用的键角弯曲为主。在结构优化过程中除了由结构扭曲引起的扭曲能稍微有所升高之外, 其余价键能和非价键能都有不同程度的降低, 其中以范德华能降低幅度最大。

图6为C18TAB/GO插层复合物优化后的键伸缩能、键角弯曲能、二面角扭曲能、反转能、范德华能及它们的拟合曲线图。由图6可以看出:在优化后的各项能量中, 键伸缩能、键角弯曲能和二面角扭曲能与C18TAB分子数具有很好的线性关系, 随着C18TAB分子数的增加呈线性增加, 反转能几乎保持不变, 而范德华能则分散不均, 这说明C18TAB/GO插层复合物优化后的键伸缩能、键角弯曲能、二面角扭曲能和范德华能来自于氧化石墨和季铵盐, 反转能主要来自于氧化石墨, 而CnTAB (n=12, 14, 16) /GO插层复合物优化后的各项能量也具有相同的变化规律。

3 结论

(1) CnTAB/GO插层复合物优化后的层间距随CnTAB分子数的增加呈阶梯状增加, 层间距的变化也受碳原子数的影响。

(2) CnTAB/GO插层复合物优化后的结构层面较初始结构发生不同程度的弯曲及倾斜, CnTAB分子在GO层间呈层状分布, 随CnTAB分子数的增加, 逐渐由一层增加至多层, 并受碳原子数的影响。

分子模拟技术 篇8

接触表面在接触滑移过程中的摩擦问题一直是摩擦学研究的热点。随着纳米尺度精密仪器的出现,在微/纳米器件的设计、制造和使用等过程中,发生在表面和界面的摩擦学、热学和化学等行为占据了主导地位[1],零件表面微观形貌对摩擦、磨损、润滑等摩擦学性能的影响至关重要。近年来,国内外学者在研究规则表面形貌提高材料的摩擦学性能方面取得了一定的成果。Menezes等[2]研究了纯铝探针在钢纹理表面上的滑动过程,发现摩擦因数主要取决于表面纹理的方向性。He等[3]研究了宏观和微观情况下聚二甲基硅氧烷(PDMS)表面纹理对摩擦因数的影响,结果表明,在微观情况下,柱状纹理表面相比光滑表面摩擦因数降低了38%。Wang等[4,5]通过无定形硅低温结晶技术修饰了不同粗糙度的表面纹理,通过实验指出,带有纹理的表面黏着力和摩擦因数显著低于非纹理表面。Xie等[6]研究了离子液体润滑条件下,三种不同表面粗糙度的硅基材料(低温氧化硅,多晶硅,Si3N4)在不同载荷和滑动速度下的摩擦特性,结果表明,具有不同表面粗糙度的摩擦副的摩擦特性相差较大,在低速滑动条件下,表面粗糙度较小的摩擦副Si3N4/Si3N4的摩擦因数较大。董云开等[7]研究了不同微米级条纹形貌硅试样的摩擦特性,结果表明,微尺度下形貌参数存在优选值以有效减小速度对摩擦因数的影响。吕文斐等[8]采用YAG纳秒激光器在45钢平板试样表面制备纹理,考察了不同润滑状态下表面纹理对材料摩擦学性能的影响,实验结果表明,在一定的润滑剂黏度、摩擦速度和载荷下,纹理表面有利于提高材料的摩擦学性能。由于实验仪器的限制,通过实验方法不能彻底了解微/纳观接触表面的深层摩擦机理,而分子动力学方法可以了解原子/分子尺度物质运动的细节行为并帮助解释实验结果,已经成为研究纳米尺度下摩擦学问题的主要方法之一[9]。

Yang等[10]运用分子动力学方法模拟了金刚石压头在铜基体上的滑动过程,研究了压头和基体之间的界面干涉以及压头形状和大小对摩擦因数的影响。Zhu等[11]研究了滑动速度在10m/s到200m/s之间时金刚石压头在银基体上的滑动过程,发现金刚石压头的剪切变形是发生黏滑现象的主要原因。佟瑞庭[12]运用分子动力学和有限元相结合的方法研究了二维情况下单晶铜刚性圆柱压头与纹理表面的滑动接触,模拟结果表明,梯形纹理的减摩效果优于矩形纹理,且在一定的凸峰高度和凸峰间距下,纹理表面可以有效地降低滑动摩擦力。

目前,对于纹理表面相关摩擦性能的研究主要是在宏观和微观下采用实验方法进行的,在纳米尺度对纹理表面摩擦学性能的相关研究较少,并且通过以往的研究发现[13,14,15],并非所有的纹理表面都有助于材料摩擦学性能的提高,因此有必要研究纳米尺度不同纹理表面的具体摩擦学行为。

在连续介质理论中,由于两弹性物体的相互接触可以等效为一个刚性体和一个弹性体的接触[16]。因此在本文的模拟中为了节约计算时间,借鉴以往分子动力学模拟方法[17],采用半球形刚性压头,并在此基础上建立了单晶铜半球形刚性压头与纹理表面黏着滑动接触的三维分子动力学模型,研究了纹理方向和纹理密度对单晶铜纹理表面黏着滑动摩擦行为的影响。模拟结果表明:在一定纹理密度范围内,表面纹理方向对滑动摩擦力有很大影响;对于一定的纹理方向,滑动摩擦力随着纹理密度的增大而增大;在一定纹理密度范围和纹理方向下,纹理表面可以有效地降低滑动摩擦力。

1 分子动力学方法与模型

模拟中采用嵌入原子势(EAM)来描述基体原子间及基体和压头原子间的相互作用。EAM势已被证明能够很好地描述金属原子间的相互作用,是目前在金属分子动力学模拟中应用最广泛的势函数。其基本方程如下:

$\begin{array}{l}E{}_{\text{t}\text{o}\text{t}}={\displaystyle \sum _{i}E}{}_i\\ E{}_i=\frac{1}{2}{\displaystyle \sum _j\text{ϕ}}(r{}_{ij})+F(\rho {}_i)\\ \rho {}_i={\displaystyle \sum _{j\ne i}f}(r{}_{ij})\end{array}$

式中,Etot为系统总能量;Ei为第i个原子的能量;ρi为第i个原子处的电子云密度,即周围原子的电子云密度在原子i处的叠加;f(ri j)为原子j在距离ri j处产生的电子云密度;F(ρi)为镶嵌能函数,代表了将原子i嵌入到密度为ρi的背景电子云中的嵌入能;ϕ(ri j)为对势函数;ri j为原子i、j间的距离。

图1为初始原子构型的俯视图和主视图,坐标轴x、y、z方向分别对应单晶铜的[1]、[0 1 0]、[0 0 1]晶向。半球形刚性压头的半径为10d(d为面心立方铜晶体的最近邻原子距离,d=2.552×10-10m),由2938个原子组成。弹性基体尺寸沿x、y、z方向分别为31a0×31a0×10a0(a0为面心立方铜晶体的晶格常数,a0=3.61×10-10m),即11.191nm×11.191nm×3.61nm。光滑表面基体由38 440个原子组成,纹理表面基体的原子数随纹理密度和纹理方向的不同而发生变化。

(a)俯视图 (b)主视图

模拟时,压头与基体初始间距控制在2.5d。在x方向和y方向施加周期性边界条件,z方向上表面为自由边界条件,z方向最底2层原子为固定原子层。原子初始速度服从Maxwell-Boltzmann分布,并按照统计力学正则系综(NVT)将基体温度控制在300K。时间积分算法采用Velocity-Verlet算法,时间步长Δt取5fs。初始计算前基体原子自由弛豫2000Δt,达到平衡状态。模拟过程分为两个阶段:第一阶段压头以恒定速度4m/s沿z轴负方向移动1.56d;第二阶段在z方向上对压头不再施加位移,压头以恒定速度4m/s沿[1]晶向滑动,滑动距离为29.73d。滑动摩擦力定义为滑动过程中基体对压头的作用力沿x方向分量总和的相反数。本文重点讨论了滑动阶段纹理密度和纹理方向对滑动摩擦力的影响。

图2为纹理表面示意图,图中,箭头所指方向为[1]晶向,即压头滑动方向,通过衡量滑动方向与纹理方向的夹角,本文设计了0°、45°和90°三种纹理方向,分别对应单晶铜的[1]、[1]和[0 1 0]晶向。图2中,a表示凸峰的宽度,b表示凹槽的宽度,凹槽深度均为2a0。通过改变a、b的值来控制纹理密度(纹理密度β定义为凸峰宽度与纹理周期的比值,即β=a/(a+b)),五组不同纹理密度表面的相关参数如表1所示。

(a)0°纹理方向 (b)45°纹理方向 (c)90°纹理方向

2 模拟结果验证及分析

2.1模型验证

光滑表面的滑动摩擦力曲线如图3所示,图中,u表示压头的滑动位移,F表示滑动摩擦力,F和u的数值均为对比量,即将EAM势中的对比单位量纲一化之后的结果(f=eV/d=6.268×10-10N为EAM势中力的对比单位,d是长度的对比单位)。文中所出现的滑动摩擦力F和滑动位移u都是量纲一化之后的数值。从图3中可见,滑动摩擦力曲线呈规则锯齿形变化。该摩擦力曲线的波动周期接近单晶铜的晶格常数,即保持在0.361nm左右。这与以往光滑表面黏滑效应的研究结果[18]相吻合,验证了本文模拟方法的正确性。

2.2纹理方向对滑动摩擦力的影响

表面纹理密度β=40%时,纹理方向分别为0°、45°、90°的滑动摩擦力曲线如图4所示,由图4可知,在平稳滑动阶段,45°方向纹理表面的滑动摩擦力振幅大于0°和90°纹理方向。纹理密度β为33%和50%时,滑动摩擦力也出现这种变化趋势。图5所示为不同纹理密度下滑动摩擦力的平均值,可以发现:对于五组纹理密度,0°方向纹理表面的滑动摩擦力始终最小。这是由于模拟中,压头沿[1]晶向滑动,基体为0°纹理方向时,压头恰好位于凹槽的上方,这样导致基体表面与压头作用的原子数减少,滑动摩擦力减小。而45°和90°纹理方向时,接触面积增大,与压头作用的基体原子数增多,黏着力增大,从而导致滑动摩擦力增大。

另外,从图5可以发现,当表面纹理密度小于50%时,45°方向纹理表面的滑动摩擦力平均值大于0°和90°方向纹理表面的滑动摩擦力平均值,这是因为单晶铜在不同晶向上原子分布密度、排列方式和原子间结合力不同。图6为纹理密度为40%时45°方向纹理表面在u/d=29.73,即滑动位移最大时的xz剖面的原子构型图,可以发现,45°方向纹理表面在(1 1 1)晶面的[1]晶向上发生晶体滑移,但0°和90°方向纹理表面没有发生晶体滑移,只在纹理表面出现部分原子迁移。当纹理密度为33%时,也出现此种情况。由于晶体滑移在原子排列稀疏的平面和晶向上阻力很大,在原子排列紧密的平面和晶向上阻力很小,因此晶体滑移一般在最密排晶面的最密排晶向上发生。对于面心立方铜晶体,晶体滑移一般发生在(1 1 1)晶面的[1]晶向上[19],这与本文的模拟结果吻合。因此当纹理密度较小时,45°方向纹理表面最易发生滑移塑性变形,这也是45°方向纹理表面滑动摩擦力较大的原因之一。

由图4可以发现,当滑动位移u/d<10时,三个纹理方向的滑动摩擦力均较小,随着滑动位移的增大,滑动摩擦力逐渐增大。纹理密度为33%、40%、50%和60%时,滑动摩擦力均出现这种变化趋势。这是由于表面纹理的引入使基体与压头之间出现原子吸附现象,随着滑动位移的增大,一些基体原子黏附在压头的下方,成为压头原子的一部分,与压头一起向前滑动,黏着作用显著增强,从而导致滑动摩擦力增大。当纹理密度为67%时,滑动摩擦力的这种前后增幅并不明显,这是由于纹理密度增大,凸峰宽度a增大,凹槽宽度b减小,凸峰原子稳定性增强,导致原子吸附作用减弱。这说明表面纹理密度较大时,原子吸附对滑动摩擦力的影响较小。

图7为光滑表面和90°方向纹理表面在不同滑动距离时的xz剖面原子构型图,其中纹理表面密度为50%。当压头与基体的间距控制在1.0d,基体为光滑表面时,滑动过程中压头与基体间只在接触点处出现突跳接触,压头滑过之后,基体原子回复初始位置,未出现原子吸附现象,而引入纹理表面后,出现了原子吸附现象。这说明纹理的引入使得纹理边界原子处于不稳定状态,所以有摩擦作用时发生原子迁移。原子吸附导致了滑动摩擦过程中原子的堆积以及位错的产生,滑动过程中在压头周围产生原子的堆积,滑动过后,纹理表面被破坏,原子重新排列,且在小范围内原子按照面心立方晶体结构重构,这与宏观晶体生长理论相符。

2.3纹理密度对滑动摩擦力的影响

在原子尺度,摩擦力与接触表面间的黏着有关,且黏着作用与接触面积密切相关,表面纹理密度直接影响到压头与基体之间的接触面积,因此纹理密度会对滑动摩擦力造成一定影响。本文在表1所示纹理密度参数模拟结果的基础上,探讨了纹理密度对滑动摩擦力的影响。图8为三组纹理密度表面在0°和45°纹理方向下的滑动摩擦力曲线,其中第1、3、5组表面的纹理密度分别为33%、50%、67%。由图8可以发现,对于0°方向纹理表面,当纹理密度相差较大时,滑动摩擦力振幅随着纹理密度的增大而增大。同时90°方向纹理表面表现出与0°方向纹理表面相同的变化规律,而45°纹理方向时,滑动摩擦力没有出现这种变化趋势,说明对于45°方向纹理表面,纹理密度对滑动摩擦力振幅的影响较小。

图9所示为三种纹理方向下不同纹理密度表面滑动摩擦力的平均值,总体变化趋势为滑动摩擦力平均值随着纹理密度的增大而增大,纹理密度相差越大,滑动摩擦力的增幅也越大。在相同的纹理方向下,纹理密度的增大直接导致滑动过程中基体表面与压头作用的原子数增加,黏着作用增强,从而使滑动摩擦力增大。说明在相同接触条件下,表面纹理密度增大是导致滑动摩擦力增大的因素之一。

如图9所示,纹理密度为60%和67%时,对于0°和90°方向纹理表面,滑动摩擦力的平均值较为接近,且滑动摩擦力平均值随着纹理密度的增大而略有减小。通过观察模拟过程的动画可以发现,纹理密度为67%时,滑动位移较大时只在压头底层原子周围吸附了少量的基体原子。纹理密度为60%时,在压头底层吸附了一层基体原子,且这层基体原子随着压头的滑动而移动。原子吸附作用的增强直接导致接触面积增大,滑动摩擦力增大,在一定程度上削弱了纹理密度对滑动摩擦力的影响。

2.4光滑表面与纹理表面滑动摩擦力的对比分析

图10对比了不同方向纹理表面与光滑表面的滑动摩擦力,其中第1、2组纹理表面密度分别为33%、40%。可以发现,对于0°方向纹理表面,当滑动位移u/d<20时,纹理表面的滑动摩擦力小于光滑表面,且90°方向纹理表面表现出与0°方向纹理表面相同的规律;而对于45°方向纹理表面,仅在滑动的初始阶段,即u/d<10时,纹理表面能够有效减小滑动摩擦力。这说明,对于较小的纹理密度,0°和90°方向纹理表面的减摩效果优于45°方向纹理表面。同时模拟结果显示:纹理密度大于50%时,仅在滑动的初始阶段,三个方向纹理表面的滑动摩擦力小于光滑表面,随着滑动距离的增大,纹理表面的滑动摩擦力大于光滑表面。而纹理密度小于50%时,0°和90°方向纹理表面的滑动摩擦力在u/d<20时均小于光滑表面。这说明纹理密度较小的表面减摩效果较好,同时表明纹理密度和纹理方向的优化组合能够更有效地减小滑动摩擦力。

对于面心立方晶体,研究证明,在纳米尺度光滑表面滑动摩擦力的变化周期等于晶体在滑动方向的晶格常数。从图10中发现,相比光滑表面,纹理表面的引入并没有改变滑动摩擦力的波动周期,这与宏观和微观纹理表面的摩擦实验结果存在差异。在文献[3]中,通过实验得到在宏观和微观情况下,聚二甲基硅氧烷(PDMS)矩形纹理表面摩擦因数的波动周期与纹理周期相等。这说明纹理表面在纳米尺度黏滑摩擦的作用机理与宏观和微观情况存在一定差异,有待进一步地研究。

3 结论

(1)当表面纹理密度小于50%时,0°和90°方向纹理表面的滑动摩擦力振幅小于45°纹理方向的滑动摩擦力振幅;对于本文所模拟的五组纹理密度,0°方向纹理表面的滑动摩擦力平均值始终最小;当纹理密度小于50%时,45°方向纹理表面最易发生滑移塑性变形。

(2)纹理表面的引入导致原子吸附现象的发生,对滑动摩擦力产生一定影响。在本文的模拟范围内,表面纹理密度为67%时,原子吸附对滑动摩擦力的影响较小。

(3)表面纹理密度影响滑动摩擦力的大小。纹理密度相差较大时,0°和90°方向纹理表面的滑动摩擦力振幅随着纹理密度的增大而增大,且表面纹理密度相差越大,滑动摩擦力的变化也越大。三种纹理方向下,滑动摩擦力的平均值随着纹理密度的增大而增大。

(4)在纳米尺度下,存在纹理方向和纹理密度的优选值,使得滑动摩擦力最小。在本文的模拟范围内,当纹理密度小于50%,在u/d<20时0°和90°方向纹理表面的滑动摩擦力小于光滑表面的滑动摩擦力。当基体表面纹理密度为33%且纹理方向为0°时,纹理表面能够最有效地减小滑动摩擦力。

(5)在纳米尺度,光滑表面滑动摩擦力的波动周期等于晶体在滑动方向的晶格常数。纹理表面对滑动摩擦力的大小产生一定影响,但没有改变滑动摩擦力的波动周期,这与微观和宏观情况存在一定差异。

分子模拟技术 篇9

由于纳米技术的发展,微机械在航空航天、医疗、通信等领域得到了广阔的应用。然而当微机械、微构件尺度减小到纳米量级,由于尺寸效应和表面效应的作用[1,2],微机械、微构件的疲劳强度、使用寿命及结构稳定性等都受到了影响。

随着切削加工技术的发展和进步,切削加工精度得到了明显的提高,加工尺度则不断减小,采用传统的试验方法来研究超精密切削加工过程十分困难,并且很多微观现象无法直接获得[3],于是计算机仿真技术被应用于微观尺度切削加工机理研究中[4]。而分子动力学(molecular dynam-ics,MD)模拟方法是一种描述微观现象的有效方法,已成为研究纳米级切削机理的一个重要工具。如Chen等[5]利用分子动力学研究了单晶硅在纳米切削过程中切削力的变化,得出在同样的切削厚度下,切削力随着刀具前角的增大而减小,以及当有切屑形成并在刀具前刀面堆积时,切削力逐渐增加等结论;罗熙淳[6]利用分子动力学方法分析不同变形区不同材料势能的变化,建议采用法向切削力和切向切削力的比值衡量初始切屑的产生。Ye等[7]、Zhang等[8]采用分子动力学方法研究了纳米切削过程中,切削速度对切削过程的影响规律。单晶铜材料是制作微机械、微构件常用的材料,具有良好的信号传输、塑性加工和抗疲劳等性能,被广泛应用于国防精密仪器、民用通信设备等领域[9]。然而,由于单晶铜宏观切削加工理论已不再适用于纳米切削,因此有必要研究纳米尺度下单晶铜表面切削特性,为微机械、微构件的制造加工提供理论依据。本文运用分子动力学模拟方法研究单晶铜材料表面纳米切削特性,采用不同的切削速度及切削厚度对单晶铜表面进行纳米切削,从其对切削力大小、切屑形成、位错及缺陷的影响等方面展开研究。

1 分子动力学建模与计算

1.1 模型建立

建立单晶铜纳米切削过程三维分子动力学模型,如图1所示。

模型基体尺寸为30a×50a×25a,a为Cu的晶格常数(a= 0.3615nm)。将单晶铜基体分成三个类型区域:固定边界区、恒温区和牛顿区。固定边界区为基体底部和左端3个原子层,在模拟过程中保持固定。恒温区为与固定边界区相邻的6个原子层,该区域使系统温度保持恒定,模拟中采用Nose-Hoover热浴法将其温度控制在293K。 牛顿区为除去固定边界区和恒温区的原子,其余的基体原子为自由原子。模拟中,采用微正则系统对模拟体系进行约束。恒温区和牛顿区原子遵从牛顿第二定律,通过运用Velocity-Verlet算法对牛顿方程进行积分,积分步长为1fs。 刀具为金刚石刀具,与参考文献[3]所建立的切削模型相比,此模型考虑了刀尖圆弧半径。另外,为了减小模拟体系中原子小于真实体系中的原子数而产生的尺寸效应,在x和y方向上施加周期性边界条件,在z方向上施加自由边界条件。

1.2 势函数选取

势函数的合理选取对模拟的精确性是非常重要的,在单晶铜纳米切削模拟过程中,包含3种不同的原子间相互作用。由于金刚石刀具被设为刚体,故可以忽略C-C原子间的相互作用。在本次模拟中,对于基体中Cu-Cu原子相互作用采用嵌入原子势(EAM势)。对EAM势,整个系统的原子势能Etot如下:

式中,φij为原子i和原子j之间的对势;rij为原子i和原子j之间的距离;Ei为原子i嵌入时发生的嵌入能;ρi为除i之外所有原子在i处产生的电子密度;ρj为原子j在原子i处产生电子密度函数。

基体原子与刀具原子(Cu-C)之间的作用采用Morse势函数来描述:

式中,Ev(rij)为原子间相互作用的势能;ED为结合能;α 为弹性模量;r0为平衡态的原子间距。

其相应参数D=0.087eV,α=0.5140nm-1,r0=0.2050nm。

1.3 分子间作用力的计算

势函数确定以后,原子之间的作用力Fij就可以通过势函数对rij求导得出,即

作用在第i个原子上的总原子力等于其周围所有原子对该原子作用力的合力[10],即

模拟中为了让系统能量趋于稳定状态,对模型进行弛豫10ps,刀具距离基体右端为1nm,然后使刀具沿着y轴负方向进行切削。为了研究不同的切削速度及切削厚度对单晶铜材料表面切削特性的影响,切削速度v分别取50m/s、100m/s、200m/s,切削厚度δ 分别为0.5nm、1nm、1.5nm,切削距离为12nm,刀尖圆弧半径为1nm。

影响材料塑性变形的重要原因是位错及晶格缺陷,但由于原子的热振动影响,中心对称参数(centro-symmetry parameter,CSP)要比滑移矢量和原子配位数法更容易辨别位错等缺陷。而相比于中心对称参数法,共近邻分析(common neighbor analysis,CNA)作为一种广泛应用于金属材料的晶体缺陷分析技术,能明确辨别缺陷的种类和位置[11],因此本文采用中心对称参数和共近邻分析来辨别位错和其他缺陷。中心对称参数描述公式[12]为

式中,Ri和Ri+6为长度相同、方向相反的近邻原子对。

2 模拟结果及讨论

2.1切削速度对单晶铜表面切削特性的影响分析

在纳米切削过程中,不同的切削速度对单晶铜材料表面切屑形貌以及基体内部的缺陷分布区域有显著影响。下面将通过对切削力、切屑形貌、缺陷分布等方面的分析来研究不同的切削速度对单晶铜材料表面切削特性的影响。图2给出切削厚度为1nm时不同切削速度下的单晶铜基体的截面图,并计算原子的中心对称参数并按其值对其着色。

对图2进行分析,可以明显看出不同的切削速度对切屑形貌原子状态及切屑体积有显著影响。切削速度越高,切屑形成的体积越大,而且切屑里原子的排列越紧密。采用共近邻分析来研究不同的切削速度对单晶铜基体内部的缺陷分布区域的影响。图3~图5所示为单晶铜材料内部在不同切削速度下的缺陷结构分布图及底部视图。

如图3~图5 所示,不同的切削速度对单晶铜基体内部位错及缺陷有显著影响。对于单晶铜塑性材料来说,在纳米切削过程中,材料的去除机理主要是由于位错形核和扩散运动引起塑性变形,塑性变形累积到一定程度后,形成切屑,实现材料的去除。当采用较低的切削速度50m/s、100m/s时,位错形核在基体内部激活的一个滑移面上运动。当采用较高的切削速度200m/s时,基体内部多个滑移面被激活,因此基体内部的缺陷分布区域较大。这是由于采用较高的切削速度导致单晶铜基体内部的位错形核及运动的时间较短而引起的[13]。

切削力反映切屑的去除过程,是解释切削现象的重要物理参数,下面根据模拟结果研究不同的切削速度对单晶铜材料表面切削特性的影响。通过对数据的提取,运用MATLAB绘图得到图6,即不同切削速度下的切削力-切削位移曲线图。

由图6可以发现,不同的切削速度下,单晶铜材料纳米切削过程中切削力均在切削初期上升,达到稳定切削状态后围绕稳定值进行波动,其主要原因是工件材料内部位错等缺陷的发生和运动引起切削作用力的波动。但可以看出,在切削初期,切削速度越大,切削力上升幅度越大;为了使不同切削速度下的切削力有一个定量直观的比较,切削力的计算取稳定切削过程4~10nm中切削力的平均值。切削速度为50m/s时,平均切削力为60.342nN;切削速度为100m/s时,平均切削力为66.976nN;切削速度为200m/s时,平均切削力为75.744nN。对不同切削速度下的平均切削力进行对比,可以发现切削速度越大,切削力越大。这一结论与Zhang等[8]的研究结论一致。其主要原因是由于采用较高的切削速度导致单晶铜基体内部位错形核及运动的时间较短,塑性变形较少;较高的切削速度导致了刀具前面形成的切屑堆积体积较大,因此也增大了切削过程的切削力。由此,可以说明切削速度对切削力的大小有显著影响。

2.2 切削厚度对单晶铜表面切削特性的影响分析

在纳米切削过程中,切削厚度对材料表面切削性能有明显影响。以下就切削厚度对单晶铜表面切削特性的影响进行分析。图7所示为切削速度200m/s时,采用三种不同的切削厚度进行切削后单晶铜基体的截面图。

对图7进行分析,可以看到,在切削过程中单晶铜材料表面发生严重的塑性变形破坏,不同的切削厚度对切屑体积大小、位错缺陷等有显著影响。同种切削速度下,切削厚度越大,在刀具前方堆积的切屑体积越大。另外,在刀具下方及左下方基体原子有不同程度的位错等缺陷,切削厚度越大,位错缺陷分布区域越大,并且有明确的位错发射。图8~图10为单晶铜材料内部在不同切削厚度下的缺陷结构分布图及底部视图。

由图8~图10可以看出,不同的切削厚度对单晶铜基体内部位错及缺陷有显著影响。随着切削厚度增大,沿切削方向基体的堆垛层错或孪晶界等缺陷原子增多,单晶铜基体塑形破坏变形严重,主要分布在刀尖前端与基体接触区域并向周围扩散。这是由于采用较小的切削厚度时,位错形核在基体内部激活的滑移面较少,并且基体内部的位错及缺陷分布区域比较规整;当采用较大的切削厚度时,基体内部多个滑移面被激活,因此基体内部的缺陷分布区域较大。

切削厚度对切削力也有显著的影响,通过对数据的提取,运用MATLAB进行绘图。图11为不同切削厚度下的切削力-切削位移曲线图。

由图11可知,不同的切削厚度下,单晶铜材料纳米切削过程中切削力均在切削初期先上升,达到稳定切削状态后围绕稳定值波动。但可以看出,在切削初期,切削厚度越大,切削力上升幅度越大。为了使不同切削厚度下的切削力有一个定量直观的比较,切削力的计算取稳定切削过程4~10nm中切削力的平均值。 切削厚度为0.5nm时,平均切削力为57.186nN;切削厚度为1nm时,平均切削力为75.744nN;切削厚度为1.5nm时,平均切削力为94.256nN。对不同切削厚度下的平均切削力进行对比,可以看出随着切削厚度的增大,切削力也增大。其原因主要是采用较大的切削厚度导致刀具前方形成的切屑体积较大,在切削运动过程中增大了摩擦阻力,导致了切削力的增大。

3 结论

(1)通过对切削过程中原子状态进行分析可知:在纳米切削过程中,不同的切削速度及切削厚度对单晶铜材料表面切屑形貌以及基体内部的缺陷分布区域有显著影响。切削速度越高,切屑堆积体积越大,切屑里原子的排列越紧密,位错缺陷分布区域越大;在同种切削速度下,切削厚度越大,在刀具前方堆积的切屑体积越大。另外,在刀具下方及左下方基体原子有不同程度的位错等缺陷,切削厚度越大,位错缺陷分布区域越大。