固体热容模型的MATLAB数值求解及可视化

MATLAB Numerical Solution and Visualization of Solid Heat Capacity Model

摘要: 固体热容是反映固体热学性质的一个重要物理参数。本文利用科学计算软件MATLAB的数值计算及绘图功能对固体热容理论中的爱因斯坦模型及德拜模型进行了数值求解，并将结果可视化。结果表明：使用MATLAB数值计算固体热容的理论模型，具有数学思想清晰、编程简单及程序执行效率高的特点；所得结果计算精度高、物理图像清晰。

Abstract: Solid heat capacity is an important physical parameter reflecting the thermal properties of solid. In this paper, the Einstein model and Debye model in solid heat capacity theory are solved numerically, and the results are visualized by using the numerical calculation and drawing function of MATLAB. The results show that the theoretical model of solid heat capacity numerically calculated by MATLAB has the characteristics of clear mathematical thought, simple programming and high program execution efficiency. The obtained physical image is clear and the calculation accuracy is high.

1. 引言

固体的热容是反映固体热学性质的一个重要物理参数，其与固体中的晶格振动、热膨胀及相变等现象密切相关。正因为如此，固体热容理论在统计物理教学中占有重要的地位，对固体热容实验结果的理论解释是统计物理学的一个重要问题 [1] [2]。

MATLAB是由美国MathWorks公司开发的一款国际公认的集数值计算、符号运算及可视化功能于一体的优秀科技应用软件，其在科学研究及工程领域有着广泛应用 [3] [4]。与其他计算机语言相比，MATLAB语言简洁，含有丰富的、可直接调用的内部函数，仅需几条简单的语句，就可以完成一大串其他计算机语言如C或Fortran才能完成的任务，且具有较高的执行效率。另外，MATLAB可视化功能强大，使用一条或几条简单的命令，就能将复杂的数学公式、数据用图形甚至动画的形式表现出来。MATLAB已成为科学研究和工程实践中必不可少的工具。

本文将MATLAB的数值计算及可视化功能用在固体热容理论模型的数值求解方面。论文首先对固体热容进行数学抽象，建立适合MATLAB程序实现的数学模型，然后利用MATLAB写出与数学模型相对应的程序并编制成m文件，最后计算了不同近似条件下固体的热容并将计算结果可视化。

2. 固体热容理论

在物理学史上，固体热容理论经历了从杜隆–珀替(Dulong-Petit)经典定律，到爱因斯坦(Einstein)理论，再到德拜(Peter Debye)理论的研究历程 [1] [2]。1820年法国科学家杜隆和珀替将气体分子的热容理论直接应用于固体，用经典统计理论处理了固体的热容 [5]。其假设固体中共有N个原子，总的自由度数为3N，根据经典的能量均分定理，每个自由度上的平均能量为k_BT，晶体总能量E即为3Nk_BT。晶体

的等容热容 $C_V={\left(\frac{\partial E}{\partial T}\right)}_V=3N{k}_B$，为一与温度无关的常数。1875年，H. F. Weber实验发现金刚石和石墨

晶体的室温热容比经典理论预测的3Nk_B小得多 [5]。高温时，该结果与实验值符合较好，在温度较低时则与实验结果(在较低温度时， $C_V\propto T^3$ )相差较大。这说明在低温时经典的能均分定理已不再适用，必须采用新的理论才能解释固体的热容现象。将晶体中N个原子的振动等效为3N个独立谐振子的振动，每个谐振动对应一个频率 $\omega$ (一种振动模式)，频率 $\omega$ 共有3N个值。

根据量子力学理论，频率为 ${\omega }_i$ 的谐振子的能量：

$E_i=\left(\frac{1}{2}+\frac{1}{{\text{e}}^{\hslash {\omega }_i/k_{B}T}-1}\right)\hslash {\omega }_i$ (1)

这3N个谐振子的总能量：

$E={\displaystyle \underset{i=1}{\overset{3N}{\sum }}{E}_i}={\displaystyle \underset{i=1}{\overset{3N}{\sum }}\left(\frac{1}{2}+\frac{1}{{\text{e}}^{\hslash {\omega }_i/k_{B}T}-1}\right)\hslash {\omega }_i}$ (2)

固体的热容为：

$C_V={\left(\frac{\partial E}{\partial T}\right)}_V={\displaystyle \underset{i=1}{\overset{3N}{\sum }}\left[k_B{\left(\frac{\hslash {\omega }_i}{k_{B}T}\right)}^2\frac{{\text{e}}^{\hslash {\omega }_i/k_{B}T}}{{\left({\text{e}}^{\hslash {\omega }_i/k_{B}T}-1\right)}^2}\right]}$ (3)

为了求解上式得到固体的热容，爱因斯坦(Einstein)在1907年对上述问题进行了简化近似，随后荷兰科学家德拜(Peter Debye)在1911年又进一步做了近似处理 [5] [6]。

3. 爱因斯坦模型求解热容–温度曲线

爱因斯坦假设晶体中所有原子都以相同的频率 ${\omega }_E$ (此频率称为爱因斯坦频率)振动，这些振动是独立的，这时晶体中N个原子的振动被等效为3N个独立的振动频率为 ${\omega }_E$ 的谐振动，由(3)式得固体的热容为：

$C_V=3N{k}_B{\left(\frac{\hslash {\omega }_E}{k_{B}T}\right)}^2\frac{{\text{e}}^{\hslash {\omega }_E/k_{B}T}}{{\left({\text{e}}^{\hslash {\omega }_E/k_{B}T}-1\right)}^2}$ (4)

于是，固体的摩尔热容：

$C_{V,m}=3R{\left(\frac{\hslash {\omega }_E}{k_{B}T}\right)}^2\frac{{\text{e}}^{\hslash {\omega }_E/k_{B}T}}{{\left({\text{e}}^{\hslash {\omega }_E/k_{B}T}-1\right)}^2}.$ (5)

式中 $R=N_A{k}_B=8.31\text{J}\cdot {\text{K}}^{-1}\cdot {\text{mol}}^{-1}$ ；N_A为阿伏伽德罗常数，其值为6.02 × 10²³ mol⁻¹；k_B为玻尔兹曼常数，值为1.38 × 10⁻²³ J·K⁻¹。

令 ${\Theta }_E=\hslash {\omega }_E/k_B$， $x=T/{\Theta }_E$，(5)式化为：

$C_{V,m}=3R{x}^{-2}\frac{{\text{e}}^{1/x}}{{\left({\text{e}}^{1/x}-1\right)}^2}$.(6)

以下使用Matlab编程绘制固体摩尔热容 $C_{V,m}$ 随x变化的曲线。

clear; clc

R=8.31;

x=0.001:0.001:1; % x为温度T与爱因斯坦温度Θ_E的比值

Cvm=3*R*x.^(−2).*exp(1./x)./(exp(1./x) −1).^2/4.2; % 此处的4.2是将单位焦耳换算为卡

plot(x, Cvm) %绘制固体摩尔热容随x变化曲线

grid on

hold on

c=xlsread(‘data1.xls’); % 读取金属Cu的热容实验数据

scatter(c(:,1), c(:,2)) % 绘制金属Cu的热容实验数据

xlabel(‘T/Θ_{E}’); ylabel(‘Cv,m4.18 J·K^{−1}·mol^{−1}’)

legend(‘爱因斯坦模型’, ‘实验数据’)

为了和已测固体的热容实验数据进行比较，绘制曲线时将热容的单位由焦耳∙开尔文⁻¹∙摩尔⁻¹(J∙K⁻¹∙mol⁻¹)换算为卡∙开尔文⁻¹∙摩尔⁻¹(4.18 J∙K⁻¹∙mol⁻¹)，所得结果见图1。需要说明的是程序中文件“data1.xls”为金刚石的热容实验数据 [7]。由图1可以看到，爱因斯坦模型与实验测得的金刚石的热容数据符合的较好。在温度很低时固体的热容逐渐趋近于零；高温时，热容趋于常数值，据此可推知高频振动对热容的贡献几乎可以忽略不计。

Figure 1. Variation of solid heat capacity with temperature obtained by Einstein model. The symbol of “O” in the figure indicates the experimental data of diamond heat capacity [7]

图1. 爱因斯坦模型得出的固体热容随温度的变化，图中的“O”为金刚石的热容实验数据 [7]

4. 德拜模型求解热容–温度曲线

德拜假设晶体为各向同性的连续介质，且晶体中不存在光学波而只有长声学波(可视为弹性波)，格波模式密度为：

$D\left(\omega \right)=\frac{V_c}{{\left(2\pi \right)}^3}{\displaystyle \underset{i=1}{\overset{3}{\sum }}{\displaystyle \int \frac{\text{d}s}{\left|{\nabla }_q\omega \right|}}}=\frac{3V_c{\omega }^2}{2{\pi }^2{v}^3}$ (7)

且有 ${\displaystyle {\int }_0^{{\omega }_D}D\left(\omega \right)\text{d}\omega }=3N$ ( ${\omega }_D$ 为最高频率)，得 ${\omega }_D^3=\frac{6{\pi }^{2}N{v}^3}{V_c}$

于是 $D\left(\omega \right)=\frac{9N{\omega }^2}{{\omega }_D^3}$。

有了模式密度后，固体热容表达式(3)可由求和变为积分：

$C_V={\displaystyle {\int }_0^{{\omega }_D}{k}_B{\left(\frac{\hslash \omega }{k_{B}T}\right)}^2\frac{{\text{e}}^{\hslash \omega /k_{B}T}}{{\left({\text{e}}^{\hslash \omega /k_{B}T}-1\right)}^2}D\left(\omega \right)\text{d}\omega }=9N{k}_B{\left(\frac{T}{{\Theta }_D}\right)}^3{\displaystyle {\int }_0^{{\Theta }_D/T}\frac{{\text{e}}^x{x}^4}{{\left({\text{e}}^x-1\right)}^2}\text{d}x}$ (8)

式中 $\frac{\hslash \omega }{k_{B}T}=x$， ${\Theta }_D=\frac{\hslash {\omega }_D}{k_B}$ 为德拜温度。

在德拜近似下固体的摩尔热容为：

$C_{V,m}=9R{\left(\frac{T}{{\Theta }_D}\right)}^3{\displaystyle {\int }_0^{{\Theta }_D/T}\frac{{\text{e}}^x{x}^4}{{\left({\text{e}}^x-1\right)}^2}\text{d}x}$. (9)

以下使用MATLAB编程绘制0~2500 K温度范围内由德拜模型计算的Pb、Cu及金刚石(德拜温度 ${\Theta }_D$ 分别为88 K、345 K及1860 K)的摩尔热容 $C_{V,m}$ 随温度T的变化曲线。

clear; clc

for n=1:3 % 此循环语句用来计算不同固体在不同温度下的摩尔热容

td=[88, 345, 1860]; %td分别为Pb、Cu及金刚石的德拜温度

td=td(1, n);

for T=0.1:10:2500 % 此循环语句用来计算并绘制给定固体的摩尔热容随温度变化曲线

t=td/T;

s=integral(@(x)x.^4.*exp(x)./(exp(x) −1).^2, 0, t);

cvm=9*8.31*(1/t)^3*s;

plot(T, cvm, ‘o’)

hold on

end

end

grid on

legend(‘-Pb’,‘-Cu’, ‘-.金刚石’)

xlabel(‘T(K)’); ylabel(‘C-{v,m}J·K^{−1}·mol^{−1}’)

Figure 2. Variation of solid heat capacity with temperature calculated by Debye model

图2. 德拜模型计算的固体热容随温度的变化

由德拜模型编程计算的固体Pb、Cu及金刚石的摩尔热容 $C_{V,m}$ 随温度T的变化曲线如图2所示。由图可知：温度较低时，三种固体的热容随温度增加而急剧增加；温度较高时，热容均趋于常数24.9 J∙K⁻¹∙mol⁻¹ (3R)，固体摩尔热容 $C_{V,m}$ 随温度T的变化趋势与固体热容实验结果相符。

5. 两种模型计算热容的比较

由爱因斯坦模型(4)式知：

$\frac{C_V}{3N{k}_B}=\frac{x^{-2}{\text{e}}^{1/x}}{{\left({\text{e}}^{1/x}-1\right)}^2}$.(10)

由德拜模型(8)式知：

$\frac{C_V}{3N{k}_B}=3{\left(\frac{T}{{\Theta }_D}\right)}^3{\displaystyle {\int }_0^{{\Theta }_D/T}\frac{{\text{e}}^x{x}^4}{{\left({\text{e}}^x-1\right)}^2}\text{d}x}$. (11)

以 $\frac{C_V}{3N{k}_B}$ 为纵坐标， $\frac{T}{\Theta }$ 为横坐标，使用MATAB编程绘制 $\frac{C_V}{3N{k}_B}~\frac{T}{\Theta }$ 的变化曲线，代码如下：

clear; clc

n=length(0.001:0.001:2);

a=zeros(1, n);

x=1;

for t=0.001:0.001:2 % 此循环语句用来求解不同T/Θ时由德拜模型计算的固体摩尔热容值

s=integral (@(x)x.^4. *exp(x)./(exp(x)−1).^2, 0, 1/t);

cv=3*t^3*s;

a(1, x)=cv;

x=x+1;

end

cvd=a;

t=0.001:0.001:2;% 以下三行代码绘制爱因斯坦模型计算的固体摩尔热容随T/Θ变化曲线

cva=t.^(−2).*exp(1./t)./(exp(1./t)−1).^2;

plot(t, cva, ‘-k’)

hold on

plot(t, cvd, ‘-k’)% 绘制德拜模型计算的固体摩尔热容随T/Θ变化曲线

c=xlsread(‘data2.xls’);

scatter(c(:,3)/345, c(:,4)/8.31/3, ‘ok’)% 绘制不同T/Θ下的金属Cu摩尔热容实验数据

xlabel(‘T/Θ’); ylabel(‘Cv/3NK_{B}’)

legend(‘爱因斯坦模型’, ‘德拜模型’, ‘实验数据’)

Figure 3. Comparison of heat capacities calculated by Einstein model and Debye Model (experimental data from reference [8] )

图3. 爱因斯坦模型与德拜模型计算热容的比较(实验数据来自文献 [8] )

两种模型计算热容的比较结果见图3，图中的“O”为金属Cu在温度较低时的实验数据 [8]，内插图为T/Θ在[0, 0.15]区间时的放大图。程序中文件“data2.xls”为金属Cu在不同T/Θ(低温)时的摩尔热容实验数据。由图可见：温度较低时，爱因斯坦模型计算的热容结果与实验数据偏差较大，而德拜模型的计算结果更接近于实验数据，这主要是由于爱因斯坦模型过于简单，忽视了各格波对热容贡献的差异 [6]。

6. 结论

本文利用MATLAB对固体热容理论的爱因斯坦热容模型和德拜热容模型进行了数值求解及结果可视化分析研究。结果表明：利用MATLAB求解固体热容理论模型，可避免冗长的数学推导及复杂的积分求解过程，具有编程简单、求解速度快等特点。本研究所得结果与固体热容理论及实验结果相符，物理图像清晰，所用方法可为相关科学研究及工程实践提供参考。

基金项目

安徽省光电信息科学与工程卓越工程师教育培养计划——省级“六卓越、一拔尖”卓越人才培养创新项目(2018zygc012)；安徽省光电物理实验实训中心——示范实验实训中心(2017sxzx15)及安徽省大学物理教学团队项目(2019jxtd046)。