2. 分子动力学基础回顾

裂纹扩展的分子动力学模拟，本质上是在原子尺度上求解“原子受什么力、如何运动、在什么约束下运动”这三个问题。本节不追求把分子动力学（MD）的全部理论铺开，而是围绕裂纹模拟真正会用到的部分，梳理出一套最小但自洽的知识集：原子间势函数决定原子之间的相互作用力，运动方程决定原子如何随时间演化，系综与边界条件决定模拟所处的物理环境，时间步长与积分算法决定这套演化能否稳定、准确地被数值求解。理解这四者的分工与耦合关系，是后续搭建裂纹模型、设置加载方式和判读结果的前提。

一、原子间势函数：力的来源

在MD中，原子不被当作连续介质，而是离散质点。原子所受的力来自势函数对位置的负梯度：F_i = −∂U/∂r_i，其中U是体系总势能，r_i是第i个原子的位置。因此势函数是整条模拟链路的起点，它决定了材料在原子尺度上的弹性、塑性、断裂等行为是否可信。

裂纹模拟中常用的势函数大致分三类。第一类是对势，如Lennard-Jones势，形式简单，适合验证算法和定性观察，但无法正确描述金属的弹性常数与空位形成能。第二类是多体势，如嵌入原子法（EAM）势，它把每个原子看作嵌入在其余原子提供的电子密度中，能较好再现金属的结合能、层错能和表面能，是金属裂纹模拟的主流选择。第三类是键级势，如Tersoff、Stillinger-Weber及反应力场ReaxFF，适合共价材料（硅、碳化硅）以及需要考虑化学键断裂与重组的情形。

选择势函数时，不能只看“能不能跑”，而要看它是否复现了与裂纹相关的关键量：弹性常数、表面能、层错能、理想断裂强度。例如用EAM模拟铜的Ⅰ型裂纹，需要确认该势给出的表面能与实验值接近，否则裂纹扩展的临界应力会明显偏离。势函数还带有截断半径r_c，超过r_c的相互作用被忽略，这直接影响计算量与近程相互作用的完整性，通常取2.5σ至5.5σ（对势）或覆盖第三近邻（EAM）。

二、运动方程与积分算法：原子如何走

原子运动遵循牛顿第二定律：m_i·d²r_i/dt² = F_i。这是一个二阶常微分方程组，MD用有限差分逐步推进求解。最常用的是Velocity-Verlet算法，其更新式为：

r(t+Δt) = r(t) + v(t)·Δt + 0.5·a(t)·Δt²
v(t+Δt) = v(t) + 0.5·[a(t) + a(t+Δt)]·Δt

该算法时间可逆、能量守恒性好，是裂纹模拟的默认选择。另一种常见形式是蛙跳法（Leap-frog），速度与位置错开半个时间步，存储量小，但速度与位置不同步，输出瞬时动能时需注意。积分算法的核心要求是辛性与稳定性：它不追求单步精度极高，而追求长时间演化下能量不漂移，否则裂纹尖端局部升温会失真。

三、系综与边界条件：模拟所处的环境

系综决定哪些宏观量被固定。裂纹模拟常用三种：NVE（粒子数、体积、能量守恒）用于观察裂纹扩展的动力学过程；NVT（粒子数、体积、温度固定）通过Nosé-Hoover或Berendsen热浴控制温度，用于研究温度对断裂的影响；NPT（粒子数、压力、温度固定）用于弛豫阶段，让体系在目标温度和压力下达到平衡。加载裂纹时，若采用应变控制，通常先NPT弛豫，再切换到NVE或NVT进行加载。

边界条件同样关键。裂纹模拟中，裂纹尖端附近必须用自由边界或固定边界来避免周期性镜像对裂纹的干扰；而远离裂纹的区域常用周期性边界条件（PBC）来模拟无限大材料。一种常见做法是：在裂纹扩展方向采用自由表面，在厚度方向采用周期性边界，在加载方向采用固定或柔性边界。若全用PBC，裂纹会与自身镜像相互作用，导致非物理的裂纹闭合或加速扩展。

四、时间步长：精度与代价的平衡

时间步长Δt必须小于体系中最快振动周期的约1/10。对金属，原子振动周期约10⁻¹³ s，因此Δt通常取1～5 fs。Δt过大，积分发散，能量爆炸；Δt过小，计算量成倍增加。裂纹尖端在扩展瞬间原子速度可能显著升高，因此建议在加载阶段用1 fs，弛豫阶段可用2～5 fs。判断Δt是否合适，可观察总能量漂移：在NVE下，若总能量在10⁵步内漂移超过0.1%，应减小Δt。

五、四要素的耦合关系

这四者不是孤立的：势函数决定力的量级和振动频率，进而约束Δt的上限；系综决定温度与压力的控制方式，影响边界条件的选取；边界条件又反过来影响势函数截断处的处理。下图给出裂纹模拟中这四要素的耦合关系与典型参数区间，可作为搭建模型时的速查参考。

六、面向裂纹模拟的最小参数清单

把上述知识落到一次具体的裂纹模拟上，可以整理成一份可直接对照的参数清单。以下示例以铜的EAM势、Ⅰ型裂纹为例，参数为教学示例值，实际使用时应根据势函数文件和目标温度调整。