简介概要

单晶镍纳米薄膜单向拉伸破坏的分子动力学模拟

来源期刊：中国有色金属学报2004年第11期

论文作者：黄丹陶伟明郭乙木

文章页码：1850 - 1850

关键词：分子动力学; 镍; 单向拉伸; 镶嵌原子法

Key words：molecular dynamics; nickel; tension; embedded atom method

摘要：应用分子动力学方法模拟了单晶镍纳米薄膜受单向拉伸破坏的过程，得出纳米尺度单晶镍薄膜的应力—应变关系、能量演化曲线和镍薄膜构型的变化及微损伤的形成和扩展过程。模拟采用原子镶嵌势描述原子间作用，得到镍单晶薄膜的弹性模量，分析了拉伸过程中系统原子能量、应力变化和外荷载的关系。结果表明：纳米薄膜的自由表面影响拉伸过程中原子的运动和薄膜整体力学性能，纳米薄膜破坏的几何特征是原子空位的连接和晶胞缺陷的扩展；单晶的断裂接近脆性断裂，模拟得到纳米薄膜的断裂强度符合Griffith脆性断裂的能量平衡理论。

Abstract: The failure process of monocrystalline nickel film under tensile stress was simulated with molecular dynamics method. The evolvement of atomic energy and arrangement of atoms in the model, the initiation and expansion of damage and the elastic modulus of monocrystalline nickel were obtained. Simulation results show that the free surfaces of single crystal take effect on the motion of atoms and mechanical properties of nano crystal. Atomic cavities and the growth of crack in crystal lead to the failure of nano film. The fracture of single crystal is similar to macro brittle rupture, and the fracture strength of the film can be explained well by Griffith's theory of rupture.

中国有色金属学报 2004,(11),1850-1855 DOI:10.19476/j.ysxb.1004.0609.2004.11.011

单晶镍纳米薄膜单向拉伸破坏的分子动力学模拟

黄丹陶伟明郭乙木

浙江大学力学系,浙江大学力学系,浙江大学力学系杭州310027 ,杭州310027 ,杭州310027

摘要：

应用分子动力学方法模拟了单晶镍纳米薄膜受单向拉伸破坏的过程,得出纳米尺度单晶镍薄膜的应力—应变关系、能量演化曲线和镍薄膜构型的变化及微损伤的形成和扩展过程。模拟采用原子镶嵌势描述原子间作用,得到镍单晶薄膜的弹性模量,分析了拉伸过程中系统原子能量、应力变化和外荷载的关系。结果表明:纳米薄膜的自由表面影响拉伸过程中原子的运动和薄膜整体力学性能,纳米薄膜破坏的几何特征是原子空位的连接和晶胞缺陷的扩展;单晶的断裂接近脆性断裂,模拟得到纳米薄膜的断裂强度符合Griffith脆性断裂的能量平衡理论。

关键词：

分子动力学;镍;单向拉伸;镶嵌原子法;

中图分类号： TB383

作者简介：黄丹(1979),男,博士研究生.通讯作者:黄　丹;电话:0571 87951769;E mail:huang524@zju.edu.cn;

收稿日期：2004-03-27

基金：国家自然科学青年基金资助项目(10302025);

Molecular dynamics simulation for failure process of monocrystalline nickel film under tensile stress

Abstract：

The failure process of monocrystalline nickel film under tensile stress was simulated with molecular dynamics method. The evolvement of atomic energy and arrangement of atoms in the model, the initiation and expansion of damage and the elastic modulus of monocrystalline nickel were obtained. Simulation results show that the free surfaces of single crystal take effect on the motion of atoms and mechanical properties of nano crystal. Atomic cavities and the growth of crack in crystal lead to the failure of nano film. The fracture of single crystal is similar to macro brittle rupture, and the fracture strength of the film can be explained well by Griffith's theory of rupture.

Keyword：

molecular dynamics; nickel; tension; embedded atom method;

Received： 2004-03-27

随着微电机系统(MEMS)研究的深入和纳米电机系统(NEMS)的发展, 纳观尺度下材料力学、物理性能的研究显得日益迫切。由于实验样品制备和实验参数控制的困难, 目前的纳米力学实验通常局限于对材料弹性模量、硬度、蠕变等参数的测试 ^[1,2,3] , 且实验结果往往差距较大, 这使计算机模拟成为研究纳观尺度下材料行为的主要方法之一。

分子动力学方法(MD)是目前纳观尺度模拟的主要手段。它直接根据原子间的相互作用模拟原子的运动过程, 找出纳米晶体结构和力学性能的本质联系。而原子间相互作用势函数的选取是分子动力学计算精确与否的关键。原子镶嵌法(EAM) ^[4] 修正了传统的对势模型对多体效应考虑不充分的缺陷, 可以更精确地描述存在界面、表面、缺陷等不规则晶格环境的原子运动。

金属镍及其合金由于耐高温、抗腐蚀、具有良好的化学稳定性, 在医疗器械、化工、国防、精密仪器等领域得到广泛应用。已有不少文献描述了镍的化学物理性质、镍基复合材料的宏观性能 ^[5,6,7] , 但对其在纳观尺度下力学、机械性能的描述 ^[8] 相对较少。

本文作者应用分子动力学理论研究金属单晶受单向拉伸损伤破坏的过程。采用原子镶嵌法描述原子间相互作用建立了单晶镍薄膜的三维计算模型, 通过原子模拟考察了薄膜原子能量、应力的变化和结构变形过程、纳米薄膜微缺陷的产生、生长以及扩展过程, 得到纳米薄膜单向静态拉伸条件下的弹性模量、应力—应变关系、能量演化曲线, 讨论了自由表面、位错运动和原子滑移对纳米材料性能和变形的影响。并从能量平衡角度用Griffith断裂理论解释了单晶镍薄膜的断裂。

1 分子动力学模拟

1.1 模拟原理与方法

镶嵌原子法适用于金属力学性能模拟。本研究采用Voter等 ^[9] 根据镶嵌原子法提出的镍的多体势函数。设原子总势能为

$E_{t o t a l} = \sum_{i} [\frac{1}{2} \sum_{j} Φ (r_{i j}) + F (ρ_{i})] ? ? ? (1)$

式中 Φ(r_ij)为相距r_ij的原子i, j之间的中心对势; F(ρ_i)为到电子云ρ_i的原子镶嵌能; ρ_i为原子i处的电子云密度。且有

Φ(r_ij)=A₁(r_c₁-r_ij)²exp(-c₁r_ij) (2)

$\begin{array}{l} F (ρ_{i}) = D ρ_{i} \ln ρ_{i} \\ ρ_{i} = \sum_{j} f (r_{i j}) ? ? ? (3) \end{array}$

其中

f(r_ij)=A₂(r_c₂-r_ij)²exp(-c₂r_ij) (4)

式(4)表示离原子i距离为r_ij的原子j对原子i处电子云的影响。 A₁、 A₂、 c₁、 c₂、 D为由材料物理性能决定的参数。 r_c₁和r_c₂分别表示计算原子对势和电子云密度的截断半径, 即在原子间距离r_ij超过截断半径时, 分别忽略两原子间的对势和电子云密度的相互影响。参照文献 [ 10] 计算面心立方金属镶嵌势函数参数的方法, 取r_c₁=1.65d, 在次近邻原子和第三近邻原子之间, r_c₂=1.95d, 在第三近邻和第四近邻原子之间。 d为最近邻原子距离, 对于面心立方晶格, $d = a_{0} / \sqrt{2} ? a_{0}$ 为晶格常数。

原子间的作用遵守牛顿第二运动定律

$F_{i} = - \frac{? E_{t o t a l}}{? r_{i}} = m_{i} ? v_{i} ? ? ? (5)$

式中 F_i、 r_i、 m_i、 v_i分别为原子力、原子位置、原子质量和原子速度。对动力学方程的求解通过Verlet算法 ^[11] 的速度形式 ^[12] 计算积分。

1.2 模型的建立

图1所示为单晶镍薄膜原子模型示意图。图1(a)所示为镍面心立方结构(fcc)单胞示意图, 图1(b)为整体模型。 x, y, z坐标轴分别对应晶胞的[100], [010], [001]晶向。采用规则晶格尺寸, 为使纳观薄膜视觉上形象化, 取6×20×20个晶胞, 共9 600个原子。模型原始尺寸为2.112 nm×7.04 nm×7.04 nm。模拟中长度单位为镍的晶格常数a₀=0.352 nm, 根据系统能量的稳定测试, 分子动力学模拟的时间步长定为t₀=3.5×10^-15 s, 能量单位为E₀=10^-19 J。将x方向的边界控制为自由表面, 在y, z方向上施加周期性边界条件, 使原子模型呈单晶薄膜结构。将系统温度初始化为绝对零度, 并在模拟中保持零温, 避免原子热激活。真实材料中原子的热运动不改变绝对零温模拟所得力学参数变化方向和趋势。

图1 单晶镍薄膜原子模型

Fig.1 Atomic model of nano nickel film

(a)—Nickel lattice; (b)—Simulation model

1.3 模拟过程

对模型进行温度初始化, 并控制温度始终为绝对零度, 避免原子热激活。首先对单晶模型原子弛豫20 000步, 使系统有充分时间达到能量最低的稳定状态。沿z方向施加平面拉伸应变0.005, 进行分子动力学模拟迭代10 000步, 时间为3.5×10^-11 s, 然后弛豫20 000步, 使系统回到平衡态。再增加拉伸应变0.005, 重复此施加应变—分子动力学模拟—弛豫过程, 使模型原子处于准静态拉伸受力状态。拉伸过程中采用Parrinello-Rahman ^[13] 方法调节x, y方向的尺寸来控制两方向的压力不变, 保持在一个大气压。持续静态加载模拟至薄膜中出现原子空位时, 降低施加荷载速率, 将应变增幅改为0.001, 重复上述过程, 最后薄膜局部破坏,破坏处和薄膜表面由于大量原子摆脱系统的作用而离开模拟空间, 模拟结束。模拟过程中每隔100步记录原子的应力矢量, 总能量, 势能, 动能, 薄膜形状、各方向横截面形状和原子座标位置。

2 模拟结果和讨论

2.1 纳米薄膜结构变化分析

图2所示为纳米镍薄膜在z向拉伸过程中x-y横截面原子构型的演化图。其中图2(a)所示为经过初始弛豫而未加外荷载的时刻, 图2(b)～2(d)分别为ε=0.10、 0.15、 0.17时刻的横截面原子排列形状。从图中可以看出: 在弛豫状态下, x方向(薄膜厚度方向)的表面原子由于配位数与内部原子不同, 在空间方向失去相邻原子而形成断键, 处于不规整晶格环境。受力不平衡导致表面原子偏离理想点阵位置并受到晶体内部原子的拉应力, 使得经过弛豫后的薄膜表面有收缩趋势(图2(a)); 随着z方向拉伸应变的增加, 表面附近的几层原子在外力和内部原子拉力的作用下开始偏离弛豫态的平衡位置, 并随着位置的变化对内部原子产生不同的作用力。当应变达到0.10时(图2(b)所示), 在薄膜的自由表面出现明显的原子运动留下的空位, 晶格破坏向薄膜内部深入几层原子, 在角落处的局部表面出现沿横截面45°方向发展的滑移面。内部几层晶胞原子排列基本不变, 仍为理想规则排列; 从ε=0.10到ε=0.15的过程中, 薄膜表面的原子空穴相互联接, 在自由表面形成裂纹, 并向薄膜内部扩展(图2(c)), 内部晶格破坏数量增加。由于加载速率为准静态, 表面裂纹尖端区域有足够的时间来发射位错, 使裂尖钝化, 故材料在此阶段表现出微弱的韧性; 外荷载继续增加时, 薄膜表面裂纹延伸, 出现连续的明显裂纹, 且有大量表面原子势能完全释放, 脱离薄膜系统原子的作用, 离开模拟空间, 薄膜内部脆弱部位的纳观孔洞生长、连接导致薄膜在拉伸方向上出现贯穿薄膜厚度的连续裂纹, 薄膜在y向也趋于断裂(如图2(d)), 模拟结束。

图2 纳米薄膜x-y截面的原子排列

Fig.2 Atomic arrangement on x-y section in nano nickel film

2.2 能量演化过程分析

图3所示为拉伸过程中原子平均总能量随应变变化的数值拟合曲线。在本研究中ε=0.03以前取模型全体原子总能量的平均值, 而此后根据模型的变形特点, 取图2中横截面位置附近的10层原子作能量平均。由图可知, 从初始弛豫自由状态到外荷载增加到ε=0.03的时刻, 系统能量基本维持在一条水平线上, 这是因为纳米薄膜比表面积大, 表面原子由于对称性被打破, 原子配位不足, 在空间方向形成断键, 而在初始弛豫阶段储存了大量的表面能。在初始加载时, 内部原子受外荷载作用运动, 动能增加, 势能也从平衡态开始上升, 单晶薄膜内部原子总能量上升, 但表面原子在外力作用下开始释放因初始受力不平衡而储存的大量表面能, 致使系统总能量能维持基本不变; 从ε=0.03开始, 系统总能量开始缓缓上升, 此后随着表面原子表面能的完全释放, 所有镍薄膜系统原子在外荷载作用下势能升高。此后, 原子活性高的表面晶胞首先开始出现纳观破坏现象, 少数原子在外荷载作用下偏离理想位置附近, 留下原子空位和孔洞。原子运动加剧, 系统能量升高速度也增加(从ε=0.08至ε=0.14)。另外, 在此阶段中局部的表面原子沿一定方向滑移, 少量的位错移动也消耗部分原子能量。在ε=0.14至ε=0.16的过程中, 模型破坏区域系统原子排列混乱, 薄膜中微缺陷的相互连接扩展使大量原子运动远离平衡位置, 原子平均能量继续上升, 纳观缺陷的连接长大导致破坏区部分自由度高的原子之间距离超过截断半径, 形成新的表面, 模型沿拉伸方向断裂, 新表面的形成使所计算部分的原子总能量上升到最高值(图3); 此后若在拉伸方向取一个部分继续加载, 可见在薄膜厚度方向也会出现贯穿的连续裂纹, 薄膜沿y方向趋于断裂。但自沿拉伸方向断裂以后, 模型中的大量原子离开模拟空间, 系统的总能量迅速下降。

图3 镍薄膜原子总能量随应变的变化

Fig.3 Variation of atomic total energy with stress in nickel film

2.3 原子应力与应变分析

应力—应变关系在材料的性能研究中有很重要的作用。图4所示为镍薄膜在外荷载作用过程中的应力演化曲线和沿拉伸方向的应力—应变关系, 曲线中应力为系统中原子的平均正应力。由于薄膜在厚度方向无法承受外荷载, 而在薄膜平面又具有二维对称性, 因此z方向的应力—应变关系和弹性模量拟合应同样适用于y方向。在模拟初始时刻(ε=0), 经过自由弛豫的纳米模型内部原子应力为零, 而自由表面原子由于空间上不对称, 受到内部原子的拉力, 表面几层原子之间的相互作用力导致系统原子平均应力不等于零。由模拟结果输出可知, 在经过弛豫而无外荷载状态, σ_x=1.945×10⁴ Pa, σ_y=σ_z=1.657×10⁵ Pa, 表明弛豫状态表面原子沿薄膜平面的面内应力超过厚度方向上的原子作用力。由图4可以看出, 在拉伸过程中, 在拉伸方向的正应力σ_z最大, 薄膜平面的另一方向应力σ_y次之, 在薄膜厚度方向上的应力σ_x最小; 而在拉伸至破坏的过程中, 在三个方向上的正应力在不同阶段以不同斜率作近似于线性的上升。在ε=0.06～0.14阶段, 由于薄膜自由表面局部滑移和在滑移线末端产生少量不全位错, 在内部晶胞产生缺陷时, 在纳观缺陷附近(如纳观裂纹的尖端)也发射位错, 位错的移动消耗外加能量, 应力曲线的斜率逐渐下降, 使纳米薄膜表现出微弱的塑性。但纳观薄膜应力—应变曲线并没有出现宏观材料拉伸实验中明显的屈服阶段, 验证了此前分子动力学模拟中晶体材料基本不存在应变强化的结论 ^[14] 。在经历原子滑移阶段后, 随着外荷载的增加, 在ε=0.148时刻, 在破坏区域与拉伸方向垂直方向的几层晶胞中微缺陷连接, 出现连续裂纹, 薄膜趋于在z方向断裂, 三个方向的正应力均达到峰值。此后若经过原子弛豫再取部分模型继续加载, 在ε=0.17时沿厚度方向也出现贯穿裂纹(图2(d)), 薄膜趋于在y方向也断裂。

图4 镍薄膜原子的应力—应变曲线

Fig.4σ—ε curves of atoms in nickel film

从初始加载到断裂的应力—应变曲线表明纳米薄膜沿拉伸方向的应力变化快, 且变化基本呈线性。由于沿厚度方向和薄膜平面内的滑移及少量位错的存在, 与拉伸方向垂直的两个方向应力变化较缓慢。单晶材料沿拉伸方向的应力—应变曲线类似于宏观脆性材料, 拉伸方向的应力最大值为σ_z=21.08 GPa, 可以视为单晶薄膜单向拉伸的断裂强度。从应力曲线上看, 在拉伸破坏模拟中纳米单晶材料没有表现出宏观多晶体的屈服阶段和塑性流动过程, 断裂过程接近脆性断裂。

2.4 弹性模量拟合

由加载过程中的应力、应变值, 可以经过最小二乘法线性拟合 ^[15] 得到纳米材料的弹性模量。在应力线性变化阶段(ε=0～0.01, 由前文分析可知此阶段薄膜自由表面原子逐渐释放表面能, 内部晶胞仍为规则排列, 原子在平衡状态附近运动, 薄膜中尚未发生永久变形或形成纳观缺陷), 将载荷增幅下降为0.001, 取模拟中所得沿拉伸方向的应力、应变值进行拟合, 得到弹性模量为186.6 GPa。此值微小于多晶镍的弹性模量(204.00～220.64 GPa), 也小于Shen ^[16] 等对无微孔隙纳米镍晶粒样品的测试结果(217 GPa)。因为模拟中采用的纳米薄膜模型中大量的表面原子运动自由度高, 原子键结合强度下降, 降低了材料的强度。值得注意的是, 本研究拟合所得的弹性模量仅为单晶镍薄膜〈100〉方向的弹性模量。

2.5 损伤演化过程与破坏分析

由图2中镍薄膜破坏发展区横截面原子排列的变化和本文中能量与应力的演化过程分析可知纳米薄膜模型中材料损伤产生、演化和破坏过程。初始弛豫时, 为达到稳定状态, 薄膜自由表面的最外几层原子受到拉力, 偏离规整晶格位置, 使薄膜表面呈收缩状, 并储存了大量表面能; 在受到薄膜面内单向拉伸应变作用时, 内部晶胞原子在外荷载作用下在平衡位置附近振动, 而活性大的表面原子则在外力作用下释放表面能, 表面晶格从收缩状态松弛; 随着外荷载的增加, 自由表面局部原子开始离开理想位置, 留下原子空位。部分表面晶胞原子沿一定方向滑移, 并在滑移线终端出现少量不全位错; 晶格破坏从表面向内部深入, 并在薄膜表面局部出现明显的滑移面(图2(b))。在应变较大时, 微缺陷(如纳观裂纹尖端)附近发射少量位错在薄膜平面和厚度方向的移动消耗外部能量, 而滑移和位错导致晶格破坏产生的微缺陷相互联接长大, 在薄膜内部产生连续裂纹, 并沿裂纹尖端扩展。表面原子空位的增多和位错发射导致薄膜沿荷载作用方向(z方向)出现类似宏观的颈缩现象(图2(c))。表面空位和薄膜内部裂纹的扩展使薄膜在y拉伸方向断裂。 z方向断裂后的薄膜经弛豫若选取部分继续施加外荷载, 表面原子脱离系统的作用逸出薄膜表面, 薄膜在y方向(垂直于荷载方向)上也出现连接、扩展到贯穿薄膜厚度的连续裂纹, 使薄膜在向趋于断裂(图2(d))。纳米薄膜的断裂并没有明显的屈服阶段, 而是在外荷载增加的过程中由于截面裂纹的生长出现突然断裂, 说明薄膜延性不高, 这是由纳米单晶与多晶材料不同的变形机制决定的。

根据Griffith关于理想脆性材料能量平衡的经典断裂理论及其适用于非理想脆性材料的修正理论, 理想晶体断裂时外荷载做功与系统应变势能和随断裂产生的自由表面能应达到平衡, 晶体断裂的最大应力应等于晶格结构可以承受的最大应力即理论强度。在本研究的模拟中薄膜沿拉伸方向断裂时, 原子最大正应力σ_z=21.08 GPa, 按照文献 [ 16] 对镍晶粒力学性能的测试, 弹性模量E=217 GPa, 本模拟中镍的断裂强度已经非常接近晶格的理论强度(一般取为0.1E), 完全符合Griffith的理想晶体能量平衡断裂理论。

3 结论

对单晶镍纳米薄膜在静态单向拉伸破坏过程的分子动力学模拟表明: 纳米薄膜由于比表面积大, 表面大量活性高、能量高的自由原子的存在, 影响了材料的能量、应力特性和变形机制。在拉伸过程中, 原子总能量持续升高, 沿拉伸方向的正应力最大, 并在拉伸方向首先产生连续裂纹, 导致薄膜断裂; 在外加应变较小时薄膜自由表面原子释放表面能, 内部晶胞的原子在理想位置附近振动。应变增加到一定程度时表面原子离开理想位置形成空位和晶格缺陷, 缺陷附近发射少量位错, 位错的移动使薄膜表现出微弱塑性。纳米镍的应力—应变曲线无明显屈服阶段, 接近宏观脆性材料的拉伸断裂曲线; 自由表面原子运动留下空穴的联接和自外向内的滑移以及微缺陷的扩展是镍薄膜断裂的几何特征; 由于大量表面原子带来的不稳定性, 模拟所得单晶镍在〈100〉方向的弹性模量E=186.6 GPa, 小于多晶镍和实验测得理想镍晶体的弹性模量; 模拟得到镍薄膜断裂强度为21.08 GPa, 几乎等于实验所测纳米镍晶粒的理论强度(约为21.7 GPa), 符合Griffith关于理想晶体断裂能量平衡的脆断经典理论。