中空铜纳米线的拉伸断裂分布与初始滑移分布的关系

刘守涛 赵健伟 王奋英 马汉杰

(1浙江理工大学信息科学与工程学院，杭州 310018)

(2嘉兴学院材料与纺织工程学院，嘉兴 314001)

(3南昌大学化学化工学院，南昌 330031)

(4浙江理工大学计算机科学与技术学院(人工智能学院)，杭州 310018)

0 引 言

金属纳米线(NW)是一种具有表面原子占比高和导电性好的材料，在微纳器件设计制造中有着潜在应用价值[1-4]，尤其是诸如内在中空微结构缺陷，不仅是器件功能化的基础，也对力学性能、电学性质以及使用寿命有着重要的影响。因此，研究金属纳米材料的中空结构与其形变断裂失效的关系成为器件设计与优化的关键基础问题。

对于金属纳米材料的应变与失效的实验研究，一般限于极慢速的近平衡态形变。例如，Oshima等[5-6]通过透射电子显微镜(TEM)成功地观察到了在超高真空条件下螺旋状多壳层金NW 的形成。Kushima 等[7]通过使用TEM 内的原位拉伸实验对全锂化硅NW 的断裂强度和塑性进行了定量测量。Li等[8]通过TEM 原位拉伸实验在单晶镍中精确测量了单个位错形成的激活参数。Wang 等[9]通过TEM 技术在原子尺度上追踪晶界的滑移，揭示了晶界滑移的微观机制。Kiani 等[10]利用TEM 对不同取向角的双晶和单晶试样进行拉伸试验，表明失效均是通过塑性塌陷发生的，并不依赖于错配角度或晶界能量，而是更多地依赖于样品的表面结构和位错活动。这些研究的体系直至断裂都保持较好的结晶状态。对于快速的准平衡态和快速的非平衡态拉伸，其实验方法需要一定时间，无法及时捕捉NW 中瞬时发生的结构变化，相比之下，分子动力学(MD)模拟实现了对各个细分时刻NW 结构变化的全过程追踪。Cheng 等[11]用MD 模拟，考察了一系列具有不同初始结构NW 的拉伸行为，发现纳米结构的局部几何特征在决定机械失效方面起着主导作用。Pang等[12]通过比较不同取向NW 的模拟，发现塑性变形模式随着NW 取向的变化而变化。Ma 等[13]采用MD模拟探究了镁单晶在不同拉伸速率下的结构转变和位错变化，结果表明拉伸速率对镁单晶的位错形成和滑移方式有显著影响，而且高速拉伸导致镁单晶发生结构转变，使其塑性变形受到限制。此外，Xing 等[14]采用MD 模拟镁双晶中[1100]对称倾斜晶界沿晶间裂纹扩展的机理，表明在低温下大多数不同错配角度的裂纹都通过脆性断裂机制扩展。

然而，基于大样本的NW 拉伸断裂模拟表明，断裂位置并非确定量，而是遵循统计分布规律[15]。该统计分布曲线中显示一个最概然的断裂分布位置，且该位置与应力波的传播密切相关。因此，断裂分布曲线表现出和温度[16]、拉伸晶向[17]、长度[18]以及缺陷[19]等条件的相关性。掌握断裂分布的统计特征，特别是明确断裂位置的影响因素，是避免金属纳米材料在载荷条件下失效的关键一步。进一步研究表明，单晶NW 以外的一些微观结构，例如孪晶界、凹凸结构[20]等均会对断裂失效分布特征带来影响。近期工作还表明NW 初期产生的滑移位置在特定条件下与断裂位置分布存在一定联系[21]，金属纳米材料的中空结构既是实现其功能化的核心，也会对其机械力学性质产生影响。因此，探讨中空结构如何决定初始滑移缺陷的产生，以及了解初始滑移面产生位置的统计分布特征与断裂失效位置分布之间的关系是至关重要的。

我们通过大样本进行统计分布分析，探讨中空结构如何影响金属NW 在拉伸过程中的2 个重要阶段，即屈服阶段和最终断裂阶段。模拟共设置了6个内部含不同中空半径大小的NW 模型，通过机器学习算法和大数据统计考察屈服点附近的初始滑移面的产生和最终断裂位置的分布，并进一步讨论两者之间的关系。

1 模型建立与计算方法

1.1 模型建立和MD计算

基于长轴Z沿[100]晶向的单晶铜NW(17a×17a×51a，a为晶格参数，0.362 nm)，移除中间给定半径的球形空间内的原子，形成中空结构，其截面如图1 所示，为便于讨论，标记模型为R1～R5，依次对应中空半径1～5个晶格常数。为了解中空结构的影响，比较了完美单晶铜NW R0，全部模型在X、Y和Z方向上均采用自由边界条件，在NW 的两端采用了各3层铜原子的固定层，以便施加匀速载荷。

图1 模型初始结构图Fig.1 Model initial architecture diagram

在MD 模拟过程中，采用Leap-Frog[22]算法求解牛顿运动方程，其中步长使用2.5 fs，同时利用Verlet和Cell 链接列表[23]构建相邻列表，将温度设定为150 K，使用校正因子法[24]对系统进行恒温处理，利用嵌入原子(EAM)势[25]来描述原子间相互作用，并采用Johnson[26-27]优化的参数，在每个模型施加载荷前进行自由弛豫以消除应力。弛豫时间平均为12 500步，每2 个样本间相隔50 步以获得不同的初始态。对应于每一个模型，分别模拟300 个不同初始态的样本，以往研究表明，该样本数足以获得NW 的统计分布特征[15-18]，R0～R5 6 个体系，共计1 800 个样本。所有模拟均使用了NanoMD[19-20,28]大规模MD 仿真软件，该软件的可靠性经过了广泛的模拟研究验证，同时模拟结果也与实验结果一致[29]。在拉伸过程中使用了0.1%·ps-1(每一端的绝对速率为19.1 m·s-1)的恒定应变率向两端进行拉伸。在每500步记录一次NW 体系中原子的类型、坐标、应力、应变和能量，以便进一步分析。

1.2 数据分析

按照位力展开公式[30]计算应力-应变曲线。径向分布函数(radial distribution function，RDF)是反映材料内部结构特征的物理量，它直接描述了材料内部原子排列的有序程度，选定一个原子为中心，半径为r、厚度为δr的球壳体积V按下式计算：

若单位体积内的粒子数为ρ，则球壳内的总粒子数为4πr2δr。

径向分布函数g(r)是距离某一原子为r时找到另一个原子的概率，ρ0为单位体积内的平均粒子数密度。若球壳内的粒子数为n(r)，那么g(r)表示如下，详细计算过程参考书籍[24]。

初始微观结构缺陷是指在屈服点产生的滑移，称为初始滑移。利用中心对称参数[31]来获取六边形紧密堆积(hcp)原子，滑移面上的原子属于hcp原子[32]。但是，在聚类过程中，考虑到初始滑移的形状和大小各不相同，并且每个样本的初始滑移数量都是随机的，因此选择了基于密度的噪声应用空间聚类(density-based spatial clustering of applications with noise，DBSCAN)算法[33]，将紧密相连的hcp原子聚类到一个滑移平面中。因为DBSCAN 算法不需要预定义群集形状或大小，因此能够识别任意形状的群集，包括稀疏的和具有复杂边界的群集，并且不需要事先指定群集的数量，这是相比K-means 不同的优点(表1)，这使得DBSCAN 在处理本数据集时更加灵活，无需事先给出聚类，其中参数(ε=0.08，MinPts=8)是必需的，它描述了一组邻域内的紧密性，其中ε是特定hcp原子的邻域距离阈值，MinPts 是特定hcp原子的邻域ε内的最小原子数。初始滑移与NW 的横向边缘相交的坐标z被定义为初始滑移位置。

1.3 DBSCAN算法思想

DBSCAN 是一种基于密度聚类算法，其主要核心思想是通过考虑数据点的密度来执行聚类，其中核心点是半径ε内至少包含最少数量MinPts个数据点的点，直接密度可达关系用于将数据点连接到核心点，而密度可达关系允许跨越稀疏区域连接簇，同时标记那些无法连接到核心点的点为噪声点。原子划分如下：

epsp={q∈H，dp,q≤eps}

其中eps表示邻域距离，epsp是点p的邻域，包含所有与点p的距离不超过esp 的点q，dp,q是点p和点q之间的距离，H是包含在应力屈服点处的所有hcp原子的数据集。

伪代码见Supporting information。

2 结果与讨论

2.1 中空NW的内表面性质与屈服点结构特征

中空结构的内表面原子所处的微观环境与均质的块体材料不同。由于受到周边原子的吸引作用，中空内表面原子存在一个指向外的表面张力，表2 汇总了各模型的结构特征和表面能(ES)。其中，NBall是去除球体中的原子数，SBall∶SNW为中空球面积与NW 侧面积之比，ES=(ED+ERB-END)/2(ED为含有中空结构NW 的总能量，ERB为移除球体的总能量，END为无缺陷的单晶铜NW的总能量)。

表2 各模型的结构特征和能量Table 2 Structural features and energies of various models

从表2 可得，随着中空半径增大，ES呈增加趋势，说明了中空结构特别是较大半径的中空结构会因其内表面张力对滑移面的产生带来影响。图2给出了在屈服点6个模型代表性样本的滑移面与中空结构的原子排布位图。从中可以看出，对于单晶铜，体系内往往会伴随着多个初始滑移的产生，且易从侧棱上开始，还可以观察到滑移面与表面呈20°～40°向上或向下发展。对于较小半径的中空结构，其滑移面的产生和发展与单晶类似，这是由于中空结构小且远离侧壁，因而对后者的影响也不显著，相比较于中空内表面原子，侧棱原子的配位数更低，能量更高，因此占据主导作用，滑移面同样产生在侧棱上。但当中空结构的半径足够大，高能内表面原子与侧面产生一定的协同作用，使得初始滑移向中空结构附近集中。例如R4 的中空截面与NW 截面比达到1∶20 时，这一特征变得较为突出。而更大的中空半径，使得最近的内表面原子距侧壁原子达到3.5a时，两者之间的影响已非常强烈，因此在模型R5 时，不仅在侧棱上产生滑移面，同样也在中空内表面上产生滑移面。这也与先前的研究一致，Sun 等[34]在研究含缺陷银NW 的形变机理时，发现内部含较大空隙无法抵抗外部横向应力时，将导致第一个位错在缺陷附近发生。这一现象说明中空结构必将影响体系拉伸后续阶段的力学行为及断裂分布。

图2 应力屈服点产生的初始滑移Fig.2 Initial slips at the yield point due to stress

2.2 拉伸形变的应力-应变关系

应力-应变曲线给出了系列模型宏观力学性质的基本描述(图3)。图中仅给出了5个代表性的中空铜NW 与无缺陷的单晶铜NW 比较。从图中变化趋势可以看出，在相同的模拟条件下，屈服点(strain≈0.042)之前，R0～R4 五个体系的应力几乎重合，但发现不同初始状态的单晶铜NW 的应力在弹性变形区域内近似线性增加[35]。此外，弹性变形过程中应力-应变曲线的斜率表示弹性模量[12]，显然弹性模量非常一致。Li等[20]在对初始变形行为的研究中也得出了类似的结果，这是因为系统没有结构性损伤。表3 给出了基于300 个样本的屈服特征的统计结果。以往的研究表明初始结构对屈服行为的影响具有多样性，同种材料的不同初始结构可以保留弹性行为，Cao[36]研究[111]晶向孪晶铜NW 的形变机理时，发现孪晶片层厚度小，孪晶NW 的屈服应力越大。Deng 等[37-38]基于类似的模型，给出了孪晶铜NW 的屈服应力和单位长度内的孪晶界个数的线性关系。上述研究均指出屈服行为与初始位错滑移的密切关系。而我们考察的前5 个体系的屈服点基本相同，同样，Wang 等[39]在研究孔洞和孪晶界对银NW形变行为时提到，屈服应力并没有随着孔径的增加而发生显著变化。而在本研究中，直至体系R5屈服应力与屈服应变才开始减小，这个结论与Sun 等[34]对含有缺陷的银NW 形变机理的研究结果相同，在银NW 达到屈服应力点时发生了相似的变化，说明中空结构的大小在一定范围内对初始滑移几乎无影响。这也意味着对塑性形变起关键作用的位错滑移并非产生于上述初始结构。

表3 不同中空结构下的屈服点处的应力-应变值Table 3 Stress-strain values at the yield point for different hollow structures

图3 不同中空结构下的应力-应变曲线Fig.3 Stress-strain curves for different hollow structures

达到屈服点之后，应力快速释放，而后会在一个较稳定的应力区间波动，维持塑性形变。该应力区间的平均高度和应变长度与中空尺寸密切相关，尺寸越大，应力越低，应变范围越窄。特别是到了半径为5 个晶格的中空体系，塑性形变区间对应的平均应力只有1.4～1.5 GPa，同时维持的应变范围也不足0.05。而R1 和R2 两个模型的塑性形变应力维持在约2.0 GPa，应变范围约0.4，说明中空结构明显降低了金属NW的延展性，其强度更低且脆性更大。在塑性形变后期，材料一般会形成颈缩并伴随着应力的快速降低。对单晶而言，塑性形变与断裂之间连续过渡，应力也是缓慢降低到0，表现了良好的塑性形变特征，随着中空半径的增加，应力下降曲线越来越陡峭，甚至表现了一定程度的指数衰减的特征，断裂应变也同样是随着中空半径的增加而降低。

图4 给出了平均的原子能量随应变的变化关系，由于中空NW 比单晶NW 多了空心的内表面，它的起点能量也随着中空半径的增加不断提升。在拉伸的弹性形变阶段，能量以近似二次方的形式增长，在屈服点附近能量也到达一个局部极值。而后，伴随着滑移的产生，不仅应力得到释放，原子平均能量也随之降低。中空结构的尺寸对释放的能量有较大的影响，中空半径越大能量释放越小，在塑性形变阶段能量恢复上升。其中单晶材料表现了极好的塑性，能量持续上升。小半径中空NW 表现了相似的特征，但在塑性形变后期上升的能量低于单晶NW，而更大的中空结构表现出不同的特征，特别是R3、R4 和R5 三个样本的能量在塑性形变中后期保持在一个较低的平台，直至断裂。图4 可以看出能量表现了与应力-应变曲线一致的变化趋势，反映出中空结构提升了原子平均能量，降低了屈服应力，减小了延展性。在图4 中虽仅绘出了5个代表性的例子，不同样本呈现一定的波动，但基本趋势符合上述论述。

图4 不同中空结构的铜NW能量随应变的变化Fig.4 Variation of energy with strain for copper NW in different hollow structures

2.3 拉伸形变过程中原子类型变化与结晶特征的统计分析

滑移面产生时，面上的原子的邻近原子配位数发生了变化，因而可以利用配位数分析[40]确定对组成滑移面原子数量的变化。图5给出了对滑移面有贡献的hcp原子数随着拉伸形变的变化关系，也同时对比了本体特征的面心立方(fcc)原子数的变化。从这6 个模型中可以看到以下的变化趋势：在屈服点之前(应变小于0.04)hcp原子数和fcc原子数无明显变化，说明在弹性区间体系保持了较完整的结晶状态。在应变大于0.05 时，体系进入了能量和应力的释放区，hcp原子数迅速升高，相应地，fcc原子数在减少，代表着持续地产生滑移面，直至0.07～0.08的应变区间。应力及能量释放完后，hcp原子数不再增加[19]，同样，fcc原子数不再减少。在随后的塑性形变区间，hcp原子数与fcc原子数都在保持振荡，表明滑移面持续产生-消失，这一过程一直持续到塑性形变的中后期。在断裂阶段，不同的中空结构表现了不同的趋势。单晶及小半径的中空结构会保留hcp原子数和fcc原子数的振荡降低，这与塑性特征一致，即在断裂前，滑移面依然在不断地产生与消失。但对于较大半径的中空NW，hcp原子数则快速降低到某一数值，同样，fcc原子数则快速增长到某一数值后，呈现出极弱的波动，直至断裂。这说明NW 在这一阶段几乎无新的滑移面产生，断裂处表现了较大的脆性。值得强调的是这6个代表性样本描绘了中空结构影响的总体趋势，就如同图4能量曲线一样，同一结构不同样本间会有一定的波动。这也进一步说明了在NW 研究中大样本统计分析的重要性，通过对多个样本的分析也明确了上述规律描述的正确性。

图5 fcc和hcp原子数对应变的变化Fig.5 Variation of atomic counts with strain for fcc and hcp

晶体NW 在拉伸过程中的近程和远程有序性可以通过RDF 来分析[24]。仅以在0.3 应变(如图3 与图4 中垂直粗线所指示)时各模型的RDF 为例，图6 分析了中空结构的影响，其中单晶R0、R1 和R2 处于塑性形变中期，伴随有大量的滑移面的产生与消失，而R3、R4 和R5 则处于塑性形变末期，NW 正形成颈缩，即将断裂。距离0.701 为第一近邻峰，峰高与峰宽随着中空半径的增加都有所增加。对于较远程，例如第四到第六近邻峰，峰高随中空半径的增加经历了先降低后升高的过程，这也与前文讨论的R0、R1 和R2 的塑性形变特征和R3、R4 和R5 的塑性形变特征一致。

图6 RDF曲线Fig.6 RDF curves

2.4 中空NW 拉伸形变的初始滑移分布与断裂位置分布之间的关系

功能性中空NW 有望在一系列器件中得以应用，但其在载荷下可能发生形变，以及断裂失效，因此了解中空结构对拉伸初始阶段产生的滑移面和断裂位置之间的因果关系对器件设计、结构优化、寿命延长等至关重要。利用DBSCAN 算法确定了初始滑移面与NW 侧面的交点，以归一化的NW 长度为参照统计了300 个彼此独立的样本，确定了初始滑移面产生位置的分布柱状图。因每个样本中产生不止一个滑移面，故统计的总数也大于样本数，前者为后者的2～4倍，意味着在屈服点附近每个样本中存在2～4个滑移面。

从图7 中可以看出，塑性形变特征明显的R0、R1 和R2 的初始滑移面的位置相对分散，既有较大概率出现于NW 两端，又有较高的可能出现于NW中间。但随空心半径增加，其影响已在R2模型中有所体现，即0.5 中间处出现概率降低，但其两侧，如0.22～0.42 和0.65～0.9 相对较高。由空心结构导致明显脆性形变的R3、R4和R5，中空结构的影响变得显著，尤其半径较大的R4 和R5 在0.3～0.45 和0.6～0.7 出现了2 个呈高斯分布函数的峰，峰高呈现一定的非对称性，这可能是NW 构筑时原子…ABCA…不对称排列顺序导致的。这种fcc金属排布的不对称性带来的影响也体现在其他方面，如Wang 等[41]对断裂分布与缺陷率关系的研究中，当应变率是1.0 ps-1且无缺陷时，断裂分布在中间位置展现出完美的高斯分布特征，当缺陷率上升至2%时，断裂位置移动到两端0.3 和0.7 处，两处表现出不对称高斯分布的特征，并且随着缺陷率的升高不对称现象更加明显。由此可知，不对称分布来源于fcc金属排布的不对称结构。

图7 不同中空结构下的初始滑移分布Fig.7 Initial slip distribution in various hollow structures

图8 给出了基于每组300 个样本的NW 断裂位置统计分布图。由图8a 可知，对于单晶NW 断裂位置的统计分布呈现靠近两端的偏态分布，这是由于本研究采用了自由边界条件，两端拉伸作用的固定层对断裂分布产生一定的影响，平均的断裂位置靠近两端(0.1和0.9)，且统计分布的直方图可以近似用偏态分布曲线拟合，峰值给出最概然断裂位置[15]，且该分布峰窄，相比之下0.2～0.8很宽范围内断裂分布稀少。图8b 中，R1 因其较小的中空半径表现出与单晶NW 相同的特征。对于图8c，中空半径更大的R2表现出由塑性形变向脆性形变的过渡特征，两端的偏态断裂分布占主导，但在中间部位也呈现出一个较弱的正态分布。而更大中空半径R3 使断裂性质由塑性形变断裂转变为脆性形变断裂，但图8d还保留了部分塑性形变断裂的特征，即NW 的两端呈偏态分布，中间呈正态分布，但两端的偏态分布相比中间的正态分布弱得多。模型R4和R5表现了完全的脆性断裂特征，仅显示了中间部分的正态断裂分布，且中空半径越大分布越窄[42]。

图8 不同中空结构下的断裂分布Fig.8 Fracture distribution in different hollow structures

对比图7 和图8 可以看到中空结构存在的条件下，初始的滑移分布与断裂分布之间的相关性。对于单晶和中空半径较小的NW，形变以塑性特征为主，断裂应变更大，初始阶段产生的结构缺陷影响不大，甚至无影响。这类NW 的断裂分布主要位于两端，但初始滑移面的分布要广的多，几乎整个NW上都有分布，故无法确定两者之间有明确的因果关系。当中空缺陷足够大时，其影响得以显现，从断裂分布看更集中在中空位置，而初始的滑移面分布也位于中空球的两侧，所以对于中空半径较大的NW，初始滑移面的位置分布与最终断裂位置分布之间有明确的因果关系，而中等大小的R2 和R3 则表现出了过渡特征。

2.5 塑性形变断裂与脆性形变断裂的微观结构分析

进一步分析微观结构，以更好地理解初始滑移分布与断裂位置分布之间的关系。中空缺陷较小时，其对NW 的拉伸断裂影响轻微，初始阶段产生的滑移面的统计分布与最终的断裂分布之间也无明显的关联，这一点也可以从不同拉伸时刻的微观原子排布结构中体现。图9a 给出了具有塑性形变断裂特征的代表性模型(R2)。因中空半径小其对屈服点处的初始滑移面无影响，从图中也可以看出多个滑移面随机在NW 侧面靠近两端处产生，进入应力释放区大量滑移面产生，其中有些新的滑移面也产生自NW的中部，这一过程对应hcp原子数的急剧升高阶段。进入到塑性形变阶段，随着拉伸的进行，部分滑移面消失，又有新的滑移面产生，因此形成了hcp原子数的波动，也对应了应力和原子平均能量的波动。在应变0.45 处，滑移面的连续产生并发展到固定层，固定层对其反射，反射后的滑移面彼此相互阻碍产生局域非晶态，加剧了塑性形变，因此产生颈缩。颈缩之外的部分应变减小且通过恒温浴的作用原子的晶态得以部分恢复，导致hcp原子数的降低，进一步拉伸导致NW 在靠近一端断裂，这也与先前研究的断裂特征一致[43-44]。

图9 不同应变下的微观原子排布结构Fig.9 Microstructural snapshots at different strains

R4 与R2 的塑性形变断裂不同，它表现了脆性的断裂特征。从图9b可以看出，较大中空结构连同侧棱诱导产生了系列滑移面，它们可以从不同的位置产生，沿{111}面发展，在应力的释放阶段，虽然NW 的其它位置也有产生新的滑移面，但主要的贡献仍来自中空结构，中空附近高密度的滑移面产生局域非晶态，并迅速形成颈缩。在应变0.35 处中空一侧已形成缝隙，但另一侧尚有粘连。在随后的短暂拉伸过程中完全断开，从原子排布结构中可以明显看出中空结构诱导产生了系列初始滑移，并在随后的应力释放阶段围绕着中空部位产生了更多的滑移面，这些都在后续的短暂拉伸中得以保留，直至形成颈缩断裂。因此，较大中空NW 脆性形变特征明显，初始滑移分布与最终断裂分布存在明确的相关性。

3 结 论

基于MD 模拟所获得的多样本数据，我们研究了具有潜在功能性的中空铜NW 在负载条件下的形变断裂的统计分布特征，并追踪了在屈服点附近初始滑移面的位置分布。对于塑性形变断裂的单晶，小中空的NW 影响断裂位置分布的因素主要集中在塑性形变区，说明初始滑移分布与断裂位置分布间无相关性，但对于脆性断裂特征的大中空NW，由中空结构诱导产生的滑移面迅速积累、产生颈缩并导致最终的断裂，因此初始的滑移面位置分布与最终的断裂位置分布之间有明确的因果关系，这在机械制造中能准确地预测和控制材料的断裂行为，以减少断裂风险或增强抗断裂性能，在应力集中处及时对结构缺陷处加强，使得应力在材料结构内尽可能地被分散，以减小材料断裂失效的风险；另外，也可在纳米器件设计时，进一步优化设计原则，为延长器件使用寿命奠定基础。我们将材料学与计算科学建立联系，与实验方法相比，计算机模拟可以降低时间成本，避免一系列的繁琐实验。个体NW 不足以总结分析断裂失效的规律，而在这种计算模拟下能生成大量数据样本，这种基于大规模样本的分析更具有统计意义。