单晶材料纳米加工的分子动力学模拟研究进展*

夏斯伟， 周 海， 徐晓明， 张春伟， 陈西府， 黄传锦， 徐彤彤

(盐城工学院 机械工程学院， 江苏 盐城 224051)

纳米加工主要是利用纳米级金刚石磨粒来实现加工，使加工形成的结构与器件的尺寸精度达到纳米量级[1]。随着电子科技、国防等行业产品的更新换代，光学晶体材料的纳米加工成为一个瓶颈问题，深入研究纳米加工的机理是解决此类问题的关键。

国内外研究纳米加工机理的方法主要有2类：(1)实验研究方法，包括纳米金刚石单点切削加工[2]、纳米金刚石研磨加工[3]、基于AFM平台的探针刻划加工、材料加工表面测试表征等；(2)计算机模拟研究方法，包括分子动力学模拟(molecular dynamics simulation，MDS)、基于光滑粒子流体的动力学模拟(SPH)等[4]。纳米加工的切削量一般控制在几个或几十个纳米量级，对实验研究的各项条件要求苛刻，且成本高昂；而计算机模拟可以辅助实验研究，还可以分析一些难以完成的特殊条件实验。如MDS模拟能够实现特定的加工环境温度、介质、压力等条件，起到沟通宏观特性与微观结构的桥梁作用，现已成为最重要的一种模拟方法[5>-6]。

1 分子动力学模拟简介

1.1 分子动力学模拟的基本原理

分子动力学模拟的基本原理是：对所研究的微观现象建立一个多体粒子系统，由量子力学确定这个系统中各粒子之间的相互作用，通过数值求解这些粒子动力学方程组，决定各个粒子在相空间的运动规律和轨迹，由统计物理原理计算出系统相应的宏观物理特性[7]。

根据模拟研究对象的不同，统计物理学将分子动力学模拟划分为平衡态分子动力学模拟EMDS和非平衡态分子动力学模拟NEMDS[8]。其中，EMDS是分子动力学模拟的基础，研究给定的系统向期望的平衡态如何演化，把合理的初始条件提供给动力学加载过程，可以预测材料在平衡态下的热力学性质；而NEMDS则主要研究在外加载荷、外加温度梯度等作用下的初始系统的动力学响应状态。一般先进行EMDS模拟，得到系统在给定条件下的平衡状态，再以此为初始条件输入进行NEMDS模拟[9]。

平衡态分子动力学模拟主要包括3方面问题：(1)建立系统模型；(2)确定系统相互作用势函数；(3)对系统进行数值计算。

1.2 经典平衡态分子动力学模拟方法

1.2.1 系统模拟模型

20世纪90年代初期，BELAK等[10]第一次提出了沿用至今的纳米加工分子动力学模拟模型。如图1，该工件的模型存在一定边界条件，工件原子被分为固定边界层、热稳定层以及牛顿层；其中，边界层原子在模拟过程中保持位置不动，忽略晶格振动，并保证对称型晶格结构；热稳定层的原子要标定以保持工件整体温度恒定；牛顿层是加工对象层，符合牛顿第二运动定律。加工刀具的模型用金刚石原子结构表示。

1.2.2 相互作用势函数

原子间的相互作用势对其坐标距离求导可以获得相互作用力，所以在分子动力学模拟计算中，势函数的选择直接决定了计算工作量，以及计算模型与真实系统的近似程度[7]。为了平衡计算精度与计算耗时，目前主要采用经验势函数方法来研究纳米加工。

经验势函数有很多，其中，单晶材料涉及的势函数可分为：对势和多体势。比如Lennard-Jones[1]、EAM[12]、Morse[11,13]等对势函数，适用于单晶铜、单晶铝晶体；Tersoff[14>-17]、ABOP[18]、Stillinger-Weber[19>-20]等多体势函数，适用于单晶硅、单晶锗。

大部分研究者将金刚石加工刀具的原子定义为刚体。近年来一些学者研究刀具磨损特性时，用于定义C-C原子间相互作用力的有ABOP[18]、Tersoff[21]等势函数。用于定义刀具与工件间相互作用力的有Tersoff[21]、Morse[22]等势函数。

图1 分子动力学模拟的模型边界条件示意图[10]

1.2.3 系统数值算法

经典分子动力学中，除固定边界层外的原子均符合牛顿第二定律，建立牛顿运动方程来描述和模拟体系中一对原子i和j，分析其在t时刻的坐标、速度及相互作用力。如公式(1)，(2)，(3)所示[7]。

rit-rjt=rijt

(1)

(2)

(3)

求解以上二阶微分方程的算法有很多，如Verlet算法[23]、Velocity-Verlet算法[24]、Beeman算法[25]、Gear预测修正算法[26]等。其中Velocity-Verlet算法具有实现位置和速度同步，占用内存低，能并行计算等优点，被广泛应用[7]，如公式(4)、(5)、(6)所示[24]，式中∂Urit为相互作用势函数。

(4)

(5)

(6)

2 单晶材料纳米加工的分子动力学模拟研究

分子动力学模拟在单晶材料纳米加工研究中应用广泛[27]，国内外很多学者针对加工过程中材料的物相转变、结构及热力学特性、介质环境对加工的影响和加工后的亚表面损伤等全过程，进行了系统深入的探索。

2.1 单晶材料纳米加工过程的特性研究

2.1.1 物相转变规律

在纳米加工单晶材料过程中，材料的脆塑性转变规律研究一直以来都备受关注。MDS通常有2类方法研究晶体材料的物相转变规律：(1)近邻分析法，包括配位数、径向分布函数、角分布函数、中心对称性参数等分析方法；(2)Vornoi图形分析法[27>-28]。

1999年起，CHEONG等[29]用MDS研究了单晶硅在纳米压痕时的物相转变规律，通过压痕过程中不同时刻单晶硅的原子图像、Si原子近邻数随时间变化曲线、键长与原子数量分布曲线等方法研究压痕过程。在加载阶段，单晶硅由立方金刚石结构逐步转变为体心四面体结构(β >-Si)，在压痕最大处的静压力为12 GPa，与实验研究获得的静压力基本相同；在卸荷阶段，体心四面体结构转化为非晶结构。对压痕的二次加载过程中，又重新生成体心四面体结构，表明单晶硅的这种物相转化是可逆的。该课题组还研究了多次纳米压痕过程中单晶硅变形特性的规律，结果表明：第一次压痕中压头与单晶硅的接触面积显著大于后续第二次、第三次压痕过程，第一次压痕与多次压痕的其他指标对比不明显[30]。2002年，该课题组通过单晶硅压痕及轴向压缩模拟，研究适合单晶硅物相转变为β>-Si的应力判据[31]，结果表明：采用最大法向应力判据最符合。轴向压缩模拟时，立方金刚石结构向β>-Si相转变为[001]晶向，最大主应力为10 GPa；纳米压痕模拟时，[001]晶向上最大法向应力等于10 GPa的区域就是单晶硅转变为β>-Si相的区域。

2001年，KOMANDURI等[32]发现纳米仿真切削单晶硅时，在加工中由于压力诱导转变，使得部分硅晶体从α>-Si体心立方结构向β-Si结构转变。结果发现，随着切削宽度与深度的比值w/d减小，单晶硅材料的侧向流动以及亚表面层损伤均迅速减少。2003年，TANAKA等[33]使用MDS研究了单晶硅超精密加工过程中的脆塑性转变机理，建立了2种尺寸的单晶硅三点弯曲模型。借助最近邻原子数分析物相转变，其中金刚石结构相的近邻原子数为4，单晶硅加工中的塑性变形区就是金刚石结构相向无定形结构相转变的位置。

2006年，KIM等[34]研究了纳米压痕中单晶硅(100)和(111)晶面的结构转变规律。模拟结果表明：单晶硅的物相转变的微观机理与其晶向、滑移体系均有直接关系。亚表面层Si>-II相转变开始于压头下方几层原子处的最大剪切应力位置处，故纳米压痕中的偏应力促进了(100)晶面上Si>-II相的转化。2008年，LIN等[35>-36]使用分子动力学纳米压痕方法，研究单晶硅的变形行为和3个晶向上的物相转变规律。借助共近邻分析法、径向分布函数分析法，监测整个过程的物相转变规律。成功地分辨出了针对具有相同近邻原子数的BC8、R8物相；同时获得了工具半径尺寸对单晶硅3个不同晶向上的物相转变规律。当压痕深度与刀尖圆角半径的比值大于0.7时，塑性变形开始发生。

由于计算机建模计算能力和原子间相互作用势函数的限制，先前的纳米切削研究没能揭示在实际未变形切屑厚度基础上的纳米切削脆塑转变规律。2015年，XIAO等[37]借助图形处理器(GPU)开发并行MDS仿真程序，研究金刚石磨粒纳米切削单晶SiC时的脆塑性转变机理。首次使用新型的VASHISHTA[38]多体势函数定义原子间作用力，分析了10～40 nm之间不同未变形切屑厚度下的材料加工、应力分布和断裂位置等特性。结果表明：随着未变形切屑厚度的增加，纳米切削经历了一个从塑性域向混合域的切削过渡，最终向纯脆性域切削模式转变的过程，切削区的拉应力也随着切屑厚度的增加而增大，直至达到临界应力发生脆性断裂，仿真获得的临界切屑厚度值为25 nm。作者同时借助构建单点金刚石>-碳化硅纳米划痕实验系统，利用获得的晶片表面划痕三维形貌与锥形槽的深度，分析纳米加工特性，实验获得的临界切屑厚度值为29 nm。实验验证了仿真分析的结果。

2.1.2 结构及力学特性

20世纪80年代末期，美国的劳伦斯利弗莫尔实验室团队率先应用分子动力学方法，由纳米压痕方法研究晶体材料，分析了晶体材料如铜、镍的显微硬度、屈服强度等力学特性[39]。90年代初期，BELAK等[10]在纳米尺度下研究材料的力学特性，对单晶铜、单晶银等开展了纳米压痕模拟，获得了载荷随压痕变化的曲线；对屈服强度、显微硬度进行了计算，获得了切削力、比能随切削深度变化的规律。在相同时期，IKAWA[40]为研究微划痕最小切削厚度，建立了单晶铜的微划痕加工模拟模型，模拟发现晶体表面未变形层厚度小于刀具刃口半径，同时获得了切削力随刀具切削距离变化的规律。

国内这方面的研究起步较晚，哈尔滨工业大学和天津大学率先开展了针对单晶硅纳米切削加工中的脆塑性转变、切屑和加工表面的机理研究[41>-43]。霍德鸿等[44]研究了单晶铝微构件的拉伸过程，分析了平衡阶段、裂纹产生以及断裂阶段的一些力学规律。段芳莉等[45]模拟了纳米颗粒与单晶硅的碰撞反弹过程，分析了粒子反弹行为、基体的弹塑性形变以及碰撞过程的能量转化过程。结果表明：碰撞中与撞击粒子相邻的基体原子立即非晶化，而非晶层以外的基体可存储弹性变形能，并为颗粒提供反弹的能量，揭示了单晶硅超精密加工的材料去除机理。

进入21世纪以后，随着计算机技术的进步，分子数值模拟技术也得到了快速发展。贾慧灵等[46]对单晶氮化钛纳米杆进行模拟轴向拉伸，通过计算得到[100]和[111] 2个晶向在不同截面尺寸、不同应变拉伸率、不同温度下的拉伸应力应变曲线，获得了屈服强度、弹性模量随以上宏观条件变化的规律，对比分析了2个晶向上的力学性能变化趋势。随着截面尺寸的增大，沿着[100]晶向拉伸，屈服强度逐渐降低；而沿着[111]晶向拉伸，屈服强度逐渐升高。随着拉伸应变率的增大，2个不同晶向上的屈服强度都逐渐增大；但相同应变率下，[100]晶向的屈服强度比[111]晶向的平均要高10.3 GPa左右。随着温度的升高，2个不同晶向的屈服强度均逐渐下降；同样，沿着[100]晶向的屈服强度要比[111]晶向的平均要高9 GPa左右。

沈忱等[47]利用MDS研究HMX晶体超精密加工过程中的力学性能，利用特殊的COMPASS力场，在278～328 K温度范围内研究其力学特性，计算出HMX晶体在2个不同晶面(011)和(010)上的杨氏模量、体积模量、剪切模量以及泊松比等力学参数。随着温度升高，HMX晶体材料的刚度、断裂强度、硬度、韧性等力学特性均呈现出先升高后下降的变化规律。

2.1.3 热力学特性

1991年，丁家强等[48]利用分子模拟技术，研究了晶体的热力学性质。在460～1400 K温度范围内，计算了α铁在不同压力下的熵值。结果表明：单晶体的熵随着温度升高而变大；不同压力下的熵值曲线表明，体系熵值受压力变化的影响很小。单组份物质晶界的表面张力由内能张力增量、熵张力增量、体积张力增量3部分组成。该课题组还计算了单晶铁的温度、压力、内能、熵、Helmholtz自由能等热力学性质[49]。

2004年，MOON等[50]提出了适用于BN晶体的一种改进经验势函数Stillinger-Weber。通过对BN晶体的结构特性、热力学特性进行分子动力学模拟，计算得到德拜温度1700 K以上的热扩散系数、比热等，热力学性质参数的模拟计算结果与理论计算、实验数据基本一致，同时得到了结合能随温度的线性变化函数。与结构性质参数一起验证了这种改进势能函数参数的可行性。

2006年，美国伊利诺伊大学香槟分校课题组[51]研究了拉伸压缩、弯曲、剪切3种应变条件对单晶硅的热力学性质的影响；计算了热扩散系数、自由能、内能、热容以及熵随温度变化曲线。结果表明：在拉伸和压缩应变状态下，单晶硅的热力学性质变化显著，在同一次拉伸和压缩状态下，变化规律不对称；在剪切应变状态下，热力学性质的变化规律呈对称分布，且不受剪切方向的影响。

2.2 单晶材料纳米加工的介质环境影响

前述文献中，纳米加工均假设在理想真空环境下，而在实际纳米加工中，研磨、抛光加工时均存在一定的介质环境。由于介质会传导应力、能量，发生相互作用，因而会对纳米加工过程产生较大的影响。国外很多研究者对介质环境下的纳米加工过程进行了分子动力学模拟，获得了一些研究成果。

2.2.1 大气介质环境的影响

1997年，INAMURA等[52]利用分子动力学模拟方法，对比了在真空、大气2种介质环境中，单晶硅晶体加工时裂纹的扩展过程。结果表明：大气环境中的氧气分子、氢气分子会与刀具尖角处的单晶硅发生应力转化，并引起化学反应。硅原子与气体中分子的化学反应在裂纹前端非常活跃，化学反应减弱了解理表面形成所需的表面能，因而促进了前端裂纹的扩展。2002年，MYLVAGANAM等[53]探索氧气对单晶硅加工过程的影响，如图2模拟了纳米压痕过程。

由图2可知，结果表明：氧气分子的初始位置及方向对氧原子压入硅基体起重要作用，单晶硅在压痕过程中会被氧化。如果压头下方存在氧气分子，那么氧化作用将损伤单晶硅的亚表面层。解释了在磨削、研磨、抛光等纳米磨削加工中，会发生氧化作用并出现亚表面层损伤的原因。

2.2.2 纯水介质环境的影响

2006年，RENTSCH等[54]利用分子动力学模拟技术，详细研究了液态水分子对材料去除过程的影响。图3为切削时间230 ps，切削长度23.0 nm时，含液体和无液体>-真空2种环境介质下的加工图像。比较2种加工介质下的温度和应力分布情况发现：2种介质对加工的影响表现在切屑尺寸和加工表面质量2个方面：含液体介质相对于无液体>-真空介质，其加工切屑更小，且已加工表面更加光滑。表明液体吸收层或反应层在加工表面起到了比较重要的作用。

图3 无液体>-真空介质下的分子模拟加工，230 ps之后(绿色线为含液体加工状态)[54]

2013年，WANG等[55]在金刚石>-单晶铜纳米压痕模拟中引入水环境介质，研究液体对纳米摩擦特性的影响，对比分析了水与真空2种介质条件下纳米压痕的力学特性，分析了载荷随压痕深度的变化、工件内部位错组织变化、显微硬度随压痕深度的变化、工件内部等效应力分布情况等规律。结果发现：相对于真空介质压痕，含水介质压痕在压入初始阶段载荷较大，但在整个压入阶段压力减小；含水介质压痕可以有效减弱刀头与工件间的黏附效应；含水介质压痕表现出正常的压痕尺寸效应；水分子的存在可以显著减小工件与刀具的摩擦力。

2015年，REN等[56]使用MDS研究了含水润滑层的金刚石>-单晶铜AFM纳米划痕加工过程。通过比较工件变形程度、划痕力、摩擦系数等指标，比较了不同水润滑深度、不同划痕深度以及划擦速度下的湿划痕与干划痕加工。仿真结果表明：含水润滑层有利于提高单晶铜划痕表面质量，对划擦力有显著影响，含水润滑划痕加工针尖的摩擦系数要显著大于干划痕过程。

2016年，WEN等[57]建立图4所示的液相环境中的Si>-SiO2截面的摩擦磨损反应模型，利用化学反应分子动力学仿真方法，研究了纳米尺度含水液相环境中Si>-SiO2界面上Si的摩擦化学磨损机理。分析Si原子的2种去除方式，结果表明：通过界面Si-O-Si桥键，2种方式都可以有效去除硅表面原子；分析了压力对化学摩擦磨损的影响，结果证明了界面间Si-O-Si桥键数目随着作用于SiO2相的表面压力的增大而增多，加速了Si原子的去除。揭示了水分子在磨损机制中的2种作用，使得硅基体表面氧化和防止界面过紧接触。

图4 含水环境中单晶硅与二氧化硅界面的化学摩擦磨损模型示意图[57]

2018年，SHI等[58]使用MDS方法研究水薄膜对金刚石纳米划擦单晶铜的影响规律。结果表明：由于水分子的加入，使得金刚石刀头部位的粒子数目增加，导致作用于金刚石磨粒的摩擦力增大，但水分子的润滑作用也减小了磨粒作用在单晶铜表面的摩擦力。去除效率研究表明：Cu原子去除量随着水薄膜层的厚度变化而产生相应的正向变化。

2.2.3 双氧水溶液介质环境的影响

为了探索单晶硅化学机械抛光原理，2017年，WEN等[59]利用化学反应分子动力学方法研究了在H2O2水溶液中，SiO2磨粒划擦加工单晶Si(100)面的工艺过程。结果表明：H2O2水溶液环境中，Si基体表面容易先被氧化，生成Si-O-Si桥键，然后由于SiO2磨粒与Si基体表面的相互划擦作用，使得基体表面的Si-Si、Si-O键受应力而断裂，从而去除Si表面的原子。相对于纯水环境的加工，在H2O2水溶液环境中，Si基体表面更容易氧化，去除Si原子效率更高。

2.3 单晶材料纳米加工后亚表面变形层研究

现今，对晶片表面加工质量及晶片使用寿命的要求越来越高，而晶片已加工表面的下方存在一个相对复杂的弹塑性变形区域，即亚表面变形层[5]。故迫切需要认识亚表面变形层的发生及扩展机理，尽可能地减少加工过程中晶体的亚表面损伤。近年来，国内外很多学者应用分子动力学模拟技术，研究了晶体材料纳米加工后的亚表面变形层特性。

2009年，张俊杰等[5]建立了AFM针尖切削单晶铜材料的分子动力学模型，通过引入原子势能与原子变形程度的匹配判定方法，分别研究了晶体材料的晶向、AFM针尖刀具的切削方向以及刀具切削速度对亚表面变形层的影响规律。结果表明：由于单晶材料的各向异性，亚表面变形层的深度受晶向的变化、刀具切削方向的变化影响很大；沿着给定的切削方向，保持切削深度和切削长度不变，切削速度在20～250 m/s范围内，亚表面变形层的深度几乎不受切削速度变化的影响。

2014年，LI等[61]研究了高速磨削单晶硅的亚表面层损伤特性。图5所示为磨削行程7 nm时获得的单晶硅截面图。通过计算共近邻原子数、分析等效应力和原子位移分布图、计算亚表面损伤层厚度(SSD)等方法，研究了塑性变形区的SSD变化机制；对比了在不同的磨削深度和磨削速度下，单晶硅的亚表面损伤特性。结果表明：随着磨削速度的增高，磨削表面粗糙度下降，亚表面损伤层厚度也降低。当磨削速度超过180 m/s时，由于超高速磨削导致磨削力和温度上升，SSD略有增加。

2016年，CHAVOSHI等[62]在单晶硅塑性加工研究中，首次在MDS中引入原子位移向量指标来评价纳米切削中材料的塑性流动特性。图6为在单晶硅切削的3个晶面的原子位移场分布中，(111)晶面出现了漩涡，表明其加工性能优于(110)、(010)晶面。同时加工区域内湍流度的变化与加工温度的增加成线性关系。

2017年，郭晓光等[63]对单晶硅多次加工时的亚表面变形损伤问题进行了研究，通过模拟对比理想无缺陷单晶硅、切削加工后有缺陷单晶硅的纳米划痕过程，探索亚表面变形层的特性。设定3种不同的划痕深度，对这2种结构的单晶硅分别进行纳米划痕模拟，然后借助共近邻分析法、配位数法来计算出亚表面变形层的厚度。结果表明；有缺陷的单晶硅再次加工时，工件表面材料去除更加容易，同时亚表面损伤层深度减少，当已有的损伤层厚度不足再次加工深度的2倍时，加工后的损伤最小。

I-已加工区，II-后刀面以下区，III-后刀面与刀尖过渡区，IV-刀尖下方区域，V-待加工区。(红色箭头标线代表停滞区域附近的材料流动方向。上方的局部放大图中θs1、θs2，h1、h2分别代表停滞区角度的边界。)

3 结语

分子动力学模拟已经同理论研究、实验研究并列，成为纳米加工领域的一项重要研究方法。本文简要回顾了经典分子动力学模拟技术的基本原理与方法，综述了分子动力学模拟方法在单晶材料纳米加工机理方面的研究成果。但还存在一些需要关注的问题：

(1)现有模拟研究大都针对的是理想无缺陷晶体材料和刚性刀具，但实际晶体均存在一定的缺陷，金刚石纳米加工刀具本身也存在磨损。

(2)分子动力学模拟局限于从皮秒到微秒的时间尺度、纳米的空间尺度，模拟的时间、空间尺度范围比较小。

(3)分子动力学模拟提供了一种预测方法，仍需要协同开展纳米加工实验研究，以相互验证。