lammps应用于热电材料
1.热传导理论
1.热导率
- 热传递机制随介质材料相的不同而改变：固体(热传导)、液体(热对流)、气体(对流和辐射)
- 热传导：热量从固体材料高温向低温部分转移的过程
- 热力学第二定律(克劳修斯表述)：“不可能把热量从低温物体传向高温物体而不引起其他变化。”
在自然条件下热量只能从高温物体向低温物体转移，而不能由低温物体自动向高温物体转移，这个转变过程是自发、不可逆的 - 热导率：又称导热系数(coefficient of thermal conductivity)，物质导热能力的量度
若物体的截面积为 1m2，相距 1m 的两个平面的温度差为 1K，则1秒内流过界面的热量就规定为该物质的热导率，标准单位为 W/(m·K) - T=0K，原子固定在平衡位置处，保持能量最低的状态
- 当T>0K时，原子在格点原子势场中热振动
2.晶格振动
- 材料各种热学性能的物理本质，均与其晶格热振动有关
- 绝缘材料中的导热现象主要由原子(晶格)的振动所引起
谐振子量子理论： H = 1 2 m x ˙ 2 + 1 2 k x 2 = 1 2 q ˙ 2 + 1 2 ω 2 q 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} H=21mx˙2+21kx2=21q˙2+21ω2q2 量子力学中，谐振子有分立的本征能量： ε n = ( n + 1 2 ) ℏ ω , n = 0 , 1 , . . . , ∞ \varepsilon_{n} =(n+\frac{1}{2})\hbar \omega,n=0,1,...,\infty εn=(n+21)ℏω,n=0,1,...,∞ - 材料热导率的差异来源于晶格原子振动对应的传递能力，声子(晶格振动的能量量子)之间的碰撞导致热量的传递
声子的能量： ℏ ω ， ω = k m \hbar\omega，\omega=\sqrt{\frac{k}{m}} ℏω，ω=mk 原子/分子间的键越强，导热能力越强 - 声子对热导的贡献：
1.长波声子具有更长的波长，易发生倒逆过程，因此对热传导性质影响更大
2.通常声学声子对热传导的贡献远大于光学声子 - 热导率公式：
κ
=
1
3
C
v
ˉ
λ
ˉ
\kappa =\frac{1}{3}C\bar{v}\bar{\lambda}
κ=31Cvˉλˉ 1.与固体的比热 C 成正比，因为比热越大，每个固体分子可带的热量也越多
2.与声子的平均速率 v ˉ \bar{v} vˉ 成正比，因为声子速度越快，热也传送得越快
3.与声子的平均自由程 λ ˉ \bar{\lambda} λˉ (与声子碰撞频率相关)成正比， λ ˉ \bar{\lambda} λˉ 越大，声子携带的热能传得越远
其中，声子的平均速度与温度变化的关联性不大，多数情况可以作为常数处理。而平均自由程比较难计算，与温度和晶格声子碰撞的方式、频率等因素有关
3.晶体热容
C V = ( ∂ U ∂ T ) V C_{V}=(\frac{\partial U}{\partial T})_{V} CV=(∂T∂U)V C V = C V l + C V e C_{V}=C_{V}^{l}+C_{V}^{e} CV=CVl+CVe 其中， C V l C_{V}^{l} CVl 为晶格振动热容， C V e C_{V}^{e} CVe 为晶体电子热容
- 固体热容主要有两部分：一个是来源于晶格振动(声子)，称为晶格热容；一个来源于电子的热运动，称为电子热容。除非在很低温，电子热运动的贡献往往很小。通常情况下， C V e < < C V l C_{V}^{e}<<C_{V}^{l} CVe<<CVl
- 实验规律：
1.在高温时，晶体的热容为： 3 N k B = 3 v R 3Nk_{B}=3vR 3NkB=3vR N 为晶体中原子的个数，kB = 1.38 × 10-23 J·K-1 为玻尔兹曼常量；v 为晶体中原子摩尔数，R = 8.31 J/K mol为普适气体常数
2.在低温时，绝缘体热容按 T 3 T^{3} T3 趋于零，导体热容按 T 趋于零 - 晶格振动的能量是量子化的，频率为 ω \omega ω 的声子能量为： E n = ( n + 1 2 ) ℏ ω E_{n}=(n+\frac{1}{2})\hbar\omega En=(n+21)ℏω 1 2 ℏ ω \frac{1}{2}\hbar\omega 21ℏω 代表零点振动能，对热容没有贡献 E n = n ℏ ω E_{n}=n\hbar\omega En=nℏω n n n 是频率为 ω \omega ω 的谐振子的平均声子数，据玻色统计理论： n ( ω ) = 1 e ℏ ω / k B T − 1 n(\omega)=\frac{1}{e^{\hbar\omega/k_{B}T}-1} n(ω)=eℏω/kBT−11 温度为 T T T 时，频率为 ω \omega ω 的振动的能量： E ˉ ( ω ) = ℏ ω e ℏ ω / k B T − 1 \bar{E}(\omega)=\frac{\hbar\omega}{e^{\hbar\omega/k_{B}T}-1} Eˉ(ω)=eℏω/kBT−1ℏω
- 爱因斯坦模型：
晶体由 N 个原子组成，每个原子有3个自由度，共有 3N 个分立的振动频率，晶体内能： U = ∑ i = 1 3 N E ˉ ( ω i ) = ∑ i = 1 3 N ℏ ω i e ℏ ω i / k B T − 1 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=i=1∑3NEˉ(ωi)=i=1∑3Neℏωi/kBT−1ℏωi 爱因斯坦模型：假设所有原子以同样的频率振动 U = ∑ i = 1 3 N E ˉ ( ω i ) = 3 N ℏ ω e ℏ ω / k B T − 1 U=\sum_{i=1}^{3N}\bar{E}(\omega_{i})=\frac{3N\hbar\omega}{e^{\hbar\omega/k_{B}T}-1} U=i=1∑3NEˉ(ωi)=eℏω/kBT−13Nℏω
高温下，比热为 3 N K b 3NK_{b} 3NKb，与经典结果一致；低温下与实验偏差较大 - 德拜模型：
若频率分布可用一个积分函数表示： ∫ 0 ω m g ( ω ) d ω = 3 N \int_{0}^{\omega_{m}}g(\omega)d\omega=3N ∫0ωmg(ω)dω=3N
g ( ω ) d ω g(\omega)d\omega g(ω)dω 表示在频率范围 ω ∼ ω + d ω \omega \sim\omega+d\omega ω∼ω+dω 可取的频率数， ω m \omega_{m} ωm为最大的频率数， q q q 和 ω \omega ω 为准连续 U = ∑ i = 1 3 N ℏ ω i e ℏ ω i / k B T − 1 = ∫ 0 ω m ℏ ω e ℏ ω / k b T − 1 g ( ω ) d ω 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 U=i=1∑3Neℏωi/kBT−1ℏωi=∫0ωmeℏω/kbT−1ℏωg(ω)dω 热容： C V = ( ∂ U ∂ T ) V = ∫ 0 ω m k B ( ℏ ω k B T ) 2 e ℏ ω / k B T g ( ω ) d ω ( e ℏ ω / k B T − 1 ) 2 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}} CV=(∂T∂U)V=∫0ωmkB(kBTℏω)2(eℏω/kBT−1)2eℏω/kBTg(ω)dω 变量代换： f D ( θ D T ) = 3 ( T θ ) 3 ∫ 0 θ D T e x x 4 ( e x − 1 ) 2 d x 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 fD(TθD)=3(θT)3∫0TθD(ex−1)2exx4dx 低温下： C v ∼ 12 π 4 R ( T / θ D ) 3 5 C_{v}\sim12\pi^{4}R\frac{(T/\theta_{D})^{3}}{5} Cv∼12π4R5(T/θD)3 - 对于大多数固体来说，Debye模型可以给出相当精确的结果
4.声子平均自由程
κ = 1 3 C v ˉ λ ˉ \kappa=\frac{1}{3}C\bar{v}\bar{\lambda} κ=31Cvˉλˉ 与声子数密度成反比： λ ˉ ∝ 1 n \bar{\lambda}\propto\frac{1}{n} λˉ∝n1 频率为 ω \omega ω 的平均声子数： n ( ω ) = 1 e ℏ ω / k B T − 1 n(\omega)=\frac{1}{e^{\hbar\omega/k_{B}T}-1} n(ω)=eℏω/kBT−11 高温下： λ ˉ ∝ T − 1 ， C 为常数， κ ∝ T − 1 \bar{\lambda}\propto T^{-1}，C为常数，\kappa\propto T^{-1} λˉ∝T−1，C为常数，κ∝T−1 低温下： λ ˉ ∝ e α T ， C ∝ T 3 ， κ ∝ e a T \bar{\lambda}\propto e^{\frac{\alpha}{T}}，C\propto T^{3}，\kappa\propto e^{\frac{a}{T}} λˉ∝eTα，C∝T3，κ∝eTa
- 热容随温度增加，声子的平均自由程 λ ˉ \bar{\lambda} λˉ 随温度减小，两者竞争使得 κ \kappa κ 出现极小值
- 低温：声子平均自由程 λ ˉ \bar{\lambda} λˉ 近似为粒径大小， C V C_{V} CV 随 T 3 T^{3} T3 增大， κ \kappa κ 值增大
- 高温：声子的平均自由程 λ ˉ \bar{\lambda} λˉ 随温度升高减小，而声子热容 C C C 趋近于常数， κ \kappa κ 值降低
- 热电材料：
塞贝克(Seebeck)效应：用于热电能量转换，如节能、特殊电源、地热发电
珀尔贴(Peltier)效应：用于热电致冷，如冰箱、空调
评价热电转化效率的热电优值： Z T = S 2 T σ κ ZT=\frac{S^{2}T\sigma}{\kappa} ZT=κS2Tσ
5.傅里叶定律
- 傅里叶定律(热量输运方程)： d Q d t = κ ⋅ A Δ T l \frac{dQ}{dt}=\kappa \cdot A\frac{\Delta T}{l} dtdQ=κ⋅AlΔT 其中， d Q d t \frac{dQ}{dt} dtdQ 为热流， κ \kappa κ 为热导率， Δ T l \frac{\Delta T}{l} lΔT 为温度梯度。描述了热量在材料内部的输运规律
- 计算出单位截面积的热流和温度梯度，就可以求得材料热导率
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 10000 compute 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 核心公式：
κ
μ
v
(
τ
m
)
=
1
Ω
k
B
T
2
∫
0
τ
m
⟨
J
ˉ
μ
(
τ
)
J
v
(
0
)
⟩
d
τ
\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
κμv(τm)=ΩkBT21∫0τm⟨Jˉμ(τ)Jv(0)⟩dτ 其中，
Ω
k
B
\Omega_{k_{B}}
ΩkB 为玻尔兹曼常数，T 为温度，
τ
m
\tau_{m}
τm 为平衡时间，J 为热通量
离散化时间步求平均： κ μ v ( τ M ) = Δ t Ω k B T 2 ∑ m = 1 M ( N − m ) − 1 ∑ m = 1 N − m J μ ( m + n ) J v ( n ) \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) κμv(τM)=ΩkBT2Δtm=1∑M(N−m)−1m=1∑N−mJμ(m+n)Jv(n) - EMD 计算热导率，核心是热通量信息的统计，而对于两体相互作用，热通量可简化为：
J
=
1
V
[
∑
i
e
i
v
i
−
∑
i
S
i
v
i
]
=
1
V
[
∑
i
e
i
v
i
+
∑
i
<
j
(
F
i
j
⋅
v
j
)
r
i
j
]
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=V1[i∑eivi−i∑Sivi]=V1[i∑eivi+i<j∑(Fij⋅vj)rij]
J
=
1
V
[
∑
i
e
i
v
i
+
1
2
∑
i
<
j
(
F
i
j
⋅
(
v
i
+
v
j
)
)
r
i
j
]
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}]
J=V1[i∑eivi+21i<j∑(Fij⋅(vi+vj))rij]
其中， e i e_{i} ei 为每原子能量(动能和势能)， S i S_{i} Si 为每原子压力张量， v i v_{i} vi 为速度矢量 - lammps可通过内置命令轻松实现能量计算，并通过能量信息计算热通量： κ = V k B T 2 ∫ 0 ∞ ⟨ J x ( 0 ) J x ( t ) ⟩ d t = V 3 k B T 2 ∫ 0 ∞ ⟨ J ( 0 ) ⋅ J ( t ) ⟩ d t \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 κ=kBT2V∫0∞⟨Jx(0)Jx(t)⟩dt=3kBT2V∫0∞⟨J(0)⋅J(t)⟩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
原文地址：https://blog.csdn.net/qq_44963682/article/details/145103821
免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！