传热学模拟大实验-二维导热温度场的数值计算

1、传热学模拟大实验-二维导热温度场的数值计算传热学模拟大实验二维导热温度场的数值计算学院：能动学院名字： 班级：专业：能源动力系统及自动化学号：一、实验题目： 有一个用砖砌成的长方形截面的冷空气空道，其截面尺寸如下图所示，假设在垂直于纸面方向上冷空气及砖墙的温度变化很小，可以近似地予以忽略。在下列两种情况下试计算：（1）砖墙横截面上的温度分布；（2）垂直于纸面方向的每米长度上通过砖墙的导热量。外矩形长为，宽为；内矩形长为，宽为。第一种情况：内外壁分别均匀地维持在0及30（等温边界条件）；第二种情况：内外表面均为第三类边界条件（对流边界条件），且已知：外壁：t1=30 , h1=10w/m内壁：t

2、2=10 , h2=4w/m砖墙的导热系数=0.53w/m；由于对称性，仅研究1/4部分即可。二、数学描述1.基本方程此问题为无内热源的二维稳态导热问题，所以可得出方程： 2tx2+2ty2=02.边界条件对边界条件进行分析，对于边界恒温的第一种情况有：边界绝热；q=0，且t1=30，t2=10对于第三类边界条件，可得：q=h×(tw-tf)三、方程离散 将上图进行区域离散，且x=y=0.1m，并建立内外节点的离散方程第一种情况：内点：ti,j=(ti-1,j+ti+1,j+ti,j-1+ti,j+1)/4边界点： ab边界（a、b除外）i=2，3，4，5； j=1 ti,1=(ti

3、-1,1+ti+1,1+2ti,2)/4 de边界（d、e除外）i=16; j=8，9，10，11t16,j=(t16,j+1+t16,j-1+2t15,j)/4 af、fe边界(包括a、f、e点）i=1；j=1，2，3，12j=12；i=1，2，3，16ti,j=30 bc、cd边界（包括b、c、d点）i=6；j=1，2，7j=7；i=6，8，16ti,j=0第二种情况：内点：ti,j=(ti-1,j+ti+1,j+ti,j-1+ti,j+1)/4边界点：ab边界（a、b除外）i=2，3，4，5； j=1 ti,1=(ti-1,1+ti+1,1+2ti,2)/4de边界（d、e除外）i=16

4、; j=8，9，10，11t16,j=(t16,j+1+t16,j-1+2t15,j)/4af、fe边界(不包括a、f、e点）i=1；j=2，3，11j=12；i=2，3，152(h1x+2)ti,j=2ti+1,j+ti,j-1+ti,j+1+2h1xt1bc、cd边界（不包括b、c、d点）i=6；j=2，3，6j=7；i=7，8，152(h2x+2)ti,j=2ti-1,j+ti,j-1+ti,j+1+2h1xt2a点：(h1x+2)t1,1=t2,1+t1,2+h1xt1b点：(h2x+2)t6,1=t5,1+t6,2+h2xt2c点：2(h2x+3)t6,7=2t5,7+t6,8+t7

5、,7+t6,6+2h1xt2d点：(h2x+2)t16,7=t15,7+t16,8+h2xt2e点：(h1x+2)t16,12=t15,12+t16,11+h1xt1f点：2(h1x+1)t1,12=t2,12+t1,11+2h1xt1四、编程方程式采用高斯-塞德尔迭代法设立初值,通过程序计算与上一次的值进行对比,设立需要的最大差值,进行多次迭代,最后得出结果.（具体编程式见后附）五、实验结果1.第一种情况:(式中k为迭代次数)k = 151s（温度场分布） = columns 1 through 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

6、 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 columns 15 through 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0每米外砖墙导热量q1=每米内砖墙导热量q2=导热量平均值 q=热平衡偏差 detq=具体程序：cleart=zeros(12,16);s=zeros(12,16);for k=1:100000 s(1:12,1)=30; s(1,2:16)=30; for i=2:5 s(i,16)=(t(i+1,16)+t(i-1,16)+2*t(i,

7、15)/4;endfor i=2:5 for j=2:15 s(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j+1)+t(i,j-1)/4; endendfor i=6:11 for j=2:5 s(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j+1)+t(i,j-1)/4; endendfor j=2:5 s(12,j)=(t(12,j-1)+t(12,j+1)+2*t(11,j)/4;endif norm(s-t)< k s break;endt=s;endif k=10000 warning('未能找到近似解');end w=zeros; f

8、or i=1:12 for j=1:16 w(13-i,j)=s(i,j); end end v=12 18 24; contour(w,v); set(gca,'xtick',1:16,'ytick',1:12); grid onq1=0; for i=2:11 q1=q1+*(30-t(i,2); end for j=2:15 q1=q1+*(30-t(2,j); end disp('每米外砖墙导热量q1：'); q1=(q1+*(30-t(2,16)+*(30-t(12,2)/q2=0; for i=6:11 q2=q2+*(t(i,5)-

9、0); end for j=6:15 q2=q2+*(t(5,j)-0); end disp('每米内砖墙导热量q2：'); q2=(q2+*(t(12,5)-0)+*(t(5,16)-0)/disp('导热量平均值q:')q=(q1+q2)/2disp('热平衡偏差detq:')detq=abs(q1-q2)2.第二种情况：(k为迭代次数)k = 263s（温度场分布） = columns 1 through 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

10、 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 columns 15 through 16 0 0 0 0 0 0 0 0 0 0 0 0每米外砖墙导热量q1=每米内砖墙导热量q2=导热量平均值 q=热平衡偏差 detq=具体程序：clear h1=10;h2=4;tf1=30;tf2=10;drxs=;x=;t=zeros(12,16);s=zeros(12,16);for k=1:100000 for i=2:11 s(i,1)=(2*t(i,2)+t(i-1,1)+t(i+1,1)+2*h1*x*tf1/drxs)/(2*(h1*x/drxs+2);endfor j=2:1

11、5 s(1,j)=(2*t(2,j)+t(1,j-1)+t(1,j+1)+2*h1*x*tf1/drxs)/(2*(h1*x/drxs+2);endfor i=2:5 s(i,16)=(t(i+1,16)+t(i-1,16)+2*t(i,15)/4;endfor j=2:5 s(12,j)=(t(12,j-1)+t(12,j+1)+2*t(11,j)/4;endfor i=7:11 s(i,6)=(2*t(i,5)+t(i-1,6)+t(i+1,6)+2*h2*x*tf2/drxs)/(2*(h2*x/drxs+2);endfor j=7:15 s(6,j)=(2*t(5,j)+t(6,j-1

12、)+t(6,j+1)+2*h2*x*tf2/drxs)/(2*(h2*x/drxs+2);endfor i=2:5 for j=2:15 s(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j+1)+t(i,j-1)/4; endendfor i=6:11 for j=2:5 s(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j+1)+t(i,j-1)/4; endends(12,1)=(t(11,1)+t(12,2)+h1*x*tf1/drxs)/(h1*x/drxs+2);s(12,6)=(t(12,5)+t(11,6)+h2*x*tf2/drxs)/(h2*x/drxs

13、+2);s(6,6)=(2*(t(5,6)+t(6,5)+t(6,7)+t(7,6)+2*h2*x*tf2/drxs)/(2*(h2*x/drxs+3);s(6,16)=(t(6,15)+t(5,16)+h2*x*tf2/drxs)/(h2*x/drxs+2);s(1,16)=(t(1,15)+t(2,16)+h1*x*tf1/drxs)/(h1*x/drxs+2);s(1,1)=(t(1,2)+t(2,1)+2*h1*x*tf1/drxs)/(2*(h1*x/drxs+1);if norm(s-t)< k s break;endt=s;endif k=10000 warning(

14、9;未能找到近似解');end w=zeros; for i=1:12 for j=1:16 w(13-i,j)=s(i,j); end end v=18 22 26; contour(w,v); set(gca,'xtick',1:16,'ytick',1:12); grid onq1=0; for i=2:11 q1=q1+h1*x*(30-t(i,1); end for j=2:15 q1=q1+h1*x*(30-t(1,j); end disp('每米外砖墙导热量q1：'); q1=(q1+*h1*x*(30-t(1,16)+*h1*x*(30-t(12,1)+h1*x*(30-t(1,1)/ q2=0; for i=7:11 q2=q2+h2*x*(t(i,6)-10); end for j=7:15 q2=q2+h2*x*(t(6,j)-10); end disp('每米内砖墙导热量q2：'); q2=(