目录

界面热导的模拟

随着现代电子器件尺寸的不断下降,我们可以看到纳米器件的复杂性会不断增加。由两种以上材料组成的器件的热导率是很难预测的。在大多数情况下,器件工程师想要利用不同材料构成界面来增大或减小热导率。

在本教程中,你将学习如何利用非平衡态分子动力学结合经验力场(势)的方法模拟热流通过界面,计算界面处的热导(热阻)。

提示

本教程使用特定版本的QuantumATK创建,因此涉及的截图和脚本参数可能与您实际使用的版本略有区别,请在学习时务必注意。

简介

通过在热源和冷源间施加温度梯度,可以模拟经硅晶体的晶界的热流。

理论

在给定温度梯度下,块体材料的热导率 $K_{bulk}$ 与热流的关系由傅立叶方程给出:

$\frac{\langle dQ/dt\rangle}{A} = -\kappa_{bulk}\langle dT/dz\rangle$

在这里, ⟨dQ/dt⟩ 代表沿着温度梯度的平均热流(也就是单位时间传输的动能),A 是垂直于热流方向的横截面积。

处在热源和冷源中间的界面通常会使温度分布产生一个 ΔT 的突变(温度的“跳跃”),因为界面会对热流产生额外的热阻 R。在这种情况下,界面处热导 G,通常被称作卡皮查(Kapitza)热导,被定义为

$G=\frac{1}{R}=\frac{\langle dQ/dt\rangle}{A\Delta T}$

方法

界面导热性质可以通过分子动力学模拟进行计算。最普遍的方法是在分子动力学模拟中引入温度梯度,也就是所谓的非平衡动量交换或反向非平衡分子动力学方法1)2)。在这种方法下,两端的两组原子被分别选择作为热源和冷源。在模拟期间,冷源中最热的原子和热源中最冷的原子根据特殊的标准进行动量交换,保证总能量和动量守恒3)。这种外部的热流则通过系统内由热源流向冷源的热流得到平衡,使体系内产生可观测的温度梯度(以及有界面时产生的温度的跳跃)。

制作晶界

这里考查一个由硅单晶的两种不同晶向形成的晶界模型。在模拟中,虽然你可以让热流沿着任意方向。但是只有沿着z方向的温度梯度可以求算。因此,构建体系时必须确保界面结构指向正确的方向。必要时可以通过旋转体系确保热流方向沿着z轴。

下面在硅的(100)表面和(110)表面之间创建界面。打开 Builder,通过 Add ‣ From database 添加 alpha-silicon,并且:

为了避免因连接两种不同晶格常量的单胞产生太大的形变,需要选择合适的重复以获得相称的单胞尺寸。单胞尺寸可以在 Bulk tools ‣ Lattice parameters 中观察并修改。

(100)面单胞的尺寸为 5.43 x 5.43 x 5.43 $Å^3$, ,(110)面 为3.84 x 5.43 x 3.84 $Å^3$。因此,(100)面的单胞应该沿着 ABC 分别进行 5x5x25 次的重复,(110)则需要 7x5x35 的重复,以使提供的超胞有着几乎相等的尺寸。你可以用 Bulk tools ‣ Repeat 插件去调节重复。设置想要的重复次数并点击Apply.

提示

一般情况下,对于非平衡热流模拟,在热流方向选择足够大超胞长度(假设是C方向)是很重要的,而结果对横截面区域的几何结构不太敏感(即,AxB)。因此,模拟所使用的超胞是典型的棍状结构。

创建一个界面可能是一件繁琐而枯燥的事情,这时可以使用 Interface Builder 帮你更简单灵活的完成这个建模过程。打开右边的 Builders > Interface 选项,将(100)超胞拖拽到第一个方框中,将 (110) 超胞放到第二个方框中。

界面的预览图应该会自动出现(如果没有则点击 Preview)。你将注意到在两个表面间有一条缝隙。你可以通过点击 Shift Surfaces 合拢该缝隙,在 z 方向填入一个负的偏移量直到两个表面间出现化学键。不仅如此,你也可以沿 y 方向调整对齐两表面。自动生成的默认值作为初始结构一般是不错的。

提示

显示出来的化学键仅仅是作为视觉辅助,没有严格的物理意义。对于界面建模来说,界面处的化学键提示了初始猜测的结构是否合适。最终的界面结构将通过分子动力学优化、退火等过程使之达到平衡结构,最大程度地接近真实情况。

在此计算中超胞在 c 方向上是非周期的,这意味着系统上下两端被两表面截断,热流仅能沿着 z 轴的正方向。横向尺寸(A和B)保留周期性。

为获得这样的设置,可以简单地增加 c 方向的长度。你可以通过打开 Bulk Tools ‣ Lattice Parameters,将 Cz 项由 268 Å 增加到 300 Å 来实现。你需要确保用的是笛卡尔坐标。

最后,片层结构应该置于超胞的中间。因此,打开 Coordinate Tools ‣ Center,不勾选 A 和 B,点击 Apply.

提示

除了这种两端终止于表面的界面模型之外,还有一种是模拟一个完整的周期性体系,使热流从热源处沿 z 轴正负双向流出4)5)。这种模型设置必须包含两个相同的界面,模拟的单胞长度也增加了一倍。

设置平衡模拟

在模拟实际热流之前,体系结构必须要达到很好的平衡结构,以减小因构造界面人为产生的缺陷,去除横向应力;最后得到一个在目标温度下平衡的体系。为了达到这个目的,你需要通过点击 Send 图标将最终确定的界面发送到Script Generator ,并添加 3 个模块:

将默认的输出文件名改为 Silicon_grain_boundary.nc,并打开 New Calculator 模块

打开 OptimizeGeometry 模块,将最大力参数调到 0.1 eV/Å 作为阈值。其余参数保留原默认值,关闭小窗口。

提示

优化会去除形成界面而出现的较大初始不稳定力,这是必要的。但是由于界面的缘故,微小的力是很可能是无法避免的。

接下来的步骤是用分子动力学手段让横向晶胞矢量得到弛豫、平衡,使之达到目标温度。整个热输运模拟过程要在平均 600K 的温度下进行。为达到这个目的,打开 MolecularDynamics 模块,将参数设置成 NPT 模拟,这是一种在常温常压下的模拟。

最后设置应该如下所示:

将脚本发送到 Editor 来增加一些额外的功能,因为你必须在结果中提取 A- 和 B-向量的平均长度。拖到底部,在脚本的最底下加上下列几行。

# Get the trajectory of cell vectors.
cells = md_trajectory.cells()
# Extract the xx- and yy-components.
vector_xx = [cell[0, 0].inUnitsOf(Angstrom) for cell in cells]
vector_yy = [cell[1, 1].inUnitsOf(Angstrom) for cell in cells]
# Calculate and print the mean value.
print "Average xx-component:  ", numpy.mean(vector_xx)
print "Average yy-component:  ", numpy.mean(vector_yy)

将脚本发送到 JobManager 并运行。

在模拟结束后,A 和 B 的平均晶格矢量将会在 log 文件中的底端被打印出来。选择 Lab floor 中的 trajectory 并在 Movie Tool 中打开它。播放运动轨迹并将其最后一帧通过 标志发送到 Builder 中。

在 Builder,打开 Bulk Tools ‣ Lattice parameters 对话框。选择 Keep Fractional ,将晶胞矢量中的 xx- 和 yy-元修改成从模拟中获得的平均值。关闭 Lattice parameters对话框。在以后的模拟中,晶格常数将会被固定。

设置非平衡动量交换模拟

现在你必须定义热流模拟过程中的热源和冷源区域,放大结构最左端,选择最外面的8个原子层(即两个单胞),如下图所示。

首先,你运行另一个处在 600K 温度下的平衡模拟,这次在恒定体积下准备一个平衡态体系。打开第一个分子动力学模块 并,

在第二个 MolecularDynamics 模块中,最终的非平衡模拟将被设置。为了实现这一目的:

提示

降低交换间隔的值将增加每次模拟的转移动能,这会造成更大的温度梯度和更明显的温度分布。这将使计算变得更容易,但是一般情况下要小心温度梯度不能太大,因为量级可能远远超实验中的温度梯度。

你可以使用 tag 里的 “heat_source” 和 “heat_sink” 确定各自热源和冷源区域,这是之前在构建模型时已经定义好的。最终的设置应该如下所示。

提示

作为演示计算,所选择的步数设置可以非常快的完成计算模拟,对于实际模拟,有必要多进行几个尝试计算并验证结果的一致性,应该考虑的是增加步数,尤其是对非平衡模拟来说。

分析结果

模拟完成后,轨迹文件将出现在 LabFloor 里。用 Movie Tool 打开非平衡轨迹文件,你会看到在整个模拟过程中平均温度是一个常数,片层没有漂移,因为求算总体的平均温度不可能体现温度的空间分布。

观察界面热导

MD Analyzer 插件打开同一个的轨迹文件,研究热导率,

注意

在这里,你应该仔细检查结果对模拟时长的收敛性。如果没有获得收敛,你应该从之前的最后一帧开始重新运行非平衡模拟(这次不需要 NVT 平衡计算了)

提示

前 100ps 尚未建立非平衡热流,因此应该被舍弃。

因此,例如你可以比较 100 ps 与 400 ps 间的分布和 400 ps 和 800 ps 间的分布。用在工具栏中的 Zoom Tool 去放大观察热导分布图中最显著的特征,或者点击 Home 返回初始视图。除了统计性的起伏外,分布应该是一样的。

温度分布显示了一些典型的特征。首先,在分布的首段和末段,可以看到一个正向和一个负向的峰值。这两个部分分别是热源和冷源。这种非线性部分来源于人为的动能转移,这是由动量交换产生的。你应该在你分析的时候去掉这一部分。

第二个显著的特征是在 z=150Å 处,温度会突然跃迁。这是界面区域,因此有额外的热阻。让我们放大该区域更进一步观察一下。

在 T1=625K 和 T2=575K 的间温度变化为 ΔT=50 K。为了计算 Kapitza 热导你还需要知道平均热流。该值可以从日志文件中的 nlprint output 中找到。在该例中,可以获得的热流值 ⟨dQ/dt⟩=0.00191eV/fs。最后你需要的是横截面积,这可以通过 A 和 B 的晶格矢量长度计算得到。根据上面的公式,结合这些数值,就可以计算 Kapitza 热导了。

鉴于单位必须要正确转换,计算热导值最简单的方法是利用 QuantumATK 的内建单位引擎。

dq_dt = 0.00191*eV/fs
A = 26.988*27.118*Angstrom**2
delta_T = 50.0*Kelvin
G = dq_dt / delta_T / A
print "Kapitza conductance in Watt/Kelvin/Meter**2:  ",G.inUnitsOf(Joule/Second/Kelvin/Meter**2)

研究晶界尺寸相关热导率

在该例中,你也可以调节热源或者冷源到边界的距离,大约为 100 Å,来研究晶界尺寸的影响。此时,你需要考虑的是热导率而不是界面热导。用 MD Analyzer 提取温度梯度,在 Zmin 和 Zmax 处输入合理的数值来选择分布的线性区域。

程序将在参数窗口中自动拟合斜率,显示温度梯度。在该例中,(100)取向晶粒的值为 -0.167 K/Å。你对右边的晶粒重复该过程,获得一条相似的斜率。左右两边的斜率不可能一样,这是由于统计起伏和不同的晶粒取向造成的。

现在我们可以在脚本中将计算得到的斜率代入到傅里叶方程中,利用单位时间平均转移的能量、横截面积计算晶粒内的热导率。

dq_dt = 0.00191*eV/fs
A = 26.988*27.118*Angstrom**2
dT_dz = 0.167*Kelvin/Angstrom
G = dq_dt / dT_dz / A
print "Grain conductivity in Watt/Kelvin/Meter:  ",G.inUnitsOf(Joule/Second/Kelvin/Meter)

结果是 25 W/K/m。对于(110)取向晶粒,热导率为 22.5 W/K/m,这些值与模拟的体块状Si的热导率相比是非常的小,在 600K 下,$\kappa_{bulk} = 150 W/K/m$,这是3) 中利用相同的 Stillinger-Weber 势的计算结果。这种结果是可以预料的,因为尺寸受限的晶粒的热导率要低于块体材料中相应的值 1) 。这是由于声子的平均自由程通常比 1 μm 要大得多,这意味着模拟的晶粒尺寸抑制了热导1),3) 。

为了研究晶粒尺寸相关性,你可以重复模拟,准备两种晶胞不同长度的结构,正如参考文献1) 所示。

块体热导率可以通过多次模拟获得,随着晶粒尺寸l的增加,绘制热导率倒数 1/κ(l) 和晶粒尺寸倒数1/l的关系,并由此外推到无限大晶粒尺寸。

参考

1)
F. Mueller-Plathe. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys., 106(14):6082–6085, 1997
2) , 3)
C. Nieto-Draghi and J. B. Avalos. Non-equilibrium momentum exchange algorithm for molecular dynamics simulation of heat flow in multicomponent systems. Mol. Phys., 101(14):2303–2307, 2003
4)
P. C. Howell. Comparison of molecular dynamics methods and interatomic potentials for calculating the thermal conductivity of silicon. The Journal of Chemical Physics, 137(22):, 2012
5)
Akbar Bagri, Sang-Pil Kim, Rodney S. Ruoff, and Vivek B. Shenoy. Thermal transport across twin grain boundaries in polycrystalline graphene from nonequilibrium molecular dynamics simulations. Nano Letters, 11(9):3917–3921, 2011