摘要:用嵌入位错线法和重合位置点阵法构建含有小角度和大角度倾斜角的双晶氧化锌纳米结构. 用非平衡分子动力学方法模拟双晶氧化锌在不同倾斜角度下的晶界能、卡皮查热阻, 并研究了样本长度和温度对卡皮查热阻和热导率的影响. 模拟结果表明, 晶界能在小角度区域随倾斜角线性增加, 而在大角度区域达到稳定, 与卡皮查热阻的变化趋势一致. 热导率随样本长度的增加而增加, 卡皮查热阻表现出相反的趋势. 然而随着温度的增加, 热导率和卡皮查热阻都减小. 通过比较含5.45°和38.94°晶界样本的声子态密度, 发现声子光学支对热传导的影响不大, 主要由声子声学支贡献, 大角度晶界对声子散射作用更强, 声学支波峰向低频率移动.
关键词: 倾斜角/
晶界能/
卡皮查热阻/
热导率

English Abstract

全文HTML

--> --> -->

氧化锌(ZnO)作为宽禁带半导体材料, 具有优异的压电性、光电性和气敏性等^[1-3]. 随着纳米材料制备技术的提高, 微纳尺度的ZnO在诸多领域有较为广泛和有效的应用^[4,5]. 由于较高的缺陷密度和表面积-体积比, ZnO纳米装置的热效率和热稳定性已经成为众多科研工作者关注的主要问题. 纳米结构ZnO中存在大量的晶界^[6], 当声子的平均自由程大于平均晶粒尺寸时, 晶界对声子的散射成为主要的散射机制. 关于热输运的许多基本问题还没有得到充分解释, 仍然需要发展不同方法来深入理解热输运细节.
目前, 国内外研究者主要对不同类型纳米结构ZnO的热性质开展研究. Zhang和Chen^[7]讨论了ZnO纳米线的直径、长度对热导率的影响, 发现其具有明显的尺寸效应, 同时引入合金、超晶格以及多孔结构等控制声子的热输运行为. El-Brolossy等^[8]采用溶剂热法制备了掺杂铝的ZnO纳米颗粒, 当颗粒尺寸减小时, 比热容明显提高, 而铝掺杂对比热容的影响较小; 由于存在几种嵌套的减热机制, 所测得的热导率比单晶ZnO小一个数量级左右. Giri等^[9]合成了氧化锌/对苯二酚(ZnO/HQ)纳米结构材料, 研究晶界密度对声子散射和热传导的影响, 结果表明晶界散射是导致热导率显著降低的主要原因. Liang和Shen^[10]探索多晶ZnO的声子和电子输运, 并确定卡皮查热阻和有效势垒高度与晶界间隙有关; 当材料中加入硫化锌(ZnS)纳米薄膜时, 卡皮查热阻增加了三倍, 势垒区域扩大了两倍. 这些研究表明纳米结构ZnO的热导率主要受到晶界对声子的散射作用的影响, 因此有必要对晶界如何影响热性质进行深层次, 多方面探讨.
声学失配模型^[11]和扩散失配模型^[12]通常被用于评估界面热阻. 然而, 这两个理论都有特定的应用范围并且不能考虑晶界的原子结构对纳米材料热传导的影响. 同时实验上需要花费大量的时间和精力得到需要的晶界结构, 例如具有特定倾斜角度的晶界, 使得定性地研究晶界对热传导的影响变得非常困难. 分子动力学(MD)模拟克服了传统的理论和实验方法的缺陷, 广泛地应用于设计不同晶界结构和研究晶界上的声子散射^[13]. 我们尝试采用这种方法去控制一定形式的晶界并得到相应的晶界性质与卡皮查热阻和热导率之间的关系.
本文采用非平衡分子动力学(NEMD)模拟研究不同倾斜角度的晶界对纳米结构ZnO的热性质的影响. 基于NEMD模拟和Frank-Bilby公式分别计算晶界能和位错密度. 选择几组小角度和大角度晶界作为对比, 研究在不同样本长度和温度下卡皮查热阻和热导率的变化. 最后利用声子态密度深入理解声子-晶界相互作用, 进一步解释倾斜角度对热阻的影响. 通过设计特定的晶界结构为获得想要的热性质提供可能.

2

-->

2.1.模　型

双晶模型通常被用于研究单一晶界的相关特征. 两个晶粒平行放置在界面两边, 将左边的晶粒沿Z轴顺时针旋转${\theta _{\rm{L}}}$角度, 将右边的晶粒逆时针旋转${\theta _{\rm{R}}}$角度, 定义晶界倾斜角${\theta _{\rm{M}}} = {\theta _{\rm{L}}} + {\theta _{\rm{R}}}$. 理论上倾斜角可以设定为任意值, 但是会产生很多原子排列混乱的缺陷晶界, 因此采用重合位置点阵方法构建合理的晶界来满足模拟要求^[14]. 在经过旋转后, 原始的两个晶粒在晶界处具有一定数量的重合位置点, 重合位置点与晶格点比例的倒数用Σ表示. 根据倾斜角, 晶界分为小角度和大角度晶界. 小角度晶界由沿[001]方向的边界位错构成, 构建了5个从5.45°到30.51°小角度晶界的双晶ZnO. 小角度晶界处的位错排列不适用于大角度晶界, 基于重合位置点阵法构建了5个从36.86°到67.38°大角度晶界的双晶ZnO. 所选晶界倾斜角与对应的Σ值见表1. 此外, 根据非线性共轭梯度法对样本进行能量最小化, 得到最稳定的结构. 选择5.45°小角度晶界和38.94°大角度晶界作为对比, 采用OVITO软件将其可视化^[15], 如图1所示. 可以清楚地看出5.45°晶界由位错线构成, 而38.94°晶界含有大量缺陷.

Tilt angle θ/(°)	Σ
5.45	221
11.42	101
20.02	149
25.98	89
30.51	65
36.86	5
38.94	9
53.13	5
60.0	3
67.38	13

表1对称倾斜角度与对应的Σ值
Table1.List of the symmetrical tilt angles and the corresponding Σ values.

图 1 含5.45°晶界(左)和38.94°晶界(右)的双晶ZnO结构
Figure1. Structures of the 5.45° grain boundary (GB) (left) and the 38.94° GB (right) in bicrystal ZnO.

2 -->

2.2.计算方法

NEMD方法的基本思想是在材料体系上施加热流. 本文采用的纳米尺度传热试样的原理模型如图2所示. 两个晶粒以不同的取向放置在模拟盒子中, 中间位置包含晶界. 为了避免边界原子的相互作用, 在热流传输方向上选择两层原子固定在模拟盒子两端. 包含五层原子的热源和热汇区域分别在固定原子层内侧. 热源区的温度要高于热汇区温度, 使体系中产生从左到右的温度梯度. 在热流传输方向采用固定边界条件, 通过设定不同的样本长度来研究卡皮查热阻和热导率的尺寸效应. 横截面方向采用周期性边界条件, 防止模拟过程中原子的运动超出边界而丢失, 保持系统中原子数量恒定, 同时消除横截面边界对原子受力的影响. 为了避免横截面尺寸对结果产生干扰, 将其大小定为10 × 10单位晶胞(46.3 ? × 46.3 ?).
图 2 NEMD模拟计算热性质的示意图
Figure2. Schematic diagram of the NEMD model for calculating the thermal properties.

所有的模拟过程都是基于LAMMPS软件进行的^[16], 利用麦克斯韦-玻尔兹曼速度分布对体系进行初始化. Buckingham势函数用于描述Zn—O和O—O的原子间相互作用, 设置参数下的模拟值与ZnO的实验结果有很好的一致性^[17], 其表达式为

${E_{\rm{R}}} = \frac{{{z_i}{z_j}{e^2}}}{r} + A\exp \left( {\frac{{ - r}}{B}} \right) - \frac{C}{{{r^6}}},$

其中对于Zn²⁺—O^2–化合键, A = 529.7 eV, B = 0.03581 nm, C = 0; 对于O^2–—O^2–化合键, A = 9547.96 eV, B = 0.021916 nm, C = 32 eV·A⁶.
首先在零温条件下对体系进行能量最小化, 优化原子位置; 然后在正则系综(NVT)下进行1 ns的平衡态模拟来达到控制的温度和体积; 最后在微正则系综(NVE)下施加温度梯度并测量热流值. 根据双晶样本的长度将其分割为一定数目的原子层. 选择热汇区中最热的原子和热源区中最冷的原子进行动能交换获得每个原子层中的平均温度. 动能交换的过程重复2 ns, 保证达到稳定热流, 其表达式为

$J = \frac{1}{{tA}}\sum\limits_n {\frac{m}{2}} \left( {v_{\rm{h}}^2 - v_{\rm{c}}^2} \right),$

式中t为模拟时间, A为样本的横截面积, n为动能交换次数, m为Zn和O原子的平均质量, v_h和v_c分别表示最热和最冷原子的速度.
在稳定状态下晶界的两侧建立线性温度梯度, 然而在体系的两端存在非线性区域. 在晶界处温度不连续, 具有非常明显的温度跳跃$\Delta T$. 卡皮查热阻R为晶界处的温度跳跃$\Delta T$与热流$J$之比,

$R = {{\Delta T}}/{J}.$

基于傅里叶定律的有效热导率${\kappa _{{\rm{eff}}}}$的表达式为

${\kappa _{{\rm{eff}}}} = - \frac{J}{{\partial T/\partial x}},$

式中$\partial T/\partial x$为温度梯度.

为了研究晶界结构对纳米ZnO热输运性质的影响, 首先计算了晶界能和位错密度, 它们是表征晶界特性的重要参数. 其中晶界能${E_{{\rm{GB}}}}$的表达式为

${E_{{\rm{GB}}}} = \frac{{{E_{{\rm{bulk}}}} - N{E_{\rm{c}}}}}{S},$

式中${E_{{\rm{bulk}}}}$为体系的总能量, N为模拟区域中的原子数量, ${E_{\rm{c}}}$为每个原子的内聚能, S为晶界面积. 晶界能随倾斜角的变化情况如图3所示. 当倾斜角在5.45°到36.86°区间时, 晶界能随倾斜角线性增加; 当倾斜角大于36.86°时, 晶界能几乎不发生变化, 达到稳定值约为2.75 J/m². 由于晶界能的变化规律在倾斜角为36.86°时发生改变, 将36.86°定义为过渡角. 当倾斜角小于过渡角时, 晶界的原子排列较疏松, 几何形状容易受到倾斜角度变化的影响, 随着倾斜角增大, 晶界的几何形状更加复杂, 表现为无定型形状; 而当倾斜角大于过渡角时, 晶界的原子排列紧密, 形成稳定的多边形几何形状, 此时倾斜角的改变对晶界能的影响几乎可以忽略. 对于小角度倾斜晶界, 可以用Read-Shockley模型计算材料的晶界能, 其表达式为
图 3 晶界能和位错密度与倾斜角的关系
Figure3. The GB energy and dislocation density as a function of misorientation angle.

${E_{{\rm{GB}}}} = {E_0}\theta (A - \ln \theta ),$

式中A为特征参数, 表示为$A = 1 + \ln (b/2\pi {r_0})$, 其中${r_0}$为位错核半径, $b$为伯格斯矢量; ${E_0}$为体系初始能量, 表示为${E_0} = Gb/4{\text{π}}(1 - v)$, 其中$G$为剪切模量, $v$为泊松比, ZnO的剪切模量$G=0.75$, 泊松比$v=0.337$. 根据此模型计算得到的晶界能满足MD的模拟值, 并且已经在Al₂O₃等材料中得到证实^[18].
位错密度$\rho $为特定倾斜角晶界单位长度的位错数量, 其表达式为^[19]

$\rho =\frac{1}{D} = \frac{{2\sin \left( {\theta /2} \right)}}{b},$

式中D为位错间距. 对于ZnO的六角纤锌矿结构, 其晶格常数为a = 3.28 ?, 伯格斯矢量的大小与晶格常数相等. 图3显示位错密度随倾斜角的变化与晶界能具有相同的趋势. 当倾斜角大于过渡角时, 近似地假设位错密度为固定值1.28 nm^–1. 这是由于在小角度区域, 每个晶界由一系列独立的位错组成. 位错之间的间隙随着倾斜角的增加而减小, 位错间距减小. 然而在大角度区域时, 位错之间相互重叠, 没有间隙, 位错间距不随倾斜角的改变而改变, 可以认为位错密度保持不变.
双晶ZnO的卡皮查热阻随倾斜角的变化关系如图4所示. 卡皮查热阻在小角度区域中几乎线性增加, 而在大角度区域趋于稳定. 这与晶界能随倾斜角的变化非常类似, 表明卡皮查热阻与晶界能之间有紧密的关联. 这一发现与之前UO₂^[20]和金刚石^[21]的研究结果一致. Read-Shockley的扩展模型能用于描述卡皮查热阻与倾斜角之间的关系, 其表达式为
图 4 卡皮查热阻与倾斜角的关系
Figure4. Kapitza resistance as a function of misorientation angle.

$R\left( \theta \right) = \left\{\!\begin{aligned}& \sin \left( {2\theta } \right)\! \left[ {R_{\rm{c}}^ < \!-\! R_{\rm{s}}^ < \ln \left( {\sin \left( {2\theta } \right)} \right)} \right], \\ & \qquad {0^\circ }\! < \!\theta \! < \!{45^\circ }, \\& R_{\rm{c}}^ < - R_{\rm{c}}^ > + \sin \left( {{{180}^\circ } - 2\theta } \right)\\ & \times\left[ {R_{\rm{c}}^ > - R_{\rm{s}}^ > \ln \left( {\sin \left( {{{180}^\circ } - 2\theta } \right)} \right)} \right],\\ & \qquad {45^\circ } < \theta < {90^\circ },\\[-10pt] \end{aligned} \right.$

式中, R_c和R_s是可调参数, 通过拟合MD数据可以确定$ R_{\rm c}^ {<} \! = $3.474 m²·K/GW, $ R_{\rm s}^ <\! = $0.348 m²·K/GW, $ R_{\rm c}^ {>} $= –0.214 m²K/GW, $ R_{\rm s}^ > $ = 0.146 m²K/GW. 图4显示拟合曲线与模拟值表现出很好的一致性, 进一步验证模拟结果的正确性.
在不同样本长度条件下研究双晶ZnO的卡皮查热阻和热导率的尺寸效应. 作为对比, 分别选择包含5.45°, 20.02°, 38.94°和60.0°晶界的样本, 样本平均长度在23.2—185.2 nm之间, 平均模拟温度设定为300 K. 从图5可以看出, 当样本长度在23.2—92.6 nm之间时, 卡皮查热阻随长度的增加急剧减小, 具有明显的尺寸效应. 结合图6可以看出, 热导率在此长度范围内显著增加, 而当长度大于92.6 nm时, 卡皮查热阻和热导率的变化趋于稳定. 通过线性外推法计算得到单晶ZnO的热导率为98.5 W/mK, 含有晶界样本的热导率的最大值仍小于单晶热导率. 声子是ZnO中主要的热载流子, 基于经典动力学理论, 晶格热导率的表达式为
图 5 300 K时卡皮查热阻与样本长度的关系
Figure5. Kapitza resistance as a function of sample length at 300 K.

图 6 热导率与样本长度的关系
Figure6. Thermal conductivity as a function of sample length.

$\kappa = ({1}/{3}) cv{\varLambda _{\rm{b}}}, $

其中$c$为比热容, $v$为声子群速度, ${\varLambda _{\rm{b}}}$为整体平均自由程. 由于材料结构变化对比热容和群速度的影响很小, 晶界结构主要影响声子的散射, 我们专注于平均自由程对尺寸的依赖性. 当样本长度远小于声子平均自由程时, 此时声子在双晶中处于弹道输运模式, 声子-晶界散射作为主要散射机制, 对声子的抑制作用明显, 散射强度随长度的增加而减小; 当样本长度远大于声子平均自由程时, 此时声子处于扩散输运模式, 声子-声子散射起主导作用, 但由于其散射强度较弱, 对声子传输几乎没有阻碍作用.
卡皮查热阻和热导率还会受到温度的影响. 选择含有5.45°和38.94°晶界的双晶样本, 设定长度为23.2 nm. 在300—700 K的温度范围内计算得到的结果如图7和图8所示. 从图7中看到, 38.94°晶界样本的卡皮查热阻随温度的增加显著减小, 而对于5.45°晶界样本其变化曲线相对平坦. 同时在不同温度下, 38.94°晶界样本的卡皮查热阻总是大于5.45°晶界样本对应值. 卡皮查热阻主要来源于声子在晶界处的散射, 当温度升高时穿过晶界的声子数量增加, 有更多的声子参与热传导过程, 导致卡皮查热阻减小. 图8显示, 对于5.45°晶界样本其热导率同样随温度的增加而减小, 而对于38.94°晶界样本, 当温度小于400 K时, 热导率呈增加趋势, 然后迅速减小. 这是由于在低温区域, 被激发的声子数量较少, 晶格的振动受到强烈限制, 热导率主要受晶界对长波声子散射的影响. 随着温度的增加, 具有长波长、低频率的声子在晶界处的穿透率大于散射率, 表现为热导率增加. 在高温区域, 材料发生热软化并且晶格剧烈振动, 激发出大量的低波长、高频率声子. 声子振动幅度随温度的增加而增加, 高频率声子之间的散射作用占据主导, 使声子的平均自由程减小, 从而导致热导率进一步降低.
图 7 5.45°和38.95°晶界结构的卡皮查热阻随温度的变化
Figure7. Temperature-dependent Kapitza resistance for the 5.45° and 38.94° GBs.

图 8 5.45°和38.95°晶界结构的热导率随温度的变化
Figure8. Temperature-dependent thermal conductivity for the 5.45° and 38.94° GBs.

晶界处的原子排列对声子模式有重要的影响. 为了深层次理解晶界对卡皮查热阻和热导率的影响, 基于速度自相关函数的傅里叶转换计算得到声子态密度, 其表达式为^[22]

$D\left( \omega \right) = \int_0^\infty {v\left( 0 \right)} \cdot v\left( t \right){{\rm{e}}^{ - {\rm{i}}\omega t}}{\rm{d}}t, $

式中$v\left( t \right)$为原子在时间t时的速度, $v\left( 0 \right)$为原子初速度, $\omega $为声子振动频率. 比较5.45°和38.94°晶界样本的声子态密度, 结果如图9所示, 设定固定长度为92.6 nm, 温度为300 K. 在ZnO中存在4个声子振动模式, 分别是低频率的横向声学分支(TA)、纵向声学分支(LA)和高频率的横向光学分支(TO)、纵向光学分支(LO). 两个样本的声子态密度在大于9 THz时几乎是重合的, 因为光学声子的群速度非常小, 对晶体内部的热传导影响可以忽略. 然而在1—7 THz频率范围, 声学模式存在显著差异. 相比于5.45°晶界样本, 38.94°晶界样本的TA和LA波峰移动到更低频率的位置, 表明具有相同频率的声子能够通过5.45°晶界, 而在38.94°晶界处发生散射, 转变为更低频率的声子. 声子模式失配现象进一步说明大角度晶界的原子排列比小角度晶界更加混乱, 对声子散射作用更强, 最终造成大角度晶界样本的卡皮查热阻大于小角度晶界的卡皮查热阻.
图 9 5.45°和38.94°晶界结构的声子态密度比较
Figure9. Phonon density of states for the 5.45° and 38.94° GBs.

本文利用NEMD方法对含有倾斜晶界的双晶ZnO的热输运性质进行模拟, 系统地分析了晶界倾斜角、样本长度以及温度对卡皮查热阻和有效热导率的影响, 研究发现:
1)在倾斜角小于36.86°时, 晶界能和位错密度随倾斜角线性增加, 当倾斜角大于36.86°时, 晶界能和位错密度几乎稳定不变;
2)卡皮查热阻随倾斜角的变化趋势与晶界能和位错密度的变化趋势相同, 证明它们之间有较强的相关性, 基于Read-Shockley的扩展模型的计算值与模拟值拟合良好;
3)卡皮查热阻和有效热导率具有明显的尺寸效应和温度效应, 在小角度区域受倾斜角的影响较大, 而在大角度区域受倾斜角的影响较小;
4)通过计算声子态密度得到光学分支几乎对热输运没有影响, 主要由声学分支贡献. 大角度晶界的声子-晶面散射作用更强, 声子频率降低, 声学分支的波峰向低频移动.