2维非稳态导热的显式格式MatLab代码
%%致谢：感谢up主 @计算传热学大叔 无敌的教学讲解视频
clc
clear
tic
format longG
rho=100;
cp=1000;
k=10;
T0=3;
TW=3;
TE=5;
TN=5;
TS=3;
s=10;
t=10000;
x=3;
y=2;
dt=10;
dx=0.1;
dy=0.1;
sizet=t/dt;
sizex=x/dx;
sizey=y/dy;
Tmatrix=zeros(sizet,sizey,sizex);%上三维了
for i=1:1
Tmatrix(i,:,:)=3;
end
for i=2:sizet
for r=1:1%第一列
for j=1:1%矩形区域左下顶点，dxw=0.5dx,dys=0.5dy
ap0=dy*dx*rho*cp/dt-k*dy/dx-2*k*dy/dx-k*dx/dy-2*k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=k*dy/dx;
aw0=2*k*dy/dx;
as0=2*k*dx/dy;
an0=k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*TW+ae0/ap1*Tmatrix(i-1,j,r+1)+as0/ap1*TS+an0/ap1*Tmatrix(i-1,j+1,r)+s*dx*dy/ap1;
end
for j=2:sizey-1%矩形左边缘,dxw=0.5dx
ap0=dy*dx*rho*cp/dt-k*dy/dx-2*k*dy/dx-k*dx/dy-k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=k*dy/dx;
aw0=2*k*dy/dx;
as0=k*dx/dy;
an0=k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*TW+ae0/ap1*Tmatrix(i-1,j,r+1)+as0/ap1*Tmatrix(i-1,j-1,r)+an0/ap1*Tmatrix(i-1,j+1,r)+s*dx*dy/ap1;
end
for j=sizey:sizey%矩形左上顶点,dxw=0.5dx,dyn=0.5dy
ap0=dy*dx*rho*cp/dt-k*dy/dx-2*k*dy/dx-2*k*dx/dy-k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=k*dy/dx;
aw0=2*k*dy/dx;
as0=k*dx/dy;
an0=2*k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*TW+ae0/ap1*Tmatrix(i-1,j,r+1)+as0/ap1*Tmatrix(i-1,j-1,r)+an0/ap1*TN+s*dx*dy/ap1;
end
end
for r=2:sizex-1
for j=1:1%矩形区域下边中部，dys=0.5dy
ap0=dy*dx*rho*cp/dt-k*dy/dx-k*dy/dx-k*dx/dy-2*k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=k*dy/dx;
aw0=k*dy/dx;
as0=2*k*dx/dy;
an0=k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*Tmatrix(i-1,j,r-1)+ae0/ap1*Tmatrix(i-1,j,r+1)+as0/ap1*TS+an0/ap1*Tmatrix(i-1,j+1,r)+s*dx*dy/ap1;
end
for j=2:sizey-1%矩形中间主体
ap0=dy*dx*rho*cp/dt-k*dy/dx-k*dy/dx-k*dx/dy-k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=k*dy/dx;
aw0=k*dy/dx;
as0=k*dx/dy;
an0=k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*Tmatrix(i-1,j,r-1)+ae0/ap1*Tmatrix(i-1,j,r+1)+as0/ap1*Tmatrix(i-1,j-1,r)+an0/ap1*Tmatrix(i-1,j+1,r)+s*dx*dy/ap1;
end
for j=sizey:sizey%矩形上边缘中部，dyn=0.5dy
ap0=dy*dx*rho*cp/dt-k*dy/dx-k*dy/dx-2*k*dx/dy-k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=k*dy/dx;
aw0=k*dy/dx;
as0=k*dx/dy;
an0=2*k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*Tmatrix(i-1,j,r-1)+ae0/ap1*Tmatrix(i-1,j,r+1)+as0/ap1*Tmatrix(i-1,j-1,r)+an0/ap1*TN+s*dx*dy/ap1;
end
end
for r=sizex:sizex%右边缘
for j=1:1%矩形区域右下角，dys=0.5dy，dxe=0.5dx
ap0=dy*dx*rho*cp/dt-2*k*dy/dx-k*dy/dx-k*dx/dy-2*k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=2*k*dy/dx;
aw0=k*dy/dx;
as0=2*k*dx/dy;
an0=k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*Tmatrix(i-1,j,r-1)+ae0/ap1*TE+as0/ap1*TS+an0/ap1*Tmatrix(i-1,j+1,r)+s*dx*dy/ap1;
end
for j=2:sizey-1%矩形右边中间，dxe=0.5dx
ap0=dy*dx*rho*cp/dt-2*k*dy/dx-k*dy/dx-k*dx/dy-k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=2*k*dy/dx;
aw0=k*dy/dx;
as0=k*dx/dy;
an0=k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*Tmatrix(i-1,j,r-1)+ae0/ap1*TE+as0/ap1*Tmatrix(i-1,j-1,r)+an0/ap1*Tmatrix(i-1,j+1,r)+s*dx*dy/ap1;
end
for j=sizey:sizey%矩形右上角，dyn=0.5dy，dxe=0.5dx
ap0=dy*dx*rho*cp/dt-2*k*dy/dx-k*dy/dx-2*k*dx/dy-k*dx/dy;
%ap0=dy*dx*rho*cp/dt-k*dy/dxe-k*dy/dxw-k*dx/dyn-k*dx/dys;
ae0=2*k*dy/dx;
aw0=k*dy/dx;
as0=k*dx/dy;
an0=2*k*dx/dy;
ap1=dx*dy*rho*cp/dt;
Tmatrix(i,j,r)=ap0/ap1*Tmatrix(i-1,j,r)+aw0/ap1*Tmatrix(i-1,j,r-1)+ae0/ap1*TE+as0/ap1*Tmatrix(i-1,j-1,r)+an0/ap1*TN+s*dx*dy/ap1;
end
end
end
Tans=zeros(sizey,sizex);
for i=1:sizey
for j=1:sizex
Tans(i,j)=Tmatrix(sizet,i,j);
end
end
Tans=Tans(sizey:-1:1,:)
% 创建网格坐标
[X,Y]=meshgrid(1:size(Tans,2),1:size(Tans,1));
% 绘制曲面图
surf(X,Y,Tans);
% 设置图形属性
title('Surface Plot');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
% 添加颜色栏
colorbar;
toc
