晶体铜微探针纳米刻划的分子动力学建模

张俊杰 孙 涛 闫永达 梁迎春 董 申

哈尔滨工业大学，哈尔滨，150001

0 引言

纳米机械加工技术如金刚石刀具切削、金刚石磨粒加工以及金刚石微探针纳米刻划等的加工精度为纳米量级，使加工形成的结构与器件的尺寸达到纳米量级，是一种重要的由上而下的纳米加工技术［1－2］。然而，缺乏对纳米机械加工机理的理解极大地制约着纳米机械加工技术的发展［3］。分子动力学（molecular dynamics，MD）模拟已被证实为研究纳米机械加工机理的一个重要工具，国内外学者对纳米机械加工的分子动力学模拟开展了大量的研究工作，并取得了很多有价值的研究结果［4－9］。然而，在他们的研究中，工件材料大多为单晶形式，对多晶材料纳米机械加工的分子动力学模拟研究较少。多晶材料发生塑性变形时，晶界不仅会阻碍位错运动，并且晶界自身也会发生变形，其变形机理远比单晶材料复杂［10］。此外，前人研究主要依据宏观统计量如切削力来分析分子动力学模拟的结果，而对工件材料内部缺陷演变的分析较少。缺陷的形核及运动是材料塑性变形的一种重要方式，并且工件材料的变形与宏观统计量直接对应。因此，为了全面、深入地理解纳米机械加工机理，需要使用分子动力学模拟并结合晶体缺陷分析技术来研究多种结构材料的纳米机械加工机理。

纳米机械加工的分子动力学模拟主要包含建立分子动力学模型、进行分子动力学模拟和处理分子动力学模拟数据。其中，建立包含具备基本物理特性的精确分子动力学模型决定着分子动力学模拟结果的正确性及精度。因此，本文针对上述不足之处，以晶体铜微探针纳米刻划为例，分别从分子动力学基本原理、晶体铜原子结构建模和晶体缺陷分析这三个方面开展了晶体铜微探针纳米刻划的分子动力学建模的研究。

1 分子动力学基本原理

图1为晶体铜微探针纳米刻划的分子动力学模型（包含晶体铜工件和金刚石微探针）示意图。晶体铜工件由具有面心立方（face centered cubic，FCC）晶格结构的铜原子组成，金刚石刀具由具有金刚石晶格结构的碳原子组成。工件在X方向、Z方向具有周期边界条件，在Y方向具有自由边界条件。同时，工件底部设置固定层原子以防止加工中工件的整体刚性移动。纳米刻划过程包含驰豫阶段和加工阶段。在加工阶段中，微探针首先从工件表面上方以恒定的速度压入工件内部直到指定的压痕深度，然后沿着X轴负方向以恒定的速度进行刻划直到指定的刻划长度。

图1 晶体铜微探针纳米刻划的分子动力学模型示意图

在建立晶体铜纳米机械加工的分子动力学模型时，需要从分子动力学基本原理角度出发设置经典牛顿运动方程的求解方法，描述模拟体系的势能函数，模拟真实环境的系综（ensemble）等参数。虽然分子动力学模拟已被证实为研究纳米机械加工机理的一个有力工具，但是分子动力学模拟结果的正确性与精度受到分子动力学模型精度的影响。在本文建立的分子动力学模型中，采用velocity－Verlet法则来求解当前分子动力学模型中原子的位置、速度和加速度，采用 Mishin等［11］提出的EAM势能函数的参数来描述晶体铜工件材料里的Cu－Cu原子相互作用，采用Morse势能函数来描述晶体铜工件与金刚石微探针之间的Cu－C原子相互作用［12］。Morse势能函数具有较大的结合能，导致工件表面与刀具之间存在一个很大的吸附力。在纳米机械加工实验中，由于工件表面很难做到原子级洁净或者存在硫醇钝化层，导致工件表面与刀具之间的吸附力显著减小。因此，为了模拟钝化层的存在对工件与刀具之间作用的影响，可以采用强排斥势能函数来模拟刚性球形金刚石微探针，以消除吸附力对微探针与工件之间作用的影响。强排斥势能函数的表达式为

式中，r为球形金刚石微探针中心与工件表面之间的距离；R为微探针的半径；d为微探针的硬度［13］。

在晶体铜微探针纳米刻划的分子动力学模拟中，需要分别对驰豫阶段和加工阶段设定相应的系综。在驰豫阶段，采用等温等压（isothermalisobaric）系综对模拟体系进行约束，要求模拟体系的温度达到设定的温度，系统压力等于或接近于0。图2所示为使用NPT系综对单晶铜微探针纳米刻划的分子动力学模拟体系在目标温度30K，目标压力为0，下驰豫40×10－12s的过程中的温度、压力变化曲线。可以看到，在初期振荡后，温度和压力均稳定在目标值，表明模拟体系驰豫完全。在随后的加工阶段中，由于工件材料的去除深度在数个或数十个原子层，因此只有很小一部分工件材料参与了加工过程，从宏观上看，纳米机械加工中工件材料体积的变化极小。在纳米尺度下，热量对材料变形的影响显著增加，有必要研究纳米机械加工中热量的变化。因此，在加工阶段可以使用微正则（microcanonical）系综对模拟体系进行约束，要求模拟体系的体积和能量保持不变。

图2 单晶铜微探针纳米刻划的分子动力学模拟体系的NPT驰豫过程

2 晶体铜建模

一个理想晶体是由晶胞在三维空间有序、周期排列而形成的。由于原子排列的周期性和疏密程度在不同的晶向上有显著的差别，因此单晶铜的导电性、导热性、弹性模量、屈服强度等物理性质均存在各向异性。纳米机械加工是一个刀具与材料高度耦合的过程，单晶铜的各向异性直接影响工件材料的塑性变形及加工结果，因此可以通过创建具有不同晶面的单晶铜工件来研究单晶铜的各向异性对单晶铜微探针纳米刻划的影响。

多晶材料发生塑性变形时，晶界阻碍位错运动，因此揭示位错－晶界交互作用是研究多晶铜微探针纳米刻划加工机理的一个重要方面。为了排除晶界变形对工件塑性变形的影响，可以通过双晶铜微探针纳米刻划的分子动力学模拟来揭示位错－晶界交互作用。双晶铜包含两个具有不同晶体取向的晶粒以及连接它们的晶界。双晶铜的原子模型可以根据重位点阵（coincidence site lattice，CSL）来创建：首先建立两个相邻的单晶铜晶粒，使它们具有相同的晶体取向。然后保持一个晶粒不动，让另外一个晶粒绕着某一旋转轴旋转一定的度数。由于发生旋转的晶粒的晶体取向不同于未发生旋转的晶粒，因此两个晶粒存在位相差，它们的晶格失配在晶粒之间形成晶界［14］。由CSL可知，不同的旋转轴和/或旋转度数得到的晶界结构也不同。按照相邻晶粒之间的位相差的大小，晶界可以分为小角度晶界（相邻晶粒的位相差小于10°）和大角度晶界（相邻晶粒的位相差大于10°）。按照旋转轴与晶界的关系，小角度晶界可以分为倾斜晶界（旋转轴平行于晶界）和扭转晶界（旋转轴垂直于晶界）。还有一类特殊的孪晶界（twin boundary，TB），晶界两侧原子的对称排列。图3所示为根据CSL创建的包含不同晶界结构的双晶铜工件的原子结构模型。

图3 基于CSL创建的双晶铜工件

多晶材料包含许多具有不同晶体取向的单晶晶粒，晶粒之间通过晶界连接。当晶粒尺寸进入到纳米尺度后，由于纳米晶体（nanocrystalline）材料的晶界体积分数（grain boundary volume fraction）显著增大，因此纳米晶体材料具有高强度、超塑性等独特的机械性能［15］。纳米尺度下，晶界对纳米晶体材料在纳米机械加工中的塑性变形及加工结果具有更为显著的影响。因此，有必要创建纳米晶体铜的原子结构模型，为开展纳米晶体铜在纳米机械加工的分子动力学模拟奠定基础。目前，常用维诺图（Voronoi diagram）来建立纳米晶体材料的原子结构模型。设定P＝｛p1，p2，…，pn｝是一组随机分布在一个平面或三维空间上的n个离散点，它将平面或三维空间分为n个Voronoi区域V（pi）。每个包含点pi的V（pi）是到pi距离最近点的集合，每个V（pi）对应一个晶粒［16］。然后在每个晶粒内根据指定或随机的晶体取向填充FCC铜原子。不同晶体取向的晶粒存在晶格失配，在相邻晶粒之间会形成晶界。当点pi分布在与笛卡儿坐标系某一坐标轴垂直的平面上时，得到的纳米晶体铜具有柱状（columnar）结构。图4所示为基于维诺图建立的纳米晶体铜的原子结构模型，其中白色原子为晶界原子，灰色原子为FCC原子。

图4 基于维诺图建立的纳米晶体铜的原子结构模型

3 晶体缺陷分析

受到晶体的生成条件、原子的热运动、对晶体进行的加工过程、掺杂及其他条件的影响，实际晶体中往往存在偏离完整性点阵结构的晶体缺陷。晶体缺陷可以分为点缺陷（替位、间隙、空位）、线缺陷（位错）和面缺陷（堆垛层错和界面）等。晶体缺陷是随着各种条件的改变而不断演变的，如材料内部的位错起源与增殖机制在外加应力作用下被激活，导致材料中的位错密度随着塑性变形的进行而增加。缺陷的形核及运动是金属材料塑性变形的一种重要形式，通过晶体缺陷分析技术来精确辨别工件材料内部形成的晶体缺陷的种类和位置，是晶体铜纳米机械加工的分子动力学模拟的一个关键之处。目前广泛使用的晶体缺陷分析技术 有 径 向 分 布 函 数 （radial distribution function，RDF）［17］、中心对称参数（centro－symmetry parameter，CSP）［18］、共近邻原子分析（common neighbor analysis，CNA）［19］、滑移向量分析（slip vector analysis，SVA）［20］、Ackland 键 角 分 布（bond angular distribution，BAD）分析［21］等。其中，共近邻分析是一种广泛应用于金属材料的晶体缺陷分析技术。

在共近邻分析中，假定有一个原子对包含原子α和β。对于一个给定的理想晶体，其结构可以用指数i、j、k、l组成的CNA图表来辨别［19］：

（1）指数i。若原子α和原子β为近邻原子，则i＝1；否则，i＝2。

（2）指数j。原子α和原子β共同享有的近邻原子（共近邻原子）的个数为j。

（3）指数k。原子α和原子β的共近邻原子之间键的个数为k。

（4）指数l。当i、j、k相同时，原子α和原子β的共近邻原子之间不同连接方式的个数为l。

如图5所示，FCC结构里所有的共近邻原子具有1421图表形式。在密排六方堆积（hexagonal close packing，HCP）结构里，共近邻原子具有1421和1422两种图表形式。体心立方（body centered cubic，BCC）结构的图表分布不稳定，3/7的共近邻原子具有1441图表形式，4/7的共近邻原子具有1661图表形式。图6为单晶铜挤压后材料内部的缺陷结构的分子动力学模拟图，原子根据计算得到的CNA值进行着色，其中，白色原子为处于表面或位错等缺陷环境的Other原子，深灰色原子为处于堆垛层错等缺陷环境HCP原子，灰色原子为处于完整晶格结构的FCC原子。为了更清楚地显示缺陷结构，可以将未发生变形的FCC消除不显示，如图6b所示。可以看到，共近邻分析可以清楚地辨别单晶铜工件中存在的缺陷结构，极大地促进了对晶体铜纳米机械加工机理的深入理解。

图5 晶体结构的CNA图表示意图［22］

图6 单晶铜挤压后工件材料内部的缺陷结构

4 结语

本文从分子动力学基本原理、晶体铜原子结构建模和晶体缺陷分析技术三方面介绍了晶体铜微探针纳米刻划的分子动力学建模。本文研究结果表明：通过合理地选择经典牛顿运动方程的积分方法，精确描述体系的势能函数，模拟真实环境的系综，可以获得到晶体铜纳米机械加工的分子动力学模拟的精确原子模型。基于晶体结构、重位点阵和维诺图等建立单晶铜、双晶铜和纳米晶体铜的原子结构模型，同时结合先进的晶体缺陷分析技术辨别纳米机械加工过程中工件内部形成的缺陷的种类和位置，可以为系统地开展晶体铜纳米机械加工的分子动力学模拟提供技术保障。

［1］Zhao X Z，Bhushan B.Material Removal Mechanisms of Single－crystal Silicon on Nanoscale and at Ultralow Loads［J］.Wear，1998，223（1/2）：66－78.

［2］Yan Y D，Sun T，Dong S，et al.Top－down Nanomechanical Machining of Three－ dimensional Nanostructures by Atomic Force Microscopy［J］.Small，2010，6（6）：724－728.

［3］陈日曜.金属切削原理［M］.2版.北京：机械工业出版社，1993.

［4］Belak J，Lucca D A，Komanduri R，et al.Molecular Dynamics Simulation of the Chip Formation Process in Single Crystal Copper and Comparison with Experimental Data［C］//Proceedings of ASPE Annual Conference.New York，1991：100－109.

［5］Fang T H，Weng C I.Molecular Dynamics Simulations of Nano－lithography Process Using Atomic Force Microscopy［J］.Surface Science，2002（1/2）：138－147.

［6］Fang F Z，Wu H，Liu Y C.Modelling and Experimental Investigation on Nanometric Cutting of Monocrystalline Silicon［J］.International Journal of Machine Tools ＆ Manufacture，2005，45（15）：1681－1686.

［7］Mulliah D，Kenny S D，Smith R，et al.Molecular Dynamic Simulations of Nanoscratching of Silver（100）［J］.Nanotechnology，2004，15（3）：243－249.

［8］罗熙淳，梁迎春，董申.单晶铝纳米切削过程分子动力学模拟技术研究［J］.中国机械工程，2000，11（8）：860－863.

［9］郭晓光，郭东明，康仁科，等.单晶硅纳米级磨削过程的理论研究［J］.中国机械工程，2008，19（23）：2847－2851.

［10］Jin M，Minor A M，Stach E A，et al.Direct Observation of Deformation－induced Grain Growth during the Nanoindentation of Ultrafine－grained Al at Room Temperature［J］.Acta Materialia，2004，52（18）：5381－5387.

［11］Mishin Y，Mehl M J，Papaconstantopoulos D A，et al.Structural Stability and Lattice Defects in Copper：Ab Initio，Tight－binding，and Embedded－Atom Calculations［J］.Physical Review B，2001，63（22）：224106.

［12］Zhu P Z，Hu Y Z，Ma T B，et al.Study of AFM－based Nanometric Cutting Process Using Molecular Dynamics［J］.Applied Surface Science，2010，256（13）：7160－7165.

［13］Ziegenhain G，Hartmaier A，Urbassek H M.Pair vs Many－body Potentials：Influence on Elastic and Plastic Behavior in Nanoindentation of FCC Metals［J］.Journal of the Mechanics and Physics of Solids，2009，57（9）：1514－1526.

［14］Tschopp M A，McDowell D L.Structures and Energies of 3Asymmetric Tilt Grain Boundaries in Copper and Aluminium［J］.Philosophical Magazine，2007，87（22）：3147－3173.

［15］Lu L，Sui M L，Lu K.Superplastic Extensibility of Nanocrystalline Copper at Room Temperature［J］.Science，2000，287（5457）：1463－1466.

［16］Du D S，Hwang F.Computing in Euclidean Geometry［M］.Singapore：World Scientific Publishing，1995.

［17］Matteoli E，Mansoori G A.A Simple Expression for Radial Distribution Functions of Pure Fluids and Mixtures［J］.Journal of Chemical Physics，1995，103（11）：4672－4677.

［18］Kelchner C L，Plimpton S J，Hamilton J C.Dislocation Nucleation and Defect Structure during Surface Indentation［J］.Physical Review B，1998，58：11085－11088.

［19］Honeycutt J D，Andersen H C.Molecular Dynamics Study of Melting and Freezing of Small Lennard－jones Clusters［J］.Journal of Physical Chemistry，1987，91（19）：4950－4963.

［20］Zimmerman J A，Kelchner C L，Klein P A，et al.Surface Step Effects on Nanoindentation［J］.Physical Review Letters，2001，87（16）：165507.

［21］Ackland G J，Jones A P.Applications of Local Crystal Structure Measures in Experiment and Simulation［J］.Physical Review B，2006，73：054104.

［22］Tsuzuki H，Branicio P S，Rino J P.Structural Characterization of Deformed Crystals by Analysis of Common Atomic Neighborhood［J］.Computer Physics Communications，2007，177（6）：518－523.