《传热学数值计算》由会员分享，可在线阅读，更多相关《传热学数值计算（41页珍藏版）》请在金锄头文库上搜索。
1、传热学数值计算作业传热学数值计算作业数值解程序： tw1=40 %三边温度 tw2=100 %一边温度正弦变化幅度 l1=40 %板长 L1:40 厘米 l2=20 %板宽 L2:20 厘米 m=41 %分划成 40*20 的网格 n=21 k=2 dx=l1/(m-1) c=ones(n,m) for i=1:m a2(i)=tw1+tw2*sin(pi*dx*(i-1)/l1) c(1,i)=tw1 ,c(n,i)=a2(i) end for j=1:n c(j,1)=tw1 c(j,m)=tw1 end while (abs(c(j,i)-k)0.0001) k=c(j,i) for i
2、=2:m-1 for j=2:n-1 c(j,i)=0.25*(c(j,i-1)+c(j,i+1)+c(j-1,i)+c(j+1,i) end endend数值解中各网格点的温度值：数值二维温度分布图像：010203040500102030406080100120140L1/cmnumerical calculation 2D temperature distributionL2/cmtemperature/celsius degree解析解程序： tw1=40 tw2=100 l1=40 l2=20 p=40 q=20 x(1)=0 for i=1:p x(i+1)=x(i)+1 end y
3、(1)=0 for j=1:q y(j+1)=y(j)+1 end for i=1:p+1 for j=1:q+1 n(j,i)=tw1+tw2*sinh(pi*y(j)/l1)*sin(pi*x(i)/l1)/sinh(pi*l2/l1) end end各网格点用解析式得到的温度值：解析二维温度分布图像：010203040500102030406080100120140L1/cmanalytical method 2D temperature distributionL2/cmtemperature/celsius degree误差分析：取 x=21，即位于板长一半处，温度随 y（宽度）的变
4、化曲线。 c1(:,1) 取自于数值解， c1(:,2) 取自于解析解 c1(:,1) c1(:,2)40.0000 40.000043.3106 43.416446.6465 46.853850.0313 50.333553.4889 53.877157.0430 57.506260.7178 61.243464.5376 65.111768.5273 69.135072.7122 73.338177.1187 77.747081.7736 82.388886.7050 87.292291.9423 92.487597.5162 98.0068103.4592 103.8840109.805
5、8 110.1555116.5925 116.8600123.8586 124.0388131.6461 131.7363140.0000 140.0000误差曲线：0510152025405060708090100110120130140L2/cmtemperature/celsius degreenumerical calculation and analytical method temperature distribution differentials at x=21 cross sectionc1(:,1) c1(:,2)由相对误差公式：d1= (c1(:,2) -c1(:,1).
6、/ c1(:,2) 可得： d1 =00.00240.00440.00600.00720.00810.00860.00880.00880.00850.00810.00750.00670.00590.00500.00410.00320.00230.00150.00070结论：数值解与解析解吻合很好。就 x=21 这一列，相对误差较小且在 1%以内，但数值解 较解析解偏小，且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。数值解程序： tw1=50 %三边温度 tw2=100 %一边温度 l1=20 %板长 20 厘米 l2=10 %板宽 10 厘米 m=41 %划分成 40*20 的
7、网格 n=21 k=1 c=zeros(n,m) c(n,1)=(tw1+tw2)/2 c(n,m)=(tw1+tw2)/2 c(1,1)=tw1 c(1,m)=tw1 for i=2:(m-1) c(1,i)=tw1 c(n,i)=tw2 end for j=2:(n-1) c(j,1)=tw1 c(j,m)=tw1 end while(abs(k-c(j,i)0.0001) k=c(j,i) for i=2:m-1 for j=2:n-1 c(j,i)=0.25*(c(j,i-1)+c(j,i+1)+c(j-1,i)+c(j+1,i) end end end数值解中各网格点的温度值：数值二
8、维温度分布图像：0102030405001020305060708090100L1/cmnumerical calculation 2D tenmperature distributionL2/cmtemperature/celsius degree解析解程序： tw1=50 tw2=100 m=40 n=20 l1=20 l2=10 tx=ones(20,40) for i=1:m+1 for j=1:n+1 x=(i-1)*0.5 y=(j-1)*0.5 k=sym(k) d=(-1)(k+1)+1)/k)*sin(k*pi*x/l1)*sinh(k*pi*y/l1)/sinh(k*pi*
9、l2/l1) h=symsum(d,k,1,200) tx(j,i)=2*h*(tw2-tw1)/pi+tw1 end end各网格点用解析式得到的温度值：解析二维温度分布图像：0102030405001020305060708090100110L1/cmanalytical method 2D temperature distributionL2/cmtemperature/celsius degree误差分析：取 x=21，即位于板长一半处，温度随 y（宽度）的变化曲线。 c1(:,1) 取自于数值解 c1(:,2) 取自于解析解 c1(:,1)= c1(:,2)=50.0000 50.0
10、00051.9799 52.088153.9731 54.185155.9906 56.299858.0435 58.440960.1418 60.616562.2951 62.834464.5116 65.101666.7987 67.424369.1622 69.807771.6064 72.255874.1337 74.771076.7447 77.354679.4377 80.005682.2090 82.721485.0526 85.497787.9602 88.327890.9214 91.203593.9244 94.115196.9554 97.0513100.0000 99.
11、8408误差分析曲线：051015202550556065707580859095100L2/cmtemperature/celsius degreenumerical calculation and analytical method temperature distribution differentials at x=21 cross sectionc1(:,1) c1(:,2)由相对误差公式：d2= abs(c1(:,2) -c1(:,1)./ c1(:,2) 可得： d2=00.00210.00390.00550.00680.00780.00860.00910.00930.00920
12、.00900.00850.00790.00710.00620.00520.00420.00310.00200.0010 0.0016 误差结论：数值解与解析解吻合很好。就 x=21 这一列，相对误差较小且在 1%以内，但数 值解较解析解普遍偏小（最后一个数除外） ，且在平板中心附近的网格点的数值解较平板边 缘数值解的相对误差大。数值解程序： t0=200; %肋基温度为 200 度 tf=0; %环境温度为 0 度 L1=20; %肋片长 20 厘米 L2=2; %肋片高 2 厘米 h=0.01; %对流换热系数 单位：w/(cm2*K) a=11; b=101; %11*101 的网格 dx
13、=L1/(b-1); k=2; %导热系数 单位：w/（cm*K） Bi=h*dx/k; ti=ones(a,b)*10; m1=ones(a,b)*3; m1(2:a-1,1)=zeros(a-2,1); m1(a,2:b-1)=ones(1,b-2); m1(1,2:b-1)=ones(1,b-2)*6; m1(2:a-1,b)=ones(a-2,1)*2; m1(1,b)=ones(1,1)*4; m1(a,b)=ones(1,1)*5; m1(1,1)=7; m1(a,1)=8; tn=ti; max1=1.0; w=0; while ( max11e-6)w=w+1;max1=0; for i=1:afor j=1:bm=m1(i,j);n=tn(i,j);switch mcase 0tn(i,j)=t0;case 1tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4*tf)/(4+2*Bi)+tf;case 2tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4*tf)/(4+2*Bi)+tf;case 3tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+
