热膨胀系数建模计算：从声子谱到温度依赖曲线的物理链路

氧化钇稳定氧化锆（YSZ）是航空发动机热障涂层（TBC）的主力材料，服役温度可以飙到1200℃以上。热障涂层失效的一个关键原因是涂层与基底（镍基高温合金）之间的热膨胀失配——基底的热膨胀系数约15×10⁻⁶ K⁻¹，YSZ在室温到1000℃范围内的热膨胀系数约10-11×10⁻⁶ K⁻¹。4-5个点的差距，在1000℃的温升下会产生约0.4-0.5%的应变失配——积累下来就是涂层剥落。

这个项目的目标是用第一性原理方法计算YSZ在0-1500K范围内的热膨胀系数曲线，并与推杆式热膨胀仪（DIL）的实测值做对标。

Grüneisen理论框架

热膨胀的驱动源是声子频率对体积的依赖性。Grüneisen参数定量描述这种依赖：

γi(q)=−∂ln⁡ωi(q)∂ln⁡Vγi​(q)=−∂lnV∂lnωi​(q)​

γ_i越大，第i个振动模式对热膨胀的贡献越大。线性热膨胀系数α_L可以通过所有振动模式的Grüneisen参数加权平均得到：

αL(T)=13B∑i,qCv(i,q,T)⋅γi(q)αL​(T)=3B1​i,q∑​Cv​(i,q,T)⋅γi​(q)

B是体积模量，C_v是每个振动模式的热容贡献。这个公式把热膨胀的问题转化成了两个可算的量：Grüneisen参数（通过声子频率的体积导数得到）和振动热容（从声子态密度的Bose-Einstein分布得到）。

YSZ的声子谱计算

YSZ不是纯ZrO₂——8mol%的Y₂O₃掺杂引入了氧空位，破坏了完美晶体的平移对称性。建模时采用了2×2×2的ZrO₂超胞（96个原子），将4个Zr替换为Y、同时去掉2个O以保持电荷中性，对应于约8%的Y₂O₃掺杂浓度。

如此大的超胞（含94个原子）做全声子谱计算在计算成本上是一个挑战。实际做法是：用VASP做结构优化（ISIF=3，残余力<0.002 eV/Å），然后在Phonopy框架下用有限位移法计算Γ点附近的声子频率。为捕捉体积依赖性，对7个体积（平衡体积的98%-104%，步长1%）重复声子计算。

YSZ的声子态密度表现出两个显著特征：低频区（<200 cm⁻¹）主要来自Zr/Y的重原子振动，高频区（400-700 cm⁻¹）来自O的振动。掺杂引入的氧空位在低频区产生了一些局域振动模，这些模的Grüneisen参数为负——意味着它们在体积膨胀时频率反而会升高，对热膨胀产生收缩贡献。

从QHA到α(T)

做完7个体积的PhDOS后，对每个温度T，拟合F(V,T) = E_0K(V) + F_vib(V,T)的最小值，得到平衡体积V_eq(T)。线性热膨胀系数从V_eq(T)的温度导数得到：

αL(T)=13Veq(T)dVeqdTαL​(T)=3Veq​(T)1​dTdVeq​​

YSZ的计算结果显示：室温下α_L约8.5×10⁻⁶ K⁻¹，随温度单调上升，在800℃达到约10.6×10⁻⁶ K⁻¹，1200℃时约11.3×10⁻⁶ K⁻¹。而DIL实验值在这三个温度点分别是9.1、10.9和11.5×10⁻⁶ K⁻¹。

全温区计算值与实验值的偏差在±0.6×10⁻⁶ K⁻¹以内（约5-6%），对于需要做涂层应力有限元分析的工程应用来说足够了。

氧空位浓度敏感性

Y₂O₃掺杂量直接影响氧空位浓度，而氧空位又通过负Grüneisen参数抑制热膨胀。做了一组对比计算：4mol%、8mol%、12mol% Y₂O₃对应的1000℃热膨胀系数分别是11.8、10.9和10.1×10⁻⁶ K⁻¹。

趋势清晰：Y含量越高，热膨胀系数越低。从热膨胀匹配的角度看，应该选低Y含量的配方（让α更接近基底），但从相稳定性和抗烧结性能看又需要足够的Y。8mol%是一个工程折中点——这在热障涂层的配方设计历史中已经被大量实验证实了，而DFT+QHA的计算框架能够从声子层面给出物理原因。

适用边界

从Grüneisen参数到α(T)的这条路，本质上是准简谐的——假定声子频率只通过体积变化受温度影响，忽略声子-声子相互作用的非简谐效应。对YSZ这种离子-共价混合晶体，QHA在1300K以下相当可靠。超过1300K（YSZ的服役温度上限约1200℃），非简谐贡献开始出现，如果要追求更高精度，需要用AIMD做直接的热膨胀模拟。

对于热障涂层的工程设计来说，这个计算框架的真正价值在于：可以将α(T)曲线作为成分的函数来扫描——换一种掺杂元素（比如Yb替代Y），不用等实验合成和热膨胀仪测量，几天计算就能给出趋势预判。