β-SiC纳米薄膜热导率的数值模拟研究

白素媛,贾孟晗，吴宝丽，丛士博，白容榕

(辽宁师范大学 物理与电子技术学院，辽宁 大连 116029)

碳化硅是第三代半导体材料之一，是一种宽禁带半导体材料，具有禁带宽度大、击穿电压高，电子迁移率高、热导率高、抗辐射能力强等优良的物理学性质，可广泛应用于人造卫星、导弹、雷达、通讯、航空航天等军事系统[1].β-SiC是立方碳化硅,具有高化学稳定性、高硬度、高热导率、低热胀系数、宽能带隙等良好的物理化学性质，在制备高级耐火材料、结构陶瓷材料、高级β-SiC制品、军工材料及镀层、涂层材料等方面具有十分广阔的市场需求[2].多晶材料建模与计算十分复杂，目前对碳化硅热导率的计算主要集中在单晶领域，对β-SiC热导率的计算方法主要有两种，分别是基于声子Boltzmann传输方程计算，以及利用分子动力学模拟计算.

声子Boltzmann传输方程计算方法是大多数宏观固体传热理论的基础，其将晶格看作准连续介质，在计算过程中引入较多的经验参数，从而限制对热导率的计算[3].分子动力学方法则比声子Boltzmann传输方程更加简单，在计算微尺度效应时非常有前景，非平衡分子动力学模拟可以直接依据原子能量的变化求得体系的温度分布关系和热流密度，利用Fourier定律计算导热系数.白素媛[4-6]利用非平衡分子动力学和实验方法研究了亚微米薄膜导热性能.Wang[7]利用非平衡分子动力学方法研究了纳米级6H-SiC法向和切向的热导率.分子动力学方法的关键是选取原子之间的势函数，势函数表征电子云重叠的量子力学作用.对多数原子来说，精确的势函数并不清楚，因此经验和半经验的势函数在分子动力学模拟中被广泛使用[8]，集成的势函数组成了力场，如DERIDING力场[9]、UNIVERSAL力场[10]、COMPASS力场[11]等被广泛使用.

本文使用COMPASS力场，采用非平衡分子动力学方法，研究计算了β-SiC纳米薄膜在室温下的热导率，此外，还研究了其热导率与薄膜厚度、热导率与系统温度的变化关系，研究内容对β-SiC工程材料应用具有一定参考价值.

1 非平衡分子动力学模拟方法(NEMD)

非平衡分子动力学方法通过对系统加入温度梯度或热流引起扰动，利用非平衡输运模型、傅立叶定律计算系统传热系数.在边界条件的选取上分为非周期性边界条件(固壁边界条件或自由边界条件)和周期性边界条件.根据学者Lukes[12]的研究成果，边界条件不会影响模拟结果，因此在本文计算中采用周期性的边界条件.周期边界条件的NEMD计算模型如图1所示.对 X、Y、Z 3个方向系统均采用了周期性边界条件，在厚度方向上(Z 方向)分成若干层.系统的整体分为热传输层、高温控制层和低温控制层.在厚度方向 1/4 和 3/4 处取相等层数的粒子分别作高温控制层和低温控制层.周期性边界条件中边界的最上端和最下端相连通，所以两个控温部分的距离等于厚度的一半.

图1 周期性边界条件计算模型

β-SiC纳米薄膜模型的建立：在系统中建立β-SiC基本单元结构模型如图2所示，整体计算模型在X、Y方向上分别取4个基本单元进行计算，在Z方向即厚度方向上取16个连续的基本单元建立超晶胞，4*4*16 β-SiC计算模型如图3所示.

图2 β-SiC基本单元结构模型 图3 β-SiC 4*4*16计算模型

2 COMPASS势函数

在进行分子动力学模拟时，选取粒子之间合适的势函数是关键，势函数对模拟准确性有很大的影响.文中对β-SiC薄膜模拟时采用的COMPASS力场，其在凝聚态物理中被广泛应用，能够预报孤立态和凝聚态分子结构、热力学性质的分子力场.通过计算分子之间键参数，利用经验方法测度范德华非键合的参数.COMPASS力场级数展开形式如下：

E=Eb+Eθ+Eφ+Ex+Ecross+Eij+Eelec=

(1)

公式(1)中，k、H、V、A、B为力常数，b为键长，θ为键角，φ为二面角，ε为能量参数，rij为两个原子i，j间距，q为电荷量.键合项、非键合项两部分组成力场能量.对角和非对角的交叉耦合项分为伸缩能Eb，键角弯曲能Eθ，键扭转能Eφ，键角面外弯曲能Ex和相互耦合能Ecross.非键合项包括范德华能Eij和库仑能Eelec，这两种能量计算原子对之间的相互作用力.

在最初始状态时，分子分布要与实际相似，例如分子的初始速度要接近真实分子速度情况.判断系统初始平衡关键要确定系统分子的速度分布是否满足平衡态速度分布规律.例如三维系统，速度要满足麦克斯韦方程：

(2)

(3)

公式(3)为公式(2)的归一化形式，其代表温度T条件下，分子速度在v周围数量占总数的百分比.采用麦克斯韦方程赋予粒子初始速度方法有其优势：可以让系统短时间内达到平衡.

3 β-SiC薄膜热导率计算

对于存在温度梯度的系统，需要判定系统的稳态是否形成，判定的方法为系统中厚度方向的温度变化曲线是否光滑.以计算模型4*4*20为例，其厚度为4.296 nm，在系统温度为300 K时，经过计算，其在50,100,150,200,250万步条件下的沿厚度方向的不同层温度为图4中所示.

图4 不同计算步数下模型温度梯度

随着计算步骤的增多，曲线逐渐光滑，在150万步以后几乎不再发生变化.与其对应的不同计算步数下的热导率见表1.

表1 不同计算步数下的β-SiC热导率

表2 不同厚度下的β-SiC热导率

由表1可得热导率随计算步数的关系如图5所示，由图5中可以看出，计算模型的热导率逐渐趋于平稳，且200万计算步数下的热导率计算结果与250万计算步数下的结果几乎相同，考虑到计算精度与计算效率的关系，最终选择计算步数200万步分别进行不同厚度薄膜热导率的计算研究.厚度方向分别选取16,20,24,28,32层基本单元，经过计算不同厚度下的β-SiC热导率见表2，绘制不同薄膜厚度下热导率曲线如图6所示，在环境温度为300 K时，β-SiC薄膜厚度在3.426～6.905 nm范围内变化时，其热导率在5.054～9.654 Wm-1K-1之间，可以看出，随着薄膜厚度的增加，其热导率几乎呈线性增加，具有明显的尺度效应.

图5 热导率与计算步数的关系 图6 热导率与薄膜厚度的关系

在β-SiC计算模型为4*4*24时，改变系统温度，研究此计算模型热导率随温度变化的关系见表3，绘制不同温度下的模型热导率曲线见图7.在温度变化范围300～1 000 K时，β-SiC薄膜热导率先增加后降低，变化逐渐变缓，温度升高时β-SiC热导率变化较小.

表3 不同温度下的β-SiC热导率

4 结论与展望

图7 不同温度下的β-SiC热导率

本文在非平衡分子动力学方法基础上,模拟了β-SiC薄膜在纳米尺度下的热导率，得出如下4点结论：

(1)基于系统中的β-SiC模型，利用Compass力场势函数计算了β-SiC薄膜热导率.

(2)在计算步数为200万时，系统内温度曲线光滑，β-SiC薄膜热导率计算更加准确.

(3)随着β-SiC薄膜厚度的增加，其热导率变大，表现出明显的尺度效应.

(4)在系统温度发生变化时，β-SiC薄膜热导率发生一定变化，随着温度升高，热导率趋于稳定.

随着科技的飞速发展，电子产品的结构与器件尺寸范围早就不局限于肉眼宏观可见，而是可从微米至亚微米甚至纳米级别.碳化硅具有抗辐射能力强、热导率高、热稳定性好等优良性质，在高温传感器件、航空航天及光电集成器件等方面应用前景广阔.在纳米碳化硅薄膜中掺杂氮元素，使其具有特殊的光电性能，还有待研究.由于目前制备碳化硅纳米级薄膜技术还不够成熟，本文的研究成果对接下来材料的研发与生产也将起到一定的帮助作用.