网站地图关于合作招聘帮助

热线：400-152-6858

测试狗科研服务

预存定制

登录

|

注册

Document

当前位置：文库百科 › 文章详情

分子动力学(LAMMPS)——实例讲解：NEMD方法计算晶格热导率

来源：科学10分钟时间：2021-05-18 17:58:03 浏览：6269次

热导率又称导热系数，反映物质的热传导能力，其定义为单位温度梯度（在1m长度内温度降低1K）在单位时间内经单位导热面所传递的热量。材料的热导率是由电子热导率和晶格热导率共同贡献的。在温度较高的情况下，材料的热导率主要由电子热导率主导；在温度较低的情况，主要由晶格热导率主导。

实验中测材料热导率是在材料两端施加热流，此时可以测得沿热流方向的横截面上的热流通量J以及温度梯度∇T ，而J和∇T 具有公式（1）的线性关系，比例系数K 就是材料的热导率。

J = -k ∇T （1）

本文我们将使用由美国能源部桑迪亚国家实验室设计开发的分子动学软件包LAMMPS[1]复现文献[2]中用分子动力学计算金属钨晶格热导率的方法。文献[2]采用的是非平衡分子动力学(NEMD) [3-6]方法计算晶格热导率。计算过程中的模拟盒子示意图如图1所示，当模拟盒子在某一温度下平衡后，将模拟盒子中间作为热端，施加热流∆ε；模拟盒子两端作为冷端对热流进行抽取，在抽取的过程当中体系中原子的动量是从热端到冷端不断交换的，此时体系处于一个非平衡的状态。因此该过程也叫非平衡分子动力学(NEMD)过程。

在施加热流后，将模拟盒子沿热流方向分成许多层，每个特定的timesteps (如10000步)将每一层的温度进行统计，当传热过程稳定后 (这里是指统计的温度保持基本稳定，一般在加热流数万步到十万步后)，对特定timesteps间隔统计的分层温度做平均并进行线性拟合得到如图2所示的温度梯度曲线（黑色）。得到温度梯度后根据公式（2）得到晶格热导率k。公式（2）中A为模拟盒子的横截面积，M为热交换的频率，t为单位时间。最后，结合公式（1）和公式（2）得出金属钨的晶格热导率。笔者取的∆ε/t =1eV/ps，A=16Å²（笔者这里设置的有些不合理但对计算结果影响不大），M =10。笔者计算得出W在300K时延<100>方向的热导率是16.7 W/KM，和文献的计算结果17.28034 W/KM基本一致[2]。

图1. NEMD计算热导率的三维模型示意图

图2. 模拟盒子分层温度统计曲线(红色),温度梯度拟合曲线(黑色)

在计算过程当中，LAMMPS的in文件中的传热过程设置脚本段如下：

unfix 1

fix 2 all nve

fix hot all heat 10 ${power} region hot

fix cold all heat 10 -${power} region cold

其中，unfix 1命令取消的是前面作特定温度下平衡的系综。其余三段为热流为传热过程的设置。

计算结束后，最关键的是对模拟盒子每一层原子温度的统计，为此我们写了一个python小脚本，可以直接运行python3 thermal_cond.py profile.txt 200 25 45 > result.txt轻松处理，python的用法，及函数的定义，脚本thermal_cond.py中都有注释。处理后得到文本文件result.txt，该文件中的信息是模拟盒子每一层温度与层中心到模拟盒子左端距离的关系。直接将result.txt拖到作图软件中，或用gnuplot等程序作图即可获得火山峰状的温度统计曲线。最后截取其中的直线段拟合即可得温度梯度∇T，本文的∇T=-0.6K/Å。

最后，介绍几个可以使结果更优的方法：

（1）将模拟盒子三个方向边长，建成W单包三个方向边长的倍数，这样模拟盒子具有完美晶体的周期性。

（2）将采集温度的模拟时间尽量增大，这样能够采集更多的样本，统计的温度就会更加准确。

该算例以及温度梯度处理脚本可以在网盘下载：

链接：https://pan.baidu.com/s/1zkpULJ2T12NFkgnQOt5KOA

提取码：CSGO

参考文献：

[1] Plimpton S 1995 J. Comput. Phys. 117 1

[2] Baoqin Fu,Wensheng Lai,Yue Yuan, et al. 2012 J. Nucl. Mater. 427 1-3

[3] P. K. Schelling, S.R. Phillpot, P. Keblinski. 2002 Phys. Rev. B 65 144306.

[4] A. Maiti, G.D. Mahan, S.T. Pantelides. 1997 Solid State Commun. 102 517.

[5] A. Maeda. 1995 Phys. Rev. E 52 34.