《MATLAB》课程论文MATLAB在热物理学中的应用姓名：田晓霞学号：12010245379专业：通信工程指导老师：汤全武学院：物理电气信息学院完成日期：2011.12.1MATLAB 在热物理学中的应用 (田晓霞 12010245379 2010级通信工程)【摘 要】 基于MATLAB 的数值计算、可视化图形处理、开放式以及可扩充体系结构的特点,并用高性能语言 MATLAB 在大学物理热物理学中的一些应用，包括在固体热容量的三种模型、理想气体定容比热回归分析和理想气体的热力学分析中的应用等对其进行数据处理。

【关键词】 MATLAB ；顺磁性固体；负温度状态；热力学；热传导；热扩散一. 问题的提出之固体热容量的三种模型热容量是热力学系统的一个重要响应函数。

经典理论曾用能量均分定理讨论了晶体在高温情况下的热容量，成功地解释了杜隆－珀替定律。

但是，经典理论不能说明低温下热容量随温度的降低而减小，以及它是系统特征量这两个实验事实。

1907年，爱因斯坦应用量子概念处理晶体振动，定性地说明了固体的热容量随温度降低而趋于零的规律。

1917年，德拜修改了爱因斯坦模型，出了3T 定律，使固体热容量理论在定量上与实验结果相符合。

1.固体热容量的经典模型－杜隆-珀替定律按照经典理论，由N 个原子或离子组成的固体可视为3N 个相互独立的经典线性谐振子的集合。

由能量均分定理，每个线性简谐振子的能量为kT ，固体的内能为U =3NkT ，热容量为3V C N k = (1)此即杜隆-珀替定律。

问题1：应用玻尔兹曼统计求经典固体的定容热容量。

(1) 解题分析经典固体可视为3N 个相互独立的经典线性谐振子的集合,每个经典线性谐振子的能量为()222212r p re m w m=+(2)其中,212rp m是两原子相对运动的动能，1212m m m m m =+为约化质量,r 是两原子间的距离，ω为振动的圆频率。

振动配分函数为dr dr p hpe zrv r ⎰⎰+-=)(2122221ωμμβ(3)求出配分函数后，再利用热力学公式13ln U NZ β∂=-∂ ， V VU C T ∂⎛⎫= ⎪∂⎝⎭ (4)可求得经典固体的热容量。

(2) Matlab 程序:syms V h beta N k T mu omiga r p; %用syms 定义10个符号变量 d=beta/2*mu;e=beta*mu*omiga^2/2; %求符号表达式的值zp=2/(d)^(1/2)*int(exp(-p^2),0,inf); %求符号表达式的值 zr=2/(e)^(1/2)*int(exp(-r^2),0,inf); %求符号表达式的值 Zv= zp*zr/h; %振动配分函数 Uv=-3*N*diff(log(Zv),beta); %求符号表达式的值 beta=1/k/T; %求符号表达式的值Uv1=eval(simplify(Uv)); %内能Cv=diff(Uv1,T); %热容量 运行结果为：Zv =2/(beta/mu)^(1/2)*pi/(beta*mu*omiga^2)^(1/2)/h Uv1 =3*N*k*T Cv =3*N*k运行结果表明，杜隆-玻替定律在固体的温度较高时与测量结果符合，但在常温和低温下与实验结果严重不符。

事实上，固体热容量是与温度和固体特性有关的量，并非该定律所描述的那样是与二者无关的常量。

杜隆-玻替定律与实验事实偏离是对经典热力学理论的严重挑战。

2.爱因斯坦模型爱因斯坦将量子观点应用于固体热容量的研究，把固体看作由3N 个频率相同的，近独立的量子线性谐振子所组成的系统，应用玻尔兹曼统计得到了固体的内能和热容量表达式，这是继普朗克辐射理论之后，利用量子理论处理问题的第二个成功范例。

问题2：应用玻尔兹曼统计求爱因斯坦固体的内能和定容热容量。

（1）解题分析量子线性谐振子的能量为12n n εω⎛⎫=+⎪⎝⎭(0,1,2,3,...n = （1） 谐振子的配分函数为1()210n Z eω∞-+=∑ （2）固体的内能和热容量分别为13ln U NZ β∂=-∂ （3）V VU C T ∂⎛⎫=⎪∂⎝⎭ （4）（2）Matlab 程序clearsyms Z1 beta n hbar w U N k T Cv; %用syms 命令定义10个符号变量 Z1=simplify(symsum(exp(-beta*hbar*w*(n+1/2)),'n',0,inf)); %应用函数规则对括号中的求和函数进行化简后得Z1U=simplify(-3*N*diff(log(Z1),'beta')); %先对其中算是求导，在对起化简 beta=1/k/T; %求beta 的表达式U1=subs(U); %应用U 的表达式求出其中的U1的值Cv=simplify(diff(U1,T)); %将U1对T 的导数求出后在进行化简得到的值 运行结果：Z1 =1/(-1+exp(beta*hbar*w))*exp(1/2*beta*hbar*w)U =3/2*N*hbar*w*(exp(beta*hbar*w)+1)/(-1+exp(beta*hbar*w)) U1 =3/2*N*hbar*w*(exp(1/k/T*hbar*w)+1)/(-1+exp(1/k/T*hbar*w))Cv =3*N*hbar^2*w^2*exp(1/k/T*hbar*w)/(-1+exp(1/k/T*hbar*w))^2/k/T^2其数学表达式为：121e e1Z βωβω=- , //3(e 1)2(e1)kT kTN U ωωω+=- ,2//2e 3(e1)kTV kT C N k kT ωωω⎛⎫= ⎪-⎝⎭ 3.德拜模型1917年,德拜完成了他的固体热容量理论，他把固体看成连续介质，认为原子的振动形成各种简正频率的弹性驻波，而把整个固体原子的微振动看作这些弹性驻波的叠加，每一个简正频率的弹性波的能量与同一频率简谐振子的能量是一样的。

而弹性波又可分为纵波和横波，并且纵波和横波的波速均为一常数。

根据这一思想，德拜从固体中原子振动的频率着手，得出固体的内能和定容热容量分别为3/3d 9()e 1D TxDTx x U N kT θθ=-⎰（1）4/32e d 9()(e 1)D xTV xDTx x C N k θθ=-⎰（2)其中，D x k TTθω==, D θ称为德拜频率。

德拜的理论在低温区与实验符合得相当好，与实验发现的低温下热容量与T 3成正比的规律相一致，因此被称为德拜T 3律。

问题3：绘制杜隆－珀替定律、爱因斯坦模型和德拜模型的固体热容量随温度变化曲线，并讨论其在高、低温两端的性质。

（1）解题分析① 杜隆-玻替定律 113V C y Nk== （1）② 爱因斯坦模型 令 E x kTTθω==，可将爱因斯坦固体热容量表达式改写为()22222e e 3e 1e 1EExT VE x TC y x N k T θθθ⎛⎫=== ⎪⎝⎭⎛⎫-- ⎪⎝⎭（2）③ 德拜模型令 D x kTTθω==，将德拜理论中热容量的表达式34/2e d 9(e 1)D x TV xD T x x C N k θθ⎛⎫= ⎪-⎝⎭⎰（3）改写为3344/322e d 1e d 333(e 1)(e 1)D yy Tx Vy yD C T y yy y y N k x θθ⎛⎫⎛⎫=== ⎪⎪--⎝⎭⎝⎭⎰⎰（4）下面，采用数值方法计算上述积分， （2）Matlab 程序clf %清图形窗口x=0:0.01:1.3; %定义一个步长为0.01，x 从0到1.3 y1=1; %杜隆－珀替定律y2=(1./x).^2.*exp(1./x)./(exp(1./x)-1).^2; %爱因斯坦模型的热容量 i=0; %以下采用循环语句计算德拜模型的数值积分for x1=0.7692:0.5:100 %用for 语句定义一个步长为0.5，初值为0.7692，终值为100的变量 i=i+1;a(i)=quadl('exp(y).*y.^4./(exp(y)-1).^2',0.001,x1); %德拜模型的热容量 y3(i)=a(i).*3./x1.^3; %给出y3的表达式 end %for 循环结束x1=0.7692:0.5:100; %定义变量x1的步长为0.5，从0到100plot(x,y1,x,y2,1./x1,y3) %用plot 函数分别作出x,y1;x,y2;x,y3坐标上的曲线 axis([0,1.3,0,1.1]), %设置坐标轴 xlabel('T/\theta'), %加x 轴说明 ylabel('Cv/3Nk') %加y 轴说明 图1 固体热容量三种理论结果的比较从图 1 可知，在高温端，爱因斯坦模型和德拜模型的曲线都趋近于杜隆-玻替定律，说明经典理论是量子理论的高温（或低频）近似。

运行结果表明，在低温端，爱因斯坦的热容量曲线比实验曲线要平缓一些,而德拜模型的热容量在低温端随温度的变化要比爱因斯坦模型来的快，与温度的三次方成比例,因此比爱因斯坦模型更符合实验结果。

二： 顺磁性固体的热力学性质顺磁性固体的理论模型是，磁性离子定域在晶体的特定格点上，认为离子间彼此相距甚远，相互作用可略去不计。

因此，顺磁性固体是由定域、近独立的磁性 离子组成的系统，遵从玻耳兹曼分布。

(1) 顺磁体的热力学性质问题4：计算顺磁体的磁化强度、内能和熵。

（1）解题分析假定磁性离子的总角动量量子数为12，磁矩大小为2e mμ=-（1）其中，μ在外场中的能量的可能值为－μB （磁矩沿外磁场方向）和μB （磁矩逆外磁方向），B 为外磁场的磁感应强度。

由此，磁性离子的能量为: B B εμμ=-+ (2)离子的配分函数为: 1eeeBBZ βεβμβμ--==+∑ (3)磁化强度：1ln N m Z B B∂=-∂ (4)内能：1ln U NZ β∂=-∂ (5)熵：11(ln ln )S Nk Z Z ββ∂=-∂ (6)（2）Matlab 程序%① 磁化强度syms Z1 beta T k mu N B %用syms 定10个符号变量Z1=exp(beta*mu*B)+exp(-beta*mu*B); %给出Z1表达式m=simplify(N./beta.*diff(log(Z1),B)) %应用函数规则对其进行化简 运行结果：m=N*mu*(exp(beta*mu*B)-exp(-beta*mu*B))/(exp(beta*mu*B)+exp(-beta*mu*B)) 数学表达式为： e e eeB B BBm N βμβμβμβμμ---=+令x = βμB ，y 1= m / Nμ，绘制x -y 1曲线。