LAMMPS计算揭示纳米晶铜拉伸变形：位错演化与晶界行为的原子尺度分析

LAMMPS计算在这个纳米晶力学项目中扮演着不可替代的角色——它需要在原子尺度上还原纳米晶铜在单轴拉伸下的塑性变形全过程。项目背景源于一个工程痛点：电沉积法制备的纳米晶铜薄膜在微电子互连中表现出与传统粗晶铜截然不同的力学行为，其强化机制到底是晶界强化还是位错强化，实验手段难以直接区分。这个项目的思路很直接——用分子动力学模拟把变形过程中每一个原子的运动轨迹记录下来，再从轨迹中提取位错的萌生、滑移和交互作用信息。

势函数的选择是这个LAMMPS计算项目里最需要审慎对待的决定。项目对比了三种广泛使用的Cu EAM势函数：Mishin（2001）、Foiles（1986）和Zhou（2004）。三者对Cu晶格常数的预测差异不到0.5%，但对层错能（SFE）的预测差异显著——Mishin势给出约45 mJ/m²，与实验值约45 mJ/m²高度一致，而Zhou势预测值偏低约30%。由于层错能直接控制不全位错之间的平衡间距，进而影响交滑移概率，项目最终认定Mishin势更适合这个场景。势函数验证过程参照了Sandia国家实验室发布的LAMMPS官方文档中推荐的测试流程（https://en.wikipedia.org/wiki/LAMMPS）。

模型构建采用Voronoi tessellation算法生成多晶结构，晶粒尺寸设为8 nm，共包含约120,000个原子，系统盒子尺寸约24×24×24 nm³。这个尺寸经过权衡——更大的系统在计算资源允许的范围内能减少周期边界对晶界网络统计性的影响，但120,000原子已是单节点128核并行运行的实际上限。弛豫阶段先在300 K下用NPT系综跑20 ps，压力耦合采用Nose-Hoover方法，时间步长1 fs，阻尼参数100倍时间步长。弛豫后系统总能量的漂移小于10⁻⁵ eV/ps，判定为达到热力学平衡。

拉伸加载沿z轴方向施加，应变率设为1×10⁹ s⁻¹。这个应变率远高于实验中的准静态拉伸（约10⁻³ s⁻¹），是MD模拟的固有局限——受限于时间步长，MD能覆盖的时间尺度通常在纳秒量级。项目在分析中明确标注了这一点，并通过对比不同应变率（5×10⁸、1×10⁹、5×10⁹ s⁻¹）下的应力-应变曲线，确认了在1×10⁹ s⁻¹下屈服应力的应变率敏感因子m < 0.05，变形机制的主导特征未发生质变。

应力-应变曲线呈现出三个特征阶段：弹性段（应变<0.02）、屈服点（应力峰值约2.1 GPa）和软化段。屈服后应力从2.1 GPa下降至1.6 GPa并进入流动平台。这个2.1 GPa的屈服强度与文献中纳米晶铜（晶粒8 nm）的实验硬度数据（约2.0 GPa，经Hertz接触模型换算）在合理范围内吻合（https://www.nature.com/articles/nmat2182）。差距不会说谎——约5%的偏差部分来源于模拟中理想化的晶界结构，实验样品中不可避免存在杂质和残余应力。

位错分析借助OVITO软件中的DXA（Dislocation Extraction Algorithm）模块完成。在应变达到0.025时，DXA首次识别出从晶界发射的Shockley不全位错，柏氏矢量1/6⟨112⟩。随着应变增加，位错密度从初始的几乎为零攀升至峰值约3.2×10¹⁶ m⁻²，随后在应变0.06左右因晶界吸收和位错反应而下降。值得注意的是，完美位错（柏氏矢量1/2⟨110⟩）在整个变形过程中占比不足15%，这与纳米晶中晶界间距短、位错来不及分解就湮灭的物理图像一致。与此同时，共格孪晶界在应变0.04时开始出现局部脱粘——原子结构快照显示孪晶界面上形成了约0.5 nm宽的无序层。这一发现印证了晶界滑动与位错滑移协同主导纳米晶铜塑性变形的机制图像。

回过头审视，LAMMPS计算揭示的核心结论是：当晶粒尺寸降至8 nm以下，位错介导的塑性逐渐让位于晶界介导的塑性，两者在过渡区间存在竞争关系。这一竞争边界的量化定位，为纳米晶金属的晶粒尺寸工程优化提供了直接的原子尺度依据。方法局限需要坦承——经典势函数无法描述电子结构变化带来的化学效应，对于含合金元素的体系需要谨慎评估适用性。此外，1×10⁹ s⁻¹的应变率虽然在该体系中未引起机制质变，但定量外推至实验条件仍需要多尺度建模的桥接工作。这些边界条件确实值得后续项目保持警醒。