2维非稳态导热显式格式数值求解MatLab程序
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,sizex,sizey);%上三维了
for i=1:1
Tmatrix(i,:,:)=3;
end
for i=2:sizet
for j=1:1
for r=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-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+1,r)+as0/ap1*TS+an0/ap1*Tmatrix(i-1,j,r+1)+s*dx*dy/ap1;
end
for r=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-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+1,r)+as0/ap1*Tmatrix(i-1,j,r-1)+an0/ap1*Tmatrix(i-1,j,r+1)+s*dx*dy/ap1;
end
for r=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-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+1,r)+as0/ap1*Tmatrix(i-1,j,r-1)+an0/ap1*TN+s*dx*dy/ap1;
end
end
for j=2:sizex-1
for r=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-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-1,r)+ae0/ap1*Tmatrix(i-1,j+1,r)+as0/ap1*TS+an0/ap1*Tmatrix(i-1,j,r+1)+s*dx*dy/ap1;
end
for r=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-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-1,r)+ae0/ap1*Tmatrix(i-1,j+1,r)+as0/ap1*Tmatrix(i-1,j,r-1)+an0/ap1*Tmatrix(i-1,j,r+1)+s*dx*dy/ap1;
end
for r=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-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-1,r)+ae0/ap1*Tmatrix(i-1,j+1,r)+as0/ap1*Tmatrix(i-1,j,r-1)+an0/ap1*TN+s*dx*dy/ap1;
end
end
for j=sizex:sizex%右边缘
for r=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-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-1,r)+ae0/ap1*TE+as0/ap1*TS+an0/ap1*Tmatrix(i-1,j,r+1)+s*dx*dy/ap1;
end
for r=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-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-1,r)+ae0/ap1*TE+as0/ap1*Tmatrix(i-1,j,r-1)+an0/ap1*Tmatrix(i-1,j,r+1)+s*dx*dy/ap1;
end
for r=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-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-1,r)+ae0/ap1*TE+as0/ap1*Tmatrix(i-1,j,r-1)+an0/ap1*TN+s*dx*dy/ap1;
end
end
Tans=zeros(sizex,sizey);
for i=1:sizex
for j=1:sizey
Tans(i,j)=Tmatrix(sizet,i,j);
end
end
Tans=Tans(sizex:-1:1,:);
end
% 创建网格坐标
[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
