热电材料利用塞贝克效应和帕尔贴效应以实现热能和电能之间的直接相互转化，在工业余热发电、热电制冷、太阳能综合利用等领域具有广泛而重要的应用前景^[1-2]。研究表明^[3-4]，材料热电转换效率的提高依赖于热电优值^[5]相关物理特性的综合提高，需要在提高电输运能力的同时降低材料的导热系数。一类新型的热电材料—“声子液体”材料^[6]内含固定的亚晶格框架，可提供良好的电输运通道，同时在晶体内部又存在其他类“液态”离子可强烈散射声子，降低材料的晶格热导率，实现了热传导和电输运的协同优化。因此，研究“声子液体”热电材料的类“液态”结构特殊性，并分析该特殊性对材料热传导性能的影响，可以为热电转换效率的提高和热传导过程机理的理解提供借鉴与指导。

“声子液体”热电材料的结构特殊性在于其中包含的类“液态”离子。Liu等^[7]发现Cu_2-xSe晶体随着温度的升高会经历结构的相变过程，即由低温的非立方相(α相)转变为高温的立方相(β相)。在高温立方相(下文称β-Cu_2-xSe)中Cu⁺不再有序排列，成为可迁移扩散的快离子，构成了类“液态”的亚晶格结构，而其具有的强烈散射声子能力使得Cu_2-xSe晶体在高温下具有极低的热导率。Ouyang等^[8]和Liang^[9]分别研究了Ag₂Te和Cu₂S两种“声子液体”热电材料，发现这些材料中存在的类“液态”离子同样可以散射声子，降低材料的热导率，进一步验证了“声子液体”类热电材料的结构特征和优异性能。此外，Kim等^[10]应用第一性原理的方法对β-Cu_2-xSe的导热过程进行机理性研究，计算材料在热传导过程中的声子传输性质，尝试解释其较低热导率产生的原因。然而，目前研究更多关注“声子液体”材料的结构性质，未能深入探讨类“液态”离子对材料导热性能的影响。此外，掺杂和空位是优化材料热电性能的常用手段，然而目前对此类材料处理手段的研究鲜见报道。

本工作采用分子动力学模拟的方法研究压强为0.1 MPa，温度为500~1000 K下β-Cu_2-xSe的热传导性能和扩散特性，通过计算材料的扩散系数D、导热系数k和声子输运特征，分析“声子液体”材料中类“液态”离子对导热系数的相关性，同时探讨掺杂和空位对材料导热系数所造成的影响，并解释其影响机制。

1 模型建构及研究方法 1.1 模拟方法及模型构建

本工作的研究对象为β-Cu_2-xSe，一种典型的“声子液体”热电材料。β-Cu_2-xSe具有反萤火石结构，Se^2-形成相对稳定面心立方亚晶格网络结构，Cu⁺则分布在硒亚晶格框架的间隙空间，同时具有随机迁移特性，其结构示于图 1。为了消除尺寸因素对导热系数计算的影响，构建2.3 nm ×2.3 nm ×4.6 nm的计算域，包含128个β-Cu₂Se晶胞，共计1536个原子。此外，通过随机增添和删减Cu⁺的手段反映掺杂(β-Cu_2.05Se)和空位(β-Cu_1.95Se)对材料在微观层面上所做出的改变。

β-Cu_2-xSe中的铜和硒均为金属原子，故粒子间相互作用势能函数采用莫尔斯势^[11]进行描述，莫尔斯势广泛应用于计算金属化合物的导热系数，具有良好的计算精度。莫尔斯势可化简为：

(1)

式中：V_ij为粒子间势能函数；D_e为势能阱深度；a为控制势阱的宽度参数；r_ij为粒子间的距离；r_e表示平衡时的粒子对键长；r_{cut_off}为截断半径。

模拟采用Namsani等^[12]拟合的莫尔斯势参数，示于表 1。

考虑到Cu_2-xSe以430 K为界，存在低温的α相和高温的β相两种不同的相态结构^[13-14]，因此选择β-Cu_2-xSe存在的温度区间500~1000 K进行研究，设定时间步长为1 ps。模拟中首先对系统进行能量最小化以消除系统内不合理的应力，然后在NPT(恒压、恒温)系综下运行50 ps，使系统达到平衡状态，并在之后的60 ns内以5 ns为周期输出材料的导热系数和扩散系数。本工作使用的控温方法为Nose-Hoover法^[15]。模拟中均采用周期性边界条件，使用分子模拟软件LAMMPS^[16]进行计算。

1.2 模型验证

为了验证莫尔斯势参数和所建立模型的正确性，计算β-Cu₂Se在压强为0.1 MPa，温度区间为500~1000 K的密度，结果示于图 2。由图 2可知，β-Cu₂Se在500 K时的密度为6.67 g/cm³，与文献[17]中实验所得高温相Cu₂Se密度(6.703 g/cm³)相一致。将计算所得500~1000 K的密度与文献[18]实验所得结果进行对比验证，相对误差在2%之内，因此模拟结果具有可靠性和准确性。此外，曲线在850 K左右出现转折，密度增加到6.573 g/cm³，这和文献[10]发现β-Cu_2-xSe在800 K左右存在的玻璃化转变相符。

1.3 数据处理方法

扩散系数通过计算系统的均方位移而得到，均方位移由式(2)所定义：

(2)

式中：MSD为粒子的均方位移；r(t)为时间t时粒子的位置；r(0)为开始时刻粒子的位置。

依据爱因斯坦扩散定律^[19]，在长时间跨度下均方位移对时间曲线斜率的1/6即为系统整体的扩散系数：

(3)

式中：D为系统中粒子的扩散系数；t为模拟时长；MSD为粒子的均方位移。

导热系数的计算采用平衡分子动力学(equilibrium molecule dynamics，EMD)模拟方法，也称为Green-Kubo法^[20-21]，该算法的基本思想是通过计算热流量自相关函数(heat flux auto correlation function，HFACF)得到材料的导热系数k，即：

(4)

式中：k_B为玻尔兹曼常数；T为热力学温度；V为系统体积；t为模拟时长；尖括号意指对时间的平均；J(0)和J(t)表示在0时刻和t时刻的热流量，并且可以由式(5)得到：

(5)

式中：J为热流量；V为体积；E_i表示第i个粒子所具有的能量，包含势能和动能；v_i为速度矢量；S_i表示作用在第i个粒子上的应力张量。

2 结果与讨论 2.1 扩散特性及导热系数相关性

“声子液体”热电材料的结构特殊性在于其中存在的类“液态”离子。本工作分析在500，700，800 K和1000 K时系统内离子的运动轨迹范围，结果示于图 3(a)。由图 3(a)可知，随着温度的升高，离子的运动范围扩大。值得注意的是，相比Se^2-，Cu⁺的移动范围更大，运动更剧烈。此外，离子“液态”性可由扩散系数来具体衡量。对于所有离子、Cu⁺和Se^2-分别应用爱因斯坦扩散定律，计算得到β-Cu_2-xSe在500~1000 K间的动态扩散系数，500~1000 K下Se^2-均未出现扩散现象，整个系统和Cu⁺的扩散系数示于图 3(b)，由图 3(b)可以看出，随着温度的升高，系统在700 K时出现扩散现象，并在800 K时出现了扩散系数的下降，这和通过密度分析发现的800 K左右玻璃化转变一致。此外，Cu⁺的扩散系数在趋势与数值上都和整体扩散系数一致。由此可见，β-Cu_2-xSe的扩散系数主要由其中的类“液态”的Cu⁺所贡献，而Se^2-构成了十分稳定的亚晶格结构。

为研究类“液态”离子的扩散和β-Cu_2-xSe导热系数的相关性，计算不同温度(500~1000 K)和不同配位数(1.95，2.00和2.05)下的导热系数和扩散系数，所得结果示于图 4。图 4(a)表明对于β-Cu_2-xSe，温度的升高和配位数的下降都会导致材料扩散系数的增加。在温度为800 K左右时出现的导热系数“逆”增长也符合β-Cu_2-xSe在800 K左右存在的玻璃化转变现象(图 4(b))。排除玻璃态转化的影响，β-Cu_2-xSe的导热系数整体随温度升高而下降。这意味着温度的升高加剧了离子的运动，导致晶格非简谐振动的增加，传热过程中晶格非简谐振动的加剧则使得声子散射更加强烈，在宏观层面上体现为导热系数的下降。

综合比较β-Cu_2-xSe的扩散系数和导热系数，β-Cu_2-xSe扩散系数的增加和导热系数的降低具有同步性。Cu⁺的扩散能力提高引起晶格非简谐振动的增加，从而导致材料导热系数的降低。以上结果表明类“液态”离子的扩散能力和材料导热系数具有很强的相关性。

2.2 掺杂和空位

为了分析掺杂和空位对β-Cu_2-xSe的导热性能的影响，计算3种不同配位数β-Cu_1.95Se，β-Cu₂Se和β-Cu_2.05Se的导热系数，结果示于图 4(b)。由图 4(b)可知，β-Cu₂Se和β-Cu_2.05Se导热系数的变化趋势和范围大致相同，均从1.56 W/(m·K)下降至0.97 W/(m·K)左右，但当温度处于500~800 K，β-Cu_1.95Se的导热系数相比前两者都更低。上述结果表明，当材料形成空位时，Cu⁺获得更多的移动空间，使得晶格存在更多的非简谐振动，从而降低了导热系数；与此相反，掺杂对导热系数的影响不明显，当晶体内未留有明显的移动空间给Cu⁺时，增加Cu⁺的数量不能明显增加晶格整体的非简谐振动，对导热系数的影响并不显著。由此可见，晶体结构的完备性对导热系数造成的影响主要体现在晶体内部是否提供足够空间给类“液态”离子进行移动迁移。

2.2.1 不同离子作用

为了进一步分析掺杂和空位对热传导过程的影响，将导热系数分解为不同离子作用贡献项以研究掺杂和空位处理后不同离子间相互作用对导热系数的影响。依据公式(6)将导热系数分为铜-铜相互作用项(Cu-Cu项)、铜-硒相互作用项(Cu-Se项)和硒-硒相互作用项(Se-Se项)^[22]，所得结果示于图 5。

(6)

式中：k为导热系数，不同下标表示对应离子间相互作用对导热系数的贡献项。

从图 5可以看出，在β-Cu_2-xSe导热系数的构成中，Cu-Se项构成了最大组成部分，贡献了约占72%的能量传递，成为最主要的热量传递方式。掺杂和空位处理后导热系数的区别也主要体现在Cu-Se项的异同。相比β-Cu₂Se和β-Cu_2.05Se，存在空位的β-Cu_1.95Se中Cu-Se项更低并且随温度的变化不明显，说明当材料内部存在足够移动空间时，类“液态”离子倾向于在内部空位间的移动，与固定框架的碰撞概率降低。此外，温度的提高对类“液态”离子移动的选位影响不大。由此可知，材料内部存在的空位使得类“液态”离子的移动集中于空间缺陷点，降低了和晶格框架的碰撞概率，有效地降低了导热系数。相比空位，掺杂的影响则不明显。

2.2.2 声子输运特征

晶格振动理论^[23]认为，晶体中的热传导主要依靠声子输运来完成。晶体内部存在温度梯度时，材料内部声子的移动把能量从高温端传递到低温端。为进一步分析空位对导热系数的影响机制，本工作应用Kong等^[24-26]提出的方法分别计算700 K下配位数2.000，1.969，1.938和1.906的声子色散曲线，结果示于图 6。由图 6可知，随着材料内部空位的增加，热传导相关的声学支声子频率下降，声学支声子的散射增强。材料内部的空位影响了声学支声子的输运，影响了晶体内部能量的传递，从而降低了材料的导热系数。

3 结论

(1)“声子液体”热电材料中类“液态”离子扩散能力的增大和材料导热系数的降低具有极强的相关性。类“液态”离子移动能力的提高会导致晶格的非简谐振动增加，从而降低材料的导热系数。

(2) 材料内部形成的空位有利于离子的移动，从而降低导热系数；材料掺杂对导热系数的影响不明显。空位的存在使得Cu⁺在移动过程中倾向于在内部空间的移动，其具体影响在于空位使得声子声学支频率降低，影响了晶体内部能量的传递，降低材料的导热系数。