分子动力学模拟是一种基于分子/原子的经典力学模型，通过数值求解体系差分方程研究材料特性的计算机模拟方法，目前已成功应用于航空发动机涡轮叶片用镍基单晶高温合金力学行为仿真，在揭示变形与损伤机理、预测力学性能等方面展现出巨大的潜力。

高压涡轮叶片是航空发动机的重要零部件，在服役过程中承受着极端热力耦合载荷，可能出现掉块、断裂或丢失等现象，导致发动机无法正常运转，甚至引起重大安全性事故。镍基单晶高温合金由于其优异的高温力学性能被广泛应用于制造涡轮叶片，对镍基单晶高温合金力学行为的深刻认识和理解是突破涡轮叶片安全性设计的关键基础。分子动力学（MD）模拟是建立在原子尺度、以原子或分子为基本单元的模拟方法，通过数值求解分子体系经典力学运动方程来计算其热力学特征与其他宏观性质。MD模拟可以从微观尺度研究镍基单晶高温合金力学行为，深入揭示材料变形与损伤的微观机理，并对试验手段难以获得的材料宏观力学性能进行外延预测，可为镍基单晶高温合金涡轮叶片安全性设计提供理论基础与基本数据。MD模拟涡轮叶片力学行为如图1所示。

MD模拟原理

MD理论认为原子在所有其他原子核和电子提供的势场作用下按牛顿定律运动，该势场可以用一个基于位置的半经验势函数进行描述^[1]。MD模拟通常使用玻恩-奥本海默（Born-Oppenheimer）近似方法简化原子核和电子的运动，并通过密度泛函理论（DFT）简化电子的分布，在此基础上结合第一性原理计算、试验数据或经验公式建立势函数。针对不同的研究对象有不同的势函数模型，目前研究金属晶体一般采用嵌入原子（EAM）势，MD模拟方法基本步骤如下。

首先，针对研究对象建立合适的原子模型，给定材料系统内所有原子的种类、位置与速度，定义系统的初始条件和边界条件，确定模拟系综（ensemble）。所谓系综是指在一定的温度、压力（压强）、载荷等宏观条件下，大量性质和结构完全相同的、处于各种运动状态的相互独立的系统的集合。在研究镍基单晶高温合金力学行为时一般采用等温等压（NPT）系综，即在模拟过程中保持总粒子数、压强与温度不变。

其次，根据原子间的势函数，基于原子的位置计算每个原子受到的合力，并基于经典力学计算原子的加速度，以此计算下一个时间步长的速度与位置。MD模拟在某一时刻原子的受力状态由势函数的导数决定，原子经过∆t的时间步长后的运动状态可以通过差分方法解牛顿第二定律的微分方程获得，如果∆t足够小，那么离散形式可以以足够的精度近似连续形式。

最后，重复上一步直到获得一系列的原子位置与速度的信息。通过分析体系中各原子的状态与相互作用，可以厘清材料的变形损伤机理，预测材料的宏观力学性质及行为。

MD模拟作为一种研究分子/原子体系结构与性质的重要方法，被广泛应用于材料科学领域，既可以研究材料的力学性质，也可以研究动能、势能、焓与体系自由能等热力学性质。因此，采用MD方法建立镍基单晶高温合金体系模型，通过计算镍原子与其他合金元素的相互作用，可以模拟镍基单晶高温合金在不同温度、载荷条件下的力学响应，支撑涡轮叶片的安全性设计与评估。

MD模拟常用工具

目前，针对金属材料的MD模拟一般采用LAMMPS软件，该软件可以根据原子间势函数和边界条件计算小至几个粒子，大到上百万甚至是上亿个粒子的体系，能够在大多数条件下完成MD模拟的建模与计算部分。LAMMPS软件的建模功能相对简单，目前普遍使用Atomsk软件进行金属材料建模，能够简便地实现位错、孪晶等缺陷的插入操作。MD模拟结果后处理一般可采用OVITO软件，可以分析MD模拟结果，对位错、相组成等重要性质进行可视化。此外，还有许多其他软件可用于MD模拟的建模、计算以及后处理。

MD模拟镍基单晶高温合金的关键技术

原子模型与边界条件

镍基单晶高温合金具有优良的蠕变强度、热疲劳强度、抗氧化性和抗热腐蚀性能，主要由γ相与γ'相构成。其中γ相是以面心立方结构的单晶Ni为主体的基体相，晶格常数为3.52Α；γ'相是以L12结构的Ni₃Al为主的强化相，以立方状沉淀相的形式周期排布在基体相中，晶格常数为3.57Α（1Α=0.1nm）。γ'相与基体γ相成立方-立方取向关系，即γ'相的晶向与γ相基体一致，但晶格常数不一致，这使得γ'相与基体γ相间存在错配度δ。错配度会使γ/γ'相界面中产生错配位错并形成位错网络，改变位错在γ通道移动的方式，提高镍基单晶高温合金的强度。

根据研究问题建立合理的镍基单晶高温合金原子模型是MD模拟的关键，镍基单晶高温合金中的γ'相分布在γ相之中，尺寸约为0.4μm，与Ni₃Al的晶格常数相差3个数量级，要想模拟一个宏观尺寸的立方体γ'相，需要模拟接近一亿个原子，对计算资源要求极高。随着算力与算法的发展，MD模拟的空间尺度限制有望突破到介观乃至于宏观，但由于计算成本限制，目前研究的原子数通常在百万量级，需要对原子模型进行简化与缩尺，并使用周期性边界条件减小模型简化引起的尺寸效应。周期性边界条件等价于在模拟区域周围放置与模拟区域完全相同的镜像，在计算边界附近原子的受力时，不仅要考虑本身模拟区域内附近的原子，还要将最近的周期性镜像区域附近的原子也考虑在内。

MD模拟镍基单晶高温合金常用的原子模型如图2所示。两相模型中Ni与Ni₃Al的晶格常数不同，在弛豫后可以自发形成位错网络。该模型已成功应用于模拟γ/γ'相界面行为，包括镍基单晶高温合金在不同应力条件下γ/γ'相界面处位错网络的损伤机制、应力应变关系以及γ/γ'相界面的错配位错运动与演化等，证明了应力状态会影响位错网络的破坏方式，从而改变镍基单晶高温合金的筏化形态和力学性能。但两相模型并未考虑γ'相的排列形式与两相尺寸差异，无法在z轴方向使用周期性边界以消除边界效应。

三明治模型在两相模型的基础上添加了一层γ相，在z轴方向上更接近真实的γ相与γ'相的排列方式，可以使用周期性边界对z轴方向的γ/γ'相界面位错网络以及力学性能进行更准确的模拟。该模型已成功应用于强化元素对镍基单晶高温合金切削行为的影响研究，证明了添加Cr、Co有助于位错缠结的产生，可以提高γ相与界面强度，提升合金的整体抗变形能力，同时证明了错配度越低，位错网络越规则，能有效防止位错运动、提高材料强度。但三明治模型只能较为准确地模拟与z轴平行的镍基单晶高温合金相界面，无法在其他方向反映γ相与γ'相的结构关系。

立方体模型一般以一个立方体γ'相与其周边的γ相基体组成，结构最接近于单晶胞元，在弛豫后可以形成六面体的位错网络，可以在三维空间内较为准确地反应镍基单晶高温合金的结构特点。立方体模型虽然也难以建模至宏观尺度，但可以使错配度与真实值基本相同，以此更为准确地研究γ/γ'相界面位错网络的性质，已经成功应用于氢元素对镍基单晶高温合金界面损伤的影响，证明了氢原子能促进界面位错段的相互反应与解离，加剧界面位错网络的后续破坏，为镍基单晶高温合金氢脆现象提供了新的认识。

加载方式

涡轮叶片长期在强时变、大梯度温度和应力场，以及高温氧化腐蚀、侵蚀环境中服役，面临疲劳、蠕变、高温腐蚀、冲击损伤等多种失效模式威胁，可能出现掉块、断裂或丢失等现象，在进行MD模拟时需要根据研究场景选择合适的加载方式。

基于MD模拟研究疲劳、蠕变等准静态力学行为时一般采用如下方式加载：首先根据加载条件令模拟区域产生变形，然后将模型中的原子映射到新的模拟区域当中，再计算各个原子的状态，如位置、势能等参数，以此模拟拉伸、剪切等力学行为。通过这种加载方式可以观察到位错延滑移系移动等塑性行为，能够更为准确地从MD角度模拟镍基单晶高温合金的宏观准静态力学行为。

MD模拟步长限制在皮秒量级，适合进行镍基单晶高温合金的高速冲击行为模拟，目前常见的加载方式分为3类，如图3所示。

小球冲击加载使用Ni或其他材料组成的小球撞击镍基单晶高温合金，可采用两相模型、“三明治”模型或立方体模型，定性研究不同晶向的冲击行为。该加载方式直观模拟冲击碰撞，但模拟冲击载荷的小球尺寸为纳米尺度，目前尚未建立纳米尺度碰撞与宏观冲击行为之间的内在关联。

动量镜（momentum mirror technique）加载使用一排具有极高初速度的边界原子作为活塞，以此模拟冲击载荷的加载，并在其余方向使用周期边界，使材料内部产生稳定推进的冲击波。此方法在卸载时直接移除活塞，使其不再提供冲击载荷。这种加载方式可以在更大的空间尺度上研究冲击响应并减少尺寸效应的影响，从微观角度验证镍基单晶高温合金的冲击行为，分析其在冲击载荷作用下的显微组织演化特征。

吸收边界法的加载方式与动量镜方法近似，区别在于卸载时在模型右端额外设置活塞，使冲击波传播至右端活塞处使该活塞获得与加载活塞相同的速度以模拟卸载过程，减小了冲击方向上的边界效应。吸收边界法已成功应用于模拟激光冲击过程，再现了弹性波与塑性波的分离现象，从微观角度验证了宏观冲击行为的特点，为镍基单晶高温合金冲击行为研究提供MD解释。

模拟结果分析

镍基单晶高温合金涡轮叶片安全性设计需要深入掌握材料在不同载荷条件下的变形损伤机理与力学响应数据。通过分析MD模拟中各原子的状态与相互作用，可以计算材料参数、总结材料变形损伤的机理、得到力学响应、预测材料的宏观力学性质及行为。

MD模拟可以根据体系中各原子运动的统计分析推测体系的各种性质：通过分析原子的位置和相互作用，可以推测体系的可能构型；通过体系的动能、位能以及体系与外界的相互作用，可以计算温度、压力和能量等参数；基于模拟得到的原子轨迹，通过统计力学计算，可以得到体系的平均能量、熵和自由能等。

MD模拟可以获得原子的位置随时间的变化轨迹，通过分析原子总体的位移、应变与位置关系，可以得到材料变形过程的应力分布与缺陷演化，通过分析晶体缺陷的演化过程，可以了解材料的损伤机理。国外学者使用改进的立方体模型模拟了1223 K下镍基单晶高温合金的拉伸与剪切行为，再现了其中位错在γ'相处钉扎、γ'析出相处弯出、γ'析出相中超位错堆积和成核等现象^[2]，揭示了微结构变化对组织演化以及宏观力学性能的影响机理，为理解位错演化规律提供了MD视角。

此外，MD模拟通过计算应变张量、应力张量和原子位移向量可以求得材料的应变分布和应力状态，得到力学响应曲线，为真实材料性质提供理论支撑。国内研究者使用MD方法，通过周期变化的载荷和温度场模拟了同相（IP）与反相（OP）作用下立方体模型镍基单晶高温合金的热机械疲劳（TMF）行为，得到了循环应力-应变曲线^[3]（见图4）。虽然此研究的模拟结果与试验数据在数值上存在一定差异，但模拟所得的曲线形状及其变化规律与实测结果呈现良好的一致性，并且通过模拟数据的变化趋势与规律，可以预测热力耦合循环载荷下镍基单晶高温合金的力学行为特征，揭示了TMF载荷作用下温度、位错演变、加载模式等因素对镍基单晶高温合金力学响应的影响。有学者使用动量镜加载模拟镍基单晶高温合金冲击加载过程，分析了带孔洞镍基单晶高温合金在冲击载荷作用下的显微组织演化特征，结合位错密度变化，证明了高温环境会激活位错爬升和交叉滑移，从微观角度验证了镍基单晶高温合金冲击行为中的高温软化机制^[4]。

结束语

MD模拟经过几十年的发展，已经成为了一个功能强大、应用广泛的研究方法，但其以原子为模拟基本单位，在模拟时受到时间尺度、空间尺度等方面的限制。随着超动力学（Hyperdynamics）与超级计算机的发展，MD方法已经可以模拟秒尺度的过程与介观尺寸的模型，再加上近几年神经网络在各个领域的发展，使用机器学习原子间势（ML-IAP）也可以大幅增加模拟时间与空间尺度的数量级，提高模拟结果的准确性。

目前，MD模拟方法已经被逐渐应用于研究镍基单晶高温合金在不同温度和载荷条件下的力学性能，成功地再现宏观试验趋势、揭示微结构变化机理，比较准确地反映镍基单晶高温合金力学行为特征，为基于现有试验结果预测极端条件下的镍基单晶高温合金力学行为提供方法和工具支撑。

（张宇宁，清华大学，硕士研究生，主要从事航空发动机高温合金本构行为研究）

参考文献

[1] VLACHAKIS D, BENCUROVA E, PAPANGELOPOULOS N, et al. Current state-of-the-art molecular dynamics methods and applications[J]. Advances in protein chemistry and structural biology, 2014,94: 269-313.

[2] YASHIRO K, NAITO M, TOMITA Y. Molecular dynamics simulation of dislocation nucleation and motion at γ/γ' interface in Ni-based superalloy[J]. International Journal of Mechanical Sciences, 2002, 44(9): 1845-1860.

[3] WU W P, DING Z J, Li Y L, et al. Molecular dynamics simulation of thermomechanical fatigue properties of Ni-based single crystal superalloys[J]. International Journal of Fatigue, 2023, 173: 107667.

[4] CHEN B, LI Y, SOPU D, et al. Molecular dynamics study of shock-induced deformation phenomena and spallation failure in Ni-based single crystal superalloys[J]. International Journal of Plasticity, 2023, 162: 103539.