分子动力学模拟拉伸 — 单轴拉伸应力-应变曲线的原子级获取

分子动力学模拟拉伸的结果和宏观拉伸实验看起来相似——都是一条先升后降的应力-应变曲线——但产生这条曲线背后的物理完全不同。宏观拉伸是准静态的（应变率约10⁻³/s）、多晶或多相的集体行为；MD拉伸是超高应变率的（10⁸-10⁹/s）、通常是单晶或纳米多晶的。这之间的应变率差约11个数量级——从MD预测的屈服强度要转换到宏观拉伸的参考值，需要倍率修正。本文把分子动力学模拟拉伸的完整参数设置和分析方法拆解清楚。

一、应变率的合理选择与实际局限

分子动力学模拟拉伸的第一个设置是应变率。MD模拟的时间窗口是纳秒级（10⁻⁹ s），你如果要在10⁹/s应变率下模拟10%的应变——跑0.1 ns（10⁻¹⁰ s）就能完成。但在实验室——10/s就是快速拉伸了，10⁻³/s是标准准静态。MD比实验快了约12个数量级——这个差距不是数值误差，是物理机制的区别。

应变率对屈服强度的影响可以用Johnson-Cook或Zerilli-Armstrong型本构来描述：σ_y = σ_y0×(ε_dot/ε_dot0)^(m)，其中ε_dot0通常取1/s作为参考应变率，m为应变率敏感指数。对于FCC金属（如铜），m约0.01-0.02——应变率从10⁸跳到10⁻³，屈服强度可能下降约15-25%。对于BCC金属（如铁），m约0.04-0.06——同样的跳跃对应约40-60%的屈服下降。BCC金属的强应变率依赖源于螺位错的非热激活滑移阻力——在高应变率下螺位错没有时间热激活跨过Peierls势垒。

实际操作中：在可行的时间窗口内选尽量低的应变率。LAMMPS中`fix deform`可以设拉伸速率。10⁸/s是一个合理的起始点（对10⁸/s、总应变20%在2 ns内完成——正常MD的产出窗口）。如果计算资源允许，可以降一两个数量级（10⁷/s）做对比，观察屈服强度的应变率收敛趋势。做3-5个不同应变率的拉伸，拟合σ_y vs log(ε_dot)的线性关系，外推到低应变率——但外推超过4个量级风险较大。

二、拉伸方向和系综

分子动力学模拟拉伸的加载方式必须匹配材料的各向异性和模拟目标。

单晶拉伸：需要沿特定晶向拉伸。FCC金属沿<100>拉伸——初始阶段弹性变形+位错在{111}面上启动，屈服后阶段滑移在多套滑移系上激活。沿<111>拉伸——早阶段位错难以启动（Schmid因子低），应力会冲得非常高，然后突然滑移→应力崩塌→新的滑移系启动→应力波动——现象丰富但对势函数的预测能力有更高要求。

系综选择：拉伸方向用”拉伸”而非NPT中的”压力控制”——因为拉伸是一种受控变形，而不是受控应力。LAMMPS中用`fix npt`在拉伸方向取消压力控制（用`NULL`），改为`fix deform`控制长度变化。非拉伸方向保持压力控制（P_yy=P_zz=1 atm）——模拟单轴拉伸的实际条件——即试样侧面不受约束。

盒子尺寸要足够大以避免周期性边界对位错行为的人工约束。经验：拉伸方向的初始尺寸至少100 nm（约300个晶格常数），截面至少20×20 nm。小于这个尺寸——位错在两个表面（周期性边界等效）之间来回弹跳——产生人为的往复强化-软化循环。

三、应力计算

分子动力学模拟拉伸中的应力来自维里应力（Virial Stress）。每个原子的应力张量σ_i^αβ = -(1/V_i)×[m×v_α×v_β + (1/2)×Σ_j r_ij^α×F_ij^β]，其中V_i是原子i的原子体积。应力σ的平均在整个模拟盒子上的统计就是宏观应力（单位GPa）。

维里应力在MD模拟中有两个处理细节：第一，应力是张量——从六个独立分量（σ_xx, σ_yy, σ_zz, σ_xy, σ_yz, σ_zx）中提取拉伸方向的法向应力分量（如单轴沿x方向拉伸取σ_xx）。第二，对于拉伸中的弹性段——应变ε_x=ln(L(t)/L_0)（真实应变或对数应变），应力σ_x=F/A(t)（真实应力，当前截面上）。

屈服应力的提取标准：0.2%偏移法（offset method）——弹性段的线性拟合线和应变偏移0.2%后的平行线与应力-应变曲线的交点。这个提取标准的物理基础是0.2%残余塑性应变作为屈服判据。

四、屈服和断裂机制的识别

分子动力学模拟拉伸中屈服后的微观机制需要用CNA和DXA等缺陷分析工具鉴别。FCC金属屈服通常对应第一个Shockley分位错（partial dislocation）在{111}面上的形核和扩展——表现为CNA中HCP原子层的出现（因为Shockley分位错将FCC堆垛变成HCP堆垛的层错）。HCP原子的数量随时间曲线与应力-应变曲线有很好的对应：屈服时HCP原子数快速上升→加工硬化阶段HCP原子数平台的建立→颈缩阶段表面原子的自由表面迁移。

断裂（fracture）在MD中是原子键的大规模破坏——通过原子间距离超过势函数的截断半径判断是否断裂。裂口的形核位置通常是：自由表面（应力集中的表面缺陷处）、晶界（晶界滑动引起的微空洞）、或位错胞壁（大量位错累积导致的局部应力释放）。在拉伸曲线中，断裂对应应力从峰值快速下降到接近零的区段。

五、从纳米到微米：尺度平移

分子动力学模拟拉伸的最终目标往往是将MD预测的屈服强度引用到宏观（微米以上尺度）结构分析中。由于MD的应变率效应和尺寸效应（纳米级的试样强度>微米级），MD直接输出的屈服强度不能作为宏观材料的屈服强度。

从MD预测到宏观的平移：做纳米多晶MD（晶粒尺寸5-50 nm），提取Hall-Petch关系σ_y=σ_0 + k_y/√d中的σ_0和k_y。当d拉到亚微米到微米级（超出MD直接模拟范围但可以通过Hall-Petch外推），可以得到对应晶粒尺寸的屈服强度。外推时应使用Hall-Petch形式（1/√d项），因为晶粒尺寸在<100 nm时，反Hall-Petch效应（位错在极细晶粒中无法容纳和增殖——变形由晶界滑移主导）在约d<10-20 nm时出现——这个非单调行为在MD中可以看到但Hall-Petch外推会忽略。

六、专业分子动力学模拟拉伸服务

 需要分子动力学模拟拉伸分析服务？

科研学术网提供专业的MD拉伸模拟分析服务：

– ✅ 全参数设置：应变率/温度/晶向/系综的定制化配置

– ✅ 微观机制解析：位错形核与演化、孪晶、相变的MD级可追溯路径

– ✅ 多温度多应变率扫描：提取温度-应变率-屈服强度的四维数据集

– ✅ 尺度平移：从纳米MD→微米晶体塑性的参数桥接

立即咨询报价 →