传热学数值计算Word文件下载.docx

1、y（j+1）=y（j）+1p+1q+1n（j,i）=tw1+tw2*sinh（pi*y（j）/l1）*sin（pi*x（i）/l1）/sinh（pi*l2/l1）各网格点用解析式得到的温度值：解析二维温度分布图像：误差分析：取x=21，即位于板长一半处，温度随y（宽度）的变化曲线。c1（:,1） 取自于数值解， c1（:,2） 取自于解析解,1） c1（:,2） 40.0000 40.0000 43.3106 43.4164 46.6465 46.8538 50.0313 50.3335 53.4889 53.8771 57.0430 57.5062 60.7178 61.2434 64.53

2、76 65.1117 68.5273 69.1350 72.7122 73.3381 77.1187 77.7470 81.7736 82.3888 86.7050 87.2922 91.9423 92.4875 97.5162 98.0068 103.4592 103.8840 109.8058 110.1555 116.5925 116.8600 123.8586 124.0388 131.6461 131.7363 140.0000 140.0000误差曲线：由相对误差公式：d1= （c1（:,2） -c1（:,1）./ c1（:,2） 可得：d1 = 0 0.0024 0.0044 0

3、.0060 0.0072 0.0081 0.0086 0.0088 0.0085 0.0075 0.0067 0.0059 0.0050 0.0041 0.0032 0.0023 0.0015 0.0007结论：数值解与解析解吻合很好。就x=21这一列，相对误差较小且在1%以内，但数值解较解析解偏小，且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。tw1=50 %三边温度tw2=100 %一边温度l1=20 %板长20厘米l2=10 %板宽10厘米m=41 %划分成40*20的网格k=1 c=zeros（n,m）c（n,1）=（tw1+tw2）/2 c（n,m）=（tw1+tw2

4、）/2 c（1,1）=tw1 c（1,m）=tw1 （m-1） c（1,i）=tw1c（n,i）=tw2（n-1）while（abs（k-c（j,i）tw1=50m=40n=20l1=20l2=10tx=ones（20,40）m+1n+1x=（i-1）*0.5y=（j-1）*0.5k=sym（k）d=（-1）（k+1）+1）/k）*sin（k*pi*x/l1）*sinh（k*pi*y/l1）/sinh（k*pi*l2/l1）h=symsum（d,k,1,200）tx（j,i）=2*h*（tw2-tw1）/pi+tw1,1） 取自于数值解 c1（:,1）= c1（:,2）= 50.0000 50

5、.0000 51.9799 52.0881 53.9731 54.1851 55.9906 56.2998 58.0435 58.4409 60.1418 60.6165 62.2951 62.8344 64.5116 65.1016 66.7987 67.4243 69.1622 69.8077 71.6064 72.2558 74.1337 74.7710 76.7447 77.3546 79.4377 80.0056 82.2090 82.7214 85.0526 85.4977 87.9602 88.3278 90.9214 91.2035 93.9244 94.1151 96.955

6、4 97.0513 100.0000 99.8408误差分析曲线：d2= abs（c1（:d2= 0.0021 0.0039 0.0055 0.0068 0.0078 0.0091 0.0093 0.0092 0.0090 0.0079 0.0071 0.0062 0.0052 0.0042 0.0031 0.0020 0.00100.0016误差结论：就x=21这一列，相对误差较小且在1%以内，但数值解较解析解普遍偏小（最后一个数除外），且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。t0=200; %肋基温度为200度tf=0; %环境温度为0度L1=20; %肋片长20厘米

7、 L2=2; %肋片高2厘米h=0.01; %对流换热系数 单位：w/（cm2*K）a=11;b=101; %11*101的网格 dx=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;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=t

8、i;max1=1.0;w=0;while （ max11e-6） w=w+1; max1=0; for i=1:a for j=1:b m=m1（i,j）; n=tn（i,j）; switch m case 0 tn（i,j）=t0; case 1 tn（i,j）=（2*tn（i-1,j）+tn（i,j-1）+tn（i,j+1）-4*tf）/（4+2*Bi）+tf; case 2 tn（i,j）=（2*tn（i,j-1）+tn（i-1,j）+tn（i+1,j）-4*tf）/（4+2*Bi）+tf; case 3 tn（i,j）=0.25*（tn（i,j-1）+tn（i,j+1）+tn（i-1,

9、j）+tn（i+1,j）; case 4 tn（i,j）=（tn（i,j-1）+tn（i+1,j）-2*tf）/（2*Bi+2）+tf; case 5 tn（i,j）=（tn（i,j-1）+tn（i-1,j）-2*tf）/（2*Bi+2）+tf; case 6 tn（i,j）=（2*tn（i+1,j）+tn（i,j-1）+tn（i,j+1）-4*tf）/（4+2*Bi）+tf; case 7 case 8 end er=abs（tn（i,j）-n）; if ermax1 max1=er; ti=tn;ti %假设垂直于纸面方向肋宽为1cmq1=0;qi=（tn（i,b）-tf）*h*0.2;q

10、1=q1+qi;q1 % x=20cm 处肋片的散热量q2=0;b-1q2j =（tn（a,j）-tf）*h*0.2;q2=q2+q2j;q2=q2*2;q2 % y=1cm与y=-1cm处肋片的散热量q=q1+q2;数值解肋片换热量：q1 =1.9022w/cm 肋端沿y方向q2 =49.4938w/cm 肋端沿x方向q = 51.3960w/cm q=q1+q2Bi=h*1/k;ta=ones（a,b）;1:if i（a+1）/2 y=-（i-（a+1）/2）*dx; else y=（a+1）/2-i）*dx;x=dx*（j-1）;ta（i,j）=（cosh（Bi0.5*（L1-x）/1）

11、+Bi0.5*sinh（Bi0.5*（L1-x）/1）*（t0-tf）/（cosh（Bi0.5*L1/1）+ Bi0.5*sinh（Bi0.5*L1/1）+tf;taqi=（ta（i,b）-tf）*h*0.2;q2j =（ta（a,j）-tf）*h*0.2;解析解二维温度图像：解析解肋片换热量：q1 =1.9006w/cm 肋端沿y方向q2 =49.5481w/cm 肋端沿x方向q = 51.4488w/cm q=q1+q2取y=0，即位于板宽一半处，温度随x（长度）的变化曲线。b1（:,1） 取自于数值解 b1（:,1）= b1（:200.0000 200.0000 197.5383 197

12、.4675 195.0883 194.9745 192.6596 192.5205 190.2587 190.1051 187.8895 187.7276 185.5543 185.3877 183.2542 183.0848 180.9896 180.8186 178.7605 178.5885 176.5667 176.3942 174.4081 174.2351 172.2842 172.1109 170.1946 170.0211 168.1389 167.9654 166.1169 165.9432 164.1280 163.9542 162.1718 161.9979 160.24

13、80 160.0741 158.3563 158.1823 156.4961 156.3222 154.6672 154.4933 152.8692 152.6953 151.1017 150.9278 149.3644 149.1906 147.6569 147.4831 145.9788 145.8052 144.3300 144.1564 142.7099 142.5365 141.1184 140.9451 139.5550 139.3818 138.0195 137.8464 136.5115 136.3386 135.0308 134.8581 133.5771 133.4046

14、132.1500 131.9777 130.7494 130.5772 129.3748 129.2028 128.0261 127.8543 126.7029 126.5314 125.4051 125.2337 124.1322 123.9611 122.8842 122.7133 121.6607 121.4901 120.4615 120.2911 119.2863 119.1162 118.1350 117.9651 117.0073 116.8376 115.9029 115.7335 114.8216 114.6525 113.7633 113.5945 112.7277 112

15、.5592 111.7147 111.5464 110.7239 110.5559 109.7552 109.5875 108.8085 108.6410 107.8835 107.7163 106.9800 106.8131 106.0979 105.9312 105.2369 105.0706 104.3970 104.2309 103.5780 103.4122 102.7796 102.6140 102.0017 101.8365 101.2442 101.0793 100.5069 100.3423 99.7897 99.6253 99.0925 98.9283 98.4150 98

16、.2511 97.7571 97.5935 97.1188 96.9555 96.4999 96.3368 95.9002 95.7374 95.3197 95.1572 94.7583 94.5960 94.2157 94.0537 93.6920 93.5302 93.1870 93.0254 92.7006 92.5393 92.2327 92.0716 91.7832 91.6224 91.3520 91.1915 90.9391 90.7788 90.5443 90.3842 90.1677 90.0078 89.8090 89.6494 89.4683 89.3088 89.145

17、4 88.9862 88.8403 88.6813 88.5530 88.3942 88.2834 88.1248 88.0314 87.8729 87.7969 87.6387 87.5800 87.4220 87.3806 87.2228 87.1987 87.0410 87.0341 86.8766 86.8869 86.7296 86.7571 86.6000 86.6446 86.4876 86.5494 86.3926曲线适当放大：d2= abs（b1（:,2） -b1（:,1）./ b1（: 0.0004 0.0006 0.0008 0.0009 0.0011 0.0012 0.

18、0013 0.0014 0.0016 0.0017 0.00180.0018肋片总散热量数值解与解析解的相对误差为q（解析解）= 51.4488w/cm q（数值解） = 51.3960w/cmt=abs（q（解析解）- q（数值解）/ q（解析解）=0.001026就y=0这行，相对误差较小且在0.2%以内，但数值解较解析解普遍偏大，且离肋基较近的网格点的数值解比离肋基较远的数值解相对误差大。解析解的肋片总散热量与数值解吻合较好，但总散热量比数值解大。tw=0; %环境温度0度L2=10; %板宽的一半10厘米c=0.5;L1=L2/c; %板长的一半20厘米p2=40; %划分成81*41

19、的网格p1=p2/c;dx=2*L2/p2;a=p2+1;b=p1+1; %导热系数 单位w/（m2*K）q=3; %单位体积的产热率w/m3ti=ones（a,b）*5;m1=ones（a,b）;m1（1,:）=zeros（1,b）;a,b）=zeros（a-1,1）;a,1）=zeros（a-1,1）;b-1）=zeros（1,b-2）;y=0;while（max1max1=0;y=y+1;m=m1（i,j）;n=tn（i,j）;switch mcase 0tn（i,j）=tw;case 1tn（i,j）=0.25*（tn（i-1,j）+tn（i+1,j）+tn（i,j-1）+tn（i,j

20、+1）+q*（dx2）/k）;er=abs（tn（i,j）-n）;if erti=tn;a1=ti（21:41,41:81）tw1=0;m=41;n=21;dx=L1/（m-1）;tx=ones（n,m）;mx=（i-1）*0.5;y=（j-1）*0.5;p=sym（p）;xi=（2*p -1）*pi/2;d=（-1）p）/（xi3）*（cosh（xi*y/L1）/cosh（xi*L2/L1）*cos（xi*x/L1）;h=symsum（d,p ,1,100）;tx（j,i）=2*q*L12/k*h+q*（L12-x2）/（2*k）+tw1;tx取x=1，即位于板长一半处，温度随y（宽度）的变化曲线。68.3121 68.3231 68.1452 68.1562 67.6444 67.6553 66.8093 66.8200 65.6393