基于分子动力学模拟的纳米压痕技术研究进展

张勤俭，李正钊，曹建国

（ 北京信息科技大学机电工程学院，北京 100192 ）

随着全球科学技术的迅猛发展，智能制造已成为大国博弈背景下建设制造强国的主攻方向。 在智能制造领域，超精密加工是重中之重。 先进的超精密加工技术是制造业高速发展的重要条件。 但超精密加工过程中的材料研究存在难点，人们无法通过实验手段来探究晶圆抛光、激光加工和超声加工等材料去除过程中的材料性能变化规律。 结合分子动力学（molecular dynamics，MD） 模拟与纳米压痕技术，将能有效地解决上述问题，进而可加速超精密加工技术的发展。 MD 模拟技术是连接理论研究与实验研究的桥梁，对于理论研究来说，MD 模拟技术可利用理论得出的数据进行复杂的数值求解；对于实验研究来说，MD 模拟技术可在原子层面展示出体系微观的结构变化。 借助计算机技术的发展，研究纳米级粒子的MD 模拟技术也逐渐完善起来。MD模拟技术的应用范围广泛，为药物设计[1-2]、材料[3-4]、化学工程[5]、凝聚态物理[6]、环境科学[7-8]等领域提供了诸多重要的理论依据，起到举足轻重的作用。

改善和应用纳米材料时，需先对该材料的基本结构有所了解。 MD 模拟技术不仅可以形象地展示材料分子的基本结构， 而且还能展示材料随着温度、压入载荷、压入深度等变化而变形的过程。 借助纳米压痕技术则可获得纳米尺度上的材料力学性能参数，可改善材料性能，比如研究合金表面氧化膜的硬度和弹性模量，有利于提高合金的耐腐蚀能力；残余压应力可以阻止材料裂纹发展，不同喷丸工艺在增强残余压应力方面的效果不同，经喷丸处理的传动机构、 大型齿轮可拥有优良力学性能，能广泛应用于航天、汽车等领域。 纳米压痕技术弥补了传统硬度测量方法在微观领域的空白， 使精密、超精密加工技术得到进一步发展，为探究微观体系的力学性能提供一种新的微观的测试技术。 因此，本文总结了分子动力学模拟和纳米压痕技术的基本原理，对采用分子动力学模拟研究纳米压痕技术的现状进行了综述。

1 分子动力学模拟技术简介

MD 模拟技术飞速发展， 通常用于研究材料在原子层面上微观变化，是利用计算机技术对微观粒子的运动状态进行模拟，经过数据分析后可了解材料在结构和微观的相互作用方面的性质，可确定材料的化学、力学和结构特性，并且可从分子层面对实验结果进行验证。 MD 模拟充当了微观长度和时间尺度之间的桥梁，弥补了实验的局限性，并可成功模拟极端温度[9]和极端压力[10]条件下的实验，也可模拟需要昂贵设备和材料的实验，可节省研究时间及经费。

1.1 分子动力学模拟的理论基础

分子动力学描述的是整个系统，需满足所有微观粒子的运动都遵循牛顿运动方程以及粒子间的相互作用遵循叠加原理这两个假设。 MD 模拟需三个步骤：①创建系统的初始状态，确定粒子的初始位置和类型；②选取粒子之间的力场，即相互作用势；③借助牛顿运动理论，进行数值求解，进而预测粒子的运动状态。

MD 模拟的基本流程[11]如图1 所示，已知粒子的初始状态，先建立分子体系模型，设置温度、粒子数、时间等条件参数；再计算系统中每个粒子的受力和加速度， 运用牛顿方程求出运动粒子的新位置，获得此时的速度和加速度，重复计算直到体系达到平衡状态；后分析体系中粒子的相轨迹。

图1 分子动力学模拟流程图

1.2 粒子间相互作用势

势函数在模拟过程中起关键作用, 选择不同势函数得到的模拟结果大相径庭。 不同的原子需设置不同的势函数，由于开展纳米压痕实验的材料通常为复合材料[12]，多采用混合力场设置，即利用两个及两个以上的势函数。 常见的势函数类型如下：

（1）Lennard-Jones 势函数

Lennard-Jones 势函数也称L-J 势，是最常见的非键结势能，多用于模拟惰性气体分子。

式中：ε 是最小势能处的能量；σ 是距离，r=σ 是体系势能为0 的位置，rc为截断半径。

（2）Morse 势函数

Morse 势与L-J 势类似，用于描述没有键连接的粒子间的作用力。

式中：D0是结合能系数；α 是势能曲线梯度系数；r0是原子间距。

（3）Tersoff 势函数

Tersoff 势函数[13-14]可描述半导体和金属等材料的结构和性质，适用于分子动力学化学计算、纳米级材料建模的研究。

式中：Vij是相互作用势，i、j 表示两个相邻的原子。

1.3 数值求解运动方程

牛顿运动方程描述了相空间中粒子的速度和位置随着时间变化的情况。 根据牛顿经典力学可知，该体系中某粒子i 的受力为总势能梯度，即受力公式为：

根据牛顿第二定律可知，此刻该粒子的加速度表达为：

式中：mi为粒子的质量。 将式（5）进行积分，得到经过驻t 后的该粒子速度和位置公式分别为：

式中：vi是某粒子i 的速度，xi是某粒子的位置，下标带有“0”的为初始值。

分子动力学方法在计算数学领域是运用有限差分法求解初值问题。 在实际体系中，用于描述相互作用关系的势函数比较复杂，牛顿运动方程的解析值难以求解，因此采用数值积分方法求解。 常见的数值求解算法有Verlet、Velocity-Verlet、Leapfrog、Beeman、Gear 及Rahman 等。

2 纳米压痕技术的基本原理和实验设备

纳米压痕技术是由DOERNER 和NIX 提出[15]，后由OLIVER 和PHARR 进行改进[16-17]。纳米压痕技术的最大优势是在压头几何形状已知的情况下，可以通过载荷-位移曲线（P-h 曲线）精确确定接触面积A。 这一优势使纳米压痕技术适合于测量小尺度材料[18]的机械性能，尤其适用于薄膜、涂层等材料。

2.1 纳米压痕技术的基本原理

Oliver-Pharr 方法（O-P 方法）[16]是研究纳米压痕最常用的分析方法。 在加载-保荷-卸载的过程中，同步记录的压痕载荷-位移曲线（P-h）见图2，建立的卸载过程前后的压痕轮廓模型见图3。 当压头压入被测试样时，被压面先产生弹性变形；随着压头继续向下移动，被压面开始产生塑性变形；卸载后，弹性变形会有一定的回复，而塑性变形无法回复，则形成压痕。

图2 载荷-位移曲线示意图[18]

图3 压痕横截面示意图[18]

在加载过程中，材料不仅发生了弹性形变还发生了塑性形变， 因此弹性接触理论不适用此阶段。在卸载过程中，由于塑性变形不可逆，压痕深度的缩小可归结于材料的弹性回复。 即P-h 曲线的卸载过程是一个纯弹性变形的过程，可采用弹性接触理论进行分析。

对于加载过程的载荷-位移关系，可表示为：

式中：P 为载荷；h 为压入深度；A 是由材料机械性能和压头几何形状所确定的常数。

对于卸载过程的载荷-位移关系， 此时压痕处存在塑性变形，因此对幂函数进行深度修正，则有：

式中：∂、m 为常数；hf为完全卸载处的最后深度。

在最大负载接触深度处，对式（9）进行微分，得到此处曲线的斜率即接触刚度S，如下：

由式（10）可见，接触刚度的精度取决于位移的精度。 测试实验是纳米级，即使是很微小的干扰如实验环境温度的变化、试验机机架的变形、试样的非弹性变形等都会引起巨大的实验结果差异，因此实验的重复性对于纳米压痕测试至关重要。

根据sneddon 公式[19]可知：

式中：hs为压头与试样接触周边处材料表面的位移量；Pmax为最大载荷；ε 与压头几何形状有关。数据分析中采用的压头参数值见表1。

表1 数据分析中不同几何形状压头的参数

利用式（11）可求出最大载荷时的接触深度为：

压痕面积A 是一个与接触深度hc有关的函数。式（13）中，第一项描述了一个完美的玻氏压头的接触面积，此时认为压头的尖端无任何曲率，而其他项的存在是为修正因尖端处钝化而导致的偏差，从而获得这一符合实际情况的表达公式。

式中：C1~C5为常数。

通过上述公式可知，被测材料的硬度为：

从卸载曲线的斜率中提取材料弹性模量为：

式中：β 是常数， 取决于压头几何形状；Er为约化弹性模量；E、v 分别为被测试样的弹性模量和泊松比，Ei、vi分别为压头材料的弹性模量和泊松比。

2.2 纳米压痕仪

近年来，基于纳米压痕技术的各种纳米压痕装置相继被研制出来，为测量材料的微观力学特性提供了更加高效、便捷的方法。 纳米压痕仪以O-P 模型为理论基础，通过设定压入载荷或压入深度这两种方式进行控制。 不同仪器在力和位移系统的范围以及数字分辨率方面的能力不同[20]，可探测多种有机材质、无机材质或软、硬度不同材质，广泛用于测量半导体材料、陶瓷材料、电池电极材料等。

纳米压痕实验装置[21]如图4 所示，纳米压痕仪的核心组成部分为机械系统、传感器系统、控制系统、数据采集和分析系统。 不同厂家生产的纳米压痕仪会有一定的差异。

图4 纳米压痕实验装置示意图[21]

在压痕实验过程中， 压头的位移由压痕仪控制。 压头根据设置的载荷逐渐接触试件，在30 s 时载荷达到最大值。 为得到更高的测试精度，在载荷最大值条件下进行保荷处理，以减少材料蠕变效应带来的实验误差。如图5 所示，在10 s 时，压头的载荷减小至最大载荷的10%，保持该状态60 s，以测量仪器的热漂移，随后进行卸载。 对于绝大多数材料来说，按照这样的时间周期进行测试，所得测试结果合理。 需要注意的是，选择最大载荷的保持时间，应重点考虑蠕变效应带来的影响。

图5 典型试验力与压痕试验时间周期[20]

3 基于MD 模拟的纳米压痕技术研究现状

3.1 测量材料表面的蠕变行为

材料的蠕变行为是指在低于屈服强度的应力作用下，材料的塑性形变随着时间的增加而增大的现象。 这一现象通常发生于金属、陶瓷、塑料和混凝土等复合材料，并且在任何温度下都会发生，尤其当约比温度>0.3 时更加明显。 升高温度[22-23]、增大压力[24-25]、细化晶粒尺寸[26]以及增加辐照[27-28]都会促进材料的蠕变过程。 其中，材料的蠕变行为受温度和压力影响较大，其次是晶格尺寸的影响。 升高温度和增加应力使晶界宽化，而较小的晶粒尺寸会加速蠕变过程， 辐照处理会使体系中产生大量空位团簇、层错四面体以及位错环等缺陷，进而影响材料的宏观性能，甚至减少使用寿命，缩短机器的服役时间。 飞秒激光加工技术[29-30]是利用辐射对靶表面进行材料去除。 这是一种非常精细的纳米级工艺，借助MD 模拟技术建立模型，为研究该工艺的复杂加工过程提供了新的可能。

MD 模拟是一种理解蠕变行为以及研究蠕变特性如何受温度、应力、晶粒尺寸等因素影响而变化的有利研究方法。 Zhao 等[31]采用了MD 模拟方法研究了不同晶粒尺寸的TiAl 合金模型（图6）。 在不同温度（923、1023、1123 K）条件下施加不同应力（1.0、1.4、1.8 GPa），基于获得的数据绘制的时间-蠕变应变曲线见图7。 可知，蠕变行为一般分为初始蠕变、加速蠕变、稳态蠕变三个阶段。 该研究发现，在同一温度条件下，应力越大，稳态蠕变速率越大、位错密度越高，同时表明在稳定蠕变阶段，受高应力作用的蠕变机制主要为位错运动机制，而在低应力作用下的蠕变机制主要为扩散蠕变机制。

图6 纳米晶TiAl 模型具有不同的晶粒尺寸[31]

图7 纳米晶TiAl 模型受温度、晶粒尺寸、应力影响的蠕变应变-时间曲线[31]

镍基单晶高温合金[32]是航空发动机涡轮叶片的关键材料。 张震男等[33]使用纳米压痕技术研究不同负载速率和外加峰值载荷对DD407 镍基高温合金蠕变行为的影响。 该研究表明：DD407 镍基高温合金的抗蠕变能力较好，但对加载速率和峰值载荷比较敏感。 此外，材料的硬度[34]也易受加载速率的影响，硬度随着加载速率的增大而降低。

改善材料的蠕变性能，对于航天、核工业、金属加工、 矿业工程和路面运输等领域具有重要意义。许多研究者通过在合金或多层膜系统中添加新物质，能有效地提高材料的抗蠕变性能。 Wu 等[35]研究了添加稀土元素Re 的镍基单晶高温合金的蠕变行为，如图8 所示，添加了Re 的合金的蠕变应变比未添加的要低。Haseeb 等[36]在多层膜系统加入超薄Ni层（70 nm），显著改善Sn-Cu 基金属间化合物的脆性。 李元等[37]对添加了石墨烯的复合钎料在高温下的蠕变行为进行了研究，发现在相同温度下，添加了石墨烯的复合钎料的蠕变深度比未添加的小。Chang 等[38]研究了纳米颗粒填充率对聚合物纳米复合材料（PNC）蠕变性能的影响，结果表明填充率越高，PNC 的蠕变柔量越小。 研究蠕变行为也是探索纳米微观结构与力学性能关系的重要切入点。 即使化学成分相同，但微观结构的差异[39]也会影响材料的蠕变变形行为。 为设计出具有优良机械性能的纳米晶粒结构材料，Zhou 等[40]通过研究发现，充氢可抑制纳米晶a-Fe 的蠕变行为， 但随着晶粒尺寸增加，抑制作用会转变为促进变形的作用。 在纳米压痕蠕变研究领域常用的连续刚度测量法[41]，是在较长时间里通过施加一定频率来得到被测材料的蠕变和应力松弛等数据结果。 该方法能直接得到材料的硬度和弹性模量随着压力深度增加而连续改变的曲线且对漂移不敏感， 较适合测量薄膜材料，也应用于研究块状材料[42]。

图8 在1 000 K 恒温、不同外加应力、有无Re 条件下的试样蠕变应变-时间曲线[35]

3.2 测量材料的硬度和弹性模量

研究薄膜和涂层等材料的微观力学性能，常采用Oliver-Pharr 方法， 并借助纳米压痕仪以及原子力显微镜[43]等实验设备，主要对材料进行三维形貌观测、成分分析等。 利用纳米压痕技术探究薄膜的力学性质时，通常要求最大压痕深度小于薄膜厚度的1/10，因而可避免基体对薄膜的力学性能产生影响。 Zak 等[44]通过纳米压痕和数值模拟研究了硅衬底上两种已知的薄膜材料体系，进一步证实了硬度测量的“10%规则”。

借助纳米压痕技术和纳米压痕测试系统，赵军等[45]测量了经过深冷处理后的铝青铜的硬度和弹性模量，发现增加压入深度会使材料硬度降低，而弹性模量不存在这种尺寸效应。 材料的硬度存在尺寸效应，原因是压头不断向下压入试样表面，使压入被测试样的压头的宽度增大，被测试样受压处的位错形核速度下降，位错密度减小；同时，随着压入深度增加， 材料的弹性模量无明显上升和下降的趋势，而是呈现波动现象。 为研究超声振动对材料微观结构和力学性能的影响，Yang 等[46]建立了压痕分子动力学模型和超声振动辅助压痕分子动力学模型，发现超声波振动会降低铝的弹性模量，并探究了声学软化现象产生的主要原因和发生在金属中的普遍性。 在金属基体弹性变形阶段，超声振动不会影响其势能。 在塑性变形阶段，金属基体产生位错吸收超声波振动的能量， 促进滑移现象产生，即发生声学软化。

王静怡等[47]对Cu-Zn 合金进行了室温扭转实验，通过剧烈的塑性变形，基于O-P 模型对该合金的力学性能进行了分析。 该研究表明：扭转圈数越大，材料的塑性变形越大、位错密度增加，并且伴随出现新的晶界和亚晶界；通过增加载荷，显微硬度和杨氏模量呈现相同的减小趋势，并且扭转圈数的增加不会对显微硬度和杨氏模量有明显影响。

纳米压痕能预测脆硬性材料[48]的纳米硬度和弹性模量， 注入离子会使材料的弹性模量和硬度降低，在对被注入离子的脆性材料进行先进制造的过程中，能有效地降低刀具的磨损以及提高机械加工性能。Wu 等[49]研究了注入不同剂量离子（表2）对单晶立方碳化硅（3C-SiC）机械性能的影响，发现随着注入量的增加，硬度和弹性模量越小，位错形成受到抑制，碳化硅材料的塑性增强（图9）。

表2 注入剂量表

图9 硬度和弹性模量变化图[49]

在切削过程中， 刀具的磨损与工件的硬度有关， 通过激光辅助加工 （laser-assisted machining，LAM）可得到更高质量的工件表面，进而能减轻刀具磨损。Liu 等[50]对LAM 过程进行研究，设置了不同激光脉冲功率（图10a），发现随着激光功率增加，具有配位数5 和6 的原子的数量降低，与配位数5 和6 对应的BCT5-Si 相和β-Si（Si-II）相结构的形成受到抑制，从而减少加工表面损伤；同时，切屑原子的数量随着激光功率的增加而增多（图10b），说明激光可降低材料硬度， 最终获得更高的材料去除率。AlCoCrFeNiCu 高熵合金可用作A301 钢的涂层，用于航空航天结构应用中的机械失效防护。Dada 等[51]通过纳米压痕技术测试了高熵合金与A301 钢基体的力学性能，并进行了比较，结果发现激光沉积的高熵合金涂层比钢基体更耐塑性变形。 NiTi 形状记忆合金多用于制作人造关节，在医疗领域得到广泛应用，由于被植入人体，长期处于人体体液环境，人们迫切需要对其特性进行了解并改善该合金表面的性能。 汪利斌等[52]探究了温度和时间对NiTi 形状记忆合金表面氧化膜的影响，发现氧化膜的纳米硬度和弹性模量会随着温度的升高先降低，在温度为600 ℃时均开始增加， 在达到峰值700 ℃之后继续降低，由此最终确定700 ℃是制备性能优良的NiTi形状记忆合金的最适宜温度。

图10 不同激光功率对单晶硅亚表面损伤的影响（T=140 ps）[50]

3.3 测量材料表面的残余应力

经过各种机械加工工艺的零部件，其内部都会产生不同程度的残余应力。 在薄膜和涂层材料中也不可避免地存在残余应力。 材料的疲劳性能[53]和力学特能[54]易受残余应力的影响，并在其失效行为中起到至关重要的作用。 因此，提高测量残余应力的准确性， 具有十分重要的研究意义和工程应用价值。近年来，残余应力的测量方法[55]逐渐从破坏性测量到无损性测量， 但每种方法都有其适用的领域，研究者应结合实际情况选取最适合的方法。

通常，借助纳米压痕技术确定表面残余应力的方法有两种。 第一种方法，基于断裂力学理论[56]，通过研究被压材料产生的裂纹长度，求出残余应力的大小和状态，相比于无应力的情况，受拉应力作用的材料的裂纹长度会增大，受压应力作用的材料的裂纹长度会缩短， 这种方法只能用于脆性材料，如陶瓷[57]、玻璃[58]等。 第二种方法[59]是通过纳米压痕实验生成的载荷-位移曲线来测定残余应力， 该方法基于残余应力对纳米压痕的接触面积、加载曲线和卸载曲线产生的影响，比如在被测材料中产生拉应力时， 压头下方的接触剪切力与拉应力的方向相同，会促进压痕过程，使卸载过程中的材料弹性回复减小，从而产生较大的残余深度，而残余应力为压应力时则效果相反。

测试材料的残余应力无固定方法，各个测试模型都有自身的优缺点。 Suresh 模型[60]是基于残余应力影响压痕接触面积大小而被提出，被测试样产生的残余应力是等双轴类型，适用于硬质材料、薄膜和涂层。 除此之外，经典测试模型[61]还有Lee 模型、Xu 模型、Swadener 模型、Wang 模型[62]等。 测量分析残余应力的方法还有X 射线衍射法[63]、拉曼光谱[64]、有限元仿真分析[65]等。

Liu 等[66]借助分子动力学仿真和纳米压痕技术研究了SiC/Al 复合材料界面微观残余应力在纳米压痕过程中的演变规律， 建立的系统模型见图11。该研究发现，在纳米压痕过程中产生的微观残余应力主要分布在压痕区域，随着压痕深度增加，压头接触区的微观残余应力分布逐渐明显， 压痕区SiC侧的微观残余应力明显大于Al 侧； 对于残余应力的释放，Al 原子通过自身的塑性变形降低了微观残余应力, 而SiC 原子则通过断裂释放了微观残余应力，SiC 原子这种通过断裂的方式会使SiC/Al 界面处的微观残余应力更集中以及组织演变更剧烈。

图11 SiC/Al 界面纳米压痕MD 模拟模型[66]

残余应力会对力学性能产生影响，特别是对载荷-位移曲线、硬度影响显著。 Zhu 等[67]利用自行设计的实验装置（图12），研究了等双轴残余应力对单晶铜的影响，验证了拉应力能降低峰值载荷、载荷曲线曲率、弹性回复深度、压痕功、堆积高度和接触面积并增加压痕残余深度。但是，压应力[68]对上述参数的影响与拉应力的情况相反。 卜家贺等[69]基于低阶应变梯度塑性理论（CMSG）[70]，提出了借助位错密度描述残余应力的理论模型，不仅发现压痕塑性影响区受残余拉应力影响，而且能准确预测压痕尺寸效应。 经过喷丸工艺处理的17Cr2Ni2MoVNb 钢[71]，引入合理的残余压应力，使材料具备了较低的摩擦系数和优良的耐磨性能。

图12 设计用于产生等双轴拉伸和压缩应力的同心弯曲装置[67]

4 结论与展望

传统的实验方法难以在微米尺度甚至纳米尺度进行测量。 而纳米压痕技术是一门在微观尺度上测量材料力学性能的新兴技术， 不仅能测试薄膜、涂层等材料的力学性能，而且可以大幅提高检测精度，有助于进一步研究微观力学特性。

通常，基于微观尺度开展的实验，其测试过程易受实验样本的表面粗糙程度、实验环境的变化和压痕点的间隔距离等多重因素影响，经常存在测量误差。 分子动力学模拟（MD 模拟）技术利用计算机技术，不仅节省了大量的人力、财力，而且可以模拟极端条件下的情况，更能形象地观察材料的微观结构变化，有助于进一步研究和提升材料的性能。 在特种加工领域，MD 模拟技术可对抛光、 超声振动、激光加工等工艺过程进行仿真，对于研究材料在加工过程中受到微观粒子的挤压、划擦以及冲击等行为而变化的力学性能机理也具有显著优势。 但是，MD 模拟技术也存在一定的局限性，比如，在模拟金属的蠕变行为时会受时间尺度的限制；在模拟纳米压痕时，会将压头钢度理想化，没有考虑压头的变形行为，等等。

未来，随着科学技术的发展，以及研究者对纳米压痕测试技术和MD 模拟技术研究的深入，上述局限有望被攻破。 可以预见，通过纳米压痕测试技术得到的硬度、弹性模量、残余应力等力学性能参数的数值会更加准确，而借助MD 模拟的纳米压痕技术也将迎来广阔的应用前景。