基于非平衡分子动力学的氧化锌薄膜热导率的模拟

白素媛，王 斌，贾孟晗，丛士博，吴宝丽

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

氧化锌薄膜相对于其它半导体材料制备所需的温度更低，制备过程更简单且制备成本低，凭借其稳定的热力学性质、良好的化学性能和光电磁等特性成为当今研究的热点[1].纳米薄膜材料是微电子器件的重要组成部分，薄膜导热性能在很大程度上决定了微电子器件和微电子系统的性能和安全性[2]，所以对薄膜材料热导率的研究具有重要的实际意义.国内外学者对氧化锌性能的开发进行了很多工作，但目前对于氧化锌薄膜的研究主要都集中在光电方面而对其传热特性还缺乏一定的认识.黄正兴等人采用瞬态热反应系统测试了溶胶凝胶法制备纤锌矿氧化锌薄膜的热导率，测试厚度范围为80～276 nm，测试结果显示薄膜导热系数随着膜厚的减小而减小[3]；Kulkarni等人运用MD方法模拟了长度为1.9～4.1nm的纤锌矿氧化锌纳米带的热物性，模拟结果远小于氧化锌块体结构的热物性数值，并且温度从500 k增大到1 500 k时，热导率下降了约52%[4]；Chantrenne和Ould-Lahoucine模拟了室温下膜厚范围在5.2～26 nm的纤锌矿氧化锌纳米薄膜热导率，发现热导率随着膜厚的增大而增大[5].基于前人运用MD方法对氧化锌热导率的相关研究，本文采用非平衡分子动力学方法和优化的Buckingham势函数，选定合理的初始条件后对厚度范围在几十纳米间的纤锌矿氧化锌纳米薄膜热导率进行数值模拟，并改变系统温度，计算分析薄膜厚度和温度的变化对薄膜热导率的不同影响.

1 非平衡态分子动力学方法

分子动力学方法是研究凝聚态物质的有效手段之一.根据是否具有温度梯度，分子动力学方法可以分为平衡分子动力学方法(EMD)和非平衡分子动力学方法(NEMD).本文采用的是非平衡分子动力学方法，其核心算法是对系统施加扰动，即通过加入温度梯度或热流的方法来建立非平衡导热过程，再根据傅立叶导热定律求得相对于整个系统的统计平均热导率.具体实现方法为：采用周期性边界条件建立起NEMD的计算模型如图1所示，系统在X、Y、Z三个方向都采用周期性边界条件，将整个模拟系统沿Z方向分成N层，保证每层中包含的原子数近似相等，再选定高温控制层和低温控制层，并控制两个温度区域保持固定的温差.选择控温层时需满足两个控温层中心的距离近似于模拟体系厚度的一半.

图1 非平衡稳态热运输模型

2 模拟过程

2.1 系综

分子动力学模拟的对象是多粒子体系，由于各种条件限制我们所模拟的系统粒子数目是有限的，但其统计物理的规律始终成立.系综意味着宏观条件下的一种约束，分子动力学模拟总是在一定系综下进行的，经常用到的平衡系综有微正则系综(NVE)、正则系综(NVT)、等温等压系综(NPT)和等温等焓系综(NPH).计算模拟时系综的选择是非常重要的，笔者选择的是在模拟过程中系统的总动量和总能量都守恒的微正则系综[6].在该系综中，体系不与外界发生能量交换，体系的粒子数和体积也不发生改变，因此想要找到可以满足系统能量的较精确的初始条件比较困难，对此可以采用速度标定的方法，即先给出一个合理的初始温度然后再不断进行能量增减的调整直到系统达到指定的能量值.

2.2 作用势

原子间作用势是用来描述原子间相互作用的函数，简称势函数.选取合理的势函数去描述原子间的相互作用力是分子模拟中至关重要的一步.势函数的选取通常是基于实验数据和相关规律的，本文根据经验选取了Buckingham势函数来描述氧化锌纳米薄膜内部原子间的相互作用力，Buckingham势[7]是一个广泛应用于氧化物材料特性研究的两体相互作用势模型，其表达式为

(1)

其中，i，j表示两种材料的参数；Eij表示i，j两种原子间的相互作用能.

模拟采用微正则系综和周期性边界条件下的NEMD方案，在周期性边界条件下，为了防止粒子间相互作用，要求模拟区域在该方向上的长度必须大于2×rc，而模拟系统的粒子数过多又不利于计算，所以在保证计算精度的前提下为尽量提高计算速度我们使系统在X和Y方向上都选取4个单元[8],Z方向上分别取112，128，144，160，176个基本单元建立超晶胞，对应的氧化锌纳米薄膜厚度分别是29.15,33.31,37.48,41.64，45.81 nm，对应的系统总原子数分别是7 168，8 192，9 216，10 240，11 264个.系统的时间步长设置为1.0 fs，每80步控制一次温度，每200步保存一次输出数据，所模拟的计算总步数为200万步.

3 结果与讨论

分子动力学方法中判断系统达到稳态非平衡状态的方法有3种，在这里我们选择通过观察温度分布曲线的光滑程度来判定模拟的非平衡过程是否进入稳态.以薄膜厚度为41.64 nm的4×4×160模拟系统为例，经过200万步的计算后得到的温度分布曲线图如图2所示，可以看到高温层和低温层之间整个温度分布曲线图是较为光滑的，系统平均温度为300 k，与高温330 k和低温270 k分别相差不到10%，控温层标记的浅色圆点表示控温区内由于有热流通过而存在着温差.以上可以判断系统模拟在进行了200万步后达到了平衡状态.

图2 氧化锌薄膜各层温度分布曲线图

3.1 纤锌矿氧化锌薄膜厚度与热导率的关系

表1 不同厚度的氧化锌薄膜热导率计算结果

表1是不同厚度的氧化锌薄膜在室温300 k下模拟热导率的计算结果，由表1可知薄膜厚度范围在29.15～45.81 nm的氧化锌薄膜的热导率在0.909 95～1.251 26 W/(m·K)之间.

热导率随膜厚的变化关系如图3所示，可以直观地观察到，在所计算的厚度范围内，薄膜材料的热导率远小于体材料的热导率(没有缺陷的氧化锌体材料热导率[9]可达到100 W/(m·K))，并随着薄膜厚度的增加而增加，接近于线性变化，表现出明显的尺寸效应.

出现上述效应的原因可做如下解释：

声子是半导体材料中主要的热载流子，声子间的相互作用甚至可以决定材料的传热性能.根据声子气动力学理论模型[10]，导出晶格的热导率表示为

(2)

其中，C为单位体积的声子比热，ν为声子平均速度，leff为材料的有效声子平均自由程，所谓平均自由程，表示的意义为声子在前后两次碰撞之间平均走过的距离.由于我们模拟的氧化锌薄膜厚度范围接近它体材料的声子平均自由程，有效声子的平均自由程近似满足

(3)

参数m由分子动力学模拟给出或者通过实验测试得出.假设声子的平均速度和比热都与薄膜厚度无关，结合以上两式可得出薄膜导热系数和膜厚的关系为

(4)

其中，α表示体材料的导热系数，b为导热系数的尺度依赖系数，所以导热系数的倒数与薄膜厚度的倒数成正比关系.影响声子平均自由程大小的因素很多[11]，有声子与声子之间的内部碰撞，晶格缺陷和杂质对声子的散射，晶界及晶体边界对声子的散射等.在微尺度传热过程中，由于微纳米结构的表面积和体积比显著上升，边界散射作用增强并成为决定热导率系数大小的最重要因素.所以当薄膜厚度从29.15 nm增大到45.81 nm时，声子的边界散射作用减小，故薄膜的热导率随之增大，表现出了明显的尺寸效应.

3.2 不同温度对纤锌矿氧化锌薄膜热导率的影响

在氧化锌薄膜热导率的模拟过程中，温度的设定是极其重要的参数，而在氧化锌薄膜的实际应用中，温度也是影响其热物性的关键因素，所以研究温度对氧化锌薄膜热导率的影响具有重要意义.以薄膜厚度为37.48 nm的4×4×144模拟系统为例，表2显示了系统温度从100 K上升到500 K时氧化锌薄膜的热导率由1.155 84 W/(m·K)下降到0.946 28W/(m·K).

表2 温度不同的氧化锌薄膜热导率计算结果

由图4可以观察到系统温度从100 K升高到500 K时，薄膜热导率减小的趋势显著，与文献[12]的结论相符.文献[12]的作者Alvarez-Quintana模拟了厚度为40，80，120，180 nm的多晶纤锌矿氧化锌薄膜以及块体纤锌矿氧化锌的热导率在温度从30 K上升到300 K时的变化情况.本文所模拟的温度范围与文献[12]有所重合，且热导率变化趋势都呈现明显的下降趋势，说明本文的计算结果是合理可信的.

图3 氧化锌薄膜热导率随膜厚的变化 图4 温度对氧化锌薄膜热导率的影响

热导率下降主要是因为随着系统温度升高，声子的振动频率增大，声子间相互作用也随之加强，导致声子偏离平衡位置的距离越来越远，引起声子的强烈散射，声子平均自由程变小，致使热导率数值下降[13].

4 结论

本文采用分子动力学模拟中的非平衡态(NEMD)方法，选取Buckingham势函数和NVE系综对纤锌矿氧化锌纳米薄膜在Z方向上热导率进行数值模拟，讨论分析了薄膜厚度和温度对其热导率的影响，并应用气动理论对研究结果进行解释分析.计算结果表明：薄膜在Z方向的热导率随着薄膜厚度的增大而增大，表现出明显的尺寸效应；随着系统温度的升高，声子散射作用加强，声子平均自由程减小，薄膜热导率也随之减小.该研究结果对人们进一步了解氧化锌纳米薄膜的传热特性提供了有益的参考价值.