说明：
1.代入计算时要特别注意量纲
2.由于代码的编写时间不同，导致同一物理量在不同代码文件中所用的符号不一致
3.如果有问题欢迎讨论，欢迎指出错误
1.比热计算
1.1 计算公式
\frac{C_\mathrm{P}}{T} = \gamma+\beta{T^2}+\sum_i\left[ A_i(\Theta_{\mathrm{E}i})^2(T^2)^{-\frac{3}{2}}\frac{\mathrm{e}^{\Theta_{\mathrm{E}i}/T}}{\left(\mathrm{e}^{\Theta_{\mathrm{E}_i}/T}-1\right)^2}\right]
\gamma = \frac{\pi^2}{3}k_\mathrm{B}NE_\mathrm{F}=1.36\times10^{-4}\times V^{2/3}{\mathrm{mol}}{n{\mathrm \gamma
}^{1/3}}\frac{m^*}{m_0}
\beta = B(12\pi Nk_{\mathrm{B}}/5)(\Theta_{\mathrm{D}})^{-3}
B=1-\sum_iA_i/3NR
R=8.314\ J\cdot mol^{-1}\cdot K^{-1}
\gamma：索末菲常数，第一项代表电子对比热的贡献；B：德拜模式对比热的贡献； N：化学式中原子数目； R：气体常数（理想气体常数）；第三项是爱因斯坦模式对比热的贡献。
1.2 计算代码
运行 Cp_fit.m 文件，func_Cp.m 是拟合的函数文件通过拟合比热与温度，func_Cp_T 同样是拟合的函数文件，拟合的函数为：Cp/T~T2
理论来讲拟合得到的参数应该是一样的，但是由于拟合的参数众多，不同的拟合函数得到的结果差距比较大，另外如果只是单纯的拟合数据可以得到非常好的拟合曲线，但是要想得到和文献中相同的参数量级需要给定合理的初值以及上下界，这一点是比较困难的。
我下面的演示也仅仅是代入的文献中的数据，虽然也进行了拟合但是效果不好，读者可以自己改进。
Cp_fit.m
% Cp的拟合
close all; clear; clc;
fileID = fopen('Cp.txt'); %获取打开文件的文件句柄(指针)
data = textscan(fileID, '%f%f','headerlines',1); %从第2行开始读取
fclose(fileID);%关闭文件句柄
%将元胞数组转换为矩阵
matrix = cell2mat(data);
%实验值
Texp = matrix(:,1);
Cpexp = matrix(:,2); %J*mol^-1*K^-1
Cp_Texp = Cpexp./Texp; %J*mol^-1*K^-1
T2exp = Texp.^2;
%用T的平方小于1000K^2的数据拟合
Cp_Texp = Cp_Texp(1:11);
T2exp = T2exp(1:11);
%拟合参数赋初值
% parameters_ini = [1,1,1,1,1,1,1,1];
parameters_ini = [6.366,1.428,1.290,22.06,6.055,47.72,17.582,103.73];
% %拟合参数的上下界
% lb = [0,0,0,0,0,0,0,0];
% ub = [100,100,100,100,100,200,100,500];
lb = [6.366,1.428,1.290,22.06,6.055,47.72,17.582,103.73];
ub = [6.366,1.428,1.290,22.06,6.055,47.72,17.582,103.73];
%用Cp~T拟合
%设置迭代次数与精度等信息, 迭代次数越大越容易找到最优解
%设置迭代参数，例如当提示Solver stopped prematurely.lsqcurvefit stopped
% because it exceeded the iteration limit,options.MaxIterations = 400 (the default value).要适当放宽最大迭代次数
options = optimoptions('lsqcurvefit','MaxFunEvals',5000,'MaxIter',10000,'TolFun',1e-20,'TolX',1e-10);
[parameters_fit,error] = lsqcurvefit(@func_Cp, parameters_ini,Texp,Cpexp,lb,ub,options);
disp(parameters_fit);
disp(error);
%0-30K
% T = 0:1:30;
%0-300K
T = 0:1:300;
figure('Name','CP')
figure(1)
hold on
plot(Texp, Cpexp,'o','Color',[1,0,0],'LineWidth',2);
plot(T,func_Cp(parameters_fit(1),T),'--','Color',[0.1020,0.9843,0.0353],'LineWidth',2);
plot(T,func_Cp(parameters_fit(1:2),T)-func_Cp(parameters_fit(1),T),'--','Color',[0,0.2157,0.9882],'LineWidth',2);
plot(T,(func_Cp(parameters_fit(1:4),T)-func_Cp(parameters_fit(1:2),T)),'Color',[0.9843,0.4941,0.0275],'LineWidth',2);
plot(T,func_Cp(parameters_fit(1:6),T)-func_Cp(parameters_fit(1:4),T),'Color',[0.9882,0.1647,0.9922],'LineWidth',2);
plot(T,func_Cp(parameters_fit,T)-func_Cp(parameters_fit(1:6),T),'Color',[ 0.4784,0.2039,0.9882],'LineWidth',2);
plot(T,func_Cp(parameters_fit,T),'Color',[0.0353,0.4863,0.0118],'LineWidth',2);
box on
ax.FontName = 'times new roman';
xlabel('\bf{T [K]}');
ylabel('\bf{C_p [J/mol\cdotK]}');
%0-30K
% set(gca, 'XTick', 0:5:30);
% xlim([0,30]);
% ylim([0,20]);
% set(gca, 'YTick', 0:4:20);
%0-300K
set(gca, 'XTick', 0:50:300);
xlim([0,300]);
ylim([0,80]);
set(gca, 'YTick', 0:10:80);
%修改坐标轴线的宽度
set(gca,'LineWidth',2);
%设置坐标轴刻度字体大小
set(gca,'fontsize',12);
%用Cp_T~T^2拟合
%设置迭代参数，例如当提示Solver stopped prematurely.lsqcurvefit stopped
%because it exceeded the iteration limit,options.MaxIterations = 400 (the default value).要适当放宽最大迭代次数
options = optimoptions('lsqcurvefit','MaxFunEvals',50000,'MaxIter',10000,'TolFun',1e-20,'TolX',1e-10);
[parameters_fit,error] = lsqcurvefit(@func_Cp_T, parameters_ini,T2exp,Cp_Texp,lb,ub,options);
disp(parameters_fit);
disp(error);
%0-30K
% T2 = 0:1:1000;
%0-300K
T2 = 0:1:90000;
figure('Name','CP/T')
figure(2)
hold on
plot(T2exp, Cp_Texp,'o','Color',[1,0,0],'LineWidth',2);
plot(T2,func_Cp_T(parameters_fit(1),T2),'--','Color',[0.1020,0.9843,0.0353],'LineWidth',2);
plot(T2,func_Cp_T(parameters_fit(1:2),T2)-func_Cp_T(parameters_fit(1),T2),'--','Color',[0,0.2157,0.9882],'LineWidth',2);
plot(T2,(func_Cp_T(parameters_fit(1:4),T2)-func_Cp_T(parameters_fit(1:2),T2)),'Color',[0.9843,0.4941,0.0275],'LineWidth',2);
plot(T2,func_Cp_T(parameters_fit(1:6),T2)-func_Cp_T(parameters_fit(1:4),T2),'Color',[0.9882,0.1647,0.9922],'LineWidth',2);
plot(T2,func_Cp_T(parameters_fit,T2)-func_Cp_T(parameters_fit(1:6),T2),'Color',[ 0.4784,0.2039,0.9882],'LineWidth',2);
plot(T2,func_Cp_T(parameters_fit,T2),'Color',[0.0353,0.4863,0.0118],'LineWidth',2);
box on
ax.FontName = 'times new roman';
xlabel('\bf{T^2 [K^2]}');
ylabel('\bf{C_p/T [J/mol\cdotK^{2}]}');
%0-30K
% set(gca, 'XTick', 0:200:1000);
% xlim([0,1000]);
%0-300K
set(gca, 'XTick', 0:10000:90000);
xlim([0,90000]);
ylim([0,1]);
set(gca, 'YTick', 0:0.1:1);
%修改坐标轴线的宽度
set(gca,'LineWidth',2);
%设置坐标轴刻度字体大小
set(gca,'fontsize',12);func_Cp.m
%Thetan(n=1,2,3)为零，赋予1是为了防止出现0/0，所得结果依然为0
function Cp = func_Cp(parameters_fit,T)
gamma = parameters_fit(1)*10^-3; %J*mol^-1*K^-2
switch length(parameters_fit)
%电子贡献
case 1
beta= 0;
A1 = 0; %J*mol^-1*K^-1
Theta1 = 1; %K
A2 =0;
Theta2 = 1;
A3 = 0;
Theta3 = 1;
%晶格贡献，Debye
case 2
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = 0; %J*mol^-1*K^-1
Theta1 = 1; %K
A2 =0;
Theta2 = 1;
A3 = 0;
Theta3 = 1;
%单共振模式
case 4
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = parameters_fit(3); %J*mol^-1*K^-1
Theta1 = parameters_fit(4); %K
A2 =0;
Theta2 = 1;
A3 = 0;
Theta3 = 1;
%双共振模式
case 6
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = parameters_fit(3); %J*mol^-1*K^-1
Theta1 = parameters_fit(4); %K
A2 =parameters_fit(5);
Theta2 = parameters_fit(6);
A3 = 0;
Theta3 = 1;
%三共振模式
case 8
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = parameters_fit(3); %J*mol^-1*K^-1
Theta1 = parameters_fit(4); %K
A2 =parameters_fit(5);
Theta2 = parameters_fit(6);
A3 = parameters_fit(7);
Theta3 = parameters_fit(8);
otherwise
disp('This is an error!')
end
Cp= gamma*T + beta*T.^3 + A1*Theta1^2*T.^(-2).*exp(Theta1./T)./(exp(Theta1./T)-1).^2 + A2*Theta2^2*T.^(-2).*exp(Theta2./T)./(exp(Theta2./T)-1).^2 +A3*Theta3^2*T.^(-2).*exp(Theta3./T)./(exp(Theta3./T)-1).^2;
endfunc_Cp_T.m
%Thetan(n=1,2,3)为零，赋予1是为了防止出现0/0，所得结果依然为0
function Cp_T= func_Cp_T(parameters_fit,T2)
gamma = parameters_fit(1)*10^-3; %J*mol^-1*K^-2
switch length(parameters_fit)
%电子贡献
case 1
beta= 0;
A1 = 0; %J*mol^-1*K^-1
Theta1 = 1; %K %实际上这一项还是为零，赋予1是为了防止出现0/0
A2 =0;
Theta2 = 1;
A3 = 0;
Theta3 = 1;
%晶格贡献，Debye
case 2
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = 0; %J*mol^-1*K^-1
Theta1 = 1; %K
A2 =0;
Theta2 = 1;
A3 = 0;
Theta3 = 1;
%单共振模式
case 4
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = parameters_fit(3); %J*mol^-1*K^-1
Theta1 = parameters_fit(4); %K
A2 =0;
Theta2 = 1;
A3 = 0;
Theta3 = 1;
%双共振模式
case 6
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = parameters_fit(3); %J*mol^-1*K^-1
Theta1 = parameters_fit(4); %K
A2 =parameters_fit(5);
Theta2 = parameters_fit(6);
A3 = 0;
Theta3 = 1;
%三共振模式
case 8
beta = parameters_fit(2)*10^-4; %J*mol^-1*K^-4
A1 = parameters_fit(3); %J*mol^-1*K^-1
Theta1 = parameters_fit(4); %K
A2 =parameters_fit(5);
Theta2 = parameters_fit(6);
A3 = parameters_fit(7);
Theta3 = parameters_fit(8);
otherwise
disp('This is an error!')
end
Cp_T= gamma + beta*T2 + A1*Theta1^2*T2.^(-3/2).*exp(Theta1./sqrt(T2))./(exp(Theta1./sqrt(T2))-1).^2 + A2*Theta2^2*T2.^(-3/2).*exp(Theta2./sqrt(T2))./(exp(Theta2./sqrt(T2))-1).^2 +A3*Theta3^2*T2.^(-3/2).*exp(Theta3./sqrt(T2))./exp(Theta3./sqrt(T2)).^2;
end1.3 计算结果
1.3.1 0-30K的拟合结果
1.3.2 0-300K的拟合结果
1.4 文献结果
代入文献中的参数，发现在0-30K温度范围内，所得到的结果与文献相符，但是当随着温度升高，基于德拜模型所得到的比热快速增大，当温度范围扩大到30-300K则严重偏离文献所给的结果，而文献的内插图中依旧是呈现非常好的拟合效果，理论上德拜模型在甚低温以及甚高温时（与德拜温度相比，Cu2Se德拜温度在292K）与实验值接近，所以在0-30K温度较低时考虑德拜模型贡献会有比较好的拟合效果，而在之后不应该还具有很好的拟合效果。单纯代入数据也可以得到德拜贡献β*T3会以T3快速增大。难道文献采用了不同的公式？
1.5 参考文献
[3]：刘丹.基于德拜模型讨论晶格比热[J].湖北教育学院学报,2007,No.122(08):10-12.
2.德拜模型下声子简谐项的比热
以下内容转载自 电声不语 微信公众号，如有侵权还请作者联系我进行删除
2.1 计算公式
C_{\mathrm{ph},\mathrm{H}}(T)=9n_{\mathrm{m}}k_{\mathrm{B}}\left(\frac{T}{\Theta_\mathrm{D}}\right)^3\int_0^{\frac{\Theta_\mathrm{D}}{T}}\frac{x^4e^x}{(e^x-1)^2}{\mathrm{d}}x
n_{\mathrm{m}}=N/V=N_\mathrm{A}/M：单位质量的原子数
N：晶胞内的原子个数；V：晶胞的体积
N_\mathrm{A}：阿伏伽德罗常数；M_\mathrm{Cu_2Se}：\mathrm{Cu_2Se}化学式中所有原子的相对原子质量之和；M=M_\mathrm{Cu_2Se}/3：\mathrm{Cu_2Se}的平均原子相对质量
\overline{V} = V/N：平均原子体积等于原胞（晶胞）除以原胞（晶胞）内原子数目；\overline{M}=M*M0/N_\mathrm{A}：平均原子质量为分子式的相对原子质量之和除以分子式中原子数目与阿伏伽德罗常数的乘积，M0：C 原子质量的1/12
k_{\mathrm{B}}：玻尔兹曼常数；x={\hbar\omega}/(k_{\mathrm{B}}T)
\Theta_{\mathrm{D}}：德拜温度
\Theta_{\mathrm{D}}=\frac{h}{k_\mathrm{B}}\left[\frac{3n}{4\pi}\left(\frac{N_\mathrm{A}\rho}{M}\right)\right]^{1/3}v_{\mathrm{s}}
h：普朗克常数；n：\mathrm{Cu_2Se}化学式中的原子数目；\rho：\mathrm{Cu_2Se}的密度；v_{\mathrm{s}}：平均声速
v_{\mathrm{s}}=\left[\frac{1}{3}\left(\frac{2}{v_\mathrm{t}^3}+\frac{1}{v_\mathrm{l}^3}\right)\right]^{-\frac{1}{3}}
v_{\mathrm{t}}：横波声速，v_{\mathrm{l}}：纵波声速
热膨胀项的贡献
C_\mathrm{D}=9Bv\alpha^2T
B：体积模量；v：比体积(单位质量的物质所占有的体积)，为密度的倒数；\alpha：线膨胀系数
B=\rho\frac{v_{\mathrm{l}}^2+2v_{\mathrm{t}}^2}{3}
等压热容
C_{\mathrm{p}}=C_{\mathrm{ph},\mathrm{H}}+C_\mathrm{e}+C_\mathrm{D}，电子比热在较高温下和较低温贡献比较大，这里计算时忽略电子比热贡献
C_\mathrm{p}=C_{\mathrm{ph},\mathrm{H}}+C_\mathrm{D}
杜隆-珀替定律(Dulong-Petit law)
固体比热：C_\mathrm{v}=3Nk_\mathrm{B}
液体比热：C_\mathrm{v}=2Nk_\mathrm{B}
2.2 计算代码
specific_heat_HCS.m
% 求材料的简谐声子热以及热膨胀系数贡献的比热，单位: J*g^-1*K^-1 HfCoSb0.8Sn0.2
clear;
clc;
Kb = 1.380649*10^-23; %玻尔兹曼常数, J/K 1J=1V*A*s 1W = 1J/s
h = 6.62608*10^-34; %普朗克常数, J*s
Na = 6.02214e23; %阿伏伽德罗常数, mol^-1
vl = 4603; %纵波声速，m/s
vt = 2689; %横波声速， m/s
d = 10.77; %密度， g/cm3
alpha = 8.9e-6; %线膨胀系数，K^-1
M = 119.73; %原子相对质量平均值为相对分子质量除以原子数目
T = 100:1000; %待计算的温度范围
T = T';
vs = ((1/3)*(1/(vl^3)+2/(vt^3)))^(-1/3); %平均声速, m/s
ThetaD = h/Kb*(3*d*1e6*Na/(4*pi*M))^(1/3)*vs; %德拜温度, K, N: 晶胞内原子数目, V: 晶胞的体积, N/V = d*Na*M
D = 9*(vl^2+2*vt^2)/3*alpha^2*1e-3; %热膨胀项的前置项
nm = Na/M; %单位质量的原子数, nm = Na/M
Cph = zeros(length(T),1);
for i = 1:length(Cph)
Cph(i) = 9*nm*Kb*(T(i)/ThetaD)^3.*integral(@(x)x.^4.*exp(x)./(exp(x)-1).^2, 0, ThetaD/T(i));
end
Cd = D.*T; %热膨胀项的比热贡献
C = Cph + Cd; %总的比热
data = [T, Cph, Cd, C];
Cdp = 3*8.314/M.*T./T; %计算Dulong-Petit值
plot(T,Cdp,'g', T,Cph,'b',T,C,'r');
specific_heat_CS.m
%求材料的简谐声子热以及热膨胀系数贡献的比热，单位: J*g^-1*K^-1 Cu2Se
clear;
clc;
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
vl = 3350; %纵波声速，m/s
vt = 2350 ; %横波声速， m/s
d = 6.7; %密度， g/cm3
M = 68.684; %原子相对质量平均值为相对分子质量除以原子数目
T = 300:773; %待计算的温度范围
T = T';
alpha = zeros(length(T),1); %线膨胀系数，K^-1
vs = ((1/3)*(1/(vl^3)+2/(vt^3)))^(-1/3); %平均声速, m/s
%vs = 2523;
ThetaD = h/Kb*(3*d*1e6*Na/(4*pi*M))^(1/3)*vs; %德拜温度, K, N: 晶胞内原子数目, V: 晶胞的体积, N/V = d*Na*M
nm = Na/M; %单位质量的原子数, nm = Na/M
Cph = zeros(length(T),1);
for i = 1:length(Cph)
Cph(i) = 9*nm*Kb*(T(i)/ThetaD)^3.*integral(@(x)x.^4.*exp(x)./(exp(x)-1).^2, 0, ThetaD/T(i));
end
D = zeros(length(T),1);
Cd = zeros(length(T),1);
for i = 1:length(T)
if T(i)<=773
alpha(i) = 2.3e-5;
D(i) = 9*(vl^2+2*vt^2)/3*alpha(i)^2*1e-3; %热膨胀项的前置项
Cd(i) = D(i)*T(i); %热膨胀项的比热贡献
else
alpha(i) = 10.7e-5;
D(i) = 9*(vl^2+2*vt^2)/3*alpha(i)^2*1e-3; %热膨胀项的前置项
Cd(i) = D(i)*T(i); %热膨胀项的比热贡献
end
end
C = Cph + Cd; %总的比热
Cdp3 = 3*8.314/M.*T./T; %计算Dulong-Petit值 固体比热的极限3NKb 3R/M R:气体常数=Na*Kb，8.314 J*mol^-1*K^-1
Cdp2 = 2*8.314/M.*T./T; %计算Dulong-Petit值 液体比热的极限2NKb 2R/M
plot(T,Cdp2,'g',T,Cdp3,'m',T,Cph,'b',T,C,'r');
legend('C_v(2NK_B)', 'C_v(3NK_B)', 'C_{ph}', 'C_{ph}+C_d');
xlabel('T [K]');
ylabel('C_p [J\cdotg^{-1}\cdotK^{-1}]');
ylim([0.2,0.4]);
xlim([300,773]);
set(gca, 'YTick', 0.2:0.05:0.4);
%set(gca ,'XTicklabel', {'300', '320', '340', '360', '380'});
2.3 计算结果
2.4 参考文献
[2]：薛丽. CuGaTe2半导体电子结构和热电性能的第一性原理研究[D].华中科技大学,2014.
3.两相复合比热计算
3.1 a(g)的m(基体)与b(g)的i(第二相)复合
V_\mathrm m = \frac{a}{\rho_\mathrm m},V_\mathrm i =\frac{b}{\rho_\mathrm i}
V = V_\mathrm m+V_\mathrm i = \frac{a\rho_\mathrm m+b\rho_\mathrm i}{\rho_\mathrm m\rho_\mathrm i}
x(\text{volume\ ratio})= \frac{\frac{b}{\rho_\mathrm i}}{V} = \frac{b\rho_\mathrm m\rho_i}{\rho_\mathrm i(a\rho_\mathrm i+b\rho_\mathrm m)}=\frac{b\rho_\mathrm m}{a\rho_i+b\rho_\mathrm m}=\frac{b}{a\frac{\rho_\mathrm i}{\rho_\mathrm m}+b}\Longrightarrow\frac{a}{b}(\text{mass\ ratio})=\frac{1-x}{x}\frac{\rho_\mathrm m}{\rho_\mathrm i}
\frac{b}{a+b}=\frac{1}{\frac{a}{b}+1}=\frac{1}{\frac{1-x}{x}\frac{\rho_\mathrm m}{\rho_\mathrm i}+1}=\frac{x}{(1-x)\frac{\rho_\mathrm m}{\rho_\mathrm i}+x}
C_\mathrm p = (1-y)C_\mathrm{pm}+yC_\mathrm{pi}
y = \frac{b}{a+b}=\frac{x}{(1-x)\frac{\rho_\mathrm m}{\rho_\mathrm i}+x},\beta=\frac{\rho_\mathrm m}{\rho_\mathrm i}
评论区