传热学课程数值计算实验报告

传热学课程数值计算实验报告 一问题重述有一个用砖砌成的长方形截面的冷气通道，其截面尺寸如下图所示，假设在垂直于纸面方向上冷空气及砖墙的温度变化很小，可以近似忽略。试计算：砖墙横截面上的温度分布；垂直于纸面方向上的每米长度上通过砖墙的导热量。墙壁内、外表面的流体的温度分别为10、30；内外表面均为第三类边界条件，且已知：=30，=10.33=10,=3.93砖墙的导热系数二离散首先考虑到整个墙体截面的对称性，对称地取其四分之一截面进行研究，可以适当简化问题。为方便与温度场电模拟实验的数据进行对比，对温度场的离散与之保持一致，即网格为0.1m的正方形，且对内、外墙表面的流体温度的模拟各设一排节点进行模拟。在实际的编程过程中，为使程序简洁且更有条理，建立一个网格矩阵，矩阵中每个元素的值根据温度场中对应位置的节点的边界情况或计算特点进行设定从而将温度场节点进行分类。网格矩阵如下。下面会结合该矩阵详细阐述迭代方程的建立过程。三方程的建立与求解在介绍迭代方程之前，需要说明的是：a、为提高迭代速度，采用Gauss-Seidel迭代方法，即总是将最新得出来的节点数据用到迭代过程中去；b、每次迭代的初始温度场记为tfi，正在迭代的新温度场记为tft；c、墙壁的内节点初值设为20。网格矩阵将温度场的节点按照边界情况或计算条件进行了分类、标记。每个不同数值含义及对应的迭代方程列举如下：0：温度场中流体的节点，迭代过程中值保持不变：1：墙壁的内节点，控制方程可有热平衡法或泰勒级数展开得出：2：外墙面边界点（角点除外），由热平衡法整理可有：3：外墙边界点（角点除外），与2相似，可由热平衡法整理得出：4：墙壁对称面，绝热边界条件：5：内边界点（角点除外），可由热平衡法整理得出:6：内墙边界点（角点除外），与5相似，由热平衡法整理得出：7:墙壁对称面，绝热边界条件：8：内墙壁角点，由于其形状位置的特殊性，方程需要单独考虑，仍然由热平衡法整理得出：9：外墙壁角点，由于其形状位置的特殊性，方程需要单独考虑，由热平衡法整理可得：值得指出的是，在实际的迭代程序编写中，对于绝热边界条件即4和7点的计算是要放在整个循环程序之外的，否则结果会出错！对于对流传热量的计算，仍然只在该四分之一墙角内计算。基本思路是，分别对内、外墙面累计加和壁面节点与周围流体的温差，而后再按照牛顿冷却公式计算出换热量。值得注意的是，对于两个绝热边界上的四个角点，它们所在的控制体积的对流换热面积是其它边界点的一半。最后可以按照实验书上的平衡误差算法计算内外壁传热的误差。下面对计算结果列举如下：在迭代精度为1e-6的条件下，最终的温度场矩阵为：迭代次数为：113次；垂直于纸面方向上的每米长度上通过砖墙内、外壁面的导热量分别为：、，平衡误差为四等温线图及温度场三维图五源程序程序用matlab语言编写如下：clear;clc;clf;h1=10.33;h2=3.93;lamda=0.53;d=0.1;%todefinesomekeyparameterstfi=ones(13,17)*20;tfi(1,:)=ones(1,17)*30;tfi(:,1)=ones(13,1)*30;tfi(8:13,8:17)=ones(6,10)*10%theinitialtemperaturefieldoftheiterationmesh1=ones(13,17);mesh1(1,:)=zeros(1,17);mesh1(:,1)=zeros(13,1);mesh1(8:13,8:17)=zeros(6,10);mesh1(3:12,2)=ones(10,1)*2;mesh1(2,3:17)=ones(1,15)*3;mesh1(2:7,17)=ones(6,1)*4;mesh1(7,8:16)=ones(1,9)*5;mesh1(8:12,7)=ones(5,1)*6;mesh1(13,2:7)=ones(1,6)*7;mesh1(7,7)=8;mesh1(2,2)=9%todefineameshmatrixofthesameshapewiththetemperaturefieldmatrixandmarkeachnodeaccordingtothespecificconditionofthecorrespondingoneintheTFMk=0;maxode=1.0;tft=tfi;while(maxode>1e-6)%precisionofiterationis1e-6maxode=0;k=k+1;fori=2:12forj=2:16m=mesh1(i,j);n=tfi(i,j);switchmcase0tft(i,j)=tfi(i,j);case1tft(i,j)=0.25*(tfi(i,j+1)+tfi(i+1,j)+tft(i-1,j)+tft(i,j-1));case2tft(i,j)=(h1*d*tft(i,j-1)+0.5*lamda*(tft(i-1,j)+tfi(i+1,j)+2*tfi(i,j+1)))/(2*lamda+h1*d);case3tft(i,j)=(h1*d*tft(i-1,j)+0.5*lamda*(tft(i,j-1)+tfi(i,j+1)+2*tfi(i+1,j)))/(2*lamda+h1*d);%case4%tft(i,j)=tfi(i,j-1);erroroccurswhencalculateadiabaticboundaryconditionintheiterationprocesscase5tft(i,j)=(h2*d*tfi(i+1,j)+0.5*lamda*(tft(i,j-1)+tfi(i,j+1)+2*tft(i-1,j)))/(2*lamda+h2*d);case6tft(i,j)=(h2*d*tfi(i,j+1)+0.5*lamda*(tft(i-1,j)+tfi(i+1,j)+2*tft(i,j-1)))/(2*lamda+h2*d);%case7%tft(i,j)=tfi(i-1,j);case8tft(7,7)=(tfi(8,8)*h2*d+0.5*lamda*(2*tft(6,7)+2*tft(7,6)+tfi(7,8)+tfi(8,7)))/(h2*d+3*lamda);otherwisetft(2,2)=(h1*d*tft(1,1)+0.5*lamda*(tfi(2,3)+tfi(3,2)))/(lamda+h1*d);endde=abs(tft(i,j)-n);if(de>maxode)maxode=de;endendendtft(2:7,17)=tft(2:7,16);tft(13,2:7)=tft(12,2:7);%thisisthecorrectpositiontocalculatethetemperatureofadiabaticboundarytfi=tft;endk%returnthetimesofiterationtf=tft%returnthetemperaturefieldsubplot(1,2,1)%visualizethetemperaturefieldmesh(tf),meshc(tf)title('temperaturefieldofaquarterofthewall','fontsize',12)zlabel('temperature/Celsiusdegree','fontsize',10)ylabel('y/0.1m','fontsize',10)xlabel('x/0.1m','fontsize',10)colorbarsubplot(1,2,2)%todelineatethethreedesirableisothermalcurvesv=[22,26,28];contour(tf,v)axis([018014],'equal')axisonholdongtext('22Celsiusdegree')gtext('26Celsiusdegree')gtext('28Celsiusdegree')x1=[2,17];y1=[2,2];x2=[2,2];y2=[2,13];x3=[7,17];y3=[7,7];x4=[7,7];y4=[7,13];y5=[13,13];x5=[2,7];y6=[2,7];x6=[17,17];plot(x1,y1,'-k',x2,y2,'-k',x3,y3,'-k',x4,y4,'-k',x5,y5,'-.k',x6,y6,'-.k')gtext('outerwallsurface')gtext('innerwallsurface')gtext('symmetryplane')gtext('symmetryplane')ylabel('y/0.1m','fontsize',10)xlabel('x/0.1m','fontsize',10)title('threeisothermalcurves','fontsize',12)atdoiw=0;atdoow=0;fori=1:12forj=1:17m=mesh1(i,j);if((m==2)||(m==3)||(m==9))%payattentiontothefixedexpression!atdoow=atdoow+30-tf(i,j);elseif((m==6)||(m==8)||(m==5))atdoiw=atdoiw+tf(i,j)-10;endendendend%toadduptemperaturedifferencesofbothouterwallandinnerwallatdoow=atdoow+0.5*(60-tf(2,17)-tf(13,2));atdoiw=atdoiw+0.5*(tf(13,7)+tf(7,17)-20);heatfluxo=atdoow*h1*d%heatfluxthroughtheouter