热电系列-德拜-卡拉威(Debye-Callaway)模型-琉璃红梅

说明：
1.代入计算时要特别注意量纲
2.由于代码的编写时间不同，导致同一物理量在不同代码文件中所用的符号不一致
3.如果lsqcurvefit函数不能运行，可能是没有安装相关工具包，可以根据提示安装
4.拟合计算的结果特别依赖初值，以及边界条件，甚至有“人为调参”的痕迹
5.这一模型应该是在低温下的，但是文献中大都选取室温到高温的数据拟合
6.如果有问题欢迎讨论，欢迎指出错误
致谢：
元素周期表的数据来自该项目源代码：TFflow，在此对作者表示感谢！

1.最小晶格热导率

1.1 基于Allen和Feldman的扩散图像

1.1.1 计算公式

\kappa_{\mathrm{min}}=\frac{n^{1/3}k_{\mathrm{B}}}{\pi}\omega_{\mathrm{avg}}

n为原子数密度, n=6.024\times10^{28}\ \mathrm{atmos}/m^3 取自ICSD的数据
\omega_{avg} = 16\ \mathrm{meV}声子振动模式的平均频率, E = \hbar \omega
通过经验发现，\omega_{avg} 通过最大的Debye频率与声速相联系, 可以对方程1做如下替换，简化为：

\kappa_{\mathrm{min}}\approx0.77n^{2/3}k_{\mathrm{B}}v_\mathrm{s}

其中

v_{\mathrm{s}}=\frac{1}{2}\left(v_\mathrm{L}+2v_\mathrm{T}\right)\approx1726\ \mathrm{m/s}

1.2 基于Cahill-Pohl方程

1.2.1 计算公式

\kappa_{\mathrm{min}}=\left(\frac{\pi}{6}\right)^{1/3}k_{\mathrm{B}}n^{2/3}\sum\limits_{i}v_{i}\left(\frac{T}{\Theta_{i}}\right)^2\int_0^{\Theta_i/T}\frac{x^3\mathrm{e}^x}{\left(\mathrm{e}^x-1\right)^2}{\mathrm{d}}x

v_i 为声速, \Theta_{i}=v_i\left(\hbar/k_\mathrm{B}\right)\left(6\pi^2n\right)^{1/3} 为截止频率, n为原子数密度。

1.2.2 计算代码

%CahillPohl25 
clc,clear;
Kb = 1.38066*10^-23;  %J/K 1J=1V*A*s 1W = 1J/s
h = 6.62608*10^-34; %J*s
hr = h/(2*pi);
T = 450; %K
n  = 6.024*10^28; %粒子数密度，单位体积内的粒子数，atoms/m^3
v = [2640,1259,1259]; %声速，m/s
D = v*(hr/Kb)*(6*pi^2*n)^(1/3)/T; % 
sum = 0;
syms x;
f = @(x)x.^3.*exp(x)./(exp(x)-1).^2;
for i = 1:length(v)
    Fn(i) = integral(@(x)f(x),0,D(i));  %求解积分
    sum = sum + v(i)*(1/D(i))^2*Fn(i);
end
kmin1 = (pi/6)^(1/3)*Kb*n^(2/3)*sum  %根据Cahill模型求得的晶格热导率的最小值
w = 16; %meV
w = w*1.60219*10^(-19)*10^(-3);%将单位换成J
w = w/hr;%将单位换成s^(-1)
kmin2 = n^(1/3)*Kb*w/pi %根据Allen和Feldman模型求得的晶格热导率的最小值,w为平均振动频率，由声子态密度求出
vs = 1726; %m/s,vs = 1/3*(vl+2vt),vl为纵波声速，vt为横波声速
kmin3 = 0.77*n^(2/3)*Kb*vs %将上式简化

% kmin1 = 0.4357；kmin2 = 0.4188；kmin3 = 0.2820

1.3 参考文献

[1]: Nunna, R., et al. (2017). “Ultrahigh thermoelectric performance in Cu₂Se-based hybrid materials with highly dispersed molecular CNTs.” Energy & Environmental Science 10(9): 1928-1935.

[2]: Zhao, K., et al. (2019). “Are Cu₂Te-Based Compounds Excellent Thermoelectric Materials?” Advanced Materials 31(49): 1903480.

2.晶格热导率

2.1 基于德拜-卡拉威模型计算晶格热导率

2.1.1 计算公式

\kappa_\mathrm{L}=\int_0^{\Theta_\mathrm{D}/T}\kappa_\mathrm{s}(x)\mathrm{d}x

其中 \kappa_\mathrm{s}(x) 为频谱晶格热导率，一般用来说明各种散射机制对不同频率声子的影响程度，x=\hbar\omega/(k_\mathrm{B}T)，\hbar 为约化普朗克常数，\omega 为声子频率，\Theta_\mathrm{D}为德拜温度。

\kappa_\mathrm{s}(x)=\frac{k_\mathrm{B}}{2\pi^2v}\left(\frac{k_\mathrm{B}}{\hbar}\right)^3T^3\frac{x^4\exp(x)}{\tau(x)^{-1}\left[\exp(x)-1\right]^2}

其中 \tau(x)为总的弛豫时间，由热电材料中存在的各种散射机制共同决定，根据马希森（Matthiessen）规则，总的弛豫时间等于各种散射机制的“和”：

{\tau}^{-1}={\tau_\mathrm{U}}^{-1}+{\tau_\mathrm{N}}^{-1}+{\tau_\mathrm{PD}}^{-1}+{\tau_\mathrm{B}}^{-1}+{\tau_\mathrm{D}}^{-1}+{\tau_\mathrm{NP}}^{-1}+\cdots

{\tau_\mathrm{D}}^{-1} = {\tau_\mathrm{DS}}^{-1}+{\tau_\mathrm{DC}}^{-1}

其中 \tau_\mathrm{U}、\tau_\mathrm{N}、\tau_\mathrm{PD}、\tau_\mathrm{B}、\tau_\mathrm{D}、\tau_\mathrm{NP} 分别为声子U过程散射、声子N过程散射、点缺陷散射、晶界散射、位错散射（位错应力散射、位错核散射）、纳米引入物（球形）散射导致的声子弛豫时间，它们的表达式如下：

{\tau_\mathrm{U}}^{-1}=\frac{\hbar\gamma^2\omega^2T}{\overline{M}v^2\Theta_\mathrm{D}}\exp\left(-\frac{\Theta_\mathrm{D}}{3T}\right)

其中 \gamma 和 \overline{M} 依次为格林艾森常数与晶体的平均原子质量，

{\tau_\mathrm{N}}^{-1}=\beta{\tau_\mathrm{U}}^{-1}

其中 \beta 为U过程散射与N过程散射之间的比例常数，

{\tau_\mathrm{PD}}^{-1}=\frac{\overline{V}(\Gamma_\mathrm{M}+\Gamma_\mathrm{S})}{4{\pi}v^3}\omega^4=A\omega^4

\Gamma=\sum_if_i\left(\frac{M_i-\overline{M}}{\overline{M}}\right)

其中 \overline{V}、\Gamma_\mathrm{M}、\Gamma_\mathrm{S} 依次为平均原子体积，质量场波动散射因子、应力场波动散射因子，

{\tau_\mathrm{B}}^{-1}=\frac{v}{L}

其中 L 为平均晶粒尺寸，

{\tau_\mathrm{DS}}^{-1}=0.6{B_\mathrm{D}}^2N_\mathrm{D}(\gamma+\Delta\gamma)^2\omega\left\{\frac{1}{2}+\frac{1}{24}\left(\frac{1-2r_\mathrm{P}}{1-r_\mathrm{P}}\right)^2\left[1+\sqrt{2}\left(\frac{v_\mathrm{l}}{v_\mathrm{t}}\right)^2\right]^2\right\}

{\tau_\mathrm{DC}}^{-1} = N_\mathrm{D}\frac{\overline V^{4/3}\omega^3}{v^2}

{\tau_\mathrm{D}}^{-1} = {\tau_\mathrm{DS}}^{-1}+{\tau_\mathrm{DC}}^{-1}

其中N_\mathrm{B}、N_\mathrm{D}、\Delta\gamma、r_\mathrm{P}、v_\mathrm{l}、v_\mathrm{t} 依次为有效伯格斯矢量、位错密度、位错应力导致的格林艾森常数的改变、泊松比、纵波声速、横波声速。

{\tau_\mathrm{NP}}^{-1}=v\left\{(2{\pi}R)^{-1}+\left[{\pi}R^2\frac{4}{9}{\left(\frac{\Delta_\mathrm{D}}{D}\right)}^2{\left(\frac{\omega R}{v}\right)}^4\right]^{-1}\right\}^{-1}N_\mathrm{P}

其中 R、\Delta_\mathrm{D}、D、N_\mathrm{P} 依次为纳米引入物的平均半径、纳米引入物的分布（数）密度。

散射机制的多样型和共存型意味着，调控晶格热导率的方式是比较广泛的，同时也可以综合引入各种散射机制共同降低晶格热导率。

2.1.2 计算代码

说明：
1.如果修改了 txt 文件名，在 setup.m 文件中也要修改相应路径，但是不要修改 txt 文件的格式。
2.放到统一文件下，并在该文件夹下执行 setup 函数即可，执行完毕会在该文件夹下生成 calData 子文件夹，用于存放计算结果。
3.整个计算均是试图复现文献2

文件夹的层级关系如下图所示

setup.m

%% 读取数据并执行计算
function setup(Tstart, Tend)
% 计算晶格热导率的起始温度：Tstart = 300; 结束温度：Tend = 850
if nargin ~= 2
    Tstart = 300;
    Tend  = 850;
end

% 导入计算材料的化学式/组分
[Elementexp, ratioexp] = loadData('./data/化学式.txt');
composition = containers.Map(Elementexp,ratioexp);

% 导入晶格热导率实验数据
[Texp, klexp] = loadData('./data/晶格热导率.txt');

% 导入拟合晶格热导率需要的参数
[DCParamsKeys, DCParamsValues] = loadData('./data/散射因子.txt');
DCParams = containers.Map(DCParamsKeys,DCParamsValues);

% 将列数据转为行数据传入参数列表进行计算
DebyeCallawayFcn(Texp', klexp', Tstart, Tend, DCParams, composition);
end

DebyeCallawayFcn.m

%% 德拜-卡拉威模型求解晶格热导率
function DebyeCallawayFcn(Texp, klexp, Tstart, Tend, DCParams, composition)
%% 加载参数与数据
% 加载化学周期表，获取相对原子质量
PTABLE = loadPeriodicTable;

% 加载常量参数
ConstParams = loadConstParams;
% 取出用到的常量参数
kB = ConstParams('kB');
h = ConstParams('h');
hbar = ConstParams('hbar');
NA = ConstParams('NA');

% 取出拟合晶格热导率所需参数
T = DCParams('T');
va = DCParams('va');
vl= DCParams('vl');
vt = DCParams('vt');
Vcell = DCParams('Vcell');
Ncell = DCParams('Ncell');
beta = DCParams('beta');
ThetaD = DCParams('ThetaD');
r = DCParams('r');
ND = DCParams('ND')*10^(15);
BD = DCParams('BD')*10^(-10);
l = DCParams('l')*10^(-9);
DeltaDdivD = DCParams('DeltaD/D');
R = DCParams('R')*10^(-9);
NP = DCParams('NP')*10^(20);
gamma = DCParams('gamma');
Gamma = DCParams('Gamma');
omegaD = ThetaD*kB/hbar;
Va = Vcell/Ncell*(10^-10)^3;

% 取出组分化学计量比参数
% 计算平均原子质量
ElementexpKeys = keys(composition);
ratioexpValues = values(composition);
M = 0;
atomCount = 0;
for i = 1:composition.Count
    M = M + ratioexpValues{i}*PTABLE(ElementexpKeys{i});
    atomCount = atomCount + ratioexpValues{i};
end
% 平均相对原子质量，Kg，组分的质量除以原子数目再除以阿伏伽德罗常数
Ma = M / (atomCount*NA)*10^-3;

omegadivomegaD = 0:0.001:1;
% 要去掉第一个点因为当omega=0是被积函数在分母上,分母不能为0,去掉第一个点对最终结果影响不大
omegadivomegaD(1) = [];
omega = omegadivomegaD*omegaD;

%% 散射函数
% U散射(phonon-phonon U)
    function y_U = func_U(T)
        y_U = (hbar*gamma^2*omega.^2*T)/(Ma*va^2*ThetaD)*exp(-ThetaD/(3*T));
    end

% N散射(phonon-phonon N)
    function y_N = func_N(T)
        y_N = beta*(hbar*gamma^2*omega.^2*T)/(Ma*va^2*ThetaD)*exp(-ThetaD/(3*T));
    end

% 点缺陷散射(Point Defects)
    function y_PD = func_PD()
        y_PD = Va*Gamma/(4*pi*va^3)*omega.^4;
    end

% 晶界散射(Boundaries)
    function y_B = func_B()
        y_B = va/l;
    end

% 位错应力场散射(Dislocation Strain fileds)
    function y_DS = func_DS()
        y_DS =  0.6*BD^2*ND*gamma^2*omega*(1/2+1/24*((1-2*r)/(1-r))^2*(1+sqrt(2)*(vl/vt)));
    end

% 位错核散射(Dislocation cores)
    function y_DC = func_DC()
        y_DC = ND*Va^(4/3)/va^2*omega.^3;
    end

% 纳米引入物(Nanoprecipates)
    function y_NP = func_NP()
        y_NP =  va*((2*pi*R)^-1+(pi*R^2*(4/9)*(DeltaDdivD)^2*(omega*R/va).^4).^-1).^-1*NP;
    end

% U+N的弛豫时间的倒数
    function y_U_N = func_U_N(T)
        y_U_N=  func_U(T) + func_N(T);
    end

% U+N+B的弛豫时间的倒数
    function y_U_N_B = func_U_N_B(T)
        y_U_N_B=  func_U(T) + func_N(T) +  func_B();
    end

% U+N+B+PD的弛豫时间的倒数
    function y_U_N_B_PD = func_U_N_B_PD(T)
        y_U_N_B_PD=  func_U(T) + func_N(T) +  func_B() + func_PD();
    end

% U+N+B+PD+DS的弛豫时间的倒数
    function y_U_N_B_PD_DS = func_U_N_B_PD_DS(T)
        y_U_N_B_PD_DS=  func_U(T) + func_N(T) +  func_B() + func_PD() + func_DS() ;
    end

% U+N+B+PD+DS+NP的弛豫时间的倒数
    function y_U_N_B_PD_DS_NP = func_U_N_B_PD_DS_NP(T)
        y_U_N_B_PD_DS_NP=  func_U(T) + func_N(T) +  func_B() + func_PD() + func_DS() +func_NP();
    end

% 总的弛豫时间的倒数
    function y_total = func_total(T)
        y_total =  func_U(T) + func_N(T) + func_PD() + func_B() + func_DS() + func_DC() + func_NP();
    end

%% 频谱晶格热导率ks
    function y_ks = func_ks(T)
        x = hbar*omega/(kB*T);
        y_ks_U_N= 4*pi*kB^4*T^3/(va*h^3).*x.^4.*exp(x)./(func_U_N(T).*(exp(x)-1).^2);
        y_ks_U_N_B= 4*pi*kB^4*T^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B(T).*(exp(x)-1).^2);
        y_ks_U_N_B_PD= 4*pi*kB^4*T^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B_PD(T).*(exp(x)-1).^2);
        y_ks_U_N_B_PD_DS= 4*pi*kB^4*T^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B_PD_DS(T).*(exp(x)-1).^2);
        y_ks_U_N_B_PD_DS_NP= 4*pi*kB^4*T^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B_PD_DS_NP(T).*(exp(x)-1).^2);
        y_ks_total= 4*pi*kB^4*T^3/(va*h^3).*x.^4.*exp(x)./(func_total(T).*(exp(x)-1).^2);
        y_ks = [y_ks_U_N; y_ks_U_N_B; y_ks_U_N_B_PD; y_ks_U_N_B_PD_DS; y_ks_U_N_B_PD_DS_NP; y_ks_total];
        y_ks = y_ks/omegaD*10^14;

        % 曲线的颜色
        strcolor = {'#000000','#7e01f5','#f40708','#1c8b8b','#4a62ff','#cf5e00'};
        % 曲线的线型
        strline = {'-','--','-','--','-','--'};
        % 在坐标区保留原有曲线，即共用坐标区
        hold(gca, "on");

        % 绘制曲线
        for i =1:length(strcolor)
            plot(omegadivomegaD,y_ks(i,:),'Color',strcolor{i},'LineStyle',strline{i},'LineWidth',2);
        end

        % 颜色填充
        for i =1:(length(strcolor)-1)
            % 获取句柄
            lineFill = fill([omegadivomegaD,fliplr(omegadivomegaD)],[y_ks(i,:),fliplr(y_ks(i+1,:))],hex2rgb(strcolor{i+1}));
            % 设置填充图形的透明度与边线的透明度
            set(lineFill, 'FaceAlpha',0.5, 'EdgeAlpha',0 );
        end

        % 为图例文本设置不同颜色
        % hex2rgb：将16进制颜色代码转为rgb
        % string：将数值型矩阵转为字符串型矩阵
        % join：将矩阵中的字符串以逗号连接
        % char：将字符串型转为字符型
        legend(strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{1})),',')),'}\bfU+N'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{2})),',')),'}\bfU+N+B'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{3})),',')),'}\bfU+N+B+PD'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{4})),',')),'}\bfU+N+B+PD+DS'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{5})),',')),'}\bfU+N+B+PD+DS+NP'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{6})),',')),'}\bfU+N+B+PD+DS+NP+DC'),...
            'FontName','times new roman','FontSize',12,'Location','South','box','off');

        % 设置横坐标刻度
        set(gca, 'XTick', 0:0.2:1);
        % 设置x,y轴的副刻度
        set(gca, 'xminortick', 'on');
        set(gca, 'yminortick', 'on');
        % 设置x轴标签
        xlabel('\bf{\omega/\omega_{D}}');
        % 设置y轴标签
        ylabel('\bf{\kappa_s [ \times 10^{-2} \cdot pW \cdot m^{-1} \cdot K^{-1} ]}');
        % 设置坐标轴边框
        box(gca, 'on');
    end

%% 不同散射机制下的晶格热导率
    function y_kl = func_kl(T,Texp,klexp)
        y_kl_U_N = zeros(1,length(T));
        y_kl_U_N_B = zeros(1,length(T));
        y_kl_U_N_B_PD = zeros(1,length(T));
        y_kl_U_N_B_PD_DS = zeros(1,length(T));
        y_kl_U_N_B_PD_DS_NP = zeros(1,length(T));
        y_kl_total = zeros(1,length(T));
        for i = 1:length(T)
            x = hbar*omega/(kB*T(i));
            y_ks_U_N = 4*pi*kB^4*T(i)^3/(va*h^3).*x.^4.*exp(x)./(func_U_N(T(i)).*(exp(x)-1).^2);
            y_ks_U_N_B = 4*pi*kB^4*T(i)^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B(T(i)).*(exp(x)-1).^2);
            y_ks_U_N_B_PD = 4*pi*kB^4*T(i)^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B_PD(T(i)).*(exp(x)-1).^2);
            y_ks_U_N_B_PD_DS = 4*pi*kB^4*T(i)^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B_PD_DS(T(i)).*(exp(x)-1).^2);
            y_ks_U_N_B_PD_DS_NP = 4*pi*kB^4*T(i)^3/(va*h^3).*x.^4.*exp(x)./(func_U_N_B_PD_DS_NP(T(i)).*(exp(x)-1).^2);
            y_ks_total = 4*pi*kB^4*T(i)^3/(va*h^3).*x.^4.*exp(x)./(func_total(T(i)).*(exp(x)-1).^2);

            y_kl_U_N(i) = trapz(x, y_ks_U_N);
            y_kl_U_N_B(i) = trapz(x, y_ks_U_N_B);
            y_kl_U_N_B_PD(i) = trapz(x, y_ks_U_N_B_PD);
            y_kl_U_N_B_PD_DS(i) = trapz(x, y_ks_U_N_B_PD_DS);
            y_kl_U_N_B_PD_DS_NP(i) = trapz(x, y_ks_U_N_B_PD_DS_NP);
            y_kl_total(i) = trapz(x, y_ks_total);
        end
        y_kl = [y_kl_U_N; y_kl_U_N_B; y_kl_U_N_B_PD; y_kl_U_N_B_PD_DS; y_kl_U_N_B_PD_DS_NP; y_kl_total];
        % 曲线的颜色
        strcolor = {'#000000','#7e01f5','#f40708','#1c8b8b','#4a62ff','#cf5e00'};
        % 曲线的线型
        strline = {'-','-','-','-','-','-'};
        % 在坐标区保留原有曲线，即共用坐标区
        hold(gca, "on");

        % 实验数据
        plot(Texp, klexp, 'mo','MarkerFaceColor','m');
        % 绘制不同散射机制下曲线
        for i =1:length(strcolor)
            plot(T,y_kl(i,:),'Color',strcolor{i},'LineStyle',strline{i},'LineWidth',2);
        end

        % 为图例文本设置不同颜色
        % hex2rgb：将16进制颜色代码转为rgb
        % string：将数值型矩阵转为字符串型矩阵
        % join：将矩阵中的字符串以逗号连接
        % char：将字符串型转为字符型
        legend('\color[rgb]{1 0 1}\bfexp',...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{1})),',')),'}\bfU+N'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{2})),',')),'}\bfU+N+B'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{3})),',')),'}\bfU+N+B+PD'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{4})),',')),'}\bfU+N+B+PD+DS'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{5})),',')),'}\bfU+N+B+PD+DS+NP'),...
            strcat('\color[rgb]{', char(join(string(hex2rgb(strcolor{6})),',')),'}\bfU+N+B+PD+DS+NP+DC'),...
            'FontName','times new roman','FontSize',12,'Location','NorthEast','box','off');

        % 设置x,y轴的副刻度
        set(gca, 'xminortick', 'on');
        set(gca, 'yminortick', 'on');
        % 设置x轴标签
        xlabel('\bf{T [K]}');
        % 设置y轴标签
        ylabel('\bf{\kappa_l [W \cdot m^{-1} \cdot K^{-1}]}');
        % 设置坐标轴边框
        box(gca, 'on');
    end

figure('Name', '频谱晶格热导率图像')
% 绘制频谱晶格热导率图像
ks = func_ks(T);
% 导出频谱晶格热导率数据
exportData('ks', [omegadivomegaD', ks']);

figure('Name', '晶格热导率图像')
% 绘制晶格热导率图像
Tcal = Tstart:10:Tend;
kl = func_kl(Tcal ,Texp, klexp);
% 导出晶格热导率数据
exportData('kl', [Tcal', kl']);
end

hex2rgb

如果因为hex2rgb函数导致报错，请将该文件放到与DebyeCallawayFcn.m同级的目录下，重新运行setup.m。

function rgb=hex2rgb(hex)
    % 16进制表
    HEXTAB=['0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'];
    % 删除开头的#字符 
    hex = hex(2:end);
    rgb = zeros(1,3);
    for i=1:3
        rgb(i) = (find(HEXTAB == upper(hex(2*i - 1)))-1)*16 + (find(HEXTAB == upper(hex(2*i)))-1);
    end
    rgb= rgb/255;
end

loadData.m

%% 加载用于计算的参数和数据文件
function [col1, col2] = loadData(filePath)
% 读取数据时跳过读取第一行的表头
% ReadVariableNames参数为false表示数据文件中没有变量名
try
    table = readtable(filePath,'HeaderLines',1);
    % 第一列数据
    col1 = table2array(table(:,1));
    % 第二列数据
    col2 = table2array(table(:,2));
catch ErrorInfo
    throw(ErrorInfo);
end
end

loadPeriodicTable.m

%% 加载元素周期表，没有形参也可以不用写小括号
function PTABLE = loadPeriodicTable
% 元素符号
chemicalSymbolSet = {
    'H' , 'He', ...
    'Li', 'Be', 'B' , 'C' , 'N' , 'O' , 'F' , 'Ne', ...
    'Na', 'Mg', 'Al', 'Si', 'P' , 'S' , 'Cl', 'Ar', ...
    'K' , 'Ca', 'Sc', 'Ti', 'V' , 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', ...
    'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', ...
    'Rb', 'Sr', 'Y' , 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', ...
    'In', 'Sn', 'Sb', 'Te', 'I' , 'Xe', ...
    'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', ...
    'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W' , 'Re', 'Os', 'Ir', 'Pt', ...
    'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', ...
    'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U' , 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', ...
    'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', ...
    'Rg', 'Cn', 'Nh', 'Fl', 'Mc', 'Lv', 'Ts', 'Og',
    };

% 相对原子质量，这里为列数据，行数据不能直接添加行内注释
atomicWeightSet = [
    1.00794;       % 1     H_  Hydrogen
    4.002602;       % 2     He  Helium
    6.941;          % 3     Li  Lithium
    9.012182;       % 4     Be  Beryllium
    10.811;         % 5     B_  Boron
    12.0107;        % 6     C_  Carbon
    14.0067;        % 7     N_  Nitrogen
    15.9994;        % 8     O_  Oxygen
    18.9984032;     % 9     F_  Fluorine
    20.1797;        % 10    Ne  Neon
    22.98976928;    % 11    Na  Sodium
    24.305;         % 12    Mg  Magnesium
    26.9815386;     % 13    Al  Aluminum
    28.0855;        % 14    Si  Silicon
    30.973762;      % 15    P_  Phosphorus
    32.065;         % 16    S_  Sulfur
    35.453;         % 17    Cl  Chlorine
    39.948;         % 18    Ar  Argon
    39.0983;        % 19    K_  Potassium
    40.078;         % 20    Ca  Calcium
    44.955912;      % 21    Sc  Scandium
    47.867;         % 22    Ti  Titanium
    50.9415;        % 23    V_  Vanadium
    51.9961;        % 24    Cr  Chromium
    54.938045;      % 25    Mn  Manganese
    55.845;         % 26    Fe  Iron
    58.933195;      % 27    Co  Cobalt
    58.6934;        % 28    Ni  Nickel
    63.546;         % 29    Cu  Copper
    65.409;         % 30    Zn  Zinc
    69.723;         % 31    Ga  Gallium
    72.64;          % 32    Ge  Germanium
    74.9216;        % 33    As  Arsenic
    78.96;          % 34    Se  Selenium
    79.904;         % 35    Br  Bromine
    83.798;         % 36    Kr  Krypton
    85.4678;        % 37    Rb  Rubidium
    87.62;          % 38    Sr  Strontium
    88.90585;       % 39    Y_  Yttrium
    91.224;         % 40    Zr  Zirconium
    92.90638;       % 41    Nb  Niobium
    95.94;          % 42    Mo  Molybdenum
    98.0;           % 43    Tc  Technetium
    101.07;         % 44    Ru  Ruthenium
    102.9055;       % 45    Rh  Rhodium
    106.42;         % 46    Pd  Palladium
    107.8682;       % 47    Ag  Silver
    112.411;        % 48    Cd  Cadmium
    114.818;        % 49    In  Indium
    118.71;         % 50    Sn  Tin
    121.76;         % 51    Sb  Antimony
    127.6;          % 52    Te  Tellurium
    126.90447;      % 53    I_  Iodine
    131.293;        % 54    Xe  Xenon
    132.9054519;    % 55    Cs  Cesium
    137.327;        % 56    Ba  Barium
    138.90547;      % 57    La  Lanthanum
    140.116;        % 58    Ce  Cerium
    140.90765;      % 59    Pr  Praseodymium
    144.242;        % 60    Nd  Neodymium
    145.0;          % 61    Pm  Promethium
    150.36;         % 62    Sm  Samarium
    151.964;        % 63    Eu  Europium
    157.25;         % 64    Gd  Gadolinium
    158.92535;      % 65    Tb  Terbium
    162.5;          % 66    Dy  Dysprosium
    164.93032;      % 67    Ho  Holmium
    167.259;        % 68    Er  Erbium
    168.93421;      % 69    Tm  Thulium
    173.04;         % 70    Yb  Ytterbium
    174.967;        % 71    Lu  Lutetium
    178.49;         % 72    Hf  Hafnium
    180.94788;      % 73    Ta  Tantalum
    183.84;         % 74    W_  Tungsten
    186.207;        % 75    Re  Rhenium
    190.23;         % 76    Os  Osmium
    192.217;        % 77    Ir  Iridium
    195.084;        % 78    Pt  Platinum
    196.966569;     % 79    Au  Gold
    200.59;         % 80    Hg  Mercury
    204.3833;       % 81    Tl  Thallium
    207.2;          % 82    Pb  Lead
    208.9804;       % 83    Bi  Bismuth
    210.0;          % 84    Po  Polonium
    210.0;          % 85    At  Astatine
    220.0;          % 86    Rn  Radon
    223.0;          % 87    Fr  Francium
    226.0;          % 88    Ra  Radium
    227.0;          % 89    Ac  Actinium
    232.03806;      % 90    Th  Thorium
    231.03588;      % 91    Pa  Protactinium
    238.02891;      % 92    U_  Uranium
    237.0;          % 93    Np  Neptunium
    244.0;          % 94    Pu  Plutonium
    243.0;          % 95    Am  Americium
    247.0;          % 96    Cm  Curium
    247.0;          % 97    Bk  Berkelium
    251.0;          % 98    Cf  Californium
    252.0;          % 99    Es  Einsteinium
    257.0;          % 100   Fm  Fermium
    258.0;          % 101   Md  Mendelevium
    259.0;          % 102   No  Nobelium
    262.0;          % 103   Lr  Lawrencium
    267.0;          % 104   Rf  Rutherfordium
    268.0;          % 105   Db  Dubnium
    269.0;          % 106   Sg  Seaborgium
    270.0;          % 107   Bh  Bohrium
    270.0;          % 108   Hs  Hassium
    278.0;          % 109   Mt  Meitnerium
    281.0;          % 110   Ds  Darmstadtium
    282.0;          % 111   Rg  Roentgenium
    285.0;          % 112   Cn  Copernicium
    286.0;          % 113   Nh  Nihonium
    289.0;          % 114   Fl  Flerovium
    290.0;          % 115   Mc  Moscovium
    293.0;          % 116   Lv  Livermorium
    294.0;          % 117   Ts  Tennessine
    294.0;          % 118   Og  Oganesson
    ];

PTABLE = containers.Map(chemicalSymbolSet, atomicWeightSet');
end

loadDCParamsDetail.m

%% 加载德拜-卡拉威模型的参数描述
function DCParamsDetail = loadDCParamsDetail
% 德拜-卡拉威模型中所使用的参数
DCParamsSet = {
    'T';
    'va';
    'vl';
    'vt';
    'Vcell';
    'Ncell';
    'Va';
    'beta';
    'ThetaD';
    'omegaD';
    'gamma';
    'Gamma';
    'r';
    'ND';
    'BD';
    'l';
    'DeltaD/D';
    'R';
    'NP';
    };
% 德拜-卡拉威模型中所使用的参数的单位以及描述
DCParamsDetailSet = {
    '温度，K';
    '平均声速，m/s';
    '纵波声速，m/s';
    '横波声速，m/s';
    '原胞/晶胞体积，Å^3';
    '原胞/晶胞内原子数目, 个';
    '平均原子体积，Å^3';
    'N过程与U过程之比，无量纲';
    '德拜温度，K';
    '德拜频率，Hz';
    '格林艾森常数，无量纲';
    '无序度散射因子，无量纲';
    '泊松比，无量纲';
    '位错密度，10^15cm^-2';
    '伯格斯矢量，Å';
    '平均晶粒尺寸，nm';
    'DeltaD/D，无量纲'
    '纳米引入物半径，nm';
    '纳米引入物数密度，10^20m^-3';
    };

DCParamsDetail = containers.Map(DCParamsSet, DCParamsDetailSet);
end

loadConstParams.m

%% 加载常量参数值
function ConstParams = loadConstParams
%{
% 玻尔兹曼常数，J/K，1J=1V*A*s
kB = 1.38066*10^(-23)
% 元电荷的电荷量，C，1C=1A*s
q = 1.60219*10^(-19)
% 普朗克常数，J*s
h = 6.62608*10^(-34)
% 约化普朗克常数，J*s
hbar = 6.62608*10^(-34)/(2*pi);
% 电子静止质量，Kg
me = 9.10953*10^(-31)
% 粒子数密度，单位体积内的粒子数，atoms/m^3
n  = 6.024*10^(28)
% 阿伏伽德罗常数常数，个/摩尔
NA = 6.02214076*10^(23)
% 统一原子质量，碳原子质量的十二分之一，Kg
AMU = 1.6606*10^(-27)
%}

ConstParamsSymbolSet = {
    'kB';
    'q';
    'h';
    'hbar';
    'me';
    'n';
    'NA';
    'AMU';
    };
% 常量参数的单位以及描述
ConstParamsValueSet = [
    1.38066*10^(-23);
    1.60219*10^(-19);
    6.62608*10^(-34);
    6.62608*10^(-34)/(2*pi);
    9.10953*10^(-31) ;
    6.024*10^(28) ;
    6.02214076*10^(23);
    1.6606*10^(-27);
    ];

ConstParams = containers.Map(ConstParamsSymbolSet, ConstParamsValueSet);
end

exportData.m

%导出计算结果
function exportData(dataType, calData)
% 导出的文件路径
filePath = '';
% 表头
headLabel = {};
% 创建文件夹
if ~exist('calData', 'dir')
    mkdir('caldata');
end
if strcmp(dataType, 'ks')
    filePath = './caldata/ks.txt';
    headLabel = {'omega/omegaD','ks_U_N','ks_U_N_B','ks_U_N_B_PD','ks_U_N_B_PD_DS','ks_U_N_B_PD_DS_NP','ks_U_N_B_PD_DS_NP_DC'};
elseif strcmp(dataType, 'kl')
    filePath = './caldata/kl.txt';
     headLabel = {'T','kl_U_N','kl_U_N_B','kl_U_N_B_PD','kl_U_N_B_PD_DS','kl_U_N_B_PD_DS_NP','kl_U_N_B_PD_DS_NP_DC'};
end

% 默认是覆写的方式
writecell(headLabel, filePath,'Delimiter','tab');
% 以追加的方式
writematrix(calData, filePath, 'WriteMode','append', 'Delimiter','tab');
end

化学式.txt

元素名称 化学比
Cu 2
Sn 0.82
In 0.18
Se 2.7
S 0.3

晶格热导率.txt

T(K)       kl(W m^-1 K^-1)
299.28636	1.30317
334.86095	1.0949
374.26664	0.92591
426.9444         	0.76788
478.11707	0.6336
528.60562	0.5249
582.10432	0.45548
633.27699	0.3961
684.58649	0.3413
734.25409	0.32668
784.74263	0.31115
810.05532	0.29654
835.23118	0.27644
860.40704	0.25178

散射因子.txt

参数 值
T 300
va 2199
vl 4105
vt 1955
Vcell 550
Ncell  24
beta 3.3
ThetaD 230
gamma  2.24
Gamma 0.07
r 0.353
ND  8
BD 5.5
l 150
DeltaD/D 0.027
R  10
NP  2

2.1.3 计算结果

2.1.4 参考文献

[1]: CALLAWAY J. Model for Lattice Thermal Conductivity at Low Temperatures. Physical Review, 1959.

[2]: MING H, ZHU G, ZHU C, Boosting Thermoelectric Performance of Cu₂SnSe₃ via Comprehensive Band Structure Regulation and Intensified Phonon Scattering by Multidimensional Defects. ACS Nano, 2021, 15(6): 10532-10541.

2.2 基于德拜-卡拉威模型拟合晶格热导率

2.2.1 计算公式

{\tau}^{-1}={\tau_\mathrm{PD}}^{-1}+{\tau_\mathrm{U}}^{-1}+{\tau_\mathrm{EP}}^{-1}+{\tau_\mathrm{B}}^{-1}

{\tau}^{-1}=A\omega^4+B\omega^2T\exp\left(-\frac{\Theta_\mathrm{D}}{3T}\right)+C\omega^2+\frac{v}{d}

{\tau}_\mathrm{B}^{-1}=\frac{v}{d}=D

PD：点缺陷散射，U：声子U散射， EP：电声散射，B：晶界散射，其中 d 为平均晶粒尺寸

2.2.2 计算代码

1.先运行variable.m生成variable.mat文件，然后运行Debye_Callaway_solve.m文件
2.ks.txt文件用于在其他绘图软件中作图，由于每调用一次func_ks1.m或func_ks2.m都会对ks.txt文件内容覆盖，所以可以采用注释的方式，依次运行与保存
3.想要加入其他散射机制或组合不同散射机制，可以对相应.m文件修改

Debye_Callaway_solve.m

%% Debye_Callaway模型求解频谱晶格热导率
close all; clear; clc;
load('variable.mat');
%频谱晶格热导率绘图
%参数值分别对应A、B、C、D, 注意量级
%下面用func_ks1
figure('Name','Spectrual Lattice Thermoal Conductivity')
%x = 0,y=0.1
hold on
parameters_fit = [20.1,2.9,2.9,2.3];%这里文献计算时B采用的是2.9*10^-19
func_ks1(parameters_fit,Theta_D,'-','b');
%x = 0,y=0.4
parameters_fit = [27.2,2.9,2.5,2.9];  
func_ks1(parameters_fit,Theta_D,'--','r');
%自定义文本代替图例
text(0.7,14,'FeV_{0.9}Nb_{0.1}Sb','FontSize',14,'Color','b')
text(0.7,12.5,'FeV_{0.6}Nb_{0.4}Sb','FontSize',14,'Color','r')
hold off

%下面用func_ks2
figure('Name','Spectrual Lattice Thermoal Conductivity')
text(.35,14,'Fe_{0.985}Co_{0.015}V_{0.6}Nb_{0.4}Sb','FontSize',14,'Color','r')
text(.35,12.5,'Fe_{0.985}Co_{0.015}V_{0.6}Nb_{0.4}Sb(BM)','FontSize',14,'Color','b')
hold on
%x = 0.015,y=0.4(BM)
parameters_fit = [29.1,2.9,6.4,10.8];
func_ks2(parameters_fit,Theta_D,'-','b');
%x = 0.015,y=0.4
parameters_fit = [28.7,2.9,5.7,2.5];
func_ks2(parameters_fit,Theta_D,'--','r');
hold off

func_ks1.m

%% 谱晶格热导率, strline为线型, strcolor为线的颜色, strname为图例的名字
function y_ks = func_ks1(parameters_fit,T,strline,strcolor)
% omega_omega_D = omega/omega_D;
global Kb h hbar v omega_D omega_omega_D 
omega_omega_D  = 0:0.001:1;
omega = omega_omega_D*omega_D;
X = hbar*omega/(Kb*T);  %防止与化学比x冲突
y_ks_U= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U(parameters_fit,T).*(exp(X)-1).^2);
y_ks_U_PD= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U_PD(parameters_fit,T).*(exp(X)-1).^2);
y_ks_U_B_EP= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U_B_EP(parameters_fit,T).*(exp(X)-1).^2);
y_ks_U_PD_B_EP= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U_PD_B_EP(parameters_fit,T).*(exp(X)-1).^2);
y_ks = [y_ks_U; y_ks_U_PD; y_ks_U_B_EP; y_ks_U_PD_B_EP];
%y_ks = y_ks/omega_D*10^14;

%打印输出谱晶格热导率
str_ks = {'ks_U','ks_U_PD','ks__U_B_EP', 'ks_U_PD_B_EP'};
fid_ks= fopen('ks.txt','wt');
fprintf(fid_ks,'%s     %s           %s                          %s                     %s               %s %s\n','omega_omega_D','',str_ks{1}, str_ks{2},str_ks{3},str_ks{4},  '(Wm^{-1}K^{-1})');
%去掉了omega_omega_D = 0的点
for i=2:length(omega_omega_D)
    fprintf(fid_ks,'    %f',omega_omega_D(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U_PD(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U_B_EP(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U_PD_B_EP(i));
    fprintf(fid_ks,'\n');
end
fclose(fid_ks);
%画图
%图窗的名字
% figure('Name','Spectrual Lattice Thermoal Conductivity')
for i =1:length(str_ks)
    hold on
    %获取句柄
    plotyks(i) = plot(omega_omega_D,y_ks(i,:),'Color',strcolor,'LineStyle',strline,'LineWidth',2);
end

% 坐标轴字体与大小
ax.FontName = 'times new roman';
xlabel('\bf{\omega/\omega_{D}}');
%pW为10^-12W
ylabel('\bf{\kappa_s [xW\cdotm^{-1}\cdotK^{-1}]}');
set(gca, 'XTick', 0:0.2:1);
%修改坐标轴线的宽度
set(gca,'LineWidth',2);
ylim([0,16]);
%设置坐标轴刻度字体大小
set(gca,'fontsize',12);
%设置坐标轴边框
box on

func_ks2.m

%% 谱晶格热导率, strline为线型, strcolor为线的颜色, strname为图例的名字
function y_ks = func_ks2(parameters_fit,T,strline,strcolor)
% omega_omega_D = omega/omega_D;
global Kb h hbar v omega_D omega_omega_D 
omega_omega_D  = 0:0.001:1;
omega = omega_omega_D*omega_D;
X = hbar*omega/(Kb*T);  %防止与化学比x冲突
y_ks_U= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U(parameters_fit,T).*(exp(X)-1).^2);
y_ks_U_B= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U_B(parameters_fit,T).*(exp(X)-1).^2);
y_ks_U_PD_EP= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U_PD_EP(parameters_fit,T).*(exp(X)-1).^2);
y_ks_U_PD_B_EP= 4*pi*Kb^4*T^3/(v*h^3).*X.^4.*exp(X)./(func_U_PD_B_EP(parameters_fit,T).*(exp(X)-1).^2);
y_ks = [y_ks_U; y_ks_U_B; y_ks_U_PD_EP; y_ks_U_PD_B_EP];
%y_ks = y_ks/omega_D*10^14;

%打印输出谱晶格热导率
str_ks = {'ks_U','ks_U_B','ks_U_PD_EP', 'ks_U_PD_B_EP'};
fid_ks= fopen('ks.txt','wt');
fprintf(fid_ks,'%s     %s           %s                          %s                     %s               %s %s\n','omega_omega_D','',str_ks{1}, str_ks{2},str_ks{3},str_ks{4},  '(Wm^{-1}K^{-1})');
%去掉了omega_omega_D = 0的点
for i=2:length(omega_omega_D)
    fprintf(fid_ks,'    %f',omega_omega_D(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U_B(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U_PD_EP(i));
    fprintf(fid_ks,'                    %f   ',y_ks_U_PD_B_EP(i));
    fprintf(fid_ks,'\n');
end
fclose(fid_ks);
%画图
%图窗的名字
% figure('Name','Spectrual Lattice Thermoal Conductivity')
for i =1:length(str_ks)
    hold on
    %获取句柄
    plotyks(i) = plot(omega_omega_D,y_ks(i,:),'Color',strcolor,'LineStyle',strline,'LineWidth',2);
end

% 坐标轴字体与大小
ax.FontName = 'times new roman';
xlabel('\bf{\omega/\omega_{D}}');
%pW为10^-12W
ylabel('\bf{\kappa_s [xW\cdotm^{-1}\cdotK^{-1}]}');
set(gca, 'XTick', 0:0.2:1);
%修改坐标轴线的宽度
set(gca,'LineWidth',2);
ylim([0,16]);
%设置坐标轴刻度字体大小
set(gca,'fontsize',12);
%设置坐标轴边框
box on

variable.m

%% 变量定义
%定义全局变量
global Kb h hbar v  Theta_D omega omega_D omega_omega_D A B C  D l 
%玻尔兹曼常数
Kb = 1.380649*10^-23;  %J/K 1J=1V*A*s 1W = 1J/s
%普朗克常数
h = 6.62608*10^-34; %J*s
%简约/约化普朗克常数
hbar = h/(2*pi); 
%平均声速
v = 3998; %m/s
%德拜温度
Theta_D = 464; %K
omega_monega_D = 0:0.001:1;
%德拜频率
omega_D = Theta_D*Kb/hbar; %s^-1
omega = omega_omega_D*omega_D;
%点缺陷散射的前置因子
A = 2.01; %A = 0.35; %*10^(-42)s^3
%U缺陷散射前置因子
B = 3.6; %B = 3.4; %*10^(-18)s/K
%电声散射前置因子
C = 0.29; %C = 0.21; %*10^(-16)
%平均晶粒尺寸
l = 1537; %nm
D = 2.3; %D = v/l/10^9; %10^9s^-1 
save variable

func_PD.m

%% 点缺陷散射(Point Defects)
function y_PD=func_PD(parameters_fit)
global omega_D A omega_omega_D
omega = omega_omega_D*omega_D;
% 如果输入的变量为空则使用在variable.mat中定义的变量的值，由于是全局变量所以为了防止变量修改，直接将后面的10^几次方
A = parameters_fit(1);
y_PD = A*10^(-43)*omega.^4;
end

func_U.m

%% U散射(phonon-phonon U)
function y_U=func_U(parameters_fit,T)
global Theta_D omega_D omega_omega_D B
omega = omega_omega_D*omega_D;
B = parameters_fit(2);
y_U= B*10^(-18)*omega.^2.*T.*exp(-Theta_D/(3*T));
end

func_EP.m

%% 电声散射
function y_EP= func_EP(parameters_fit)
global  omega_omega_D omega_D C
omega = omega_omega_D*omega_D;
C = parameters_fit(3);
y_EP = C*10^(-17)*omega.^2;
end

func_B.m

%% 晶界散射(Boundaries)
function y_B=func_B(parameters_fit)
global D
D = parameters_fit(4);
y_B = D*10^9;
end

func_U_B.m

%% U+B的弛豫时间的倒数
function y_U_B = func_U_B(parameters_fit,T)
y_U_B=  func_U(parameters_fit,T)+  func_B(parameters_fit);
end

func_U_PD.m

%% U+PD的弛豫时间的倒数
function y_U_PD = func_U_PD(parameters_fit,T)
y_U_PD=  func_U(parameters_fit,T) + func_PD(parameters_fit);
end

func_U_B_EP.m

%% U+B+EP的弛豫时间的倒数
function y_U_B_EP = func_U_B_EP(parameters_fit,T)
y_U_B_EP=  func_U(parameters_fit,T) +  func_B(parameters_fit) + func_EP(parameters_fit) ;
end

func_U_PD_EP.m

%% U+PD+EP的弛豫时间的倒数
function y_U_PD_EP = func_U_PD_EP(parameters_fit,T)
y_U_PD_EP=  func_U(parameters_fit,T) + func_PD(parameters_fit) + func_EP(parameters_fit) ;
end

func_U_PD_B_EP.m

%% U+PD+B+EP的弛豫时间的倒数
function y_U_PD_B_EP = func_U_PD_B_EP(parameters_fit,T)
y_U_PD_B_EP=  func_U(parameters_fit,T) + func_PD(parameters_fit) +  func_B(parameters_fit) + func_EP(parameters_fit) ;
end

2.2.3 计算结果

从结果来看，FVNS 组分与文献完全符合，而 FCVNS 组分的U+PD+EP在 w=0 时与文献有出入，可能是文献中计算错误。

2.2.4 参考文献

[1]: Zhu T J , Fu C G , Xie H H , et al. Lattice thermal conductivity and spectral phonon scattering in FeVSb-based half-Heusler compounds[J]. EPL (Europhysics Letters), 2013, 104(4):46003.

[2]: Chenguang Fu, Hanhui Xie, T. J. Zhu, Jian Xie, and X. B. Zhao , Enhanced phonon scattering by mass and strain field fluctuations in Nb substituted FeVSb half-Heusler thermoelectric materials, Journal of Applied Physics 112, 124915 (2012).

目录CONTENT

热电系列-德拜-卡拉威(Debye-Callaway)模型