matlab在热学中的应用
§4-5 固体的热力学性质
本节利用MATLAB 来处理固体热容量的三种模型、顺磁性固体及负温度状态。
4.5.1 固体热容量的三种模型
热容量是热力学系统的一个重要响应函数。经典理论曾用能量均分定理讨论了晶体在高温情况下的热容量,成功地解释了杜隆-珀替定律。但是,经典理论不能说明低温下热容量随温度的降低而减小,以及它是系统特征量这两个实验事实。1907年,爱因斯坦应用量子概念处理晶体振动,定性地说明了固体的热容量随温度降低而趋于零的规律。1917年,德拜修改了爱因斯坦模型,导出了T 定律,使固体热容量理论在定量上与实验结果相符合。
(1)固体热容量的经典模型-杜隆-珀替定律
按照经典理论,由N 个原子或离子组成的固体可视为3N 个相互独立的经典线性谐振子的集合。由能量均分定理,每个线性简谐振子的能量为kT ,固体的内能为U =3NkT ,热容量为
3
C V 3Nk
此即杜隆-珀替定律。
● 题目(ex4511)
应用玻尔兹曼统计求经典固体的定容热容量。 ● 解题分析
经典固体可视为3N 个相互独立的经典线性谐振子的集合, 每个经典线性谐振子的能量为
e =
1
p r 2+m 2w 2r 2) (2m
其中,
m 1m 212
p r 是两原子相对运动的动能,m =为约化质量, r 是两原子间的2m m 1+m 2
距离,ω为振动的圆频率。振动配分函数为
1v
Z 1=
h
蝌e
-
b 2
(p r +m 2w 2r 2) 2m
d p r d r
bm ×w 22
r 2
+?
-1骣çç蝌=e çh ç桫-?
b 2
p r
2m
骣-÷ççd p r ÷e ÷ç÷ç÷桫
÷
d r ÷÷÷÷
求出配分函数后,再利用热力学公式
U =-3N
∂⎛∂U ⎫
ln Z 1 , C V = ⎪∂β⎝∂T ⎭V
可求得经典固体的热容量。
● 程序(ex4511)
syms V h beta N k T mu omiga r p; 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 个频率相同的,近独立的量子线性谐振子所组成的系统,应用玻尔兹曼统计得到了固体的内能和热容量表达式,这是继普朗克辐射理论之后,利用量子理论处理问题的第二个成功范例。
● 题目(ex4512)
应用玻尔兹曼统计求爱因斯坦固体的内能和定容热容量。 ● 解题分析
量子线性谐振子的能量为
εn = +n ⎪ ω (n =0,1,2,3,.......)
⎝2⎭
谐振子的配分函数为 Z 1=
∞0
1
-(+n ) ω2
⎛1⎫
∑e
固体的内能和热容量分别为
U =-3N
∂⎛∂U ⎫ln Z 1, C V = ⎪∂β⎝∂T ⎭V
● 程序(ex4512)
clear
syms Z1 beta n hbar w U N k T Cv; Z1=simplify(symsum(exp(-beta*hbar*w*(n+1/2)),'n',0,inf)) U=simplify(-3*N*diff(log(Z1),'beta')) beta=1/k/T; U1=subs(U)
Cv=simplify(diff(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
e 3N ω(e ω/kT +1) e ω/kT ⎛ ω⎫
即 Z 1=β ω, U =, C V =3Nk ⎪ ω/kT ω/kT 2
e -12(e-1) -1) ⎝kT ⎭(e
(3)德拜模型
1917年, 德拜完成了他的固体热容量理论,他把固体看成连续介质,认为原子的振动形成各种简正频率的弹性驻波,而把整个固体原子的微振动看作这些弹性驻波
1
β ω2
2
的叠加,每一个简正频率的弹性波的能量与同一频率简谐振子的能量是一样的。而弹性波又可分为纵波和横波,并且纵波和横波的波速均为一常数。根据这一思想,德拜从固体中原子振动的频率着手,得出固体的内能和定容热容量分别为
U =9NkT (
x =其中,
T
θD
)
3
⎰
θD /T
x 3d x T 3θD /T e x x 4d x
; C V =9Nk () ⎰。 x x 20e -1θD (e-1)
ωθD
=, θD 称为德拜频率。德拜的理论在低温区与实验符合得相当好,kT T
与实验发现的低温下热容量与T 3成正比的规律相一致,因此被称为德拜T 3律。
● 题目(ex4513)
绘制杜隆-珀替定律、爱因斯坦模型和德拜模型的固体热容量随温度变化曲线,并讨论其在高、低温两端的性质。
● 解题分析
① 杜隆-玻替定律 y 1=② 爱因斯坦模型 令 x =
C V
=1 3Nk
ωθE
=,可将爱因斯坦固体热容量表达式改写为 kT T
C V e T e x ⎛θE ⎫2
y 2== ⎪=x 22x 3Nk ⎝T ⎭⎛θE
⎫e -1() e T -1⎪
⎝⎭
2
θE
③ 德拜模型 令 x =
ωθD
= k T T
3
将德拜理论中热容量的表达式
⎛T ⎫
C V =9Nk ⎪
θ⎝D ⎭
⎰
θD /T
e x x 4d x
(ex -1) 2
3
改写为
⎛T ⎫C
y 3=V =3 ⎪
3Nk ⎝θD ⎭
3
⎰
θD /T
e y y 4d y ⎛1⎫
=3 ⎪
(ey -1) 2⎝x ⎭
⎰
x
e y y 4d y
y 2
(e-1)
下面,采用数值方法计算上述积分,
● 程序(ex4513) clear,clf
x=0:0.01: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 i=i+1;
a(i)=quadl('exp(y).*y.^4./(exp(y)-1).^2',0.001,x1); %德拜模型的热容量 y3(i)=a(i).*3./x1.^3;
end
x1=0.7692:0.5:100;
plot(x,y1,'k-',x,y2,'.r-',1./x1,y3,'-bo') axis([0,1.3,0,1.1]), xlabel('T/\theta'), ylabel('Cv/3Nk')
legend('杜隆-珀替定律',' 爱因斯坦模型',' 德拜模型')
从图4-5-1 可知,在高温端,爱因斯坦模型和德拜模型的曲线都趋近于杜隆-玻替定律,说明经典理论是量子理论的高温(或低频)近似。 实验表明,在低温端,爱因斯坦的热容量曲线比实验曲线要平缓一些, 而德拜模型的热容量在低温端随温度的变化要比爱因斯坦模型来的快,与温度的三次方成比例, 因此比爱因斯坦模型更符合实验结果。
4.5.2 顺磁性固体的热力学性质
顺磁性固体的理论模型是,磁性离子定域在晶体的特定格点上,认为离子间彼此相距甚远,相互作用可略去不计。因此,顺磁性固体是由定域、近独立的磁性
离子组成的系统,遵从玻耳兹曼分布。
(1) 顺磁体的热力学性质 ● 题目(ex4521)
计算顺磁体的磁化强度、内能和熵。 ● 解题分析
假定磁性离子的总角动量量子数为,磁矩大小为
12
μ=-
e 2m
其中,μ在外场中的能量的可能值为-μB (磁矩沿外磁场方向)和μB (磁矩逆外磁场方向),B 为外磁场的磁感应强度。由此,
磁性离子的能量为 ε=-μB +μB
离子的配分函数为 Z 1=
∑e βε=e βμ
-
B
-βμB
+e
磁化强度:m =-
N ∂∂
ln Z 1 内能:U =-N ln Z 1 B ∂B ∂β
熵:S =Nk (lnZ 1-β
∂
ln Z 1) ∂β
● 程序(ex4521)
%① (ex45211) 磁化强度
syms Z1 beta T k mu N B
Z1=exp(beta*mu*B)+exp(-beta*mu*B);
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 βμB -e -βμB
即 m =N μβμB -βμB
e +e
令x = βμB ,y 1= m / N μ,绘制x -y 1曲线。 clf
x=-3:0.01:3;
y1=(exp(x)-exp(-x))./(exp(x)+exp(-x));
plot(x,y1,'r-'), xlabel('\beta\muB'),ylabel('m/\muN')
grid on 运行结果如图 4-5-2所示。 % ②(ex45212) 内能
Z1=exp(beta*mu*B)+exp(-beta*mu*B); U=-N.*diff(log(Z1),beta) 运行结果:
U=-N*(mu*B*exp(beta*mu*B)-mu*B*exp(-beta*mu*B))/(exp(beta*mu*B) +exp(-beta*mu*B)) 即
e βμB -e -βμB
U =-N μB βμB -βμB
e +e
令x = βμB ,y 2= U / NkT ,绘制x -y 2曲线。
syms x y; x=-3:0.01:3;
y2=-x.*(exp(x)-exp(-x))./(exp(x)+exp(-x));
plot(x,y2,'r-') xlabel('vB/kT')
ylabel('U/NkT')
grid on
运行结果如图 4-5-3所示。 % ③ (ex45213) 熵 syms N k beta mu B
z1=exp(beia*mu*B)+exp(-beia*mu*B);
S=N*k*(log(z1)-beta*diff(log(z1),beta))
S=N*k*(log(exp(beta*mu*B)+exp(-beta*mu*B))-beta*(mu*B*exp(beta*mu*B)- mu*B*exp(-beta*mu*B))/(exp(beta*mu*B)+exp(-beta*mu*B))) 即 S =Nk ln(e
βμB
+e
-βμB
e βμB -e -βμB
) -Nk βμB βμB -βμB
e +e
令x = βμB ,y 3= S / Nk ,绘制x -y 3曲线。
x=-3:0.01:3;
y3=(log(exp(x)+exp(-x))-x.* (exp(x)-exp(-x))./(exp(x)+exp(-x)));
plot(x,y3)
xlabel('vB/kT') ylabel('S/Nk') grid
读者可根据玻耳兹曼关系来分析图4-5-4的物理意义。 4.5.3 负温度状态
● 题目(ex4531)
设核自旋量子数为, 在外磁场B 下由于磁矩可与外磁场逆向或同向,其能量有两个可能值±
1
2
Be
,简记为±ε。以N 表示系统所含的总核磁矩数,x 和y 分别表2M
示能量为+ε和-ε的核磁矩数。求系统的熵和内能之间的关系并绘制U -S 曲线,由此讨论负温度状态及其意义。
● 解题分析
由热力学知下述关系成立
1⎛∂S ⎫= ⎪ T ⎝∂U ⎭x
在一般系统中熵随内能单调增加,因此T 为正。但是,也存在一些系统,当内能增加时熵反而减小,此时系统处在负温度状态,本题目就是这样的例子。
由题知,总核磁矩数为 N =x +y ,总能量为 U =(x -y ) ε。对该问题,利用玻耳兹曼关系,有 S =k ln
N !
x ! y !
再利用近似公式 l n m ≈m l n m ,可求解。
● 程序(ex4531)
% (ex45311) 熵和内能的关系 clear
syms N e x y U k;
[x,y]=solve('N=x+y','U=x*e-y*e');
s=k*(N*log(N)-x*log(x)-y*log(y)) 运行结果:
s=k*(N*log(N)-1/2*(U+e*N)/e*log(1/2*(U+e*N)/e)-1/2*(-U+e*N)/e*log(1/2*(-U+e*N)/e))
1U 1U 1U 1U
即 S =k (N l n N -+N ) (+N ) (+N ) ((+N ) )
2ε2ε2ε2ε
%(ex45312) 绘制S -U 图形
%下面令yy=S/N/k,xx=U/N/e。 clear
xx=-1:0.01:1;
yy=log(2)-1/2.*(1+xx).*log(1+xx)-1/2.*(1-xx).*log(1-xx);
plot(xx,yy),grid on 运行结果如图4-5-5所示。
讨论: clear syms N e x y U k T; s=k*(N*log(N)-1/2*(U+e*N)/e*lo
g(1/2*(U+e*N)/e)-1/2*(-U+e*N)/e*log(1/2*(-U+e*N)/e));
y=1./T;
y=simplify(diff(s,U))
y=-1/2*k*(log((U+e*N)/e) -log((-U+e*N)/e))/e
可以看出,在U 0时(曲线的右半部分),y =来说明其物理图象。
1
为正, 系统处在正温系统;在T
1
为负, 系统处在负温系统。请读者结合核自旋系统T
§4-6 理想气体的热力学性质
理想气体是近独立粒子体系的典型理想模型之一, 也是热力学和统计物理学理论的一个重要应用对象。在统计物理学中,我们把理想气体分子当作可分辨粒子来看待,用玻耳兹曼分布进行处理。
用玻耳兹曼统计处理问题的一般方法是:根据分子的能量函数计算出配分函数及其对数,然后代入统计热力学公式求气体的内能、压强、熵和其它热力学函数。本节利用MATLAB 的符号运算功能来推导单原子和双原子理想气体的配分函数以及有关热力学函数。
4.6.1 单原子理想气体的热力学性质
● 题目(ex4611)
由玻耳兹曼分布研究单原子理想气体的热力学性质。 ● 解题分析
理想气体分子可视为三维自由粒子,其能量函数为
ε=
122
(p x +p y +p z 2) 2m
-b 22
(p x +p 2y +p z ) 2m
利用配分函数的表达式可计算出理想气体分子的配分函数Z 1
1
Z 1=3
h
蝌蝌蝌e
ゥ
d x d y d z d p x d p y d p z
-b 2
p z 2m
b 2-p x V 2m
=3蝌e d p x
h -?
¥
e
?
-
b 2
p y 2m
d p y e
d p z
b 2-p V
=3(2òe 2m d p ) 3 h 0
将得到的配分函数取对数后代入下列统计热力学公式
U =-N N ∂∂∂ln Z 1, S =Nk (lnZ 1-βln Z 1) , ln Z 1, p =β∂V ∂β∂β
H =U +PV , F =U -TS , G =U -TS +PV
便可求出理想气体的内能、状态方程、熵、焓、自由能和吉布斯函数等热力学性质。
● 程序(ex4611)
syms V h a beta m N k T pi;
Zt=V/h^3*(2/(beta/2/m)^(1/2)*(sym('pi'))^(1/2)/2*erf(inf))^3
% 这里使用了计算误差函数的指令erf
U=-N*diff(log(Zt),beta);
S=N*k*(log(Zt)-beta*diff(log(Zt),beta));
beta=1/k/T;
U1=eval(U)
Cv=diff(U1,T)
P1=N/beta*diff(log(Zt),V)
S1=vpa(eval(S),3)
F=vpa(simplify(U1-T*S1),3)
H=vpa(simplify(U1+P1*V),3)
G=vpa(simplify(U1-T*S1+P1*V),3)
运行结果:
Zt =2*V/h^3*2^(1/2)/(beta/m)^(3/2)*pi^(3/2)
U1 = 3/2*N*k*T
Cv = 3/2*N*k
P1 = N*k*T/V
S1 = N*k*(log(15.7*V/h^3/(1/k/T/m)^(3/2))+1.50)
F = -2.75*N*k*T-1.*N*k*T*log(V/h^3*k*T*m/(1/k/T/m)^(1/2))
H = 2.50*N*k*T
G = -1.75*N*k*T-1.*N*k*T*log(V/h^3*k*T*m/(1/k/T/m)^(1/2))
4. 6.2 双原子理想气体的热力学性质
● 题目(ex4621)
由玻耳兹曼分布研究双原子理想气体的热力学性质。
● 解题分析
双原子理想气体分子可分为弹簧振子模型和刚性哑铃模型两种情况,我们讨论前一种情况。双原子理想气体分子的能量的经典表达式为
e =11骣1122222÷2222çp +p +p +p +p +p +m w r ) () (÷çx y z q j ÷r 2ç2m 2I 桫sin q 2m
其中,第一项为平动动能,第二项为转动动能,第三项为振动动能。利用玻耳兹曼分布配分函数的表达式可计算出双原子理想气体分子的配分函数Z 1为
Z 1=Z 1Z 1r Z 1
其中,平动配分函数为 t r v
Z 1t =1
h 3-βεt e ⎰⎰⎰⎰⎰⎰d x d y d z d p x d p y d p z
β2β2β2-p x -p y -p z V =3⎰e 2m d p x ⎰e 2m d p y ⎰e 2m d p z h -∞-∞-∞∞∞∞
β2-p V 2m =3(2⎰e d p ) 3 h 0∞
转动配分函数为
1Z 1r =2h ⎰⎰⎰⎰e
2p -β2I 2(p θ+1sin θ2p ϕ) d p θd p ϕd θd ϕ -b
sin q p j 21=2h
振动配分函数为 p +? d j 蝌00d q 蝌e -? -b 2p q 2I d p q e d p j
-1Z =e 蝌h v
1b 2(p r +m 2w 2r 2) 2m d p r d r
-bm ×w 22r 21=h +? e 蝌
-? -b 2p r 2m d p r e d r
然后代入统计热力学公式 U =-N N ∂∂⎫ ∂ln Z 1,S =Nk ⎛ln Z 1,P = ln Z -βln Z 11⎪ ⎪β∂V ∂β∂β⎝⎭
求双原子分子理想气体的内能、压强和熵。
● 程序(ex4621)
syms V h beta m N k T I theta phi mu omiga pi;
a=beta/2/m;
Zt= vpa(V/h^3*(2/(a)^(1/2)* pi^(1/2)/2*erf(inf))^3,3); % 平动配分函数
b=beta/2/I/ (sin(theta))^2; c=beta/2/I;
zpth= 2/(b)^(1/2)*pi^(1/2)/2*erf(inf);
zpph=2/c^(1/2)*pi^(1/2)/2*erf(inf);
Zr=simplify(int(int(zpth*zpph,theta,0,pi),phi,0,2*pi)/h^2); % 转动配分函数
d=beta/2/mu; e=beta*mu*omiga^2/2;
zp=2/(d)^(1/2)*pi^(1/2)/2*erf(inf);
zr=2/(e)^(1/2)*pi^(1/2)/2*erf(inf);
Zv= zp*zr/h; %振动配分函数
Z1=simple(Zr*Zt*Zv) %双原子分子配分函数
U=simple(-N*diff(log(Z1),beta)) %双原子分子气体内能
P=simple(N/beta*diff(log(Z1),V)) %双原子分子气体的压强
S=simple(N*k*(log(Z1)-beta*diff(log(Z1),beta))) %双原子分子气体的熵
运行结果为
Z1 =32/beta^3*I*pi^(9/2)/h^6*V*2^(1/2)*m^2/(beta*m)^(1/2)/omiga
U =7/2*N/beta
P =N/beta/V
S =N*k*(log(32/beta^3*I/h^6*V*m^2/omiga*pi^4*2^(1/2)/(beta*m/pi)^(1/2))+7/2)
§4-7 热传导过程
热传导过程的物理规律一般用偏微分方程来描述,除非特殊的边界条件和规则形状外,一般是难以给出解析解的。 在大多数情况下,我们可以通过数值计算来处理这些问题。如§2-4所介绍的,MATLAB 为我们提供了一个求解偏微分方程问题
的有力工具——PDETOOL ,利用这个工具,在用户图形界面(GUI)上可以方便地求解较复杂问题的数值解,并给出动态图形。本节介绍用PDETOOL 处理热传导过程的两个典型实例。
4.7.1 受热金属块的热传导
● 题目(ex4711)
一受热的有矩形裂纹的金属块,金属块的左侧被加热到100℃,在右侧热量以恒定速率降低到周围空气的温度,所有其它边界都是独立的。假定金属块的初始温度为0℃。
● 解题分析
该题目所列热传导方程为常见的抛物线型方程,即
d ∂u -∆u =0 ∂t
其中,u 为温度,t 为时间,d 为热能系数。由题意可列出下列边界条件
① u =100 左侧(Dirichlet 条件)
∂u =-10 右侧(Neumann 条件) ∂n
∂u =0 其他边界(Neumann 条件) ③ ∂n ②
● 用GUI 求解
在命令窗中键入pdetool ,打开图形用户界面,画CSG 模型。
首先画一个矩形(R1),其四角的坐标为x=[-0.5 0.5 0.5 –0.5], y=[-0.8 –0.8 0.8 0.8]。然后画另一个矩形(R2),代表矩形裂纹,它的四角坐标为x=[-0.05 0.05 0.05 –0.05], y=[-0.4 –0.4 0.4 0.4]。选择Options 选项,打开Grid Spacing对话框,在-0.05和0.05处输入X 轴的附加短线,以帮助画出代表裂纹的矩形。然后显示网格和“Snap-to-grid ”风格,这样就很容易画出矩形裂纹了。
金属块的CSG 模型用简单公式来表示:R1-R2。 单击按钮,输入Boundary 模式,选择边界并指定边界条件。 单击 按钮,打开PDE Specification对话框,输入PDE 系数。
需要解决的PDE 表达式为
d ∂u -∇⋅(c ∇u ) +au =f ∂t
初值u 0=u (t 0) ,计算的时间保存在数组tlrst 中。本例中取d=1,c=1, a=0, f=0
。单击
按钮,将网格初始化,再单击按钮改进网格。依次点击Solve\Parameters,打开Solve Parameters对话框,输入初值u0=0,时间[0:0.5:5]。 单击按钮,在11个不同的时间段内求解热传导方程,默认时显示最后时间解的interpolated 图。
若要动态显示热传导的过程,可在Plot Selection 对话框中选择Animation 核选框,选择colormap hot ,单击Plot 按钮,在单独的图形窗口中开始解的图形记录,然后会重复照亮5次。注意到金属块的温度升高很快,为改进照明效果并注重第一秒钟,可试着改变时间列表的表达式logspace(-2,0.5,20)。也可试着改变热能系数d 和右边界的热流量来研究它们是如何影响热扩散的。
● 用命令行函数求解
① 首先建立几何模型和边界条件的M 文件
金属块的几何模型在名为crackg.m 的M 文件中,它给出本题目PDE 模型的几何数据。
NE=CRACKG %给定边界段的数目。
D=CRACKG(BS) %对每个边界段给定一列矩阵。
在 BS 中指定:
Row 1 包含初始参数值
Row 2 包含最终参数值
Row 3 包含左边界区域数量
Row 4 包含右边界区域数量
[X,Y]=CRACKG(BS,S)
% 给定边界点的坐标。BS 指定边界段,S 指定相应参数值。BS 可以是标量。
边界条件保存在名为crackb.m 的M 文件中,本题边界条件数据的函数式M 文件如下:
function [q,g,h,r]=crackb(p,e,u,time)
bl=[1 1 1 1 1 1 1 1;0 0 0 0 0 1 0 0;1 1 1 1 1 1 1 1;3 1 1 1 1 1 1 1;48 48 48 48 48 1 48 48;45 48 48 48 48 3 48 48;49 48 48 48 48 48 48 48;48 48 48 48 48 48 48 48;49 49 49 49 49 49 49 49;48 48 48 48 48 49 48 48;0 0 0 0 0 48 0 0;0 0 0 0 0 48 0 0];
if any(size(u))
[q,g,h,r]=pdeexpd(p,e,u,time,bl);
else
[q,g,h,r]=pdeexpd(p,e,time,bl);
end
② 创建初始网格
调用函数为 initmesh 。
[p,e,t]=initmesh('crackg');
③ 求解热传导方程
u=parabolic(0,0:0.5:5,'crackb',p,e,t,1,0,0,1);
% 用函数Parabolic 求解抛物线型热传导方程, 得到的解u 是一个11列的矩阵
% 其中每一列对应于11个时间点上的解。
④ 由函数生成不同时间点上解的图形
pdeplot(p,e,t,'xydata',u(:,11),'mesh','off',
'colormap','hot')
%用pdeplot 函数生成不同时间点上
%解的图形,如t=5.0时
运行结果:
156 successful steps
0 failed attempts
314 function evaluations
1 partial derivatives
29 LU decompositions 313 solutions
of linear systems
绘制的图形如图4-7-1所示。 4.7.2 放射性棒的热扩散
● 题目(ex4721)
考察一根放射性棒,在其左端持续不断地加以热源,右端恒温。在外边界条件上,棒体与外界发生热交换。同时,由于放射性过程,整个棒体均匀地产生热量。设初始温度为零,求解的问题为
r c ¶u -炎(k ? u ) ¶t f
其中,ρ为密度,c 为棒的热容量,k 为导热系数,f 为放射性热源的强度。该热传导问题是一个三维抛物型方程问题,在柱坐标下可被简化为二维问题。
设棒的密度为7800 kg/m3,热容量为500 Ws/kg ·℃, 导热系数为40 W/m2·℃, 热源强度为为20000 W/m3,右端温度为100℃,外边界周边温度为100℃,热交换系数为50 W/m2·℃, 左端的热流强度为5000 W/m2。
● 解题分析
这是一个圆柱问题,故需将方程改为由r 、z 和θ所确定的柱坐标形式,考虑到对称性,其解与θ无关,方程可转换为
r r c
边界条件为: 抖u -抖t r (kr 抖u ) -抖r z (kr u ) =fr z
(c u ) +qu =g 。对 ① 棒的左端为Neumann 条件。广义Neumann 条件为n 籽
(c u ) =5000r 。 于本问题,c =kr ,边界条件可表达为:z = -1.5, n 籽
② 棒的右端为Dirichlet 条件 z =1.5, u =100
(k u ) =50(100-u ) , ③ 外边界为广义Neumann 条件,表达式为 n 籽对于本
(c u ) +50r ? u 题可表为:r =0.2, n 籽50r 100。
④ 在原问题中柱轴r =0不是边界,但现在降维成二维问题以后它是边界,必须
(c u ) =0。 给出人为边界条件: r =0, n 籽
设初始温度为 u (t=0)=0。
● 解题步骤
利用GUI 按下列步骤求解:
① 打开PDE 图形用户界面,建立几何模型。将棒体画成基线平行于X 轴的矩形,令X 轴为z 方向,Y 轴为r 方向,矩形四角的坐标分别设为(-1.5,0),(1.5,0),(1.5,0.2)(-1.5,0.2),则棒的长度为3,半径为0.2。
② 双击边界,打开Boundary Condition对话框,在对话框中输入边界条件。对于左边界,输入Neumann 条件:q =0,y =5000*y 。右端用Dirichlet 条件:h =1,r =100。对外边界,用Neumann 条件:q =50*y ,g =50*y *100,坐标轴用Neumann 条件:q =0,g =0。
③ 依次点击PDE\PDE Specification , 打开PDE Specification 对话框。PDE ;类型选为抛物型 "Paraboic" ;输入系数:c =40*y ,a =0,d =7800*5000*y ,f =20000*y 。 ④依次点击Solve\Parameters,打开“Solve Parameters ”对话框,取Time=0:100:2000,u (t 0)=0, 其余取默认值。最后点击“OK ”确定。
⑤ 依次点击Plot\Parameters,打开“Plot Selection”对话框,依次点击“Color ”、“Contour ”和“Plot ”,得到图4-7-2所示的图形。如果在该对话框中选择“Animation ”,
则可以观察热量在右侧和外边界上是怎样传播的。图4-7-3 给出了动态演示中的一帧画面。
思考与练习
4.1 参考程序(ex4122) ,求解迭特里奇状态方程
p (v -b ) =RT e -a
RvT
的临界参量。(提示:可先求出v c 和T c ,最后代入方程求出p c )
4.2 关于用MATLAB 处理分形问题的实例,读者可通过网站“萝卜驿站(http://luobo.yculblog.com)进一步学习,该网站给出了不少有趣的例子。
4.3 利用jacobian 和det 函数,从理想气体分子麦氏速度分布导出速率分布。
4.3 试编写一个M 文件,求解下列问题:分子从器壁的小孔射出,求在射出的分子束中,分子的平均速率、方均根速率和平均能量。
4.3 试编写一个M 文件,求解下列问题: 一容积为1L 的容器,盛有温度为300 K、压强为3.0×104 Pa 的氩气,氩的摩尔质量为0.040 kg。若器壁上有一面积为1.0×10-3 cm 2的小孔,氩气将通过小孔从容器内逸出,经过多长时间容器里的原子数减少为原有原子数的?
4.4 自选题目并利用程序(ex4411) 计算理想气体在等值过程的热量、功和内能。
4.4 迭特里奇状态方程的形式为 1e
(p +a )(v -b ) =RT 。 v n
试利用程序(ex4421) 求等温过程中外界对系统所做的功、系统从外界吸收的热量和内能增量的表达式。
4.5 程序(ex4512) 给出了利用循环语句计算数值积分的方法,读者可仿照此方法处理统计物理中的类似问题。
4.6 试编写一个M 文件,求解下列问题:已知粒子遵从经典玻耳兹曼分布,其能量表达式为
e =1222(p x +p y +p Z ) +ax 2+bx 2m
其中,a , b 是常数,求粒子的平均能量。
4.7 试编写一个M 文件,求解下列问题:气柱的高度为H ,截面为S ,处在重力场中。试证明此气柱的内能和热容量为
NmgH , U =U 0+NkT -(emgH /kT -1)
N (mgH ) e mgH /kT 10。 C v =C v +Nk -mgH /kT 22(e-1) kT 2