热障涂层系统(Thermal barrier coating system，TBCs)一般由起隔热作用的表面陶瓷层(Top coat，TC)、粘结层(Bond coat，BC)和服役过程中因氧化产生的热生长氧化层(Thermally grown oxide，TGO)组成。TBCs的热导率很低，可以为其涂覆的热端部件提供热防护，使热端部件能满足更高的工作温度要求。此外还可以保护热端部件在高温下工作时免受腐蚀和氧化，因此目前已被广泛应用在燃气轮机和航空发动机上^[1-2]。由公开报道可知，6mol%~8mol% YSZ具有高熔点、低热导率、与基体相匹配的热膨胀系数、高透氧性及较好的高温热稳定性等特点，是目前应用最广泛的TBCs表面材料^[3]。TBCs服役状态下的热导率大小是判定其性能好坏的重要指标，而YSZ中离子扩散会对其热导率造成一定的影响，同时也会影响TGO的生长速率，进而影响TBCs的服役寿命^[4-5]。

由于在高温工作环境下，传统实验方法很难获取TBCs的实时服役状态，而且目前的实验手段对涉及TBCs微观尺度层面的研究十分困难，很难通过实验方法对服役状态TBCs中离子扩散情况和微观热输运机制进行研究。分子动力学(Molecular dynamics，MD)方法可以模拟高温体系中原子的运动，解析体系在高温下宏观现象表现出的微观机制，因此采用MD方法可以克服上述问题^[6-8]。MD是指对于原子核和电子所构成的多体系统，用计算机模拟原子核的运动过程，从而对系统的结构和性质进行运算，其中每一原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动^[9]。为了深入了解YSZ在服役温度下热导率变化及其离子扩散的内在机制，指导设计开发出性能更加优异的TBCs，采用MD方法进行研究是一种行之有效的方法。

对于分子动力学模拟来说，模拟体系的大小对热导率的结果有一定的影响。如果模拟体系偏小，存在的尺寸效应将无法保证体系内的全部声子-声子间散射，但过大的体系会大大增加模拟的运行时间，为此本文分析了模拟尺寸对热导率的影响规律。

虽然已有文献采用MD方法对Y₂O₃稳定ZrO₂ (Yttria-stabilized zirconia，YSZ)热导率模拟计算进行了相关研究^[10-14]，但在服役温度1273~1473 K条件下，微观模型尺寸对YSZ离子扩散的影响、温度对YSZ热输运特性及离子扩散的影响尚未明确，特别是以温度对YSZ热输运特性及离子扩散影响为根据来阐明1273~1473 K下YSZ热输运机制的研究还未见报道。

本文通过非平衡分子动力学模拟的方法对立方8YSZ(c-8YSZ)材料在1273~1473 K服役温度下的热输运特性、离子扩散情况进行研究，从微观层次探讨尺寸效应、温度对YSZ热导率及离子扩散的影响，揭示1273-1473 K下YSZ热输运机制。

1.   模拟步骤及过程

采用LAMMPS MD模拟的方法，本文选用计算准确性上更具优势的Langvin局部热浴法进行YSZ材料热导率的计算^[15]。图1为YSZ模型系统热浴示意图，模型沿热流方向的长度为L_x。

图  1  8mol%Y₂O₃稳定ZrO₂(8 YSZ)模型系统热浴示意图

Figure  1.  Model system heat bath diagram on 8mol%Y₂O₃-stabilized ZrO₂ (8 YSZ)

L_x—Model length; z—Direction of heat flux

下载: 全尺寸图片 幻灯片

根据傅里叶定律计算热导率，如下式所示：

其中：J是沿z方向的热流；A是模型垂直热输运方向横截面积；ΔT是拟合得到的温度梯度。

模型的构建过程如下，首先根据晶格常数构建立方ZrO₂单胞，然后在x、y、z 坐标方向进行扩胞，从而建立立方ZrO₂超胞模型，之后根据8YSZ中各成分的掺杂百分比计算出Y³⁺与Zr⁴⁺在超胞中的比例，按比例用Y³⁺离子随机取代Zr⁴⁺离子进行掺杂，且在掺杂过程中每掺杂两个Y³⁺离子产生一个氧空位，来保持模型系统的电中性，最后完成c-8YSZ模型构建。温度的变化会影响单胞的晶格常数a(nm)，它们的关系如下式所示^[16]：

模拟过程中x、y、z方向的边界条件均设置为周期性边界条件，采用Velocity-Verlet 算法求解运动方程，各离子的初始速度按Maxwell-Boltzmann分布给定，用Langevin热浴方法对模型系统进行控温。选取Buckingham势函数和库伦势来描述YSZ离子间的相互作用，其中Buckingham势描述短程相互作用，库伦势和峰值脉冲功率(PPPM)描述长程相互作用^[17]：

其中：r_ij代表i、j两个离子间的间距；q_i和q_j表示有效离子电荷；A_ij、ρ_ij和C_ij是Buckingham势的参数。模拟中采用的势函数参数如表1所示^[18-19]。对Zr、Y和O粒子，分别赋予有效离子电荷+4.0、+3.0和−2.0。

表  1  分子动力学中所采用的势函数参数

Table  1.  Potential functions used in molecular dynamics

Interatomic interactions/i-j	Aij/eV	ρij/nm	Cij/(eV·nm6)
Zr4+-O2−	985.869	0.0376	0
Y3+-O2−	1325.60	0.0349	0
O2−-O2−	22764.3	0.0149	2.789×10−5
Notes: Aij, ρij and Cij—Buckingham potential parameters.

下载: 导出CSV 

| 显示表格

模拟中采用的时间步长为1 fs，首先采用共轭梯度法对模型系统进行静力学平衡使初始结构达到能量最小化；然后依次在等温等压系综(NPT)和正则系综(NVT)下分别驰豫200 ps，进行动力学平衡，使系统在给定温度下达到稳定结构；接着转到微正则系综(NVE)在Langvin热浴下平衡300 ps，对系统进行松弛，最后通过900000步的计算收集数据。

2.   结果与讨论

2.1   尺寸效应对YSZ材料热导率的影响

对于块体材料的热导率计算，模拟体系的大小对热导率的结果有一定的影响，如果模拟体系偏小，无法保证体系内的全部声子-声子间散射，但过大的模拟体系会大大增加模拟的运行时间，为此分析了模拟尺寸效应对热导率的影响规律。以c-8YSZ作为研究材料，在1473 K温度下，截面积A为4a×4a、5a×5a、6a×6a，模型长度L_x分别取5a、25a、45a、60a进行热导率分析计算，结果如图2(a)所示。可知模型截面大小对YSZ热导率的计算几乎没有影响，而热输运方向结构长度L_x对YSZ热导率模拟计算影响较大。为进一步研究模型长度与热导率的关系，后续模型截面选取5a×5a，L_x分别取5a、15a、25a、35a、45a、50a、55a、60a、65a进行模拟计算。从图2(a)看出，由于5a×5a在L_x=45a以后的热导率变化趋势开始平缓，几乎没有什么变化，因此选取5.70 W/(m·K)作为热导率标准值，定义Δ为热导率计算值与标准值的相对误差，计算结果如图2(b)所示。可知，热输运方向模型长度对热导率计算结果影响较大，热输运方向热导率随模型长度L_x增大而增大；且当L_x≥50a时，Δ趋近于0%，所得热导率基本趋于稳定(变化率＜5%)，没有表现出明显的尺寸效应，这说明YSZ中依旧存在自由程大于25 nm的声子影响着热导率，但对热导率的影响很小，可忽略不计。故在本文后面的模拟中，近似使用L_x= 60a的模拟体系计算得到的热导率结果作为材料热导率的计算结果。

图  2  8YSZ热导率随模型结构长度的变化关系(a为晶格常数)

Figure  2.  Effect of model structure length on calculated thermal conductivity of 8YSZ (a is the lattice constant)

Δ—Relative error between calculated and standard values of thermal conductivity

下载: 全尺寸图片 幻灯片

2.2   尺寸效应对YSZ材料离子扩散的影响

由于YSZ中阳离子扩散相较于O²⁻的扩散较慢，且差距较大，研究意义不大，因此本文主要对O²⁻扩散进行研究^[20]。当截面尺寸为5a×5a，L_x分别为5a、10a、15a、25a、35a、45a、50a、55a、60a、65a时，在1473 K温度下对O²⁻的扩散情况进行模拟，8YSZ体系中不同模型尺寸的O²⁻均方位移(Mean square displacement，MSD)如图3所示。可知，O²⁻的MSD和时间呈线性关系，用O²⁻在模拟体系中的迁移量来判断O²⁻在体系中的振动强度。随着热输运方向模型尺寸L_x的增大，8YSZ系统在1473 K时的斜率逐渐增大，扩散越快，振动越强烈。当热输运方向模型尺寸小于25a时，O²⁻扩散受模型尺寸的影响，模型尺寸越大，扩散越快；当热输运方向模型尺寸大于25a时，离子扩散与模型尺寸无关。

图  3  8YSZ在1473 K下不同模型中O²⁻的均方位移(MSD)

Figure  3.  Mean square displacement (MSD) of O²⁻ ions in different models at 1473 K for 8YSZ

下载: 全尺寸图片 幻灯片

图4分别为L_x取5a、10a、15a、25a、60a时，1473 K温度下x、y、z方向及总的O²⁻的MSD模拟计算结果。可知，当热输运方向模型尺寸小于15a时，z方向(即热输运方向)的氧离子迁移速率小于x、y方向，且随着模型尺寸的增大3个方向的迁移速率趋于相同；当热输运方向模型尺寸大于等于15a时，x、y、z方向的氧离子迁移速率基本相同，当L_x为60a时，3个方向的氧离子迁移速率完全重合，因此认为热输运方向结构长度L_x的影响已经完全消除。

图  4  8YSZ在1473 K下不同温度中O²⁻的MSD

Figure  4.  MSD of O²⁻ ions in different temperatures at 1473 K for 8YSZ

下载: 全尺寸图片 幻灯片

2.3   1273 ~1473 K温度下YSZ热输运机制

取5a×5a×60a模型为研究样本，图5为在1273~1473 K范围内温度对c-8YSZ中O²⁻的MSD的影响。可知O²⁻扩散受温度影响较大，且随温度升高，MSD斜率随之增大，O²⁻扩散变快。

图  5  8YSZ在不同温度下O²⁻的MSD

Figure  5.  MSD of 8YSZ O²⁻ ions at different temperatures

下载: 全尺寸图片 幻灯片

图  6  高温下立方8YSZ(c-8YSZ)模拟的热导率、实验和文献中四方YSZ (t-YSZ)模拟的热导率

Figure  6.  Thermal conductivity of cubic 8YSZ (c-8YSZ) simulations at high temperatures, experimental and tetragonal YSZ (t-YSZ) simulations in the literature

下载: 全尺寸图片 幻灯片

图6为在1273~1473 K范围内温度对c-8YSZ的热导率的影响模拟计算结果及其他研究人员的实验和计算结果^[21-22]。可以看出c-8YSZ的实验导热系数在1273 K开始，随着温度的增加热导率数值几乎没有发生变化。而文献[22]报道的四方YSZ(t-YSZ)的模拟热导率随温度升高而升高。通常声子数与温度成正比，声子数与温度的关系可以用下式表示^[23]：

其中，$n$、${k}_{\mathrm{B}}$和$T$分别表示声子数、玻尔兹曼常数和温度。而声子在碰撞前的移动距离(平均自由程)决定了热导率，如下式所示。声子的平均自由程可以表示为^[24]

其中，${l}_{i}(\omega ,T)$、${l}_{p}\left(\omega \right)$ 和 ${l}_{b}$分别代表声子-声子散射、点缺陷散射和晶界散射。根据式(5)可知，点缺陷、晶界和位错等缺陷的存在导致声子的数量的增加，进而提高了声子碰撞的概率，碰撞进一步增强了声子的散射，改变了声子的初始传热方向并释放其能量。由于碰撞引起能量损失，声子的传输发生减弱或终止(声子的平均自由程降低或保持不变)。

根据式(4)~(5)可以推测出t-YSZ热导率随着温度的升高，决定热导率的声子数增加，声子的平均自由程随温度的升高而增加所导致的结果。

而本研究所模拟c-8YSZ与实验测得的热导率的变化趋势一致，即在高温服役状态下，c-8YSZ热导率基本趋于稳定，原因是在高温环境下，声子的平均自由程已经达到极值，使影响热导率的因素保持稳定。而t-8YSZ与c-8YSZ热导率在高温下存在不同的趋势可能是由于其结构不同导致，如较小的晶格参数可能会导致较窄的原子间距而阻碍声子的传输。

在该研究中每个测试温度下，理论和实验结果的差异是由于在MD模拟中使用的理想模型是没有晶界散射的。由于本文的MD模拟忽略了所有晶界和部分缺陷，只存在氧空位，因此c-8YSZ的理论导热系数大于实验导热系数。

上述结果说明，高温服役状态下O²⁻扩散对热输运的影响不大，YSZ热输运主要依靠声子进行传递，因此增加声子散射是有效抑制YSZ服役状态下的热输运能力的方法。可通过氧化物掺杂如CeO₂、Sc₂O₃、Al₂O₃和提高孔隙率等方法来增强高温下YSZ基热障涂层的声子散射，进而降低涂层材料的热导率^[25-27]。

3.   结 论

(1) 横截面积的大小对热导率几乎无影响，而热输运方向长度尺寸对热导率影响较大；当热输运方向模型尺寸小于50晶格常数(a)时，计算得到的热导率随热输运方向模型尺寸增大而增大，而当热输运方向模型尺寸大于50a时，计算得到的热导率趋于稳定。

(2) 当热输运方向模型尺寸小于25a时，O²⁻扩散受模型尺寸的影响，且随模型尺寸的增加，O²⁻扩散越强；当热输运方向模型尺寸大于25a时，O²⁻扩散与模型尺寸无关。

(3) 当热输运方向模型尺寸小于15a时，z方向的O²⁻扩散速率小于x、y方向，且随着模型尺寸的增大3个方向的扩散速率趋于相同；当热输运方向模型尺寸大于等于15a时，x、y、z方向的O²⁻扩散速率基本相同。

(4) 在1273~1473 K的温度范围内，立方8YSZ(c-8YSZ)热导率基本趋于稳定，而氧离子扩散受温度影响较大，且随温度升高离子扩散变快，因此在高温服役状态下，O²⁻扩散对YSZ热输运的影响不大，YSZ热输运主要靠声子进行，增加声子散射是有效抑制YSZ服役状态下的热输运能力的方法。