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

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

2、 / m2C砖墙的导热系数0.53W / m Ch2tw2：数学描述该结构的导热问题可以作为二维问题处理，并且其截面如图 1-1 所示，由于 对称性，仅研究其 1/4 部分即可。其网络节点划分如图b上述问题为二维矩形域内的稳态、 无内热源、 常物性的导热问题， 对于这样的 物理问题，我们知道，描写其的微分方程即控制方程，就是导热微分方程：2tx2第一类边界条件：内外壁分布均匀地维持在 0 C 及 30 C ；t w1 =30 Ctw2 =0 C砖墙的导热系数第三类边界条件：内外表面均为第三类边界条件，且已知：t130 C, h110 W / m2 Ct210 C, h24 W / m2 C0.

3、53W / m C三：方程的离散 如上图所示，用一系列与坐标轴平行的网络线把求解区域划分成许多子区 域，以网格线的交点作为需要确定温度值的空间位置， 即节点， 节点的位置已该 点在两个方向上的标号 m 、n 来表示。每一个节点都可以看成是以它为中心的小 区域的代表，如上（ m ，n ）：对于（ m ，n ）为内节点时：由热平衡法可以得到，当 x = y 时：t 1 (t t t t )tm,n (t m 1,n tm 1,n tm,n 1 t m,n 1) 4 对于（ m ，n）为边界节点时：恒温边界只需特殊考虑位于绝热平直边界上的节点：t m,n( t m 1,n4t m ,n 1对流边界分

4、为角点、绝热边界点和对流边界点11.绝热边界点： tm,n( tm 1,n42. 对流边界点：hx22 t m,nt m, n 12tm 1,ntm,n 1 )2h x2t tm 1, nm, n 1tf3. 外角点：hxm, n2hxt m 1 , nt m, n 1tf4. 内角点：hx3tm,n2 t m, n 1t m 1, n2h xt m, n 1t m 1,ntf四：编程思路及流程图开始输入已知参数说明边界条件取定初始试探值TA(i, j)=0计算新的内节点和边界点温度T(i, j)T=TA计算内外边界散热量及热平衡偏差输出内外边界散热量及热平衡偏差MATLAB 程序如下 fun

5、ction =wendu() t=zeros(12,16); tf=zeros(12,16);Q1x=0;Q1y=0;Q1=0;Q2x=0;Q2y=0;Q2=0; n=0;for i=1:12t(i,1)=30;endt(1,j)=30;endt0=t;t=diedai1(t);for i=2:5for j=2:16while t(i,j)-t0(i,j)>=0.00001;t0=t;t=diedai1(t);endendendfor i=6:12for j=2:5while t(i,j)-t0(i,j)>=0.00001;t0=t;t=diedai1(t);endendendt0

6、=tfor i=2:11Q1x=Q1x+0.53*(t(i,1)-t(i,2); end Q1x=Q1x+0.53*(t(12,1)-t(12,2)/2; for j=2:15Q1y=Q1y+0.53*(t(1,j)-t(2,j); end Q1y=Q1y+0.53*(t(1,16)-t(2,16)/2; Q1=(Q1x+Q1y)*4 for i=6:11Q2x=Q2x+0.53*(t(i,5)-t(i,6); end Q2x=Q2x+0.53*(t(12,5)-t(12,6)/2; for j=6:15Q2y=Q2y+0.53*(t(5,j)-t(6,j); end Q2y=Q2y+0.53

7、*(t(5,16)-t(6,16)/2; Q2=(Q2x+Q2y)*4 n=2*abs(Q1-Q2)/(Q1+Q2)t0=tf;tf=diedai2(tf);for i=2:5for j=2:16while tf(i,j)-t0(i,j)>=0.00001;t0=tf;tf=diedai2(tf);endendendfor i=6:12for j=2:5while tf(i,j)-t0(i,j)>=0.00001;t0=tf;tf=diedai2(tf);endendendt0=tfQ1x=0;Q1y=0;Q1=0;Q2x=0;Q2y=0;Q2=0;n=0;for i=1:11Q1

8、x=Q1x+10*0.1*(30-tf(i,1);endQ1x=Q1x+10*0.05*(30-tf(12,1);for j=2:15Q1y=Q1y+10*0.1*(30-tf(1,j);endQ1y=Q1y+10*0.05*(30-tf(1,16);for i=6:11Q2x=Q2x+4*0.1*(tf(i,6)-10); end Q2x=Q2x+4*0.05*(tf(12,6)-10); for j=7:15Q2y=Q2y+4*0.1*(tf(6,j)-10); end Q2y=Q2y+4*0.05*(tf(6,16)-10); Q1=(Q1x+Q1y)*4; Q2=(Q2x+Q2y)*4

9、; n=2*abs(Q1-Q2)/(Q1+Q2)for i=7:12tf(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)colorbarfunction t1=diedai1(t)for j=2:1

10、5t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4;endt(i,16)=(2*t(i,15)+t(i-1,16)+t(i+1,16)/4;endfor i=6:11for j=2:5t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4;endendfor j=2:5t(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)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53

11、+1);for j=2:15t(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); end t(1,16)=(2*t(2,16)+2*t(1,15)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);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); end t(12,1)=(2*t(12,2)+2*t(11,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);

12、for i=2:5for j=2:15t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4;endt(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:15t(6,j)=(2*t(5,j)+t(6,j-1)+t(6,j+1)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+2); end t(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:11t(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); end t(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:11for j=2:5t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4;endend for j=2:5t(12,j)=(