热力学模拟计算

我做过一个合金体系的热力学模拟计算，客户想知道他们的三元合金在不同温度下哪个相更稳定，以及相变时热容会出现什么信号。我这个项目让我把吉布斯自由能 G(T) 和热容 Cp 的关联重新讲了一遍。

项目背景是某镍基合金的 α 固溶体、β 金属间化合物和液相之间的竞争。左侧图就是三条 G-T 曲线：α(solid) 蓝色、β(solid) 绿色、Liquid 红色。低温区 α 相 G 最低，是稳定相；随着温度升高，β 相 G 下降更快，在 T_α-β ≈ 300 K 附近两条固相线相交，β 相成为稳定相；继续升温，液相 G 降到最低，在 T_m ≈ 300 K 附近发生熔化。这里要注意，300 K 是为了教学演示取的典型值，真实合金的相变温度可能是 1200 K 以上，但分析方法完全一样。

吉布斯自由能的定义是 G = H − TS。温度升高时，熵项 −TS 越来越重要，高熵相（通常是液相）在高温区占优。计算上，我用了第一性原理热力学：先在 0 K 用 DFT 算各相的总能 E0，再通过准谐波近似（QHA）把振动熵 S_vib 加进来，得到 G(T)。振动熵对固相影响大，对液相影响小，所以固-液 G 差会随温度缩小。对于液相，我还额外用 CALPHAD 数据库做了外推，因为纯 DFT 处理无序液相比较棘手。

右侧图是热容信号 Cp(T)。可以看到低温段三条曲线的基线不同，反映了声子贡献的差异；在相变温度附近出现一个尖锐的潜热峰。这个峰的面积就是相变潜热 ΔH，峰的宽度与转变类型有关：一级相变峰尖锐，二级相变峰平缓。我算了 α→β 和 β→Liquid 两个峰，后者高度明显更高，说明熔化潜热大于固-固相变潜热。

我这个项目踩过的坑主要是振动熵收敛。QHA 需要不同体积下的声子谱，一开始我用的体积点太少，导致热膨胀系数为负。后来我在平衡体积附近取了 7 个体积点，每个体积做 PHONOPY 声子计算，3×3×3 超胞，才得到合理的 α(T)。另一个坑是磁性贡献：镍基合金有磁性，G 里还要减去磁熵。我用公式 S_mag = R ln(2S+1)，其中 S 是平均原子自旋，这一步对 α/β 的相对稳定性影响约 5 meV/atom，不能忽略。

工程经验上，我建议做热力学模拟时，先确认相图实验点或文献数据，用 DFT+CALPHAD 混合方案最稳。纯 DFT 适合做 0 K 能量排序和趋势预测，但定量相变温度通常偏高或偏低 50-100 K。交付时我给了客户一个温度-相稳定性表，告诉他 300 K 以下用 α，300-1200 K 用 β，1200 K 以上考虑液相，同时标注了 ±50 K 的理论不确定性。

关于相图热力学和 CALPHAD 的结合方法，我整理在 [热力学模拟计算](https://www.keyanxueshu.com/category/xtjs/xcfjs/) 栏目里。[热力学模拟计算](https://www.keyanxueshu.com/)

再看一个细节：右侧 Cp 峰的面积等于相变潜热 ΔH，可以用数值积分得到。我们这张图里 α→β 的峰比较平缓，β→Liquid 的峰又高又尖，说明固-固相变是弱一级或二级相变，熔化是强一级相变。做热力学模拟时，Cp 峰的位置往往比 G-T 交点略偏，因为相变区有滞后。我的经验是，把 G-T 交点作为平衡相变温度 T_eq，把 Cp 峰温作为实验可观测的 T_peak，两者差通常在 5-20 K。交付给客户时，我会同时给出这两个温度，并说明哪一个更适合和 DSC 曲线对比。如果峰很宽，可能是模拟盒子太小或相变过程中出现亚稳相，需要加大超胞重跑。