高通量应力状态图谱揭示单晶铁bcc-to-hcp转变的加载路径依赖性

《Computational Materials Science》：A high-throughput stress-state map of the bcc-to-hcp transformation in single-crystal iron from molecular dynamics

　　

采用分子动力学（MD）方法，在连续应变分配网格上对单晶铁（Fe）中bcc-to-hcp转变进行了高通量应力状态图谱研究。研究构建了221条LAMMPS轨迹，覆盖[100]和[110]两个晶体取向。结果表明，操作屈服压力（P_Y）保持在较窄区间（15.3–26.2 GPa），而操作转变压力（P_T）跨度超过4倍（23.8–106.9 GPa）。P_Y对应此前在该势函数上识别的弹性-塑性 onset，塑性先于转变发生；在整个网格上其变化远小于P_T。转变事件处的微观结构和偏应力（σ_eq）随加载路径变化。基于转变体积（V_T）和转变时刻等效应力（σ_eq|T）的梯度提升代理模型可以R²=0.99（RMSE 2.1 GPa）恢复P_T；高斯过程（GP）模型在相同输入下给出校准的不确定度。有限速率下铁中的转变取决于应力状态和取向，而非单一压力。文献静态平衡参考值约为13 GPa，而本目录中的操作P_T标签跨度为23.8–106.9 GPa。该区间下端与斜坡压缩窗口（15–38 GPa）重叠，并延伸至强过压的单轴路径。该网格上的近静水压情形在10⁹ s^?1应变率下使用修正Ackland势给出P_T≈38 GPa。

**研究背景与问题**

铁在高压下从体心立方（bcc，α相）到六方密排（hcp，ε相）的转变是凝聚态物理、地球物理和冲击动力学领域的经典问题。自1956年Bancroft等人在冲击加载铁中首次识别多形性事件以来，该转变已被静态压缩、冲击加载和原子模拟广泛研究。静态X射线衍射确定室温平衡边界约在13 GPa附近，并绘制了其对温度、应力状态和合金化的敏感性。该转变为马氏体相变——无扩散、具有明确晶格对应关系和原子面 shuffle。

然而，动态加载下的转变过程远比静态相图复杂。冲击铁的波剖面测量显示显著的动力学和滞后特征；斜坡压缩实验报告α→ε onset应力约15–38 GPa，取决于应变率和加载时间，相对平衡边界存在系统过压。超快X射线衍射证实了纳秒冲击加载中的结构变化，并分辨了激光冲击铁中的亚纳秒三波剖面和卸载路径。接近~10⁹ s^?1速率的实验以及斜坡和动态金刚石压砧动力学进一步强调了速率和路径效应。实践中，压力诱导相变取决于加载几何和速率，其程度不亚于热力学平衡。

原子尺度分子动力学提供了上述唯象学的微观对应。Kadau等人用数百万原子模拟表明冲击诱导密排晶粒在皮秒时间尺度形核，阈值以上出现分裂的弹性-转变波结构，且微观结构和转变压力强烈依赖晶体学冲击方向。相关多晶MD将模拟衍射与实验联系，表明晶界塑性常先于α–ε事件。平行研究强调塑性与转变不是先后关系而是耦合过程：孪生和位错活动可根据加载协议抑制或促进形核。纳米晶和缺陷样品进一步表明预先存在的晶界和空洞同时移动屈服和转变阈值。

定量MD比较需要原子间势能同时再现平衡α–ε体积塌缩和高应变率塑性响应。本研究使用拟合到13.75 GPa附近bcc–hcp转变的修正Ackland嵌入原子势。该势下纯Fe的静水压压缩在10 K和10⁹ s^?1下已将第一压力峰解析为26.7 GPa的弹性-塑性极限，随后在39.3 GPa发生α→ε事件。位错提取显示密排界面处存在部分位错，故转变前有塑性先行。跨实验和模拟，onset应力与应力空间中加载路径的关联至少与单一平衡压力同等强烈。

这一路径依赖性正是本研究的动机。静水压在主应力轴置换下不变且对偏量分量不敏感；而马氏体相变具有晶格对应关系和类剪切畸变，直接耦合偏应力和晶体取向。研究人员在连续双参数网格上生成一族单晶MD模拟，对[100]和[110]取向重复扫描，并用相同规则从每条轨迹提取事件。在该数据集中，P_Y变化因子为1.7而P_T跨度超过4倍。转变 onset对加载路径和晶体学取向的依赖性此前已确立，但本文基于网格的分析突显了具有相应转变标签的事件定量分离于单一压力-体积分支。该纯Fe目录可推广至合金系统的程度及其局限性在第5.6节讨论。

**关键技术方法**

研究采用100×100×100常规晶胞（200万原子）单晶BCC铁，a₀=0.28665 nm，初始温度10 K，时间步长0.001 ps，目标体积工程压缩ε_v=0.235，周期性边界条件，LAMMPS结合修正Ackland嵌入原子势。两个单晶取向族（[100]族36例，[110]族185例）共享同一(θ,φ)应变网格但模拟盒轴嵌入晶体方式不同。压缩通过Cartesian轴施加，各轴应变率正比于应变分配比。每例以共轭梯度能量最小化加零压各向同性盒弛豫初始化，NPH（P=0）平衡后重启文件复用于同取向所有压缩案例以消除噪声。各向异性压缩在10 K下按每轴应变率施加。局部结构在OVITO中用共同邻位分析（CNA）后处理，原子分为BCC、HCP、FCC-like和"其他"四类。应力描述符将Cauchy应力张量分解为平均压力p和偏应力s_ij，报告P_T、P_Y、V和σ_eq≡√J₂。事件从压力-体积轨迹用两个非重叠体积窗口读取：屈服窗口7.0–6.4 cm³ mol^?1，转变窗口6.4–5.6 cm³ mol^?1；在每窗口内定位Savitzky–Golay平滑（宽度21，三阶）曲线上dP/dV最负特征，取紧邻其前的原始压力峰为事件压力。

**研究结果**

**轨迹原型**：近静水压验证案例（θ=0.955，φ=0.785）轨迹呈最简单特征——平滑单调压缩至约V=6.5 cm³ mol^?1处斜率短暂变化（屈服标签P_Y≈26.2 GPa），继续压缩后在V=6.2 cm³ mol^?1附近压力陡降标记为P_T。对应P_T对[100]为38.0 GPa、[110]为38.1 GPa，一致到0.1 GPa，符合立方晶体体积压缩预期。单轴案例则截然不同且两取向间有别：[100]单轴案例P_Y≈15.4 GPa后攀入高压亚稳臂直至P_T≈82.4 GPa结构塌缩；[110]案例θ=0（沿[001]单轴压缩）P_Y≈15.3 GPa几乎相同，P_T≈86.4 GPa且转变体积更小（V_T≈5.82 vs 5.86 cm³ mol^?1）。

**应变分配球面图谱**：将每个事件映射到应变分配球面，屈服压力球面在两取向上几乎单色，仅在静水纬度附近有轻微隆起——这是同势下静水压缩识别的弹性-塑性极限的应力状态扩展。转变压力球面则强烈织构化：亮区标记高P_T（长亚稳BCC路径），暗区标记早转变。两取向不同：[110]球面赤道附近有高P_T带达106.9 GPa，而[100]球面在[001]极附近有高P_T帽约89 GPa。该模式令人联想到Kadau等人的冲击方向图，但此处是在连续三轴网格而非三个离散冲击轴上采样。

**三轴加载路径下的屈服样和转变压力**：全部221条轨迹中操作屈服压力跨15.33–26.22 GPa（因子1.7），操作转变压力跨23.81–106.91 GPa（因子4.5）。[100]中位P_Y=17.20 GPa、[110]中位17.80 GPa；而[100]中位P_T=64.29 GPa、[110]中位41.21 GPa：屈服在两取向间近乎不变，转变则强烈取向依赖。屈服样压力为屈服体积窗口内首次强弛豫；鉴于同势先前工作证明塑性先于α→ε事件（静水弹性-塑性极限26.7 GPa），将其解释为塑性变形 onset。转变标签则与晶格对应关系绑定：BCC母相须通过协调shuffle找到剪切兼容的HCP产物；加载促进还是阻碍该shuffle取决于偏应力和加载轴相对晶体的取向。部分路径在温和压力下早达转变窗口，部分将BCC结构载至近107 GPa才让位。P_Y对这些路径差异相对不敏感，P_T则放大之。

**压力-体积事件与微观结构**：在P–V坐标中屈服事件坍缩为V≈6.6–6.7 cm³/mol、P_Y∈15–26 GPa的紧凑簇，两取向族间无可视分离。转变事件则描绘单一单调类状态方程分支从(V,P)≈(6.4, 24 GPa)至(5.7, 107 GPa)——等价于事件分支上P_T∈24–107 GPa。该分支为两取向共享但采样不均：[110]密网格填充整个高压尾，而较稀疏的[100]网格集中于较大θ，主要占据同一曲线较低压较大体积部分。该覆盖差异解释了为何P_T的取向效应在轨迹极端最显著而非族间系统偏移。最终微观结构亦路径依赖：最终f_HCP和f_FCC接近f_HCP+f_FCC≈1，符合产物局部密排但含层错、孪晶和界面被CNA归类为FCC-like的预期。散布从近纯HCP（f_HCP≈0.96）至重层错产物（f_HCP≈0.20，f_FCC≈0.78）。[110]案例以较低中位压力转变（41.2 vs 64.3 GPa）但终以较低中位f_HCP（0.607 vs 0.886）：沿[110]对应的早转变在该数据集中并未产生更干净的馬氏体产物。

**按加载类别的微观结构结果**：静水压点是最清晰验证基准。两族共享近各向同性分配，若两细胞设置和提取管线对齐则力学响应应重合：P_Y=26.22 GPa（[100]）和26.21 GPa（[110]）；P_T=37.98和38.11 GPa（差0.13 GPa）；V_T=6.218 cm³ mol^?1一致到三位小数。静水加载下(e)(f)面板不坍缩（f_HCP=0.76 vs 0.61），这是预期的：相同平均压力不保证相同产物织构，一旦路径和取向进入形核 pathway。远离静水锚点，加载对称性以与补充材料一致的顺序组织力学标签。P_Y对类别不敏感：每类基本占据同一15.3–26.2 GPa带，静水案例在上端（26.2 GPa）单轴在下端（15.3–15.8 GPa）。P_T和V_T则强分离类别：静水转变近38 GPa和V_T≈6.22 cm³ mol^?1，而[100]单轴和双轴参考达82和79 GPa at V_T≈5.86 cm³ mol^?1；[110]单轴和双轴条更宽（P_T从26至86 GPa）因该族采样四个不同单轴方向和三个双轴对而非各一个代表顶点。一般三轴加载填充这些极限间区间（P_T 23.8–106.9 GPa）并承载大部分取向对比：[100]三轴案例平均P_T=62.5 GPa vs [110]的49.2 GPa。微观结构遵循类似粗模式但类内仍取向依赖：[100]单轴和双轴参考最高f_HCP最低f_FCC（f_HCP≈0.85），[110]双轴态中间（平均f_HCP≈0.74），静水[110]产物分数最低（f_HCP≈0.61）。三轴条跨0.19≤f_HCP≤0.94（[100]）和0.27≤f_HCP≤0.96（[110]），故加载类别是微观结构的有用第一坐标但不能替代完整(θ,φ)网格。

**物理信息代理模型**：仅用转变时刻定义的量（无加载方向输入）拟合紧凑代理模型。梯度提升回归器（GBR）为scikit-learn GradientBoostingRegressor，500树、最大深度3、学习率0.05、子采样1.0、随机种子0，在合并221例上用5折交叉验证（打乱、固定种子）评估；多项式P_T(V_T)映射用二次阶同折。三个映射为P_T=f(V_T, σ_eq|T)、P_T=g(V_T)（二次EOS）、V_T=h(P_T, σ_eq|T)。P_T从V_T和σ_eq|T映射R²=0.99、RMSE 2.1 GPa（5折交叉验证）：一旦轨迹进入转变窗口，P_T跟踪BCC分支体积V_T，σ_eq|T解释剩余路径依赖。二次P_T(V_T)关系R²=0.96、RMSE 4.3 GPa。V_T从P_T和σ_eq|T预测R²=0.98、RMSE 0.025 cm³ mol^?1：事件处压力、体积和σ_eq一旦知其二则近乎冗余。

**不确定度量化**：在相同特征上添加Matérn-ν=5/2高斯过程回归器：常数幅度核加白噪声，超参数每训练折最大似然优化（四次随机重启），输入每折标准化，目标归一化。交叉验证精度与上述GBR匹配（面板(a) R²=0.996、RMSE 1.5 GPa；面板(b) P_T(V_T) R²=0.96；面板(c) V_T R²=0.99）。归一化残差(y?μ)/σ在采样噪声内服从N(0,1)：经验±2σ覆盖率(a)–(c)为96%、93%、93%（标称95%），标准化残差RMS为0.99、1.01、0.98。这些表明一旦轨迹达转变窗口，标签(P_T, V_T, σ_eq|T)位于紧密曲面：事件级GBR以交叉验证RMSE 2.1 GPa（R²=0.99；GP RMSE 1.5 GPa）恢复P_T，约比P_T的加载方向跨度（23.8–106.9 GPa，因子4.5）小四十倍。

**讨论与结论**

P_T的4.5倍散布不能简单归为冷（10 K）快（10⁹ s^?1）MD的动力学伪影：该解释会预测P_T跨所有加载路径大致均匀偏移——温度、速率、细胞尺寸和势作为各向同性因子。然而P_T随取向、应变分配和事件处偏应力强烈变化（而P_Y在相同路径上变化很小），故散布跟踪加载路径而非仅速率或温度。有限速率压缩下P_T解释为BCC母相失去亚稳性的压力而非平衡相边界。

虽然α→ε onset对加载路径、偏应力和晶体学取向的依赖性此前已由实验和MD确立，本目录不声称该依赖性为新发现；而是提供数百个在单一一致势、温度、速率、细胞尺寸和事件规则下模拟的系统比较。除确认已知依赖性外，网格分析突显两个关键特征：第一，屈服样压力和转变压力不同步缩放——同一网格上P_Y变1.7倍而P_T变4.5倍；第二，两取向的转变事件落在单一压力-体积分支上，表明宽区间代表母相沿该特定分支失去亚稳性的压力而非无关 onset 的集合。

事件级拟合以仅两个转变时刻输入在近五倍压力范围上达R²≈0.99。主导控制是V_T——相变 onset 时BCC分支上的体积——而非加载描述符；同一时刻的σ_eq|T拾取体积单独遗漏的部分。二次P_T(V_T)映射和逆映射V_T(P_T, σ_eq|T)显示同样紧密耦合而无对加载方向的过拟合。GP给出留出折上±2σ覆盖率在标称95%的几个百分点内；一旦知V_T和σ_eq|T，P_T典型后验半宽为几GPa。

对冲击和斜坡实验，onset压力无加载几何难以解释。对MD势，拟合到受控应力状态（通常静水）比匹配路径依赖可观测量的文献平均值更安全。与冲击和斜坡数据及后续Fe–X模拟的比较需完整应力张量而非仅静水压力。

与早期工作的直接比较需分离平衡参考和路径依赖操作标签。静态X射线研究将室温α–ε边界置于约13 GPa附近；本文所用修正Ackland EAM拟合到约13.75 GPa的静水MD值。本文近静水案例给出P_T≈38 GPa，这不是平衡测量：高应变率（和冷温度）将 onset 移至静态参考之上，如其他快MD压缩铁中所示。箔斜坡压缩实验报告α→ε onset 15–38 GPa取决于应变率；这些值位于本文P_T分布（23.8–106.9 GPa）较低部分，与早转变加载路径一致。单轴和强三轴案例延伸至107 GPa，即偏应力阻碍BCC–HCP shuffle时预期的过压。冲击和原位X射线工作确认沿[001]的相同晶格对应；本文静水分支在可比操作压力下匹配该对应，而连续(θ,φ)网格分辨单方向冲击MD或离散三轴研究仅在孤立点采样的织构。近期三轴MD同样报告阈值和HCP/FCC形貌随三轴度偏移；本文散布遵循整个单纯形上相同趋势。

对后续Fe–X工作——行星核Fe–Ni、与超级地球内部相关的Fe–Si和Fe–Ni–Si、以及压力下Fe–O/Fe–C——本目录是参考框架：相同(θ,φ)网格、相同事件定义、纯Fe宿主，合金轨迹可与之比较。溶质可能改变的不止转变压力：它们可改变相稳定性、转变势垒、局部应力场、形核路径以及P_Y和P_T的取向依赖性。计划扩展包括同一(θ,φ)网格上的置换和间隙溶质（Ni、Si、O）、选定单纯形点的缺陷模板以及有限温度和磁性。

本文在连续应变分配网格上以两个取向（221条LAMMPS轨迹）绘制了单晶铁中α-Fe→ε-Fe转变。P_Y保持在窄带（15.3–26.2 GPa）而P_T跨度超过4倍（23.8–106.9 GPa）。P_Y是此前在该势上识别的弹性-塑性 onset，塑性先于转变；在整个网格上其变化远小于P_T。转变事件处的微观结构和偏应力随加载路径变化。基于V_T和σ_eq|T的梯度提升代理模型以R²=0.99（RMSE 2.1 GPa）恢复P_T（事件级V_T+σ_eq|T模型）；相同输入上的高斯过程模型给出校准的不确定度。有限速率下铁中的转变取决于应力状态和取向而非单一压力。文献静态平衡参考约为13 GPa，而本目录中的操作P_T标签跨度为23.8–106.9 GPa。该区间下端与斜坡压缩窗口（15–38 GPa）重叠并延伸至强过压单轴路径。该网格上近静水案例在10⁹ s^?1下用修正Ackland势给出P_T≈38 GPa。