当前位置：首页 > news >正文

lammps应用于热电材料

news 来源：原创 2025/6/11 9:28:06

文章目录

- 1.热传导理论
- - 1.热导率
  - 2.晶格振动
  - 3.晶体热容
  - 4.声子平均自由程
  - 5.傅里叶定律
- 2.lammps计算Ar热导率
- 3.lammps模拟SiGe热电材料
- 4.平衡分子动力学(EMD)

1.热传导理论

1.热导率

热传递机制随介质材料相的不同而改变：固体(热传导)、液体(热对流)、气体(对流和辐射)
热传导：热量从固体材料高温向低温部分转移的过程
热力学第二定律(克劳修斯表述)：“不可能把热量从低温物体传向高温物体而不引起其他变化。”
在自然条件下热量只能从高温物体向低温物体转移，而不能由低温物体自动向高温物体转移，这个转变过程是自发、不可逆的
热导率：又称导热系数(coefficient of thermal conductivity)，物质导热能力的量度
若物体的截面积为 1m²，相距 1m 的两个平面的温度差为 1K，则1秒内流过界面的热量就规定为该物质的热导率，标准单位为 W/(m·K)
T=0K，原子固定在平衡位置处，保持能量最低的状态
当T>0K时，原子在格点原子势场中热振动

2.晶格振动

材料各种热学性能的物理本质，均与其晶格热振动有关
绝缘材料中的导热现象主要由原子(晶格)的振动所引起
谐振子量子理论： $H=\frac{1}{2}m\dot{x}^{2}+\frac{1}{2}kx^{2}=\frac{1}{2}\dot{q}^{2}+\frac{1}{2}\omega^{2}q^{2}$ 量子力学中，谐振子有分立的本征能量： $\varepsilon_{n} =(n+\frac{1}{2})\hbar \omega,n=0,1,...,\infty$
材料热导率的差异来源于晶格原子振动对应的传递能力，声子(晶格振动的能量量子)之间的碰撞导致热量的传递
声子的能量： $\hbar\omega，\omega=\sqrt{\frac{k}{m}}$ 原子/分子间的键越强，导热能力越强
声子对热导的贡献：
1.长波声子具有更长的波长，易发生倒逆过程，因此对热传导性质影响更大
2.通常声学声子对热传导的贡献远大于光学声子
热导率公式： $\kappa =\frac{1}{3}C\bar{v}\bar{\lambda}$ 1.与固体的比热 C 成正比，因为比热越大，每个固体分子可带的热量也越多
2.与声子的平均速率 $\bar{v}$ 成正比，因为声子速度越快，热也传送得越快
3.与声子的平均自由程 $\bar{\lambda}$ (与声子碰撞频率相关)成正比， $\bar{\lambda}$ 越大，声子携带的热能传得越远
其中，声子的平均速度与温度变化的关联性不大，多数情况可以作为常数处理。而平均自由程比较难计算，与温度和晶格声子碰撞的方式、频率等因素有关

3.晶体热容

$C_{V}=(\frac{\partial U}{\partial T})_{V}$ $C_{V}=C_{V}^{l}+C_{V}^{e}$ 其中， $C_{V}^{l}$ 为晶格振动热容， $C_{V}^{e}$ 为晶体电子热容

固体热容主要有两部分：一个是来源于晶格振动(声子)，称为晶格热容；一个来源于电子的热运动，称为电子热容。除非在很低温，电子热运动的贡献往往很小。通常情况下， $C_{V}^{e}<<C_{V}^{l}$
实验规律：
1.在高温时，晶体的热容为： $3Nk_{B}=3vR$ N 为晶体中原子的个数，k_B = 1.38 × 10^-23 J·K^-1 为玻尔兹曼常量；v 为晶体中原子摩尔数，R = 8.31 J/K mol为普适气体常数
2.在低温时，绝缘体热容按 $T^{3}$ 趋于零，导体热容按 T 趋于零
晶格振动的能量是量子化的，频率为 $\omega$ 的声子能量为： $E_{n}=(n+\frac{1}{2})\hbar\omega$ $\frac{1}{2}\hbar\omega$ 代表零点振动能，对热容没有贡献 $E_{n}=n\hbar\omega$ $n$ 是频率为 $\omega$ 的谐振子的平均声子数，据玻色统计理论： $n(\omega)=\frac{1}{e^{\hbar\omega/k_{B}T}-1}$ 温度为 $T$ 时，频率为 $\omega$ 的振动的能量： $\bar{E}(\omega)=\frac{\hbar\omega}{e^{\hbar\omega/k_{B}T}-1}$
爱因斯坦模型：
晶体由 N 个原子组成，每个原子有3个自由度，共有 3N 个分立的振动频率，晶体内能： $U=\sum_{i=1}^{3N}\bar{E}(\omega_{i})=\sum_{i=1}^{3N}\frac{\hbar\omega_{i}}{e^{\hbar\omega_{i}/k_{B}T}-1}$ 爱因斯坦模型：假设所有原子以同样的频率振动 $U=\sum_{i=1}^{3N}\bar{E}(\omega_{i})=\frac{3N\hbar\omega}{e^{\hbar\omega/k_{B}T}-1}$
高温下，比热为 $3NK_{b}$ ，与经典结果一致；低温下与实验偏差较大
德拜模型：
若频率分布可用一个积分函数表示： $\int_{0}^{\omega_{m}}g(\omega)d\omega=3N$
$g(\omega)d\omega$ 表示在频率范围 $\omega \sim\omega+d\omega$ 可取的频率数， $\omega_{m}$ 为最大的频率数， $q$ 和 $\omega$ 为准连续 $U=\sum_{i=1}^{3N}\frac{\hbar\omega_{i}}{e^{\hbar\omega_{i}/k_{B}T}-1}=\int_{0}^{\omega_{m}}\frac{\hbar\omega}{e^{\hbar\omega/k_{b}T}-1}g(\omega)d\omega$ 热容： $C_{V}=(\frac{\partial U}{\partial T} )_{V}=\int_{0}^{\omega_{m}}k_{B}(\frac{\hbar\omega}{k_{B}T})^{2}\frac{e^{\hbar\omega/k_{B}T}g(\omega)d\omega}{(e^{\hbar\omega/k_{B}T}-1)^{2}}$ 变量代换： $f_{D}(\frac{\theta_{D}}{T})=3(\frac{T}{\theta})^{3}\int_{0}^{\frac{\theta_{D}}{T}}\frac{e^{x}x^{4}}{(e^{x}-1)^{2}}dx$ 低温下： $C_{v}\sim12\pi^{4}R\frac{(T/\theta_{D})^{3}}{5}$
对于大多数固体来说，Debye模型可以给出相当精确的结果

4.声子平均自由程

$\kappa=\frac{1}{3}C\bar{v}\bar{\lambda}$ 与声子数密度成反比： $\bar{\lambda}\propto\frac{1}{n}$ 频率为 $\omega$ 的平均声子数： $n(\omega)=\frac{1}{e^{\hbar\omega/k_{B}T}-1}$ 高温下： $\bar{\lambda}\propto T^{-1}，C为常数，\kappa\propto T^{-1}$ 低温下： $\bar{\lambda}\propto e^{\frac{\alpha}{T}}，C\propto T^{3}，\kappa\propto e^{\frac{a}{T}}$

热容随温度增加，声子的平均自由程 $\bar{\lambda}$ 随温度减小，两者竞争使得 $\kappa$ 出现极小值
低温：声子平均自由程 $\bar{\lambda}$ 近似为粒径大小， $C_{V}$ 随 $T^{3}$ 增大， $\kappa$ 值增大
高温：声子的平均自由程 $\bar{\lambda}$ 随温度升高减小，而声子热容 $C$ 趋近于常数， $\kappa$ 值降低
热电材料：
塞贝克(Seebeck)效应：用于热电能量转换，如节能、特殊电源、地热发电
珀尔贴(Peltier)效应：用于热电致冷，如冰箱、空调
评价热电转化效率的热电优值： $ZT=\frac{S^{2}T\sigma}{\kappa}$

5.傅里叶定律

傅里叶定律(热量输运方程)： $\frac{dQ}{dt}=\kappa \cdot A\frac{\Delta T}{l}$ 其中， $\frac{dQ}{dt}$ 为热流， $\kappa$ 为热导率， $\frac{\Delta T}{l}$ 为温度梯度。描述了热量在材料内部的输运规律
计算出单位截面积的热流和温度梯度，就可以求得材料热导率
1.按照特定的方法对材料某一长度方向施加热流，计算热功率
2.经过一定时间后材料内部温度达到稳态，统计温度梯度
3.结合截面积和长度参数，计算热导率

2.lammps计算Ar热导率

模拟对象：FCC密排的Ar结构
模拟过程：使用傅里叶定律计算Ar的热导率，引入热流在Ar内部产生稳定的温度梯度

lammps代码：

units   lj
atom_style   atomic
lattice   fcc 0.6
region   box block 0 10 0 10 0 20
create_box   1 box
create_atoms   1 box
mass   1 1.0
velocity   all create 1.35 46423
pair_style   lj/cut 2.5
pair_coeff   1 1 1.0 1.0
neighbor   0.3 bin
# 定义0-1为热区，10-11为冷区 
region   hot block INF INF INF INF 0 1                       
region   cold block INF INF INF INF 10 11
# 计算冷热两区的温度
compute   Thot all temp/region hot
compute   Tcold all temp/region cold
# 第一次跑平衡
fix   1 all nvt temp 1.35 1.35 0.5
thermo   100
run   1000
velocity   all scale 1.35
unfix   1
# 第二次跑平衡
fix   1 all nve
fix   hot all heat 1 100.0 region hot                      # 常数热流交换法，对热区注入热量 ，每1时间步向hot区注入100能量单位的热量
fix   cold all heat 1 -100.0 region cold                       # 同时冷区抽出能量，每1时间步从cold区抽出100能量单位的热量
thermo_style   custom step temp c_Thot c_Tcold
thermo   1000
run   10000compute   ke all ke/atom             
variable   temp atom c_ke/1.5                 # 定义温度变量，方便计算各层温度
compute   layers all chunk/atom bin/1d z lower 0.05 units reduced                           # 沿Z方向分20层，选低起始点，units reduced 将z方向无量纲看作1，每层厚度0.05
fix   2 all ave/chunk 10 100 1000 layers v_temp file profile.heat         # 计算各层的温度并输出到文件
variable   tdiff equal f_2[11][3]-f_2[1][3]                   # 计算第11层和第1层的温差
# 计算冷/热区温差，统计平均并输出 
fix   ave all ave/time 1 1 1000 v_tdiff ave running start 13000           # 从13000步开始，对上述温差进行时间平均
thermo_style   custom step temp c_Thot c_Tcold v_tdiff f_ave
run   20000

fix ehex的原理、用法和fix heat一致，相比于fix heat，改进了算法，可以使用非对称的能量交换，改善了由于注入/抽出能量带来的能量漂移问题

Muller-Plathe(速度交换)法：fix thermal/conductivity

units   lj
atom_style   atomic
lattice   fcc 0.6
region   box block 0 10 0 10 0 20
create_box   1 box
create_atoms   1 box
mass   1 1.0
velocity   all create 1.35 46423
pair_style   lj/cut 2.5
pair_coeff   1 1 1.0 1.0
neighbor   0.3 bin
# 第一次跑平衡
fix   1 all nvt temp 1.35 1.35 0.5
thermo   100
run   1000
velocity   all scale 1.35
unfix   1
# 第二次跑平衡
compute   ke all ke/atom             
variable   temp atom c_ke/1.5                 # 定义温度变量，方便计算各层温度
fix   1 all nve
compute   layers all chunk/atom bin/1d z lower 0.05 units reduced                           # 沿Z方向分20层，选低起始点，units reduced 将z方向无量纲看作1，每层厚度0.05
fix   2 all ave/chunk 10 100 1000 layers v_temp file profile.mp       # 计算各层的温度并输出到文件
fix   3 all thermal/conductivity 10 z 20 swap 2                 # z方向分为20层，每10步，第1层原子速度最大的两个原子与第11层原子速度最小的两个交换速度，该命令计算交换的能量值
variable   tdiff equal f_2[11][3]-f_2[1][3]                   # 计算第11层和第1层的温差
thermo_style   custom step temp epair etotal f_3 v_tdiff
thermo   1000
run   20000# 计算热导率
fix   ave all ave/time 1 1 1000 v_tdiff ave running
thermo_style   custom step temp epair etotal f_3 v_tdiff f_ave
run   20000

3.lammps模拟SiGe热电材料

SiGe纳米线建模：
构建半径为1 nm，长度周期为 10 nm的圆柱形 SiGe(元素比1：1)纳米线

# 初始化
units   metal
dimension   3
boundary   f f p
neighbor   2.0 bin
neigh_modify   every 10 delay 0 check yes
atoms_style   atomic
lattice   diamond 5.43
region   box block -10 10 -10 10 0 20
create_box   2 box
region   SiGe cylinder z 0 0 2 0 20 units lattice
create_atoms   1 region SiGe
set   type 1 type/ratio 2 0.5 13658
mass   1 14
mass   2 32
write_data   SiGe.data

热传递过程模拟：

# 初始化
units   metal
timestep   0.001
atom_style   atomic
read_data   SiGe.data
pair_style   tersoff
pair_coeff   * * SiCGe.tersoff Si(D) Ge
mass   1 14
mass   2 32
minimize   1e-4 1e-6 100 1000
# 在100K下热平衡
dump   1 all atom 100 dump.atom
fix   1 all nvt temp 100 100 0.5
run   20000
unfix   1
# MP方法产生温度梯度
compute   ke all ke/atom
variable   KB equal 8.625e-5
variable   temp atom c_ke/1.5/${KB}
fix   1 all nve
compute   layers all chunk/atom bin/1d z lower 0.05 units reduced
fix   2 all ave/chunk 10 100 1000 layers v_temp file profile.tmp
fix   3 all thermal/conductivity 50 z 20
variable   tdiff equal f_2[11][3]-f_2[1][3]               # 第1层、11层温度差
thermo_style   custom step temp epair etotal f_3 v_tdiff
# 计算热导率
fix   ave all ave/time 1 1 1000  v_tdiff ave running
variable   dQdt equal f_3
variable   DeltaTemp equal f_ave
fix   extra all print 1000 "${dQdt} ${DeltaTemp}" file thermal_conductivity.dat screen no
thermo_style   custom step temp epair etotal f_3 v_tdiff f_ave
run   100000

注意事项：
1.利用傅里叶定律计算热导率的关键和难点是要在系统中产生稳定、线性的能量梯度分布
2.无法获得温度梯度时，需要调整系综、热流交换速率和时间步等参数
3.温度梯度不宜太大，一般需要保持高温区和低温区的温差在100K以内，平均值在给定的系综温度附近，不宜有太大漂移

4.平衡分子动力学(EMD)

Fourier 定律计算热导率：
1.模拟过程中对象各部分能量处于非平衡状态，常数热流交换(与外界有人为的能量交换)、MP方法(各组原子之间有人为的能量交换)
2.经过一段时间后，扩散效率与人为交换能量相当时，体系各组的温度趋向稳定
3.基于Fourier定律计算热导率的方法通常被称为非平衡分子动力学模拟(NEMD)
注：各向异性材料的热导率是一个张量，NEMD一次模拟只能获得一个方向的热导率，适合低维材料热导率的模拟
EMD 特点：
1.不需要引入温度梯度，可同时计算其他性质(粘度等)，效率高
2.尺寸效应相对小，理论上可计算无限大体系
3.要求体系必须达到平衡，弛豫时间较长
EMD 实现方式：
1.在微正则系综(NVE)下进行长时间弛豫
2.记录体系热通量，直到自相关积分收敛
3.平衡后，通过Green-Kubo公式计算热导率
基于Green-Kubo线性响应理论，材料的热导率与平衡态下涨落-耗散的热流相关，可以通过平衡态的热流自相关函数来计算热导率
EMD 核心公式： $\kappa_{\mu v}(\tau_{m})=\frac{1}{\Omega_{k_{B}}T^{2}}\int_{0}^{\tau_{m}}\left \langle \bar{J}_{\mu}(\tau)J_{v}(0)\right \rangle d\tau$ 其中， $\Omega_{k_{B}}$ 为玻尔兹曼常数，T 为温度， $\tau_{m}$ 为平衡时间，J 为热通量
离散化时间步求平均： $\kappa_{\mu v}(\tau_{M})=\frac{\Delta t}{\Omega_{k_{B}}T^{2}}\sum_{m=1}^{M}(N-m)^{-1}\sum_{m=1}^{N-m}J_{\mu}(m+n)J_{v}(n)$
EMD 计算热导率，核心是热通量信息的统计，而对于两体相互作用，热通量可简化为： $J=\frac{1}{V}[\sum_{i}e_{i}v_{i}-\sum_{i}S_{i}v_{i}]=\frac{1}{V}[\sum_{i}e_{i}v_{i}+\sum_{i<j}(F_{ij}\cdot v_{j})r_{ij}]$ $J=\frac{1}{V}[\sum_{i}e_{i}v_{i}+\frac{1}{2}\sum_{i<j}(F_{ij}\cdot(v_{i}+v_{j}))r_{ij}]$
其中， $e_{i}$ 为每原子能量(动能和势能)， $S_{i}$ 为每原子压力张量， $v_{i}$ 为速度矢量
lammps可通过内置命令轻松实现能量计算，并通过能量信息计算热通量： $\kappa = \frac{V}{k_{B}T^{2}}\int_{0}^{\infty }\left \langle J_{x}(0)J_{x}(t) \right \rangle dt=\frac{V}{3k_{B}T^{2}}\int_{0}^{\infty}\left \langle J(0)\cdot J(t)\right \rangle dt$ 最后使用简化后的Green-Kubo方程计算热导率

EMD 计算热导率关键命令：

# 计算每个原子的动能、势能以及应力
compute   myStress all stress/atom NULL virial
compute   flux all heat/flux myKE myPE myStress
# 对计算的热流分量的自相关函数进行采样平均
fix   1 all ave/correlate 100 10 100 c_flux[1] c_flux[2] c_flux[3]& type auto file profile.heatflux ave running