Matlab模拟2D非稳态导热数值解
clc;clear all;close all;
%定义几何尺寸、物性参数、网格大小、时间步长
xlength=1; %x方向长度
ylength=1; %y方向长度(高度）
a=1e-4; %导温系数
nx=15;%x方向网格数目
ny=15;%y方向网格数目
delta_x=xlength/nx;%x方向单个网格长度
delta_y=delta_x;
t=1600;%总时间
nt=150;%时间步的个数
delta_t=t/nt; %时间步
%给网格点编号，定义边界条件和初试条件
n=((xlength/delta_x)+1)^2; %点的数目
m=sqrt(n);%行上点数目或列上点数目
r=(t/delta_t)+1;%时间步数目;
Fo=a*delta_t/delta_x^2;%网格Fo数
if Fo<0.25
disp('[有稳定解]');
else
disp('[解发生震荡]');
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({['非稳态/瞬态导热'];['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
