二维导热物体温度场的数值模拟

1、二维导热物体温度场的数值模拟作者： 学 号： 学院(系)：能源与动力工程学院专业：能源动力系统及自动化班级： 二维导热物体温度场的数值模拟 一：物理问题 有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸和示意图如图1-1所示，假设在垂直纸面方向上冷空气及砖墙的温度变化很小，可以近似地予以忽略。在下列两种情况下试计算：（1）砖墙横截面上的温度分布；（2）垂直于纸面方向的每米长度上通过砖墙的导热量。 第一种情况：内外壁分布均匀地维持在0及30； 第二种情况：内外表面均为第三类边界条件，且已知： 砖墙的导热系数 二：数学描述 该结构的导热问题可以作为二维问题处理，并且其截面如图1-1所示，由于对称

2、性，仅研究其1/4部分即可。 其网络节点划分如图 f a c （m，n） b = n e m d 上述问题为二维矩形域内的稳态、无内热源、常物性的导热问题，对于这样的物理问题，我们知道，描写其的微分方程即控制方程，就是导热微分方程： 第一类边界条件：内外壁分布均匀地维持在0及30；=30=0 第三类边界条件：内外表面均为第三类边界条件，且已知： 砖墙的导热系数三：方程的离散如上图所示，用一系列与坐标轴平行的网络线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，即节点，节点的位置已该点在两个方向上的标号m、n来表示。每一个节点都可以看成是以它为中心的小区域的代表，如上（m

3、，n）：对于（m，n）为内节点时：由热平衡法可以得到，当=时： 1 对于（m，n）为边界节点时：l 恒温边界只需特殊考虑位于绝热平直边界上的节点：l 对流边界分为角点、绝热边界点和对流边界点。1.绝热边界点：2.对流边界点：3.外角点：4.内角点：四：编程思路及流程图 开始 输入已知参数说明边界条件取定初始试探值ta（i, j）=0计算新的内节点和边界点温度t(i, j)t=ta比较所有节点|ta(i, j)- t(i, j)|<epsilo ？否是计算内外边界散热量及热平衡偏差输出内外边界散热量及热平衡偏差结束画出温度场模拟图 matlab程序如下function=wendu()t=z

4、eros(12,16);tf=zeros(12,16);q1x=0;q1y=0;q1=0;q2x=0;q2y=0;q2=0;n=0;for i=1:12 t(i,1)=30;endfor j=1:16 t(1,j)=30;endt0=t;t=diedai1(t);for i=2:5 for j=2:16 while t(i,j)-t0(i,j)>=0.00001; t0=t; t=diedai1(t); end endendfor i=6:12 for j=2:5 while t(i,j)-t0(i,j)>=0.00001; t0=t; t=diedai1(t); end ende

5、ndt0=t for i=2:11 q1x=q1x+0.53*(t(i,1)-t(i,2);endq1x=q1x+0.53*(t(12,1)-t(12,2)/2;for j=2:15 q1y=q1y+0.53*(t(1,j)-t(2,j);endq1y=q1y+0.53*(t(1,16)-t(2,16)/2;q1=(q1x+q1y)*4for i=6:11 q2x=q2x+0.53*(t(i,5)-t(i,6);endq2x=q2x+0.53*(t(12,5)-t(12,6)/2;for j=6:15 q2y=q2y+0.53*(t(5,j)-t(6,j);endq2y=q2y+0.53*(t

6、(5,16)-t(6,16)/2;q2=(q2x+q2y)*4n=2*abs(q1-q2)/(q1+q2) t0=tf;tf=diedai2(tf);for i=2:5 for j=2:16 while tf(i,j)-t0(i,j)>=0.00001; t0=tf; tf=diedai2(tf); end endendfor i=6:12 for j=2:5 while tf(i,j)-t0(i,j)>=0.00001; t0=tf; tf=diedai2(tf); end endendt0=tf q1x=0;q1y=0;q1=0;q2x=0;q2y=0;q2=0;n=0;for

7、 i=1:11 q1x=q1x+10*0.1*(30-tf(i,1);endq1x=q1x+10*0.05*(30-tf(12,1);for j=2:15 q1y=q1y+10*0.1*(30-tf(1,j);endq1y=q1y+10*0.05*(30-tf(1,16);for i=6:11 q2x=q2x+4*0.1*(tf(i,6)-10);endq2x=q2x+4*0.05*(tf(12,6)-10);for j=7:15 q2y=q2y+4*0.1*(tf(6,j)-10);endq2y=q2y+4*0.05*(tf(6,16)-10);q1=(q1x+q1y)*4;q2=(q2x+

8、q2y)*4;n=2*abs(q1-q2)/(q1+q2) for i=7:12 for j=7:16 tf(i,j)=10; endendsubplot(211);pcolor(t)shading interp;colormap(hot)hold oncontour(t,3,'k')hold offcolorbark=caxis;subplot(212);pcolor(tf)shading interp;colormap(hot)hold oncontour(tf,3,'k')hold offcaxis(k)colorbar function t1=dieda

9、i1(t) for i=2:5 for j=2:15 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; end t(i,16)=(2*t(i,15)+t(i-1,16)+t(i+1,16)/4;endfor i=6:11 for j=2:5 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; endendfor j=2:5 t(12,j)=(2*t(11,j)+t(12,j-1)+t(12,j+1)/4;endt1=t; function t1=diedai2(t)t(1,1)=(t(1,2)+t(2,1)+

10、2*10*0.1*30/0.53)/(2*(10*0.1/0.53+1);for j=2:15 t(1,j)=(2*t(2,j)+t(1,j-1)+t(1,j+1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);endt(1,16)=(2*t(2,16)+2*t(1,15)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);for i=2:11 t(i,1)=(2*t(i,2)+t(i-1,1)+t(i+1,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);endt(12,1)=(2*t(12,2)+2*t(11

11、,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);for i=2:5 for j=2:15 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; end t(i,16)=(2*t(i,15)+t(i-1,16)+t(i+1,16)/4;endt(6,6)=(2*(t(6,5)+t(5,6)+t(7,6)+t(6,7)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+3);for j=7:15 t(6,j)=(2*t(5,j)+t(6,j-1)+t(6,j+1)+2*4*0.1*10/0.53)/(2*(4*

12、0.1/0.53+2);endt(6,16)=(2*t(5,16)+2*t(6,15)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+2);for i=7:11 t(i,6)=(2*t(i,5)+t(i-1,6)+t(i+1,6)+2*10*0.1*4/0.53)/(2*(4*0.1/0.53+2);endt(12,6)=(2*t(12,5)+2*t(11,6)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+2);for i=6:11 for j=2:5 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; endendfor j=2:5 t(12,j)=(2*t(11