热传导与热辐射大作业_动视

热传导与热辐射大作业

2025-09-25 14:08:41 责编:小OO

点击下载本文 文档为doc格式

热传导与热辐射大作业报告1作业问题 (3)

2问题求解 (4)

2.1第一部分求解分析解 (4)

2.1.1第一步问题的数学模型： (4)

2.1.2第二步问题分解 (4)

3计算结果 (9)

3.1温度计算结果 (9)

3.2热流密度计算结果 (9)

3.3等温线计算结果 (12)

3.4不同内热源时温度的计算结果 (15)

4误差的分析 (17)

5个人体会和小结 (18)

附录1 (19)

分析解计算程序 (19)

附录2 ............................................................................................................. 错误！未定义书签。

数值解计算程序 (20)

附录3 ............................................................................................................. 错误！未定义书签。

热流密度分析解程序 (22)

1作业问题

一矩形平板a x ≤≤0, b y ≤≤0，内有均匀恒定热源0g ，在0=x 及0=y 处绝热，在a x =及b y =处保持温度1T ，初始时刻温度为0T ，如右图1所示：

1、求0>t 时，矩形区域内的温度分布()t y x T ,,的解析表达式；

2、若m a 18=，m b 12=，301m W g =，K 6T 001=，K T 2000=，热传导系数 K m W k ⋅=0.1，热扩散系数2

0.8m s α=。请根据1中所求温度分布用

MATLAB 软件绘出下列结果，加以详细物理比较和分析：

(a) 300s 内，在同一图中画出点)4,0(、)8,0(、()0,6、)0,12(、)6,9(（单位：m ）温度随时间的变化曲线；

(b) 200s 内，画出点)4,18(、)8,18(、()12

,6、)12,12(、)6,9(（单位：m ）处热流密度值随时间的变化曲线；

(d) 300s 内，在同一图中画出点()0,9（单位：m ）在0g 分别等于31m W ，32m W ，

33m W 情况下的温度变化曲线；

3、运用有限差分法计算2中(b)和(d)，与解析结果进行比较，且需将数值解与解

析解之间相对误差减小到1‰； 4、附上源程序和个人体会；

以报告形式整理上述结果，用A4纸打印上交。

2问题求解

2.1第一部分求解分析解

分析：本问题是二维非齐次方程，非齐次边界条件问题。因此属于第四类问题，温度T （x,y,t ）可以分解为非齐次稳态的解Ts （x,y ）和齐次非稳态的解Th(x,y,t)的和。

2.1.1第一步问题的数学模型：

控制方程：2220

221g T T T x y k t

α∂∂∂++=∂∂∂

边界条件：1100,0;,;0,0;,;,0;

d T

x d x

T T x a d T

y d x

T T y b T T t ⎧==⎪⎪==⎪⎪⎪

==⎨

⎪==⎪⎪==⎪⎪⎩

2.1.2第二步问题分解

根据分析可以将问题分析为两部分，一个是非齐次稳态和一个齐次非稳态问题，如下

第一个非齐次稳态问题

控制方程：220220s s

T T g x y k

∂∂++=∂∂

边界条件：110,0;,;0,0;,;

s s s s d T x d x T T x a d T y d y T T y b ⎧==⎪⎪=

=⎪⎨

⎪==⎪⎪=

=⎩

设方程的解()()(),,,S

T x y p x y x y θ

=+，其中通过查表2-4可以得出特解()20,2g p x y x k =-

，则可以设方程的解()()2

0,,2S

g T x y x y x A k

θ=-+,将其带入控制方程：2222

0x y θθ

∂∂+=∂∂ 边界条件：2

012

010,0;,;20,0;,;2d x d x g T a A x a k

d y d y g T x A y b k

θθθθ⎧==⎪⎪⎪=-

-=⎪⎪

⎨

⎪==⎪⎪⎪=-

-=⎪⎩

可以看出方程的边界条件有两非齐次条件，可以令参数2

012g A T a k

=-,则方程的非齐次个数可以减少到一个，此时的边界条件变为：

()2

00,0;0,;0,0;,;

2d x d x x a d y d y g x a y

b k

θθθθ⎧==⎪⎪==⎪⎪⎨==⎪⎪⎪=-=⎪⎩

此问题转换为齐次方程一个非齐次边界条件问题，可以用分离变量法，设

()()()

,x y X x Yy θ=，将上式则方程可以分离为以下两个更简单的问题：

1、控制方程：222

0m X

X x

β∂+=∂ 边界条件：0,0;0,;

d X

x d x X x a ⎧==⎪

⎨⎪==⎩

2、控制方程：2220m Y

Y y

β∂-=∂

边界条件：

0,0;dY

y dy

== 其中方程1的解可以查表2-2，得()()()2

,c o s ;1/;m m m X x x N a

βββ==特征方程为：()c o s 0m a β=，所以可以得()c o s ()m

X x x β=，对于方程2，可以设其解为()()c o s h m Yy y β=。则方程的解设为()()()1

,c o s c o s h m m m

m x y C x y θββ∞

==∑，再利用边界条件y=b 时，则可以得到方程：()()()2

201

c o s c o s h 2m m m

m g x a C x b k ββ∞=-=∑，方程两边同时利用

()0

cos a

x dx

β⎰算子，利用三角函数的正交性可以求出

()()

()()22

001c o s c o s h 2a m m m m g C x a x d x N b k βββ=-⎰带入已知量可以得到()()03

2s i n c o s h m

m m m

g a C k a b βββ=-，所以可以求出

()()()

()()

2s i n ,c

o s c o s h c o s h m

m m

m g a x y x y k a b

βθββββ∞

==-∑，则非齐次稳态问题部分的解为：()()()()()()()0222200

113

12s i n ,c o s c o s h 22c o s h m s m m m m m

g a g g T x y a x Ta x T x y k k k a b βθ

ββ

∞

==+-+=-+-∑ 第二个齐次非稳态问题

设()()(),,,,,h s

T x y t T x y t T x y =-，并将其带入原问题的控制方程可以得

控制方程：2222

21h h h

T T T x y t

α∂∂∂+=∂∂∂ 边界条件：00,0;0,;0,0;0,;,0;

h h

h h s d T x d x

T x a d T y d y

T y b T T T t ⎧==⎪⎪==⎪⎪⎪

==⎨

⎪⎪==⎪=-=⎪⎪⎩

利用分离变量法设()()()h T X x Yy t =Γ，将其带入到上述控制方程和边界条件中可以得到三个分离方程。

1、控制方程：222

0n X X x

β∂+=∂ 边界条件：0,0;0,;

d X

x d x X x a ⎧==⎪

⎨⎪==⎩

2、控制方程：2220v Y

Y y

β∂+=∂

边界条件：0,0;0,;d Y

y d y Y y b ⎧==⎪

⎨⎪==⎩

3、时间变量方程：

20d dt

αλΓ

+Γ= 同理可以查表2-2得到方程1和方程2的解，()()()2,c o s ;1/;n n n

X x x N a

βββ==和()()()2,c o s ;1/;v v v

Y y y N b

βββ==方程3的解为()()22

n v t t e αββ-+Γ=，所以齐次非稳态问题的通解为()()()()()()2211

c o s c o s n v

t h n v n v

n v T X x Y yt C x y e α

ββββ∞∞

-+===Γ=∑∑ 利用初始条件，可以得到方程()()011

c o s c o s s n v

n v T T C x y ββ∞∞

==-=

∑∑

，利用正交算子()0

cos a

x dx

β⎰和

()0

cos b

y dy

β⎰，则可以求出

()()()20001

22224s i n s i n n v v

n v nv n n n v a b g g C T T a b k k βββββββββ⎡⎤⎛⎫⎢⎥=--+ ⎪+⎢⎥⎝⎭⎣⎦

，所以齐次非稳态的解可以

写

为：

()()()()()()22200012

22211s i n s i n 4c o s c o s n v t n v v

h n v n v n v n n n v a b g g T T T x y e a b k k αββββ

βββ

ββ

ββββ

∞∞-+==⎡⎤⎛⎫⎢⎥=--+ ⎪+⎢⎥⎝⎭⎣⎦

∑∑

将非齐次稳态解和齐次非稳态解结合起来，可以得出此问题的最终为： ()()(),,,,,s h

T x y t T x yT x y t =+，带入各个量可以得到最终的结果为 ()()()()()()()()()()()()2202

20131200012222

112s i n ,,c o s c o s h 2c o s h s i n s i n 4c o s c o s n v m m m

m m

m t n v v

n v n v nv n n n v g a g T x y t a x T x y k k a b a b g g TT x y e a b k k αββββββββββββββββββ∞

=∞∞

-+===-+-+⎡⎤⎛⎫⎢⎥--+ ⎪+⎢⎥⎝⎭⎣

⎦∑∑∑

到此问题的分析解就完成了，下面就是利用Matlab 对问题的具体描绘，以图形显示。

3计算结果

3.1温度计算结果

（a ）要求在同一幅图中给出（0，4）、（0，8）、（6，0）、（12，0）、（9，6）几个点在300秒内温度随时间变化的关系。附录1中给出了，分析解的程序，运行程序可以求出给点在300秒内的温度变化。取出T 中相应的值画图，下图给出了结果。

100

150200

250300

200250300350400450500550600650700时间t

温度︒C

300秒内不同点温度随时间变化曲线（分析解）

（0,4）（0,8）（6,0）（12,0）（9,6）

3.2热流密度计算结果

（b ）由傅里叶导热定律可知X 方向的热流密度为：x T

q k

∂=-∂，Y 方向的热流密度为y T

q k

∂=-∂，附录2中给出了求给点热流密度值的程序，运行程序可以得到给点不同时间的热流密度值，画出图形如下。数值解的程序：运行数值温度程序

Point=[9,6];%给定点的坐标

xx=point(1)/dx+1;%将给定点转换为网格点 yy=point(2)/dx+1;%将给定点转换为网格点 t=1:200;

zzs(t)=T(xx-1,yy,round(t./dt)); zzs2(t)=T(xx,yy,round(t./dt)); zzs4(t)=T(xx,yy-1,round(t./dt)); qx=k*(zzs2-zzs)/dx;%x 方向的热流密度值 qy=k*(zzs4-zzs2)/dx;%y 方向的热流密度值

20406080

100120140160180200

-300-250

-200

-150

-100

-500

时间t

x 方向热流密度值

200秒内不同点x 方向热流密度值（分析解）

（18,4）（18,8）（6,12）（12,12）（9,6）

20406080

100120140160180200

-300-250

-200

-150

-100

-500

时间t

y 方向的热流密度值

200秒内不同点y 方向热流密度值（分析解）

（18,4）（18,8）（6,12）（12,12）（9,6）

20406080

100120140160180200

-300-250

-200

-150

-100

-50

时间t

热流密度值

不同点在200秒内x 方向的热流密度值（数值解）

（9,6）（18,8）（18,4）（6,12）（12,12）

20406080

100120140160180200

-300-250

-200

-150

-100

-50

时间t

热流密度值

不同点在200秒内y 方向热流密度值随时间变化（数值解）

（12,12）（9,6）（6,12）（18,8）（18,4）

3.3等温线计算结果

（c ）不同时间的温度随时间变化的关系图，附录1给出了求各点温度的计算程序，运行程序，提取相应的温度值，化等温线如下。

程序命令： load Txf

dwx = Txf(:,:,50);%改变时间值可以画出不同时间的等温线 zhb=contour(dwx); clabel(zhb)

420

440

460

480

500

520

540

560

580

600

y 方向

x 方向

1224681012

1618

500

510

520

530

540

550

560

570

580

590

600

600600

600

600600y 方向

x 方向

75秒时物体内部等温线

246

810

12141618

560

565

570

575

580

585

590

595

600

600600

y 方向

x 方向

246810121416

598

599

600

600600

600

600600

600601

601602

602603

603

604

605

606

607

y 方向

x 方向

125秒时物体内部等温线

68101214

1618

600

600600

600600602604

606

608610612

614616

618

620

y 方向

x 方向

150秒时物体内部的等温线

1012141618

3.4不同内热源时温度的计算结果

（d ）300s 内不同内热源时，点（9,0）处的温度变化曲线如下。运行程序，改变其中的内热源可以画出不同内热源时的，给定点随温度的变化。

50100150

200250300350

200300

400

500

600

700

800

time

T e m p e r t a u r e

点（9，0）在不同内热源时的温度随时间变化图

g=1g=2g=3

050100150

200250300350

200

300

400

500

600

700

800

时间t

温度 C

（9,0）在不同内热源时温度随时间的变化（数值解）

4误差的分析

利用公式||

zz zzs w zz

计算数值解和分析解的相对误差，其中zz 表示不同内热源时分析解在（9,0）时不同时刻是的温度值，zzs 表示不同内热源时数值解在（9,0）时不同时刻是的温度值。画出误差图如下，可以看出满足千分之一的要求。

50100150

200250300350

012345678x 10

-4

时间

相对误差

g0不同时数值计算与分析解的相对误差

g0=1g0=2g0=3

虽然大作业做完了，但是它给我带来的“痛苦”远没有结束。通过这次的大作业使我学习到了很多的知识，一方面是matlab软件，以前有一点基础，因此编程还比较容易，但是往往是自己有优势的方面容易出问题。在解析解和数值的的精度上，我就吃了“；”的大亏，程序都写好了，但是分析解和数值解就是相差3%左右，就是达不到要求的千分之一，我为这个问题，苦苦想了好几天，一只进程都没有，最后，从新检查程序才发现循环中一个地方多了一个“；”，导致没有进入循环。我的天啊，看到这个小问题，我当时就要崩溃了，去掉了这个bug后，程序表现非常好。还有，通过这次大作业也激励起我学习的热情，原本没有学习的动力，但是，有了这个大作业，同学们都奋勇向前，努力的最先做出来，我也不甘落后，每天都苦思，虽然结果不一定是最好的，但是这其中的过的让我很享受。

该死的程序都附上了，还做了详细的说明很明了。我不想谈编程的感受了，编程从来就没有好受。我只是想说一下，数值解和分析解的理解。已经我总是认为，有一台电脑，就可以算任何的传热问题，只要遇到一点复杂的问题，首先就想到能不能用电脑去模拟一下，去算一下。通过这次的学习，我认识到有时候，数值解并不是最适合的，对于一些问题，如果可以从分析角度去解，最好是通过分析解好。一方面分析解有很好的物理意义，可以让我们更进一步的了解问题的本质，了解解决问题的方法，简化问题的方法；另一方面，有时分析解的运算量有时就比数值解小多了。这点，我深有体会。有时，分析解输入电脑，一会就出结果，但是数值模拟可能要运行很长的时间，并且随着精度的增加有时电脑的内存都可能吃不消。我运行这个大作业的程序产生了5G的数据量，多得吓死人。

好了，不说了废话。最后希望老师能将这样的活动继续办下去，同时增加点趣味，不要以为的用课本上的题目，可以有点创新。还有可以分组，大家多交流，能更好的到达学习的目的。

附录

分析解计算程序

clear all; %清除系统中原有的变量

clc; %清除屏幕

a=18; %x方向的长度

b=12; %y方向的长度

g0=3; %内热源

T1=600; %边界条件温度

T0=200; %初始温度

k=1.0; %导热系数

af=0.8; %热扩散系数

kk=15; %循环变量，取15时误差小于0.2*10-4，足够精确

for i=1:19 %matlab循环从1开始

x=i-1; %x的长度

for j=1:13

y=j-1;

Tsp=0;%初始化计算变量

for m=1:kk

bm=(2*m-1)/(2*a)*pi; %非齐次稳态问题的特征值

Tsp=Tsp+2*g0*sin(bm*a)/(k*a*bm^3*cosh(bm*b))*cosh(bm*y)*cos(bm*x) ;

end

Ts = T1+g0/(2*k)*(a^2-x^2) - Tsp;%非齐次稳态问题的解

for t = 1:350 %显示的时间范围

Th = 0;

for p=1:kk

for q =1:kk

rp=(2*p-1)/(2*a)*pi;%齐次非稳态问题的x方向的特征值

kq=(2*q-1)/(2*b)*pi;%齐次非稳态问题的y方向的特征值

Th=Th+4/(a*b)*sin(rp*a)*sin(kq*b)/(rp*kq)*((T0-T1-g0/(k*rp^2))+g0 *kq^2/(k*(kq^2+rp^2)*rp^2))*cos(rp*x)*cos(kq*y)*exp(-af*(rp^2+kq^ 2)*t);%齐次非稳态的温度解

end

Txf(i,j,t) = Ts + Th;%最终的温度解

end

endsave Txf;

数值解计算程序

clear all;%清除系统内存里的变量

clc;%清除屏幕

point=[9,0];%给定点的坐标

g = 1; %内热源

sj=350; %给定时间为350秒

M=55; %x方向的网格数量

N=2/3*(M-1)+1; %保证x和y方向的网格长度一样

dx=18/(M-1); %网格的长度

dt=(dx)^2/6; %保证Fo数的值小于0.25

P=ceil(sj/dt); %时间步长

T=zeros(M,N,P);%初始化温度矩阵，一维表示x方向，二维表示y方向，三维表示时间变量

T0=200;%初始条件温度值

T1=600;%边界条件温度值

a=0.8;%热扩散系数

k=1;%热传导系数

Fo=a*dt/(dx)^2;%网格傅里叶数

for i=1:M %给定初始化条件，matlab中循环从1开始

for j=1:N

T(i,j,1)=T0;%所有点在初始条件是的温度为T0

end

for kk=2:P %给定边界条件

for i=1:N

T(M,i,kk)=T1;%（x=18的边界温度恒定）

end

for j=1:M

T(j,N,kk)=T1;%（y=12的边界温度恒定）

end

for p=1:P

for i=1:M-1

if i==1

for j=1:N-1

if j==1

T(i,j,p+1)=2*Fo*(T(i+1,j,p)+T(i,j+1,p))+(1-4*Fo)*T(i,j,p)+a*dt/k* g;%x=0，y=0的处温度的计算式

elseT(i,j,p+1)=Fo*(T(i,j+1,p)+T(i,j-1,p)+2*T(i+1,j,p))+(1-4*Fo)*T(i,j ,p)+a*dt/k*g;%x=0的边界点温度计算式

end

else

for j=1:N-1

if j==1

T(i,j,p+1)=Fo*(T(i+1,j,p)+T(i-1,j,p)+2*T(i,j+1,p))+(1-4*Fo)*T(i,j ,p)+a*dt/k*g;%y=0的边界温度计算式

else

T(i,j,p+1)=Fo*(T(i,j+1,p)+T(i,j-1,p)+T(i+1,j,p)+T(i-1,j,p))+(1-4* Fo)*T(i,j,p)+a*dt/k*g;%中间点的计算公式

end

xx=point(1)/dx+1;%将给定点转换为网格点

yy=point(2)/dx+1;%将给定点转换为网格点

save T热流密度分析解程序

clear all;%清除内存内变量

clc;%清屏

a=18;

b=12;

g0=1;

T1=600;

T0=200;

k=1.0;%热传导系数

af=0.8;%热扩散系数

kk=15;%循环变量

for i=1:19

x=i-1;

for j=1:13

y=j-1;

qxs=0;%初始化热流量中间变量

for m=1:kk

bm=(2*m-1)/(2*a)*pi;

qxs=qxs+2*g0*sin(bm*a)/(a*bm^2*cosh(bm*b))*cosh(bm*y)*sin(bm*x); end

qx1= g0*x - qxs;%稳态热流量

for t = 1:200

qx2 = 0;

for p=1:kk

for q =1:kk

rp=(2*p-1)/(2*a)*pi;

kq=(2*q-1)/(2*b)*pi;

qx2=qx2+4*k*rp/(a*b)*sin(rp*a)*sin(kq*b)/(rp*kq)*((T0-T1-g0/(k*rp ^2))+g0*kq^2/(k*(kq^2+rp^2)*rp^2))*sin(rp*x)*cos(kq*y)*exp(-af*(r p^2+kq^2)*t);

end

qx(i,j,t) = qx1 + qx2;%x方向的热流密度矩阵

end

endend

for i=1:19

x=i-1;

for j=1:13

y=j-1;

qy1=0;

for m=1:kk

bm=(2*m-1)/(2*a)*pi;

qy1=qy1+2*g0*sin(bm*a)/(a*bm^2*cosh(bm*b))*sinh(bm*y)*cos(bm*x); end

for t = 1:200

qy2 = 0;

for p=1:kk

for q =1:kk

rp=(2*p-1)/(2*a)*pi;

kq=(2*q-1)/(2*b)*pi;

qy2=qy2+4*k*kq/(a*b)*sin(rp*a)*sin(kq*b)/(rp*kq)*((T0-T1-g0/(k*rp ^2))+g0*kq^2/(k*(kq^2+rp^2)*rp^2))*cos(rp*x)*sin(kq*y)*exp(-af*(r p^2+kq^2)*t);

end

qy(i,j,t) = qy1 + qy2;%y方向的热流密度矩阵

end

save qx

save qy下载本文

显示全文

全部频道