传热学大作业

传热学第四章大作业二维稳态导热问题的数值解法姓名：班级：学号：第一题：如图所示，一个无限长矩形柱体，其横截面的边长分别为和，常物性。该问题可视为二维稳态导热问题，边界条件如图中所示，其中=0.6m，=0.4m，=60℃，=20℃，。编写程序求解二维导热方程。绘制x=/2和y=/2处的温度场，并与解析解进行比较。已知矩形内的温度场的解析解为：解：(1)建立控制方程及定解条件控制方程：定解条件：（2）区域离散化（确立节点）将矩形区域分为M*N个网格，其中x方向上的步长;y方向的步长。设节点为（m,n）。（3）建立节点离散方程对节点（m,n）有：内节点:化简得边界节点：（4）编程求解，程序见附录取M=N=50得到矩形区域各节点温度（见附件一），为方便在这里仅给出M=N=10时温度分布数据，如下表：606060606060606060606060.4486560.8287961.0919261.213461.1883861.0274160.7531460.39773606060.9277461.7095862.2468162.491962.4371862.1051761.5423660.81427606061.4716862.6990863.5315463.9036663.8092363.285562.4050261.26912606062.1241863.8641265.0208665.5227665.3720864.6243963.381261.78318606062.9474565.2883666.8012267.4319167.2002766.1835164.5150962.37953606064.040367.0806368.9736469.7243169.377468.0322465.856663.08431606065.5775969.3828171.6537872.5045271.998270.2500867.4632863.92776606067.9025872.3668974.9651375.8878875.1709172.9295569.4027464.945536066.1803471.7557176.1803479.021138079.0211376.1803471.7557166.1803460M=N=50时矩形区域节点温度分布图如图：M=N=50时x=/2处数值解与解析解温度场分布如下：M=N=10时将数值解与理论解进行比较，如下表所示：y00.4/90.4/9*20.4/9*30.4/9*40.4/9*50.4/9*60.4/9*70.4/9*80.4数值解6061.213462.491963.9036665.5227667.4319169.7243172.5045275.8878880理论解6061.1744762.4128263.7824465.3578167.2246469.4844972.2602975.7030380相对误差00.064%0.127%0.190%0.252%0.308%0.345%0.338%0.244%0M=N=50时y=/2处数值解与解析解温度场分布如下图：M=N=10时将数值解与理论解进行比较，如下表所示：x00.6/90.6/9*20.6/9*30.6/9*40.6/9*50.6/9*60.6/9*70.6/9*80.6数值解6062.1241863.8641265.0208665.5227665.3720864.6243963.381261.7831860理论解6062.1372464.016765.4116966.1539566.1539565.4116964.016762.1372460相对误差00.021%0.238%0.598%0.954%1.182%1.204%0.993%0.569%0第二题：将第一题中处的边界条件变为，其他条件不变。编写程序求解二维导热方程并计算从y=0处导入的热量。当时，该二维导热问题可简化为一维导热问题。在一维的近似下，试计算从y=0处导入的热量，并比较不同L2/L1下的比值。由该问题的解析解可知：L2/L10.0070.010.050.080.10.99870.99120.9560.930.912解：编程得到M=N=50温度分布数据见附件二，得到温度分布图如下：这里也仅给出M=N=10时温度分布数据，见下表：606060606060606060606058.9248458.0307557.4169257.1105157.1105157.4169258.0307558.92484606057.7692255.9369354.697254.0848354.0848354.697255.9369357.76922606056.4364853.5797551.6986550.78750.78751.6986553.5797556.43648606054.7896250.7889648.2692447.0839947.0839948.2692450.7889654.78962606052.605147.3399744.2467242.854242.854244.2467247.3399752.6051606049.4740242.925739.4683438.0055138.0055139.4683442.925749.47402606044.5750837.1376833.8034832.5066732.5066733.8034837.1376844.57508606036.1261529.5260127.2331326.4314826.4314827.2331329.5260136.126156020202020202020202020从y=0处导入的热量可以近似看作从y=0向y=0.4/9处传递的热量，y=0.4/9处的温度分别为：x00.6/90.6/9*20.6/9*30.6/9*40.6/9*50.6/9*60.6/9*70.6/9*80.6T6058.9248458.0307557.4169257.1105157.1105157.4169258.0307558.9248460则整理得当时，该二维导热问题可简化为一维导热问题。当时，令M=N=100，得=1129571.429看作一维导热时不同下比值如下表：L2/L10.0070.010.050.080.11129571.42978954015281292890729081142857.143800000160000100000800000.98840.98690.95510.92890.91140.99870.99120.9560.930.912附录Matlab程序：第一题：functionchuanredazuoyeclearclcL1=0.6;%矩形长度L2=0.4;%矩形宽度Tw1=60;Tw2=20;globalMN%设置网格数M=input('请输入将区间[0,L1]等分的个数M:');N=input('请输入将区间[0,L2]等分的个数N:');XDIF=L1/M;%x方向上的步长YDIF=L2/N;%y方向上的步长axis([0,M,0,N]);grid%设置网格U=initial(M,N,Tw1,Tw2,L1,XDIF);AIM=YDIF/XDIF;AIP=YDIF/XDIF;AJM=XDIF/YDIF;AJP=XDIF/YDIF;AP=2*YDIF/XDIF+2*XDIF/YDIF;CON=0;%离散方程系数T=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP);A=flipud(T);mesh(A);title('温度分布')xlabel('x')ylabel('y')zlabel('T')functionU=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP)while1temp=U;fori=2:M-1forj=2:N-1U(i,j)=(AIM*U(i,j-1)+AIP*U(i,j+1)+AJM*U(i-1,j)+AJP*U(i+1,j))/AP;%设置迭代条件endendeps=abs(U-temp);ifmax(max(eps))<1e-8break;%限制迭代次数endendfunctionU=initial(M,N,Tw1,Tw2,L1,XDIF)U=zeros(M,N);%赋温度矩阵初值U(:,1)=Tw1;U(:,M)=Tw1;U(N,:)=Tw1;fori=1:MU(1,i)=Tw1+Tw2*sin(pi*i*XDIF/L1);end%初始和边界条件的设定第二题：functiongaibianjietiaojianclearclcL1=0.6;%矩形长度L2=0.4;%矩形宽度Tw1=60;Tw2=20;globalMN%设置网格数M=input('请输入将区间[0,L1]等分的个数M:');N=input('请输入将区间[0,L2]等分的个数N:');XDIF=L1/M;%x方向上的步长YDIF=L2/N;%y方向上的步长axis([0,M,0,N]);grid%设置网格U=initial(M,N,Tw1,Tw2);AIM=YDIF/XDIF;AIP=YDIF/XDIF;AJM=XDIF/YDIF;AJP=XDIF/YDIF;AP=2*YDIF/XDIF+2*XDIF/YDIF;CON=0;%离散方程系数T=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP);A=flipud(T);mesh(A);title('改变边界条件后温度分布')xlabel('x')ylabel('y')zlabel('T')functionU=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP)while1temp=U;fori=2:M-1forj=2:N-1U(i,j)=(AIM*U(i,j-1)+AIP*U(i,j+1)+AJM*U(i-1,j)+AJP*U(i+1,j))/AP;%设置迭代条件endendeps=abs(U-t