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

二维导热物体温度场的数值模拟作者：学号：学院(系)：能源与动力工程学院专业：能源动力系统及自动化班级：二维导热物体温度场的数值模拟一：物理问题有一个用砖砌成的长方形截面的冷空气通道，其截面尺寸和示意图如图1-1所示，假设在垂直纸面方向上冷空气及砖墙的温度变化很小，可以近似地予以忽略。在下列两种情况下试计算：（1）砖墙横截面上的温度分布；（2）垂直于纸面方向的每米长度上通过砖墙的导热量。第一种情况：内外壁分布均匀地维持在0C及30C；第二种情况：内外表面均为第三类边界条件，且已知：10,3011hCtCmW2/4,1022hCtCmW2/砖墙的导热系数CmW/3.5011ht，1wt22ht2wt二：数学描述该结构的导热问题可以作为二维问题处理，并且其截面如图1-1所示，由于对称性，仅研究其1/4部分即可。其网络节点划分如图fac（m，n）bxx=ynyemd上述问题为二维矩形域内的稳态、无内热源、常物性的导热问题，对于这样的物理问题，我们知道，描写其的微分方程即控制方程，就是导热微分方程：02222ytxt第一类边界条件：内外壁分布均匀地维持在0C及30C；1wt=30C2wt=0C第三类边界条件：内外表面均为第三类边界条件，且已知：10,3011hCtCmW2/4,1022hCtCmW2/砖墙的导热系数CmW/3.50三：方程的离散如上图所示，用一系列与坐标轴平行的网络线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，即节点，节点的位置已该点在两个方向上的标号m、n来表示。每一个节点都可以看成是以它为中心的小区域的代表，如上（m，n）：对于（m，n）为内节点时：由热平衡法可以得到，当x=y时：)tttt(41t1,1,,1,1,nmnmnmnmnm①对于（m，n）为边界节点时：恒温边界只需特殊考虑位于绝热平直边界上的节点：)tt2t(41t1,,1,1,nmnmnmnm对流边界分为角点、绝热边界点和对流边界点。1.绝热边界点：)tt2t(41t1,,1,1,nmnmnmnm2.对流边界点：f1,,11,,tx2htt2tt2xh2nmnmnmnm3.外角点：f1,,1,tx2httt1xh2nmnmnm4.内角点：f,11,,11,,tx2htttt2t3xh2nmnmnmnmnm四：编程思路及流程图开始输入已知参数说明边界条件取定初始试探值TA（i,j）=0计算新的内节点和边界点温度T(i,j)比较所有节点|TA(i,j)-T(i,j)|EPSILO？计算内外边界散热量及热平衡偏差输出内外边界散热量及热平衡偏差T=TA是否画出温度场模拟图结束MATLAB程序如下function[]=wendu()t=zeros(12,16);tf=zeros(12,16);Q1x=0;Q1y=0;Q1=0;Q2x=0;Q2y=0;Q2=0;n=0;fori=1:12t(i,1)=30;endforj=1:16t(1,j)=30;endt0=t;t=diedai1(t);fori=2:5forj=2:16whilet(i,j)-t0(i,j)=0.00001;t0=t;t=diedai1(t);endendendfori=6:12forj=2:5whilet(i,j)-t0(i,j)=0.00001;t0=t;t=diedai1(t);endendendt0=tfori=2:11Q1x=Q1x+0.53*(t(i,1)-t(i,2));endQ1x=Q1x+0.53*(t(12,1)-t(12,2))/2;forj=2:15Q1y=Q1y+0.53*(t(1,j)-t(2,j));endQ1y=Q1y+0.53*(t(1,16)-t(2,16))/2;Q1=(Q1x+Q1y)*4fori=6:11Q2x=Q2x+0.53*(t(i,5)-t(i,6));endQ2x=Q2x+0.53*(t(12,5)-t(12,6))/2;forj=6:15Q2y=Q2y+0.53*(t(5,j)-t(6,j));endQ2y=Q2y+0.53*(t(5,16)-t(6,16))/2;Q2=(Q2x+Q2y)*4n=2*abs(Q1-Q2)/(Q1+Q2)t0=tf;tf=diedai2(tf);fori=2:5forj=2:16whiletf(i,j)-t0(i,j)=0.00001;t0=tf;tf=diedai2(tf);endendendfori=6:12forj=2:5whiletf(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;fori=1:11Q1x=Q1x+10*0.1*(30-tf(i,1));endQ1x=Q1x+10*0.05*(30-tf(12,1));forj=2:15Q1y=Q1y+10*0.1*(30-tf(1,j));endQ1y=Q1y+10*0.05*(30-tf(1,16));fori=6:11Q2x=Q2x+4*0.1*(tf(i,6)-10);endQ2x=Q2x+4*0.05*(tf(12,6)-10);forj=7:15Q2y=Q2y+4*0.1*(tf(6,j)-10);endQ2y=Q2y+4*0.05*(tf(6,16)-10);Q1=(Q1x+Q1y)*4;Q2=(Q2x+Q2y)*4;n=2*abs(Q1-Q2)/(Q1+Q2)fori=7:12forj=7:16tf(i,j)=10;endendsubplot(211);pcolor(t)shadinginterp;colormap(hot)holdoncontour(t,3,'k')holdoffcolorbark=caxis;subplot(212);pcolor(tf)shadinginterp;colormap(hot)holdoncontour(tf,3,'k')holdoffcaxis(k)colorbarfunctiont1=diedai1(t)fori=2:5forj=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;endfori=6:11forj=2:5t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j))/4;endendforj=2:5t(12,j)=(2*t(11,j)+t(12,j-1)+t(12,j+1))/4;endt1=t;functiont1=diedai2(t)t(1,1)=(t(1,2)+t(2,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+1));forj=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));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));fori=2:11t(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,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2));fori=2:5forj=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));forj=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));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));fori=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));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));fori=6:11forj=2:5t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j))/4;endendforj=2:5t(12,j)=(2*t(11,j)+t(12,j-1)+t(12,j+1))/4;endt1=t;输出结果如下：恒温边界误差对流边界两个温度场如下五、问题讨论与电模拟实验相比，数值模拟的误差更小，精度更高。但两种方法都在误差允许的范围之内，都可以使用。本次实验让我对传热学的概念有了更深刻的理解，同时对于传热学数值研究有了深刻的体会，对于迭代的方法能够熟练掌握，编程的思路更加连贯，总的来说，受益匪浅。