1.无序因子
无序度散射因子 \Gamma 是以 Slack 和 Abeles 提出的模型为基础得到的：
\Gamma=\Gamma_{\mathrm{M}}+\Gamma_{\mathrm{S}}
\Gamma_{\mathrm{M}} 为质量波动散射因子，\Gamma_{\mathrm{S}} 为应力场波动散射因子
\Gamma_{\mathrm{M}}=\frac{\sum\limits_{i=1}^{n} c_{i}\left(\frac{\overline{M_{i}}}{\overline{\overline{M}}}\right)^{2} \Gamma_{\mathrm{M}}^{i}}{\sum\limits_{i=1}^{n} c_{i}}
\Gamma_{\mathrm{S}}=\frac{\sum\limits_{i=1}^{n} c_{i}\left(\frac{\overline{M_{i}}}{\overline{M}}\right)^{2} \Gamma_{\mathrm{S}}^{i}}{\sum\limits_{i=1}^{n} c_{i}}
其中, c_{i} 表示各晶格位置的原子总数 (例如对于 SnTe 来说，共两种晶格位置而Sn位和Te位的 c 均为1 ), 该化合物的平均原子质量为：
\overline{\overline{M}}=\frac{\sum\limits_{i=1}^{n} c_{i} \overline{M_{i}}}{\sum\limits_{i=1}^{n} c_{i}}
一般来说，经常会有几种不同的原子占据同一晶格位置，假设第 i 个晶格位置上的第 k 个原子的质量为 M_{i}^{k} , 半径为 r_{i}^{k} , 占比为 f_{i}^{k} 。则第 i 个晶格位置上的平均质量、平均原子半径，以及质量场和应力场波动散射因子为：
\overline{M_{i}}=\sum_{k} f_{i}^{k} M_{i}^{k}
\overline{r_{i}}=\sum_{k} f_{i}^{k} r_{i}^{k}
\Gamma_{\mathrm{M}}^i=\sum_k f_i^k\left(1-\frac{M_i^k}{\overline{M_i}}\right)^{2}
\Gamma_{\mathrm{S}}^k=\sum_k f_i^k \varepsilon\left(1-\frac{r_i^k}{\overline{r_i}}\right)^{2}
\varepsilon 为可调弹性因子（phenomenological adiust elastic properties factor）
\varepsilon=\frac{2}{9}\left[(G+6.4 \gamma) \frac{1+v_{\mathrm{p}}}{1-v_{\mathrm{p}}}\right]^{2}
v_{\mathrm{p}}=\frac{1-2\left(v_{\mathrm{t}} / v_{\mathrm{l}}\right)^{2}}{2-2\left(v_{\mathrm{t}} / v_{\mathrm{l}}\right)^{2}}
\gamma=\frac{3 \alpha_{1} B}{C_{\mathrm{p}} \rho}
E=\frac{\rho v_{\mathrm{t}}^{2}\left(3 v_{1}^{2}-4 v_{\mathrm{t}}^{2}\right)}{v_{1}^{2}-v_{\mathrm{t}}^{2}}
B=\frac{E}{3\left(1-2 v_{\mathrm{p}}\right)}
G 是体模量相对变化和键长的比例，对与Ⅳ-Ⅵ族化合物，G 通常为3，v_{\mathrm{p}} 为泊松比，约为0.23，平均横波声速 v_{\mathrm{t}} 和平均纵波声速 v_{\mathrm{l}}, \gamma为格林艾森常数(\mathrm{Gr\ddot{u}neisen}), \alpha_{\mathrm{l}} 为线膨胀系数，约为 2.07*10^{-5} K^{-1}，B 为体模量，C_{\mathrm{p}} 为比热，\rho为密度，E 为杨氏模量。
2.合金模型计算晶格热导率
参考自文献1和3
U^{2} =\left(\frac{\pi^{2} \Theta_{\mathrm{D}} \Omega}{h v^{2}}\right) \kappa_{\mathrm{L}}^{\text {pure }} \Gamma
\Theta_{\mathrm{D}}=\frac{h}{k_\mathrm{B}}\left[\frac{3n}{4\pi}\left(\frac{N_A\rho}{M}\right)\right]^{1/3}v_{\mathrm{s}}
v=\left[\frac{1}{3}\left(\frac{2}{v_\mathrm{t}^3}+\frac{1}{v_\mathrm{l}^3}\right)\right]^{-\frac{1}{3}}
\kappa_{\mathrm{L}}^{\text {cal }} =U^{-1} \tan ^{-1} (U) \kappa_{\mathrm{L}}^{\text {pure}}
其中 U 是无序参数，\Theta_\mathrm{D} 是由声速计算得到的德拜温度，\Omega 是平均原子体积，V 是原胞的体积，N_\mathrm{A} 是阿伏伽德罗常数，M 是原胞内原子的平均相对原子质量和，h 是普朗克常数，v_{\rm{t}} 是横波声速，v_{\rm{l}} 是纵波声速，v 是平均声速，n 是原胞中的原子数目；\rho 是样品密度，\kappa_{\mathrm{L}}^{\text{pure}} 是基体的晶格热导率，\kappa_{\mathrm{L}}^{\text{cal}} 则是计算得到的理论晶格热导率。
2.1 计算示例
数据来文献1
2.1.1 数据以及代码文件
data.txt
Samples 横波声速（m/s） 纵波声速（m/s） 平均声速（m/s） 杨氏模量（GPa） 体模量（GPa） 格林艾森常数-gamma 可调参数-varepsilon 晶格热导率（W/m/K）
0.00 1871 3068 2066 52.5 32.4 1.54 93.9 1.44811
0.02 1857 3341 2068 52.5 32.4 1.53 92.8 1.31547
0.04 1853 3343 2064 51.9 32.0 1.52 91.7 1.35364
0.06 1837 3091 2034 50.3 31.0 1.51 90.6 1.05717
0.08 1770 2846 1950 48.5 29.9 1.50 89.6 0.89148
0.10 1828 3015 2020 48.9 30.2 1.48 88.5 1.06229Slack_Abeles.m
%因为Te为并未掺杂所以无论质量场波动散射因子还是应力场散射因子都为0，计算时可不考虑
%Sn0.85-xSb0.15MgxTe
clc,clear,close all;
%实验数据导入
data_import;
Kb = 1.38066*10^-23; %玻尔兹曼常数, J/K 1J=1V*A*s 1W = 1J/s
h = 6.62608*10^-34; %普朗克常数, J*s
Na = 6.02214e23; %阿伏伽德罗常数, mol^-1
M0 = 1.661*10^(-27); %C原子质量的1/12, Kg
kappaL_pure = data{9}(1); %Sn0.85Sb0.15Te在300K下的晶格热导率，W/m/K
%相对原子质量
M_Sn = 118.71;
M_Sb = 121.75;
M_Mg = 24.305;
M_Te = 127.60;
%原子（离子）半径，单位：埃，计算结果对离子半径比较敏感
r_Sn = 0.93;
r_Sb = 0.92;
r_Mg = 0.86;
r_Te = 0.94;
Gamma_M = zeros(6,1);
Gamma_S = zeros(6,1);
Gamma = zeros(6,1);
kappaL_cal = zeros(6,1);
rho = zeros(6,1);
V = 252.8*10^(-30); % Å^3，晶胞体积
for i=1:6
%Sn位晶格的Sn原子占比
f_Sn_Sn = (0.85-data{1}(i));
%Sn位晶格的Sb原子占比
f_Sn_Sb = 0.15;
%Sn位晶格的Mg原子占比
f_Sn_Mg = data{1}(i);
%Sn位的平均原子质量
Mave_Sn = f_Sn_Sn*M_Sn+f_Sn_Sb*M_Sb+f_Sn_Mg*M_Mg;
%Sn位的平均原子（离子）半径
rave_Sn = f_Sn_Sn*r_Sn+f_Sn_Sb*r_Sb+f_Sn_Mg*r_Mg;
%Te位晶格的Te原子占比
f_Te_Te = 1;
%Te位的平均原子质量
Mave_Te = f_Te_Te*M_Te;
%Te位的平均原子（离子）半径
rave_Te = f_Te_Te*r_Te;
%Sn位的晶格原子总数
c_Sn = 1;
%Te位的晶格原子总数
c_Te = 1;
%化合物的平均原子质量
Mave_formula = (c_Sn*Mave_Sn+c_Te*Mave_Te)/(c_Sn+c_Te);
%Sn位的质量场波动散射因子
Gamma_M_Sn = f_Sn_Sn*(1-M_Sn/Mave_Sn)^2+f_Sn_Sb*(1-M_Sb/Mave_Sn)^2+f_Sn_Mg*(1-M_Mg/Mave_Sn)^2;
%Te位的质量场波动散射因子
Gamma_M_Te = f_Te_Te*(1-M_Te/Mave_Te)^2;
%总质量场波动散射因子
Gamma_M(i) = (c_Sn*(Mave_Sn/Mave_formula)^2*Gamma_M_Sn+c_Te*(Mave_Te/Mave_formula)^2*Gamma_M_Te)/(c_Sn+c_Te);
%Sn位的应力场波动散射因子
varepsilon = data{8}(i);
Gamma_S_Sn = varepsilon*(f_Sn_Sn*(1-r_Sn/rave_Sn)^2+f_Sn_Sb*(1-r_Sb/rave_Sn)^2+f_Sn_Mg*(1-r_Mg/rave_Sn)^2);
%Te位的应力场波动散射因子
Gamma_S_Te = f_Te_Te*varepsilon*(1-r_Te/rave_Te)^2;
%总应力场波动散射因子
Gamma_S(i) = (c_Sn*(Mave_Sn/Mave_formula)^2*Gamma_S_Sn+c_Te*(Mave_Te/Mave_formula)^2*Gamma_S_Te)/(c_Sn+c_Te);
%总波动散射因子
Gamma(i) = Gamma_M(i) + Gamma_S(i);
n = 2;
%原胞内原子的平均相对原子质量
M = Mave_Sn + Mave_Te;
M = M/n*10^(-3); %g/mol化为Kg/mol 摩尔质量
%平均原子质量 M/n
vt = data{2}(i);
vl = data{3}(i);
% v = data{4}(i);
%平均声速
v = ((1/3)*(2/vt^3+1/vl^3))^(-1/3);
%杨氏模量 1Gpa = 10^9 Pa = 10^9 N/m^2 = 10^9 Kg*m/s
E = data{5}(i)*10^9;
rho(i) = E*(vl^2-vt^2)/(vt^2*(3*vl^2-4*vt^2));
%德拜温度
Theta_D = (h/Kb)*((3*n)/(4*pi)*(Na*rho(i)/M))^(1/3)*v;
%平均原子体积，用原胞（晶胞）体积/原胞（晶胞）内原子数，此晶胞内有8个原子
Omega = V/(n*4);
%无序参数U,U2=U^2
U2 = (pi^2*Theta_D*Omega/(h*v^2))*kappaL_pure*Gamma(i);
kappaL_cal(i) = (sqrt(U2))^-1*atan(sqrt(U2))*kappaL_pure;
end
figure(1)
plot(data{1}, Gamma_M,'r-*', data{1}, Gamma_S, 'g-*', data{1}, Gamma,'b-*');
legend({'\bf{\Gamma_M}','\bf{\Gamma_S}','\bf{\Gamma}'},'Location','northwest');
xlabel('x')
figure(2)
plot(data{1}, kappaL_cal,'r--');
hold on
scatter(data{1},data{9},50*ones(1,6),'r','filled');
xlabel('x')
ylabel('\bf\kappa_L (W\cdotm^{-1}\cdotK^{-1})')
box on
xlim([-0.01,0.11])
ylim([0.4,2.0]);data_import.m
%% 实验数据的导入，这是我从文献中提取的
clc,clear
%打开文件句柄
file_data = fopen('data.txt','r');
%以空格分隔，%'EmptyValue', 0，将重复的分隔符视为一个分隔符，从第2行开始读取，忽略第一列
data = textscan(file_data , '%f %f %f %f %f %f %f %f %f', 'Delimiter', ' ', 'MultipleDelimsAsOne',1, 'headerlines',1);
%关闭文件句柄
fclose(file_data );
% matrix_data= cell2mat(data);2.1.2 计算结果
不出意外和文献计算结果还是有差异，计算结果对离子半径较为敏感。
3.参考文献
[1]: 傅铁铮. P型SnTe基热电材料的电声输运及性能优化[D].浙江大学,2019.
评论区