T1(i,j)=M1*(T(i-1,j)+T(i+1,j))+M2*(T(i,j-1)+T(i,j+1))+(1-2*M1-2*M2)*T(i,j); end end

for i=2:8 for j=2:8

T(i,j)=T1(i,j); end end if(k==5) l=l+1;

Tc(l)=T(5,5); k=0; end

end

i=1:9;j=1:9;

[x,y]=meshgrid(i); figure(1);

subplot(1,2,1);

mesh(x,y,T(i,j))% 画出10min后的温度场 axis tight;

xlabel('x','FontSize',14);ylabel('y','FontSize',14);zlabel('T/℃','FontSize',14) title('1min后二维温度场模拟图','FontSize',18) subplot(1,2,2);

[C,H]=contour(x,y,T(i,j)); clabel(C,H);axis square;

xlabel('x','FontSize',14);ylabel('y','FontSize',14); title('1min后模拟等温线图','FontSize',18) figure(2); xx=1:600;

plot(xx,Tc,'k-','linewidth',2)

xlabel('时间/s','FontSize',14);ylabel('温度/℃','FontSize',14);title('中心点的冷却曲线','FontSize',18)

else disp('Error!') % 如果不满足稳定性条件，显示“Error！” end

实验结果：

中心点的冷却曲线500450400350温度/℃3002502001501000100200300400500600时间/s

1min后二维温度场模拟图100.06100.05100.04℃/T100.03100.02100.0110098768546342y12x5min后二维温度场模拟图106105104℃/T10310210110098768546342y12x10min后二维温度场模拟图100.06100.05100.04℃/T100.03100.02100.0110098768546342y12x1min后模拟等温线图9100.018100.02100.037410100.100.01.00201000..61050100014.100.0651000.000.010010.030.0y542110041.0016010.0.0003.00.05013101010.0100.042100.0320.2001100.011123456789x5min后模拟等温线图9101820111031070132041011051165010110310604460y50111024102011103105311040301121012011123456789x10min后模拟等温线图9100.01820.0100.031100017.100.040.002011100004.06015.050.100010.0.00.10604y5103100.1000.00124100.160.000.0200.01010130.3100100.50.0400.031210100.01200.011123456789x

10min后二维温度场模拟图(不满足稳定性条件)x 102121T/℃0-1-286426428y21x

0.5中心点的冷却曲线(不满足稳定性条件)x 100-0.5温度/℃-1-1.5-2-2.50100200300400500600时间/s