传热学数值计算作业数值解程序：tw1=40 %三边温度tw2=100 %一边温度正弦变化幅度l1=40 %板长L1:40厘米l2=20 %板宽L2:20厘米m=41 %分划成40*20的网格n=21k=2dx=l1/(m-1)c=ones(n,m)for i=1:ma2(i)=tw1+tw2*sin(pi*dx*(i-1)/l1)c(1,i)=tw1 ,c(n,i)=a2(i)endfor j=1:nc(j,1)=tw1c(j,m)=tw1endwhile (abs(c(j,i)-k)>0.0001)k=c(j,i)for i=2:m-1for j=2:n-1c(j,i)=0.25*(c(j,i-1)+c(j,i+1)+c(j-1,i)+c(j+1,i)) endendend数值解中各网格点的温度值：数值二维温度分布图像：解析解程序： tw1=40 tw2=100 l1=40 l2=20 p=40 q=20 x(1)=0 for i=1:px(i+1)=x(i)+1 end y(1)=0 for j=1:qy(j+1)=y(j)+1 endfor i=1:p+1 for j=1:q+1n(j,i)=tw1+tw2*sinh(pi*y(j)/l1)*sin(pi*x(i)/l1)/sinh(pi*l2/l1) end end各网格点用解析式得到的温度值：50L1/cmnumerical calculation 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e解析二维温度分布图像：误差分析：取x=21，即位于板长一半处，温度随y （宽度）的变化曲线。

c1(:,1) 取自于数值解， c1(:,2) 取自于解析解 c1(:,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.5376 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.000050L1/cmanalytical method 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e误差曲线：由相对误差公式：d1= (c1(:,2) -c1(:,1))./ c1(:,2) 可得： d1 = 0 0.0024 0.0044 0.00600.00720.0081 0.0086 0.0088 0.0088 0.0085 0.0081 0.0075 0.0067 0.0059 0.0050 0.0041 0.0032 0.0023 0.0015 0.0007 0结论：数值解与解析解吻合很好。

就x=21这一列，相对误差较小且在1%以内，但数值解较解析解偏小，且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。

510152025L2/cmt e m p e r a t u r e /c e l s i u s d e g r e enumerical calculation and analytical method temperature distribution differentials at x=21 cross section数值解程序：tw1=50 %三边温度tw2=100 %一边温度l1=20 %板长20厘米l2=10 %板宽10厘米m=41 %划分成40*20的网格n=21k=1c=zeros(n,m)c(n,1)=(tw1+tw2)/2c(n,m)=(tw1+tw2)/2c(1,1)=tw1c(1,m)=tw1for i=2:(m-1)c(1,i)=tw1c(n,i)=tw2endfor j=2:(n-1)c(j,1)=tw1c(j,m)=tw1endwhile(abs(k-c(j,i))>0.0001)k=c(j,i)for i=2:m-1for j=2:n-1c(j,i)=0.25*(c(j,i-1)+c(j,i+1)+c(j-1,i)+c(j+1,i)) endendend数值解中各网格点的温度值：数值二维温度分布图像：解析解程序： tw1=50 tw2=100 m=40 n=20 l1=20 l2=10tx=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*l2/l1) h=symsum(d,k,1,200) tx(j,i)=2*h*(tw2-tw1)/pi+tw1 end end各网格点用解析式得到的温度值：50L1/cmnumerical calculation 2D tenmperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e解析二维温度分布图像：误差分析：取x=21，即位于板长一半处，温度随y （宽度）的变化曲线。

c1(:,1) 取自于数值解 c1(:,2) 取自于解析解 c1(:,1)= c1(:,2)=50.0000 50.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.9554 97.0513 100.0000 99.840850L1/cmanalytical method 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e误差分析曲线：由相对误差公式：d2= abs(c1(:,2) -c1(:,1))./ c1(:,2) 可得： d2= 0 0.0021 0.0039 0.0055 0.00680.00780.0086 0.0091 0.0093 0.0092 0.0090 0.0085 0.0079 0.0071 0.0062 0.0052 0.0042 0.0031 0.0020 0.00100.0016误差结论：数值解与解析解吻合很好。

就x=21这一列，相对误差较小且在1%以内，但数值解较解析解普遍偏小（最后一个数除外），且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。

510152025L2/cmt e m p e r a t u r e /c e l s i u s d e g r e enumerical calculation and analytical method temperature distribution differentials at x=21 cross section数值解程序：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=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 ( max1>1e-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+1)+tn(i-1,j)+tn(i+1,j));case 4tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2*tf)/(2*Bi+2)+tf;case 5tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2*tf)/(2*Bi+2)+tf;case 6tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4*tf)/(4+2*Bi)+tf;case 7tn(i,j)=t0;case 8tn(i,j)=t0;ender=abs(tn(i,j)-n);if er>max1max1=er;endendendti=tn;endti%假设垂直于纸面方向肋宽为1cmq1=0;for i=1:aqi=(tn(i,b)-tf)*h*0.2;q1=q1+qi;endq1 % x=20cm 处肋片的散热量q2=0;for j=1:b-1q2j =(tn(a,j)-tf)*h*0.2;q2=q2+q2j;endq2=q2*2;q2 % y=1cm与y=-1cm处肋片的散热量q=q1+q2;q数值解中各网格点的温度值：数值二维温度分布图像：数值解肋片换热量：q1 =1.9022w/cm 肋端沿y 方向 q2 =49.4938w/cm 肋端沿x 方向 q = 51.3960w/cmq=q1+q2150L1/cmnumerical calculation 2D temperature distributionL2/cmt e m p e r a t u r e /c e l s i u s d e g r e e解析解程序：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=L1/(b-1);k=2; %导热系数单位：w/（cm*K）Bi=h*1/k;ta=ones(a,b);for i=1:1:afor j=1:1:bif i>(a+1)/2y=-(i-(a+1)/2)*dx;else y=((a+1)/2-i)*dx;endx=dx*(j-1);ta(i,j)=(cosh(Bi^0.5*(L1-x)/1)+Bi^0.5*sinh(Bi^0.5*(L1-x)/1))*(t0-tf)/(cosh(Bi^0.5*L1/1)+ Bi^0.5*sinh(Bi^0.5*L1/1))+tf;endendta%假设垂直于纸面方向肋宽为1cmq1=0;for i=1:1:aqi=(ta(i,b)-tf)*h*0.2;q1=q1+qi;endq1 % x=20cm 处肋片的散热量q2=0;for j=1:1:b-1q2j =(ta(a,j)-tf)*h*0.2;q2=q2+q2j;endq2=q2*2;q2 % y=1cm与y=-1cm处肋片的散热量q=q1+q2;q各网格点用解析式得到的温度值：解析解二维温度图像：解析解肋片换热量：q1 =1.9006w/cm 肋端沿y 方向 q2 =49.5481w/cm 肋端沿x 方向 q = 51.4488w/cmq=q1+q2150L1/cmanalytical method 2D temperature distributionL2/cmt e m p e r a t u r e /c l e s i u s d e g r e e误差分析：取y=0，即位于板宽一半处，温度随x（长度）的变化曲线。