%%致謝：感謝up主 @計算傳熱學(xué)大叔 無敵的教學(xué)講解視頻

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%矩形區(qū)域左下頂點，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%矩形區(qū)域下邊中部，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%矩形區(qū)域右下角，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,:)

% 創(chuàng)建網(wǎng)格坐標(biāo)

[X,Y]=meshgrid(1:size(Tans,2),1:size(Tans,1));

% 繪制曲面圖

surf(X,Y,Tans);

% 設(shè)置圖形屬性

title('Surface Plot');

xlabel('X-axis');

ylabel('Y-axis');

zlabel('Z-axis');

% 添加顏色欄

colorbar;

toc