clc;clear all;close all;

%定義幾何尺寸、物性參數(shù)、網(wǎng)格大小、時(shí)間步長

xlength=1; %x方向長度

ylength=1; %y方向長度(高度）

a=1e-4; %導(dǎo)溫系數(shù)

nx=15;%x方向網(wǎng)格數(shù)目

ny=15;%y方向網(wǎng)格數(shù)目

delta_x=xlength/nx;%x方向單個(gè)網(wǎng)格長度

delta_y=delta_x;

t=1600;%總時(shí)間

nt=150;%時(shí)間步的個(gè)數(shù)

delta_t=t/nt; %時(shí)間步

%給網(wǎng)格點(diǎn)編號(hào)，定義邊界條件和初試條件

n=((xlength/delta_x)+1)^2; %點(diǎn)的數(shù)目

m=sqrt(n);%行上點(diǎn)數(shù)目或列上點(diǎn)數(shù)目

r=(t/delta_t)+1;%時(shí)間步數(shù)目;

Fo=a*delta_t/delta_x^2;%網(wǎng)格Fo數(shù)

if Fo<0.25

disp('[有穩(wěn)定解]');

else

disp('[解發(fā)生震蕩]');

end

T=zeros(m,m,r);

Ti=200;%初始溫度

Ttop=100;

Tleft=100;

Tright=100;

Tbottom=100;

Tmax=max([Ttop,Tleft,Tright,Tbottom,Ti]);

for k=1:r

for i=1:m

for j=1:m

if(i==1)&&(j==1)

T(i,j,k)=(Tbottom+Tleft)/2;

elseif (i==1)&&(j==m)

T(i,j,k)=(Ttop+Tleft)/2;

elseif (i==m)&&(j==m)

T(i,j,k)=(Ttop+Tright)/2;

elseif (i==m)&&(j==1)

T(i,j,k)=(Tbottom+Tright)/2;

elseif (i==1)&&(j>1&&j<m)

T(i,j,k)=Tleft;

elseif (i==m)&&(j>1&&j<m)

T(i,j,k)=Tright;

elseif (j==m)&&(i>1&&i<m)

T(i,j,k)=Ttop;

elseif (j==1)&&(i>1&&i<m)

T(i,j,k)=Tbottom;

else

T(i,j,k)=Ti;

end

end

end

end

for k=1:r-1

for i=2:m-1

for j=2:m-1

T(i,j,k+1)=T(i,j,k)+Fo*(T(i+1,j,k)+T(i-1,j,k)+T(i,j+1,k)+T(i,j-1,k)-4.*T(i,j,k));

end

end

end

figure(1)

Tinistial=imagesc(T(:,:,1));colorbar

title(['Temperature Profile','time(t)=',num2str(0),'s',' Ttop=',num2str(Ttop),'^oC'])

set(gca,'xtick',[]);

xlabel(['Tbottom=',num2str(Tbottom),'^oC'])

yyaxis left

set(gca,'ytick',[]);

ylabel(['Tleft=',num2str(Tleft),'^oC'])

yyaxis right

set(gca,'ytick',[]);

ylabel(['Tright=',num2str(Tright),'^oC'])

figure(2)

%Tfinal=imagesc(T(:,:,r));colorbar;shading interp;

Tfinal=pcolor(T(:,:,r));colorbar;shading interp;

title(['Temperature Profile','time(t)=',num2str(t),'s',' Ttop=',num2str(Ttop),'^oC'])

set(gca,'xtick',[]);

xlabel(['Tbottom=',num2str(Tbottom),'^oC'])

yyaxis left

set(gca,'ytick',[]);

ylabel(['Tleft=',num2str(Tleft),'^oC'])

yyaxis right

set(gca,'ytick',[]);

ylabel(['Tright=',num2str(Tright),'^oC']);

figure(3)

x=1:m;

y=1:m;

for k=1:r

h=surf(x,y,T(:,:,1));colormap(jet);shading interp;colorbar;

axis([0 m 0 m 0 Tmax]);

title({['非穩(wěn)態(tài)/瞬態(tài)導(dǎo)熱'];['time(t)=',num2str((k-1)*delta_t),'s']});colorbar;

drawnow;pause(0.001);

refreshdata(h);

if k~=r

T(:,:,1)=T(:,:,k+1);

else

break;

end

end